Применение обратного метода расчета трехмерного пограничного слоя к задаче обтекания крыла с учетом влияния вязкости Текст научной статьи по специальности «Физика»

УЧЕНЫЕ ЗАПИСКИ Ц А Г И Т о м XX* 198 9

№ 5

УДК 629.735.33.015.3.025.1 : 532.526 532.526—3

ПРИМЕНЕНИЕ ОБРАТНОГО МЕТОДА РАСЧЕТА ТРЕХМЕРНОГО ПОГРАНИЧНОГО СЛОЯ К ЗАДАЧЕ ОБТЕКАНИЯ КРЫЛА С УЧЕТОМ влияния вязкости

О. В. Карась, В. Е. Ковалев

Разработан обратный метод расчета трехмерного сжимаемого ламинарного и турбулентного пограничного слоя, используя численную схему предиктор — корректор и схему Келлера. На основе полуобратного метода расчета вязко-невязкого взаимодействия разработан алгоритм расчета обтекания крыла при трансзвуковых скоростях полета с учетом влияния вязкости, включая тонкие зоны отрыва. Приведен пример расчета изолированного крыла. Показаны полученные в расчете на соответствующих режимах полета срывные зоны в окрестности задней кромки и в области скачка уплотнения.

При анализе современных аэродинамических компоновок крыла представляет значительный интерес расчет как безотрывных течений, так и сложных трехмерных отрывных течений. На реальных конфигурациях отрывные течения могут существовать в окрестности задней кромки крыла на верхней поверхности, замкнутые срывные зоны вблизи скачка уплотнения или в области подрезки на нижней поверхности при малых углах атаки. Характерные замкнутые срывные зоны возникают при трубных числах Яе в окрестности передней кромки в областях ламинарного отрыва с последующим переходом и присоединением.

Значительный прогресс в развитии методов расчета таких течений достигнут в связи с использованием обратных методов в теории пограничного слоя, которые позволяют в процессе расчета проходить точку нулевого тренш?. В этом случае в качестве граничного условия (кроме условия прилипания) задается либо толщина вытеснения, либо коэффициент трения на стенке, в отличие от классической задачи, где задается скорость на внешней границе пограничного слоя.

Основной прогресс при расчете обтекания профилей вязким трансзвуковым потоком газа был достигнут с помощью так называемого полуобратного метода расчета [1—3]. Для реализации этого метода требуется обратный алгоритм расчета пограничного слоя и прямой метод расчета невязкого течения.

В настоящее время известно ограниченное число работ, где предложены обратные методы расчета трехмерного пограничного слоя. В работах [4, 5] предложены обратные методы расчета пограничного слоя на скользящем крыле. Общий обратный метод расчета трехмерного пограничного слоя изложен в работе [6]. В работе [5] с помощью полуобратного метода расчета вязко-невязкого взаимодействия решена задача о ламинарной короткой срывной зоне в окрестности передней кромки на скользящем крыле.

В работе предложен конечно-разностный обратный метод расчета трехмерного пограничного слоя и дано приложение этого алгоритма к расчету крыла конечного размаха при трансзвуковых скоростях полета методом вязко-невязкого взаимодействия, включая слабые зоны отрыва. Изложенный ниже обратный метод расчета построен на основе метода, разработанного в работе [7], где приведен алгоритм расчета двумерных отрывных течений путем решения обратной задачи для ламинарного несжимаемого пограничного слоя с заданной толщиной вытеснения.

1. Конечно-разностный метод расчета трехмерного пограничного слоя. Исходная система дифференциальных уравнений, описывающая течение газа в трехмерном пограничном слое, имеет вид:

¿(p«A2sln6)+ ¿ (рда/ijsin 6) + A (pvAjft* sin 6) = 0;

и du . w du , — du , ci t. oi л f ,

P л7 Ж + Ж + +Pcsc6^'ffi' +

, csc2 6 dp . ctg 0 ese G dp . d ! du —

. +pkauw = —HT-di + —hT~-5; + ^\*'dj~pB *)

и dw , w dw i — dw , с, , , f., „ ,

РлТ^ + РлТи + Р^ - P ctS ^ "Ь P esc bklu + j (i)

