Том XXXIII
УЧЕНЫЕ ЗАПИСКИ ЦАГИ 20 0 2
№ 1—2
УДК 533.6.011.5:629.7.025.1 629.735.33.015.3.025.1
КРЫЛО МАКСИМАЛЬНОГО КАЧЕСТВА В СВЕРХЗВУКОВОМ ПОТОКЕ НА РЕЖИМЕ ВЯЗКО-НЕВЯЗКОГО ВЗАИМОДЕЙСТВИЯ
С. Д. ЖИВОТОВ, В. С. НИКОЛАЕВ
Решена вариационная задача об оптимальной форме срединной поверхности трапециевидного крыла на режиме сверхзвукового вязко-невязкого взаимодействия.
Предложен полуэмпирический интегральный метод расчета турбулентного пограничного слоя как модификация метода Лойцянского — Лапина. Проведены расчеты оптимальных форм и их аэродинамических характеристик при разных относительных толщинах крыла и радиусах затупления передней кромки, различных числах Маха и Рейнольдса при ламинарном и турбулентном режимах течения в пограничном слое.
При многопараметрических расчетах аэродинамических характеристик и решении задач оптимизации эффективны приближенные методы, основанные на методе полос и гипотезе локальности, когда местные коэффициенты давления и трения считаются зависящими лишь от ориентации элемента поверхности к вектору скорости набегающего потока, местного числа Рейнольдса и режима течения в пограничном слое (ламинарный, турбулентный).
В данной статье исследован режим вязко-невязкого сверхзвукового взаимодействия, при котором справедлива схема разделения течения на невязкую возмущенную зону и пограничный слой, а силовое воздействие потока на тело складывается из сил давления без учета вязкостных эффектов, добавочного давления, индуцированного вытесняющим воздействием пограничного слоя, и сил поверхностного трения. В рамках этих предположений в [1] предложен расчетный алгоритм для режима слабого взаимодействия, который был усовершенствован и распространен на режим умеренного взаимодействия в случае ламинарного пограничного слоя в [2].
Методы расчета сверхзвукового обтекания элементов летательных аппаратов на режиме турбулентного пограничного слоя исследовались в [3] — [9]. Численные методы расчета основаны на непосредственном решении уравнений в частных производных турбулентного пограничного слоя, их погрешность определяется в основном погрешностью гипотез, используемых для замыкания системы уравнений. Интегральные методы опираются на использование интегральных соотношений импульсов и энергии, модели турбулентности, лежащие в их основе, носят полуэмпирический характер, при их реализации прибегают к приемам, упрощающим расчеты. Однако они значительно проще и требуют меньше временных затрат, а лучшие интегральные методы не уступают по точности численным дифференциальным методам, особенно в случае малых градиентов давления.
В данной статье предложена модификация метода [5], [8], позволяющая уточнить формулу связи коэффициента трения с числом Рейнольдса при небольших значениях параметра трения.
С учетом влияния трения и толщины вытеснения в рамках метода локальной пластины и теории полос прямым численно-аналитическим методом Ритца получены систематические данные об оптимальных формах тонких крыльев как для ламинарного, так и для турбулентного пограничных слоев.
1. Интегральный полуэмпирический метод расчета турбулентного пограничного слоя.
Основополагающие подходы к приближенному расчету интегральных характеристик
сжимаемого турбулентного пограничного слоя сформулированы и разработаны в рабое Франкля
— Войшеля [3]. Учет сжимаемости в [3] ограничен случаем небольших чисел Маха и малого отличия температуры стенки от температуры набегающего потока: при расчетах зависимости коэффициента трения на плоской пластине от числа Яе использовано разложение по двум малым параметрам. Лойцянский и Лапин [5], [8], используя идеи [3], развили полуэмпирический интегральный метод расчета сжимаемого турбулентного пограничного слоя на пластине в случае произвольного числа М и произвольной температуры стенки. При расчете интегралов для толщины потери импульса и толщины вытеснения они использовали предположение о большой величине параметра трения й1=уе/ут, где ут — динамическая скорость, индекс «е» относится к величинам на внешней границе пограничного слоя. Интегралы для толщины потери импульса 5** и толщины вытеснения 5* (точнее для Яе** = реие5 **/ц,е и Яе* = реие5 */ц,е) они представили в виде асимптотических рядов и нашли два первых члена разложения. Они были использованы авторами для определения формпараметра Н* = 5 */ 5 ** и в законе сопротивления
— связи Яе** и ^ (или с у), С = . I/ , где индекс <т» относится к условиям на стенке,
V / сУ ре
2т
су=—^ — локальный коэффициент поверхностного трения. Однако при построении
Рвие
необходимой для расчета аэродинамических сил зависимости су от Яех = реиех/Це авторы [5],
[8] использовали лишь один главный член асимптотического ряда для Яе**, а для компенсации так подобрали константу интегрирования, чтобы в случае несжимаемой жидкости получить согласование с эмпирической формулой Кармана, полученной на основе опытных данных Кемпфа:
0,242 =0,41 + 1в (Яе х с у ) .
Последнее упрощение вносит в методику [5], [8] дополнительные предположения эмпирического характера, которые представляются необязательными и устранимыми в рамках идей [5], [8]. С этой целью в данной статье предпринята модификация метода [5], [8].
При расчете интегралов в выражениях для Яе** и Яе* в настоящей статье, как и в [8], использована линейная зависимость полной энергии от скорости (интеграл Крокко), которая при числе Прандтля Рг ^ 1 основана на приближенном выражении для температуры теплоизолированной стенки Т :
Г да
Т„г = Те (1+г ^ м2),
где г = Рг13 — коэффициент восстановления, у — отношение удельных теплоемкостей. В расчетах настоящей статьи г = 0,9; у = 1,4, а значения универсальных эмпирических констант турбулентности к = 0,4, а = 11,5. Приведем зависимости температуры и плотности от скорости
Т ГУ гу^1 м2 1 + гу^1 ме
Т = ^ = 1 _юи _р[ и |; р= 2 ‘ ; ш= 1 _ ^ е
V ие У
сИ+—м21 |1+—м2
Здесь ^ = Тк / Т0 — температурный фактор, Т0 — температура торможения набегающего потока. Далее в формулах для Яе** и Яе* используется некоторый интеграл I, зависящий от распределения температуры по пограничному слою, см. [8]:
2В + ю arcsin^ - arcsin
ю
ю
>/4р
ю
Приведем полученные в [8] и используемые в настоящей статье выражения для Яе**, Яе* и формпараметра Н*. Эти выражения определены с относительной ошибкой порядка 1/^2, т. е. удержаны два первых члена асимптотических рядов:
Re** = e_________а Mw e<-i
к Me
Re* = е-Каа^*ек^ 1 ( 1+Р
1 -
2 - 1,5ю-р 1 -ю-р
Л
к Me
1-ю-р
С ( ю^
р 3-1,5ю-р+—
1 Ч 2р.
кС(1+ру 1-ю-р
H * =
1+р
1 -р-ю
к(1 + р)С
Для получения зависимости Rex от с/ (или С) используется интеграл
Re**
Re. = -^ } С 2d Re
**
Р
w 0
В настоящей статье в отличие от [8] при последнем интегрировании удержан не один, а два члена асимптотического представления, что позволило получить выражение для Яех с той же относительной погрешностью 1/^2, что и приведенные выше выражения для Яе** и Яе*. В итоге имеем
Rex =
_ e а Mw pe екС'1 С2
к Me Pw
1-
Г2-1,5ю-р 2 ^ к+
к^ 1-ш-р к-I
В [8] выражение в квадратных скобках принято равным 1 и вводится «компенсирующий» множитель сь выбираемый из условия согласования зависимости су от Яех с эмпирической формулой Кармана для несжимаемой жидкости.
Подставляя в последнюю формулу численные значения констант турбулентности, выражение для параметра трения ^ и логарифмируя, получим
, /т-. ч 0,2457 • I»/1 - ш - р 1 / \ 1
lg(Re x ) = -----------f=--------р - lg(cf ) + lg
y/Cf
Me
- 0,2381 -
0,7677,Jc~
f2 - 1,5ю-р 2^
В соответствующей формуле [8] отсутствует последнее слагаемое правой части, а вместо предпоследнего слагаемого (-0,2381) стоит (-0,41), как следствие упомянутой «компенсации» (кроме того, имеется небольшое различие в множителе первого слагаемого, в [8] он равен 0,242). Последнее соотношение весьма удобно для итерационного процесса ввиду того, что каждое из слагаемых есть монотонно убывающая функция су.
Зависимость Яе** от су (закон сопротивления), используемая обычно при сопоставлении моделей турбулентности, совпадает с формулой [8], кроме небольшой разницы в коэффициенте первого слагаемого
, , _ 0,2457• IJ1 -ю-р ,
lg(Re* *) =-------------р=--------+ lg
4cf
Vw^
Me
0,7677JCf
- 0,5391-------3LL (2-1,5ю-р).
1-ю-р
1
В [8] при расчете интеграла для Яе** был опущен первый пристеночный участок интегрирования (вязкий подслой), а второй (турбулентное ядро) — продолжен до стенки. В настоящей статье проведено асимптотическое исследование возможности такого подхода. В итоге получена следующая асимптотическая оценка для относительной погрешности при вычислении Яе**, вносимой вышеупомянутым упрощающим приемом:
ак ка ка ,
-------Є + Є -1
При принятых в статье значениях а и к эта погрешность равна —130е
-0,4-/^
и стремится к
нулю при стремлении ^ к да. Далее в работе универсальные константы турбулентности а и к не используются, а а применяется для обозначения угла атаки.
Интегральное соотношение импульсов при нулевом градиенте давления и отсутствии вдува дает следующую связь между толщиной потери импульса и коэффициентом трения:
7 гч
а 5 а х
Рємє2
і£
2
где г«, — напряжение трения на стенке. При расчете на режиме вязко-невязкого взаимодействия, когда учитывается добавочное давление, индуцированное вытесняющим воздействием пограничного слоя, необходимо знать величину а5*Шх. Так как 5*=5**Н*, имеем:
а о *
йх йх
а 5 ** н * +5 ** ан *
а х
Из асимптотического представления величин 5*, 5**, Н* следует, что отношение второго слагаемого к первому имеет порядок 1/^2, и в рамках метода этим слагаемым можно пренебречь. В итоге имеем:
а 5
йх
а 5 н * _ н * с' /
ах 2
н * _-
1+р
1 -ш-р
Г, , 3,536^(1 -ш-р) I
||1 + 1+р V' •
Толщина вытеснения 5* в турбулентном пограничном слое при сверхзвуковой скорости в невязкой возмущенной зоне численно намного больше толщины потери импульса 5**. Величина формпараметра Н*=5 */5** возрастает с ростом температурного фактора и особенно с ростом числа Ме. Так, при Ме = 4 величина Н* в зависимости от температурного фактора меняется от 4 до 7, а при Ме = 8 — от 15 до 25.
Коэффициент добавочного индуцированного давления содного порядка с коэффициентом трения. Так, в рамках линейной теории
ао * н * с
/
-1ах л/М
-1
В рамках гиперзвуковой стабилизации (Мдаа >> 1, а — местный геометрический угол атаки) из формулы метода касательных клиньев следует
а 5 *
с ^ = 2(у +1)а--= (у +1)аН*сг .
ах
Таким образом, при определении влияния вязких эффектов на локальные аэродинамические характеристики трение и индуцированное давление следует учитывать одновременно: как следует из двух последних формул, срі и су близки и по порядку, и численно. На режиме слабого
X
взаимодействия дополнительный учет индуцированного давления не представляет каких-либо вычислительных трудностей и лишних временных затрат ввиду непосредственной аналитической связи между а5 */йХ и с/. В настоящей статье оба эффекта учитываются при расчетах аэродинамических характеристик тонких крыльев.
В статье проведено тестирование предложенного метода сравнением с другими расчетными методами в широком диапазоне чисел Рейнольдса. На рис. 1 представлены данные расчета коэф-
Рис. 1. Зависимость коэффициента трения на пластине в несжимаемом газе от числа Рейнольдса:
— — настоящий метод, • — формула Кармана, V — закон 1/7
фициента трения на пластине в несжимаемом газе по методу настоящей работы, формуле Кармана и по степенной зависимости:
с/ _■
0,0263
Яе-
1/7
На рис. 2 представлены данные расчета су в сжимаемом газе по методу настоящей статьи и по полуэмпирической теории Ван Дриста II [4], достаточно хорошо согласующейся с результатами многочисленных экспериментов [6] и рекомендованной для тестирования моделей турбулентности по результатам исследований, проводимых в рамках специальной международной программы под эгидой Станфордского университета [9]. Последние сравнения проведены при числах Ме=4, 6 и 8 и при значениях / 1У1Г (ґкг _ Ткг/То ), равных 0,2; 0,6 и 1. При расчетах показатель степени в степенной зависимости коэффициента вязкости от температуры принимался равным т = 0,76.
567В9 10 56789 10 56789
а) б) в)
Рис. 2. Зависимости коэффициента трения на пластине в сжимаемом газе от числа Рейнольдса:
----- / /^ = 0,2,-------/^ / /^ = 0,6 , ----- / / / = 1 (расчет по настоящему методу);
• — ^ ^ = 0,2, ■ — /^ / /^ = 0,6, А — / / / = 1 (расчет методом Ван Дриста)
Соответствующие расчеты коэффициентов трения Су и индуцированного давления с . при
ламинарном режиме течения в пограничном слое проводились по методике [1], [2].
2. Оптимальное изолированное крыло максимального качества. В качестве объекта исследования выбрано трапециевидное крыло заданной формы в плане: угол стреловидности по передней кромке равен % = 55°, задняя кромка перпендикулярна плоскости симметрии, отношение концевой хорды к корневой Ь1/Ь0 = 0,25. Рассматривались крылья с острыми передними кромками при параболическом законе распределения толщин профилей по потоку и при затупленных передних кромках с постоянным радиусом затупления по передней кромке, что обеспечивало постоянство локального потока тепла к кромкам. Температура поверхности тела принималась постоянной (предположение абсолютно теплопроводного тела). Отношение радиуса затупления к длине корневой хорды Ь0 однозначно связано с относительной толщиной крыла г (отношение толщины корневого сечения к корневой хорде):
Распределение толщины крыла при варьировании формы его срединной поверхности считалось неизменным. В связанной системе координат координаты у верхней и нижней поверхностей выражаются через половину местной толщины у и координату срединной поверхности ys: у = ± у1. Для затупленного крыла
где г — координата г концевой хорды. Относительные толщины профилей по потоку были одинаковы для тупых и острых кромок.
Были подсчитаны аэродинамические характеристики таких крыльев с плоской (недеформированной) формой срединной поверхности и при том же распределении толщин — с деформированной формой срединной поверхности. Эта деформация подбиралась так, чтобы величина аэродинамического качества была максимальной при некотором угле атаки, т. е. помимо параметров, характеризующих форму срединной поверхности, в качестве подбираемых параметров был и угол атаки.
Решение данной вариационной задачи проводилось прямым методом Ритца. В качестве минимизирующей последовательности выбиралась последовательность полиномов от двух переменных (х и г) специального вида, обеспечивающего выполнение граничных условий, а именно: передняя кромка и корневая хорда лежат в одной плоскости (х, г). От этой плоскости отсчитывается угол атаки крыла а, который в процессе решения оптимизационной задачи варьировался вместе с коэффициентами полинома (в данной постановке ограничились 10 коэффициентами, т. е. суммарное число варьируемых параметров — 11). В итоге решения поставленной задачи получим крыло, обладающее круткой, т. е. профили крыла по потоку расположены под разными углами атаки, а средние линии этих профилей имеют ненулевую кривизну.
Сравнение величины максимального аэродинамического качества Ктах крыла оптимальной формы с максимальным качеством недеформированного плоского крыла показывает, что величина выигрыша в Ктах зависит от многих факторов: геометрических (острая или тупая кромка, относительная толщина крыла) и аэродинамических (числа Маха и Рейнольдса, температурный фактор).
Диапазон рассмотренных чисел Рейнольдса был разным для ламинарного и турбулентного пограничного слоя, однако для сравнения были проведены расчеты и при одинаковых числах
27г2 ¿0 32008 %
Рейнольдса. Число Рейнольдса крыла Яе0 определялось по скорости и плотности набегающего потока, коэффициенту вязкости при температуре торможения набегающего потока и по длине корневой хорды. При определении безразмерных аэродинамических коэффициентов соответствующие силы относились к скоростному напору набегающего потока и площади крыла в плане.
Рассмотрено три варианта сочетания газодинамических параметров (числа Маха набегающего потока Мш, температуры торможения набегающего потока Т0 и температурного фактора
^ — отношения температуры поверхности к температуре торможения Т0):
1. Мш = 5, То = 1400 К, = 0,7;
2. Мш=10, Т = 5000 К, ^ = 0,2;
3. Мш =20, Т = 20 000 К, К = 0,05.
В дальнейшем при обозначении вариантов ограничимся одной меткой — указанием на число Мш. Получены данные для четырех чисел Яе0=104, 105, 106, 107 в случае ламинарного пограничного слоя (при неучете вязких эффектов формально Яе0=<х>) и чисел Яе0, равных 106, 107, 108, 109, в случае турбулентного пограничного слоя. Таким образом, при числах Рейнольдса, равных 106 и 107, имеются расчеты и для ламинарного, и для турбулентного слоев.
Основным расчетным вариантом был следующий: Ъх /Ъ0 _ 0,25, х_ 0,05, тупые передние
кромки, _ 10 . Кроме того, проводился анализ влияния сужения (Ъх /Ъ0 _ 0,1 и 0,5) при сохранении угла стреловидности % _ 55°, относительной толщины (х _ 0,02 и 0,08) и заострения кромки в случае ламинарного слоя. Помимо задачи о максимуме качества Ктах рассматривалась задача о минимуме сопротивления при заданной подъемной силе (суа _ 0,1 и 0,25).
Некоторые результаты расчетов представлены на рис. 3 — 6. На рис. 3 приведены зависимости аэродинамического качества от угла атаки при различных числах Рейнольдса для оптимальной деформации срединной поверхности крыла и для недеформированного крыла (число Мш = 10, Ъ\1Ъ0 = 0,25, х = 0,05, затупленные кромки, ламинарный слой). На рис. 4 приведены зависимости аэродинамического качества от угла атаки для оптимального и недеформированного крыльев в случаях ламинарного и турбулентного режимов течения в пограничном слое (Яе0=106, Мш=10, Ъ1/Ъ0 = 0,25, х = 0,05, затупленные кромки). Зависимость максимального качества оптимального и недеформированного крыльев от числа Рейнольдса для ламинарного и турбулентного режимов течения в пограничном слое представлена на рис. 5 (Мш=10, Ъ\/Ъ0 = 0,25, х=0,05, затупленные кромки). Зависимость от числа Рейнольдса относительного выигрыша в величине максимального качества оптимального крыла по сравнению с недеформированным крылом для различных чисел Маха дана на рис. 6 (Ъ1/Ъ0=0,25, х=0,05, затупленные кромки, ламинарный и турбулентный режимы). С ростом числа М относительный выигрыш заметно возрастает. Выигрыш за-
м»=ю
Рис. 3. Сравнение зависимостей аэродинамического качества от угла атаки для оптимального крыла и недеформированного крыла при разных значениях чисел Рейнольдса:
— оптимальное крыло,------недеформированное крыло;
1 — Яе0 = 104, 2 — Яе0 = 105, 3 — Яе0 = 106, 4 — невязкое обтекание
Рис. 5. Зависимость максимального качества оптимального крыла и недеформированного крыла от числа Рейнольдса при разных режимах течения в пограничном слое:
— оптимальное крыло,-------недеформированное крыло;
1 — ламинарный слой, 2 — турбулентный слой, 3 — невязкое обтекание
Рис. 4. Сравнение зависимостей аэродинамического качества от угла атаки для оптимального крыла и недеформированного крыла при разных режимах течения в пограничном слое:
— оптимальное крыло,------недеформированное крыло;
1 — ламинарный слой, 2 - турбулентный слой
О----------------1 ! І І І III--------1 І І І І II!-----------1 І I ( IIII------------1 ...........І'-------1 і ИІ-------------------1 ’ I П ГП
Ю4 Ю5 1 06 1 07 1 0е 109
Рис. 6. Зависимость относительного выигрыша аэродинамического качества за счет деформаций от числа
Рейнольдса:
— ламинарный слой,--------турбулентный слой; 1 — Мад = 5,
2 — Мад = 10, 3 — Мад = 20, 4 — невязкое обтекание
метно снижается с уменьшением числа Яе, в чем сказывается меньшая чувствительность аэродинамических характеристик к изменению формы при увеличении влияния вязких эффектов. Относительный выигрыш в величине максимального качества также возрастает при увеличении относительной толщины. Так, при значениях г, равных 0,02; 0,05; 0,08, относительный выигрыш ЛК равен соответственно 1,5%, 7,6%, 12,8% (затупленные кромки, Мда=10, Ь1/Ь0 = 0,25, Яе0 =да).
Обратная тенденция по числам Рейнольдса имеет место для оптимального крыла с острыми
кромками: увеличение значения AK при Re0=104 равно 13%, а при Re0=<x> оно равно 5% (Мш=10, Ъ\/Ъ0 = 0,25, х = 0,05, ламинарный слой).
На рис. 7, 8 представлены некоторые данные расчета формы оптимальных крыльев (крутка, кривизна профилей) при числе Мш=10 и относительной толщине х = 0,05 (bi/b0 = 0,25, затупленные кромки, Т0 = 5000 К, tw = 0,2). На рис. 7 приведена зависимость крутки профилей крыла по
Рис. 7. Зависимость крутки профилей крыла по размаху при разных режимах течения в пограничном слое:
турбулентный слой, Re0 = 10 ;
ламинарный
слой, Re0 = 10 ;....невязкое обтекание
Рис. 8. Зависимость кривизны профилей крыла по размаху при разных режимах течения в пограничном слое:
— турбулентный слой, Reo = 107; - - - ламинарный
слой, Reo = 105;..невязкое обтекание
размаху, т. е. (асеч-а) как функция а^. На рис. 8 приведена зависимость от а/а1 кривизны профилей крыла. За меру кривизны профилей срединной поверхности у, принималась следующая
величина. При каждом значении а находился тах| у — улин |, где улин — хорда профиля срединной
*
поверхности. Пусть этот максимум достигается при значении х . Тогда за меру кривизны
(У5 УДЫМ ~ ~
принимаем величину----------------, где хМос и ххв — координаты передней и задней кромок
Ххв — Хмос
профиля (положительная кривизна соответствует выпуклости вверх). Как следует из анализа полученных данных, оптимальная крутка в рассчитанных случаях весьма мала (порядка 1°). В то же время оптимальная кривизна вполне соизмерима с относительной толщиной профиля (здесь половина толщины 0,025) и даже превышает ее.
В рассматриваемых случаях подобная деформация привела к следующему увеличению максимального качества по сравнению с недеформированным крылом. При турбулентном слое (Яс0 = 107) — ДК=0,16 (5,7%), при ламинарном слое (Яс0 = 105) — ДК=0,19 (6,9%), при невязком обтекании (Ясо=<х>) — ДК=0,22 (7,6%).
ЛИТЕРАТУРА
1. Александров В. Ю., Галкин В. С., Нерсесов Г. Г., Николаев В. С. Приближенный метод аэродинамического расчета летательных аппаратов при больших сверхзвуковых скоростях полета//Труды ЦАГИ.— 1990. Вып. 2492.
2. Животов С. Д., Николаев В. С. Метод расчета локальных аэротермодинамических характеристик плоских тел в сверхзвуковом потоке с учетом влияния вязкости и граничных условий скольжения//Ученые записки ЦАГИ.— 1998. T.XXIX, №1—2.
3. Франкль Ф., Войшель В. Трение в турбулентном пограничном слое около пластины в плоскопараллельном потоке сжимаемого газа при больших скоростях//Труды ЦАГИ.— 1937. Вып. 321.
4. Van Driest E. The problem of aerodynamic heating//Aeron. Eng. Review.— 1956. Vol. 15, N 10.
5. Лойцянский Л. Г., Лапин Ю. В. Применение метода Кармана к расчету турбулентного пограничного слоя на пластине в газовом потоке//Труды ЛПИ им. М. Н. Калинина.— М.— Л.: Машгиз.— 1961, № 217.
6. Hopkins E. J., Inouye M. An evaluation of theories for predicting turbulent skin friction and heat transfer on flat plate at supersonic and hypersonic Mach number//AIAA Journal.— 1971. Vol. 9, N 6.
7. Ре пик Е. У., Романов В. Н. Оценка эффективности расчетных методов определения коэффициента турбулентного поверхностного трения при сверхзвуковых и гиперзвуковых числах М//Труды ЦАГИ. —1980. Вып. 2084.
8. Лапин Ю. В. Турбулентный пограничный слой в сверхзвуковых потоках газа.— М.: Наука.— 1982.
9. Bradshaw P., Launder B. E., Lumley J. L. Collaborative testing of turbulence models//AIAA Paper.— 1991, N 91-0215.
Рукопись поступила 9/VII2000 г.