, ctg 6 csc 0 dp csc20 dp . d l dw ——,\

+ pftall№, = -*4— ----+ J '

и dH , w dH , —dH P 1ц ~дх ~1~ P h2 dz ?V dy

где pv = pv + p' v'; ut = (и2 -4- w2 + 2uw cos 8)t/2.

Координата у направлена по нормали к поверхности крыла, переменные х, z определяют систему неортогональных координатных линий с углом Q(x, z) между ними на поверхности, и, v, w — составляющие вектора скорости по координатам х, у, г, р — плотность, р — дав-

Í. ds, , dst

ление, Ц — коэффициент динамическои ВЯЗКОСТИ, «1 = gj- , п2 = -fe----------

метрические коэффициенты, Sj, s2 — длина дуги вдоль координатных линий X и Z.

Параметры kit k2, ka, hi, характеризующие кривизны координатных линий, выражаются в виде:

1 [ d ,, dhi'

kt =

[¿~(Л2С08в)-^];

[■¿(A» cos 6)-gf];

/гх Л2 sin

1 Г d ... оч dh

hi Ä2 sin I

Граничные условия имеют вид:

у = 0, и = т = V = 0, дН\д у = О у = 8 и=ие = (х, г), т = чюе{х,г),

3 — толщина пограничного слоя.

Для замыкания системы (1) используется концепция турбулентной вязкости, согласно которой связь между касательными напряжениями Рейнольдса и профилем скорости осредненного течения выражается следующим образом:

В работе используется двухслойная модель турбулентной вязкости Себечи — Смита, модифицированная для расчета отрывных течений.

В качестве начальных условий для расчета параметров течения на передней кромке используется система уравнений, описывающая течение на линии торможения [8]. Для расчета начальных условий в корневом сечении крыла могут быть использованы система уравнений для плоскости симметрии или для бесконечного скользящего крыла [8].

В настоящее время методы расчета трехмерного пограничного слоя, построенные на основе конечно-разностных шаблонов Келлера — Себечи, получили широкое распространение [8—И]. В данной работе используется метод, построенный на основе схем предиктор — корректор и схемы Келлера. Указанный метод имеет второй порядок точности и сочетает безытерационность метода предиктор — корректор и простоту реализации схемы Келлера.

Первоначально общая система уравнений (1) (без уравнения энергии), используя преобразование:

^2 —'^2^1 ^*12-----^12^-'» ^*21---^21

где Ие = Мга ¿/V,

и определение двухкомпонентного векторного потенциала:

Рик2 вШ 6 = дР/дУ; Р1УЬ18тв = дв/д V;

удовлетворяющего уравнению неразрывности, приводится к системе двух нелинейных дифференциальных уравнений третьего порядка.

Согласно метода Келлера [12] эта система уравнений преобразуется к системе из шести дифференциальных уравнений первого порядка:

Р-9и% + р2№д£-рАдр

дХ

+ §)У+Л~Р„& +

Р,&$!!+Р,

__ р др

-- г 27

+ л8-Ц-+(^)';

Г' — Р19 р/7; и' = I/; ;1РдЛ-.Р.Г£.л.Ш

РЧдХ

ч

дР

+ й)т+м =

дХ ^2а дг + (^),5

С' = Р297\Г;

№' = Т,

(3)

где

£ = - Рм р[/2 + Р1Б р"1Г2 + Р16 р£/ №; М = - Я24 Р^2 + Р28 ?и2 + Р26 Р^;

коэффициенты Р\ — Ргэ — зависят от параметров системы координат /г1( /г2, &1, /гг, &12, й21, 0 и могут быть легко определены.

— штрих (’) —дифференцирование по У.

Граничные условия (2) запишем в виде:

у=о, и= 1г = /7 = о = о,

У = Уе1, и = ие (х, 2:)/Иое, № = (х, г)/«со-

(4)

В этом разделе уравнение энергии не решается, а распределение температуры определяется из интеграла Крокко для теплоизолированной стенки [13]:

Т1Те=>\+0,2гМ1[1-(щ1и(еУ),

где г = 0. 84 — для ламинарных течений; г = 0. 89 — для турбулентных течений.

Система нелинейных дифференциальных уравнений (3) решается численным методом второго порядка, построенным на основе метода предиктор — корректор и схемы Келлера. На рис. 1 изображен конечно-разностный шаблон, применяемый для разностной аппрокси-

Ша 6л он предиктор ,п-1,т

Шаблон корректор ,п-1,т

п,т-1

п,т-1

Рис. 1

□ производные по X о >’ 1

х * » У

/ расчетная точка

мации системы (3). Функции и, V, б, Ш и Т аппроксимируются в узлах сетки X», гт, У{. На этапе предиктор система (3) аппроксимируется следующим образом:

[Г]

/ --о-

= р

19

] -

„_±,т

г,п — 1, т I г Г1П — 1, т Гд№П 4 пп — 1, т гтчл — Ь 1 4

1рС/|/ч _я

_ Л

г/лгмГ 2>т_Го о о7,тш тд°Лп-1'т

[(ЬТ) ^ — |Р27 - Р28 ^ Р2 р У7 ^ + Р3 Т д2^_х_ ~

(5)

_Мп~*•».

1

' 2

На этом шаге неизвестными являются сеточные функции Л 2 ,

которые определяются из решения линейной системы алгебраических уравнений.

На шаге корректор система (3) после конечно-разностной аппроксимации имеет вид:

Г С*1 \п> т пЯ. ш Гг 71Я» Ш

{гх1± = рп ;

1 2 1 2

1 1111 п — , т п~ тг • т гт,-,п—гг » т п — ~,т п — 7Г,тг,-.__п —- , т

г) 2’ г/71 2 ГаиП 2 п 2 ГХ1 2 Г^Т 2

[р ]/-1 И/-4 3 [л/-4 к]у-4

_Л

\thT\n 2’т \о дР О дР 0-\Г/д]ХГ , П 7-^1" 2

[(¿>7) ] , —\Р^дх Р2 2дг Р2р1Г д2 + Р3 ^¿2 1

/ л I. J/ л

1

п —— , т

■М. г ■

(6)

На этапе корректор неизвестными являются значения сеточных функций т, ___ Первоначально реализуется шаг предиктор на всей

линии х = сопз1 по размаху. Затем, таким же образом проводится шаг корректор.

При решении классической (прямой) задачи в результате разностной аппроксимации системы уравнений (3) и граничных условий (4) на этапах прогноза и коррекции получаем 6/4-6 (/—число интервалов по У) линейных алгебраических уравненийчотносительно 6 7 + 6 неизвестных (6/ уравнений дает система (5) или (6) и 6 уравнений следуют из граничных условий).

Полученная система алгебраических уравнений, как и в методе Келлера — Себечи, имеет трехдиагональную блочную структуру и решается блочным методом исключения Гаусса.

Необходимо отметить, что схема Келлера и ее модификации приводят к системам нелинейных алгебраических уравнений, которые для решения требуют привлечения итерационных алгоритмов.

Разработанная схема расчета трехмерного пограничного слоя, как и схема Келлера — Себечи, имеет второй порядок точности по всем трем координатам. Однако по сравнению с методом Келлера — Себечи предагаемый алгоритм расчета имеет ряд преимуществ. Во-первых, он приводит к системе линейных алгебраических уравнений и, как следствие, к безытерационному алгоритму, во-вторых, расчет коэффициентов этой системы менее трудоемок.

2. Обратные задачи трехмерного пограничного слоя. При решении обратной задачи компоненты скорости, а, следовательно, и градиенты давления дР/дХ и дР/дЪ на внешней границе пограничного слоя — неизвестны. Поэтому в этом случае в качестве неизвестных добавляются еще две величины дР/дХ и дР/д! и, соответственно, для замыкания системы уравнений необходимо добавить еще два уравнения, связывающие скорость на внешней границе с градиентами давления (здесь используются уравнения Бернулли в предположении изэнтро-пичности течения).

Соответствующая система линейных алгебраических уравнений, (5) или (6) может быть записана в компактной матрично-векторной форме:

[Я] + (7)

Матрица [/?] имеет блочную трехдиагональную структуру.

Вектор А (решение системы уравнений (7)) может быть представлен в виде суммы трех векторов [7]:

+ + (8>

У, X, Z -решения систем уравнений, отличающихся только правыми частями:

[/>] К = /?; [Я] *=/?,; [0]г = Яг.

При решении задачи с заданным коэффициентом трения, величи-дР дР

ны дХ и Ш 0ПРеДеляются из граничных условии:

С/яУЯё = К0;

С/2 УТе « 2рт Т0

и выражений для У0 и Т0 из (8). Затем из (8) получим значения искомых сеточных функций.

Этот алгоритм может быть применен как для расчета безотрывных пограничных слоев, так и для расчета пограничных слоев с отрывом и присоединением. Для расчета отрывных течений применяется аппроксимация [14], которая в настоящее время используется для расчета двумерных областей отрыва, где скорости возвратного течения

относительно невелики 0,11 . Смысл этого приближения состоит

в исключении конвективного члена в уравнении импульса, что приводит к численной устойчивости решения при интегрировании против потока в зоне возвратного течения при сравнительно невысокой погрешности аппроксимации.

При решении обратной задачи с заданной толщиной вытеснения выражения для дР/дХ и дР/дЪ определяются из граничных условий:

8 ;/^=увРвг/в-^/р19;

ъ:/П~е=Уере\Уе-Ое!Рм (9)

и выражений для ве, ие и №е из (9).

Отметим, что в данном методе можно задавать в качестве граничных условий одну компоненту С/* или е; и значение поперечной скорости №е■ Эта постановка используется в задаче вязко-невязкого взаимодействия, изложенной в следующем разделе.

3. Расчет вязко-невязкого взаимодействия на крыле конечного размаха. Обратный метод расчета смешанного сжимаемого трехмерного пограничного слоя используется в методе расчета трансзвукового обтекания крыла с учетом влияния вязкости. В данной работе алгоритм расчета крыла построен на основе так называемого полуобрат-ного метода [1]. В отличие от классической схемы, используемый метод позволяет рассчитывать отрывные течения. Для практической реализации этого метода необходимы три основные компоненты:

1. Метод расчета трехмерного потенциального обтекания.

2. Обратный метод расчета пространственного пограничного слоя.

3. Метод сращивания потенциального решения с решением в пограничном слое.

В данном случае для расчета потенциального обтекания применяется модификация метода [15]. Принятое в исходной версии метода преобразование физической области в расчетную ведет к существенному уменьшению числа расчетных узлов на поверхности крыла по мере движения от корня к законцовке для крыльев с сужением, что крайне неудобно при использовании этого метода в задаче вязко-невязкого взаимодействия. Поэтому в модифицированной версии метода принято новое преобразование, которое приводит к приемлемой координатной сетке на поверхности крыла, где расчетные узлы совпадают с линиями равных процентов. Расчет потенциального течения осуществляется на основе решения полного уравнения для потенциала ф в квазилинейной форме. Итерационное решение конечно-разностного аналога уравнения для потенциала осуществляется методом релаксации вдоль линии в' криволинейных координатах, соответствующих отображению пространства вне крыла на внутренность прямоугольного параллелепипеда.

Схема расчета вязкого обтекания крыла по полуоб-ратному методу показана на рис. 2. Расчет проводится в следующей последова-

Рис. 2

тельности. На первом этапе рассчитывается потенциальное течение и определяется распределение скоростей на поверхности крыла. Затем рассчитывается пограничный слой на большей части поверхности прямым методом, а в зонах, близких к отрыву и в областях отрыва — обратным методом с заданным коэффициентом трения (компонентой С,х). В результате этого расчета получаем распределение скоростей на внешней границе пограничного слоя. Далее уточняется потенциальный расчет с модифицированным граничным условием на поверхности крыла.

Соответствующее граничное условие определяется из условия равенства нормальных скоростей на поверхности, разделяющей внешнее (потенциальное) и внутреннее (пограничный слой) решения. Это условие приводит либо к концепции поверхности вытеснения, либо к изменению граничного условия на поверхности крыла (уФО) . В первом случае поверхность, где выполняется условие непротекания, совмещается с поверхностью вытеснения. Чтобы найти вытесняющую поверхность в трехмерном случае, необходимо решить дополнительное дифференциальное уравнение [16]. В этом случае возникает необходимость в пересчете геометрии крыла в итерационном процессе вязко-невязкого взаимодействия, что ведет к снижению быстродействия метода. Учитывая, что расчет потенциального обтекания- занимает большую часть времени, становятся очевидными преимущества второго пути сращивания, при котором нет необходимости определять поверхность вытеснения и, соответственно, изменять исходную геометрию в итерационном процессе вязко-невязкого взаимодействия. При этом потенциальное обтекание модифицируется посредством изменения граничного условия на поверхности крыла [16]:

1

V =

ре /г, h2 sin (

д (ие fe h2bxsinb) д (We ре h-y 5г sin 0) дх дг

где 8* и 8* — толщины вытеснения по направлениям z = const и х = = const.

На следующем этапе задаваемая компонента коэффициента трения модифицируется по разности полей скоростей (компоненты ие), полученных из расчета потенциального течения и пограничного слоя обратным методом:

(C/,1/Re)(i+1)= (C^|/Re)(i) — kv{uve-Ue) . '

Значение коэффициента Kv определяет сходимость итерационного процесса и его значение для турбулентного пограничного слоя должно находиться в пределах 5-ь 15.

Учет эффектов, связанных с перепадом давления поперек пограничного слоя в окрестности скачка и задней кромки и нормальной скорости в следе за крылом, в данной работе не производился.

4. Примеры расчета. На основе изложенного алгоритма составлена программа расчета вязкого обтекания крыла при трансзвуковых скоростях полета.

На рис. 3 представлены результаты расчета вязкого обтекания стреловидного крыла с профилем Piercy (с = 0-12) в сечении при числах М = О и Re = 0,4-108. Этот случай интересен тем, что крыло имеет малое удлинение Я = 3, большую стреловидность х = 45° и течение около него существенно трехмерно. Линия перехода ламинарного пограничного слоя в турбулентный в расчете задавалась в окрестности мини-

Верхняя поверхность, а=8°

Внешние линии с тока у

Предельные линии тока

сс = 8°

Потенциальный расчет

О

0,5

Рис. 3

Верхняя поверхность крыли

—■— линия перехода

-------Внешние линии тока

-------пристеночные линии тока

Замкнутая область отрыВа за скачком

ъ=П,8

Рис. 4

мума давления. Дано сравнение рассчитанных зависимостей су(и) и Су сеч (2) с соответствующими экспериментальными данными работы [17]. Согласование расчетных и экспериментальных зависимостей су Сеч {г) при а = 2° и 4° хорошее. Некоторое расхождение расчетной зависимости сеч (г) с экспериментальной при а=8° по-видимому связано с развитием протяженного (по хорде) концевого срыва, обычно имеющем место на крыльях такого класса, начиная с некоторого угла атаки. Отмеченное явление усиливается с увеличением угла атаки и находится за пределами применимости предлагаемого метода расчета.

Приведенный пример показывает, что положенные в основу алгоритма численные методы позволяют расчетным путем получить основные аэродинамические характеристики на безотрывных режимах с достаточной для практических приложений точностью.

На рис. 4, 5 приведен пример расчета изолированного крыла (удлинением ?ь=10, с переменной толщиной с = 0,10—0,15 и с заданной круткой по размаху) пассажирского самолета. Расчет проведен при М = = 0,84, Re=2,5-106, а=4,8°. На верхней поверхности крыла на таком режиме существуют две зоны отрыва: одна замкнутая — в области скачка уплотнения, вторая — в окрестности задней кромки.

На рис. 4 приведено распределение компонент коэффициента трения (cfx и CfZ) и распределение давления в сечении 2=0,8 этого крыла. Приведен результат потенциального расчета, выполненный при одном и том же с вязким расчетом = 0,535. Смещение скачка к передней кромке в вязком потоке по сравнению с потенциальным достигает ~8% на консольной части крыла.

На рис. 5 показан вид линий F = const в сечении z = 0,8 на верхней поверхности крыла. При значении х = 0,6 видна характерная замкнутая срывная зона в окрестности скачка уплотнения, а в области задней кромки — зона диффузорного отрыва.

Последний пример иллюстрирует возможности предлагаемого метода для расчета сложных отрывных течений. Однако необходимо отметить, что авторы не претендуют на расчет взаимодействия скачка уплотнения с пограничным слоем в данной постановке, а предлагаемый алгоритм в основном направлен на расчет диффузорного отрыва, часто возникающего на крейсерских режимах при трубных числах Re. Пределы применимости данного метода и его возможности могут быть определены в процессе его использования в практических приложениях и сравнением с экспериментальными данными в широком диапазоне режимов полета.

Расчет последнего примера включает 20 вязких итераций. Время расчета около 3 часов на ЭВМ. с быстродействием ~ 106 оп/с. (Потенциальный расчет выполняется на сетке 192x32x24. При расчете пограничного слоя используется сетка 60x21 на одной поверхности крыла и 40 узлов по нормали к поверхности.)

ЛИТЕРАТУРА

1. Le Baleur J. С. Strong matching method for computing transóme viscous flows including wakes and separation. — Lifting airfoils. La Recherche Aerospat'iale, 1981, N 3.

2. Va ts а V. N. and Verdón J. M. Viscid/inviscid interaction análisis of separating traling-edge flows. — AIAA Journal, 1985, vol. 23, N 4.

3. V a n D a 1 s e m W. P. and S t e g e r J. L. Finite-differense simula-tion of transonic separated flow using a full potential-boundary layer ap-

proach. — AIAA-83-1689, AIAA 16 fluid and plasma dynamics conference, July, 1984.

4. Radwan S. F. and Lekoudis S. G. Boundary layer calculatios in the inverse mode for incompressible flows over infinite swept wings. — AIAA J., 1984, vol. 22, № 6.

5. Davis R. L., С a r t e r J. E. and Reshotko E. Analysis of transitional separation bubbles on infinite swept wings. — AIAA-85-1685, AIAA 18 fluid and plasma dynamics and lasers conference, July 1985.

6. D e 1 e г у J. M. and Formery M. J. A finite difference method for inverse solutions of 3-D turbulent boundary-layer flow. AIAA 83-0301, AIAA 21-st aerospase sciences meeting, Jan. 10—13, 1983/Reno, Nevada.

7. Cebeci T. Keller H. B. and Williams P. G. Separating boundary layer flow calculations. — Journal of Computational Phisics, 1979, vol. 31.

8. Cebeci Т., К a u p s K. A., Ramsey J. A. A general method for calculating three-dimensional compressible laminare and turbulent boundary layers on arbitrary wings. — NASA CR 2777, 1977.

9. С e b e с i T. Calculation of three-dimensional boundary layers. —

II. Three dimensional flows in cartesian coordinates. — AIAA J., 1975,

vol. 13, N 8.

10. L inch F. T. Resent applications of advanced computational methods in the aerodynamics design of transport aircraft configuration. — Aeronautical J., 1978, N 19.

11. Смит П. Д. Численный расчет трехмерных пограничных слоев.— В кн.: «Трехмерные турбулентные пограничные слои»./Под ред. X. Ферн-хольца и Е. Краузе. — М.: Мир, 1985.

12. Keller Н. В. Numerical methods in boundary-layer theory.— Ann. Rev. Fluid Mech., 1978, vol. 10.

13. Шлихтинг Г. Теория пограничного слоя. — М.: Наука, 1974.

14. R е у h n е г Т. A., F 1 u g g е-I о t z I. The integration of a chock vawe with a laminar boundary layer. — Int. Journal of Non-Linear Mech.— 1968, vol. 3, N 2, June.

15. J a m e s о n А., С a и g h e у D. A. Numerical calculation of the transonic flow past a swept wing. — N.-Y. University. ERDA report COO-3077-140.

16. Bradshow P., Cebeci T. Engineering calculation methods of turbulent flows. — Hemispere-McGraw Hill, Washington, 1983.

17. Falkner V. М., Lehrian D. E. Low-speed measurements of the pressure distribution at the surfase of a swept-back wing. — ARC RM 2741, 1953.

Рукопись поступила 21/VJI ]988