Расчет искривленной пристеночной турбулентной струи Текст научной статьи по специальности «Физика»

Том XVI

УЧЕНЫЕ ЗАПИСКИ Ц А Г И

19 85

№ 3

УДК 532.525.6

РАСЧЕТ ИСКРИВЛЕННОЙ ПРИСТЕНОЧНОЙ ТУРБУЛЕНТНОЙ СТРУИ

В. М. Сыч

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

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

В расчетных исследованиях плоских пристеночных струй можно выделить три основных направления. Аналитические методы [1—4] применялись главным образом при расчетах ламинарных пристеночных струй несжимаемой жидкости. Рассчитывался основной участок струи (после смыкания вязких слоев), где течение близко к автомодельному. Интегральные методы позволяют рассчитывать турбулентные пристеночные струи (достаточно полная библиография приведена в работе [5]). К недостаткам этих методов относится необходимость задания эмпирических профилей газодинамических переменных и напряжения трения. Жидкость обычно считается несжимаемой. Кривизна поверхности в работе [5] по существу не учитывается. Наиболее перспективными в настоящее время представляются конечно-разностные методы, позволяющие без особых затруднений учитывать сжимаемость жидкости и использовать различные модели турбулентности.

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

газодинамических величин. Это связано с присутствием в поперечном профиле струи как участков со значительной турбулентной вязкостью (пристеночный пограничный слой, слой смешения), так и участков, где течение близко к невязкому (ядро струи). Поэтому большинство численных методов при расчетах струйных течений не дают удовлетворительных результатов. В работах [6—11] используется разностный метод Патанкара—Сполдинга [12], однако подробности расчетной схемы не приведены и неясно, был ли сделан расчет начального участка струи с относительно тонкими вязкими слоями. Проведенные нами расчеты по программе [12] в этом случае показали, что возникают значительные счетные осцилляции [8], искажающие решение. Предпринятая попытка раздельного расчета вязких слоев с итерационной стыковкой через ядро струи вела к значительному усложнению расчетной схемы и не обеспечивала достаточную точность при расчете невязкого ядра.

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

1. Рассмотрим полные уравнения Навье—Стокса в плоской ортогональной системе координат (х, у), связанной с криволинейной поверхностью (рис. 1), приведенные, например, в работе [13]. Оценивая

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

В рассматриваемой задаче в качестве малых параметров целесообразно взять величины в^ = VE|um и &2 =УЕ1К- Здесь уЕ, Уе-—координата и поперечная скорость на внешней границе струи, ит — максимальная продольная составляющая скорости в поперечном сечении, Я — радиус кривизны поверхности. В дальнейшем рассматривается изобарический режим истечения струи из сопла (точнее, квазиизобарический, с неизобаричностью из-за кривизны), поэтому основной малый параметр 81 имеет тот же порядок, как, например, в автомодельном турбулентном слое смешения (е^0,03-^-0,04). Этот подход позволяет, в частности, обосновать расчетное приближение для достаточно толстых струй.

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

дх и ди

1 + у1Я дх

1__д_

(1+у/Л)» ду

ди . ¡IV V—- --

ду Я + у,

+

1

др _

1 + у/Я дх

\ 1

др_

ду

дк . „ дН.

—- + V-

у/Я дх ду

ри5

Я + у

I 1 + >

др

+

1 + У/Я ду

эфф

дТ ду

у/Я дх

ду I

Я + у

(2)

(3)

(4)

К, этим уравнениям добавляются уравнение состояния, определение удельной энтальпии к и соотношения для эффективной и турбулентной вязкости.

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

дх\у

дх\

ри

ду

а<>

У

; (]> = ш2 фЕ = | рийу.

Тогда уравнение (1) удовлетворяется тождественно и вместо уравнений (2) — (4) получим:

0++('+*)

тЕ<& ди

^Е ди

+

2тЕи>3

?я

+

+

0+-П

2ш др

р и дх

д*>

ри

дю

Я I 2«ф| д<*

«м

я

о дН .

2(0 ~дГ +

тЕ дН

ди>

Я + у '

(5)

(6)

р и

дН

Рг У эфф 2шф|. ди>

]+

+ ■

д / и» |Г / и \ 1 ?(1 + я) ди Иэфф — \|^ФФ-\"р7]эфф ] 2Чя "лГ - —

да> \ Фе

Здесь Я — полная энтальпия.

Эффективные параметры определяются формулами

¡^эфф

(-^Л = + Рг =5; 0,7; Рг ^ Рг рг ^ Ргт '

/

.0,9.

В дальнейшем удобно перейти в уравнениях к безразмерным переменным. В качестве характерных величин берутся значение ит0 — продольной скорости в начальном сечении ядра струи и параметры р0, Та, с при нормальных атмосферных условиях. Тогда составляю-

ра

щие скорости и, V относятся к ит0, энтальпия Н— к Фтй, температура Т и плотность р —к Та и рв, коэффициенты вязкости р., Рт — к IV давление р — к линейные размеры — к ра/(ра ит 0),

функции тока <р, фя —к [Аа, массовый расход через внешнюю границу струи оте— к ра ит0- В этом случае уравнения будут иметь прежний вид, кроме соотношений

2, р — ЕирТ,

с* = Е —

и 1 .

Ср Та

Ра

'т О

2 '

Ра "то

7 = Ср!Съ-

Вместе с уравнениями (5) — (7) используем следующие определе-

ния:

у О)

ф СУс) = (й = (^)1/2' У =

Р"

о о

Массовый расход через внешнюю границу струи тЕ находится из соотношения, аналогичного [12]:

тЕ = -

с1х

ду

■у/Я)

(8)

ди

/? )~ду

Эта формула была получена из условия предельного перехода на внешней границе продольного уравнения импульсов в соотношение Бернулли. В общем случае такой переход некорректен. Однако в случае автомодельного течения в слое смешения и на основном участке струи справедливость формулы для тЕ может быть доказана строго. Поэтому следует ожидать достаточную точность соотношения (8) в условиях локальной автомодельное™, имеющей место в изобарических струях.

Для замыкания задачи остается привести выражения для ¡д, и цт-Зависимость коэффициента ламинарной вязкости от температуры считается степенной (х==Р'76. При расчете коэффициента турбулентной вязкости часто используется локальная модель турбулентности Пранд-тля [12]. Эта модель служит в рассматриваемой задаче для задания начального профиля |ят. Длина пути смешения вычисляется по характерным толщинам слоя смешения и пограничного слоя. У стенки используются соотношения Ван-Дриста с поправками Себеси и Смита [15]. В дальнейшем для расчетов применяется дифференциальная модель Коважного—Секундова [16]. В переменных (х, со) соответствующее уравнение имеет вид:

2св

дх

+

(О)'

д и>

д_(с дчТ <М ' ди,

С,=

1 I У + а1 Р" Я / 2О>^2е

(1) =0, 7т = [Хт/р,

а3 + 1

2ы

и

ди

(9)

В уравнения (5)—(9) поперечная скорость V явно не входит. После решения системы она может быть вычислена из уравнения (1) по формуле

V =

1 + У/Я

' дх )|ш

т

Значения V могут быть использованы при расчете течения вне струи и решения укороченных уравнений в следующем приближении ПО 81.

2. При численных расчетах наибольшую сложность представляет совместное решение уравнений сохранения импульса, поскольку в (5)

др а

величина ^"неизвестна и может быть найдена только из уравнения

(6). Используя численный метод [12], можно построить как итерационную, так и безытерационную (более экономную) схему решения уравнений (5) и (6). Однако численные эксперименты показывают появление счетных осцилляций в ядре струи на свободной границе, когда вязкие слои составляют малую часть поперечного профиля переменных (начальный участок). Это вызывает накапливание погрешности, существенно искажающей решение, и значительно затрудняет автоматизацию сеточного построения в расчетной области.

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

Исключим уравнения (5) с помощью уравнения (6) следую-

щим образом. Интегрируя (6) в пределах от со до 1, найдем выражение для р, которое затем дифференцируем по х с учетом (5). В результате получим интегродифференциальное уравнение

др

дх

дРЕ , (РЕ~Р) йЯ

+ Я йх + Я

дх

йу

1 +У1Я

(10)

где с точностью до членов —8162

_ Р-эфф Л , У\щди.

Мш 2

. 1 + уе!К 1+у/Я")

тельно

Дифференцируя уравнение (10) по у и решая линейное относи-дР

дх

др_ дх

уравнение, окончательно получим

(1 + уЕт *рван

(1 + у№ дх

я ах

и> /

Тогда вместо уравнения (5) имеем

~Ь 1*3 — ¿4 — 1^2 — ¿5,

(И)

где

¿i=l +

L3 =

2(о

Jü

yj 2(» —

R дх

1 4 Ув) dpE

' R ) dx

г ди

Г =

+

(1 +у IR)

4-1L ' R

[m^+{>+iff]

2<W- dR

da>

k

du>

i +ylR

k = \í

эфф

1 + -*-

R2 ¿x

Pu

í

ll(od<o

R I

Составляя для уравнения (11) монотонную разностную схему, аналогично [17] получим, что разностный оператор

¿4/

и,

1+а/

аппроксимирует дифференциальный оператор Ь5 локально со вторым порядком на неравномерной по со сетке в наиболее важных областях

с достаточной вязкостью. Там же, где функции меняются не-

существенно и приемлем первый порядок аппроксимации. Здесь

[ Г1 I Д(0 — Г1 (сг+1 - 2шг + <ог_г>

а1=~ i ~~~

2

2 k¡

Да)

[átf И,+1 — (gti + grü uí + 8.Ü ui-i\>

Д<о

ig? + gT)"i

i +

Дсо ;

sti = sr

u¿ + i — о>г-ь

У<+1 Я

1 + JÜ-R

ki + k; +1 2Kr иг)'

+

& ui-i

1 +

Уг—i

Si

R

kl + k¡_ j

(I gr 1 - gr) ¿o

4Ai ¡

yi + y¡+l} 1

IR

áV=M 1 +

yi

S\ i Si д

R da>

2 (<ог — (0._j) (I gr I + gr)

4*1^1

У." + Уг-1 ü/?

Г 1

L p(i+y/«)2J'

t — номер узла в поперечном направлении.

Оператор L6 аппроксимируется с первым порядком. Аппроксимация Lt и L3 проводится стандартным образом. Для аппроксимации в продольном направлении используется неявная схема с «весами». Второй порядок соответствует схеме Кранка—Никольсона.

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

Интегрирование уравнений для определения р [из уравнения (6)], и у проводится, как и в монографии [12], по правилу трапеций. Однако численные эксперименты показали, что предположение о линейном изменении функций между узлами не обеспечивает достаточную точность у внешней границы струи (где вычисляется тЕ). В дальнейшем

в этой области при интегрировании используется квадратичная интерполяция между узлами. Возникающая при этом сложность из-за необратимости формул интегрирования для фи у [18] преодолевалась ли-

у>-+1

у1+1

неаризацией Дф|+1 = | рийу по приращениям б(ры), и бУг (б/ =

у7'+1

I

= /'+1—/') при вычислении у! + 1 на новом слое /+1 по х. В результате для определения приращений поперечной координаты получается система двух линейных алгебраических уравнений. Причем для внутренних узлов это приращение, как и Д-ф,, определяется двумя способами из-за выбора при квадратичной интерполяции дополнительного узла как выше, так и ниже рассматриваемого интервала (г, ¿+1). В качестве окончательного значения берется среднеарифметическое. Указанный метод интегрирования уравнения неразрывности применяется в нескольких узлах у свободной границы и мало влияет на общую экономичность численной схемы, повышая вместе с тем устойчивость расчетов.

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

Решение уравнений осуществляется сразу во всем поперечном сечении струи. Оптимальное разбиение он (г/г) сохраняется на одном или нескольких шагах в маршевом направлении (по х), затем находится новое оптимальное разбиение, на которое профили газодинамических функций перестраиваются с использованием сплайн-интерполяции [19].

3. Численные расчеты турбулентной струи около криволинейной поверхности проводились при ширине щелевого сопла с1~уЕо^ 10 мм, что соответствует характерному числу Рейнольдса Ией—Ю5. Скорость на выходе из сопла задавалась не непосредственно, а находилась из изоэнтропических формул при задании параметров а = роф/ра и Т0ф = = 300 К в форкамере и расчетном истечении с учетом поправки на неизобаричность из-за кривизны. Начальный профиль скорости в слое смешения брался из автомодельного решения. В пограничном слое профиль степенной, гладко сопряженный с линейным в ламинарном подслое.

Наклон профиля у стенки и толщина пограничного слоя определяются числом Рейнольдса. На стенке использовались условия прилипания и изотермичности (ГЮ = 300 К). На внешней границе струи в общем случае задавалась скорость спутного потока иЕ и ТЕ =293 К. Число узлов в поперечном сечении в большинстве расчетов принималось N=50. При N=100 и 200 результаты менялись в пределах 1%. Продольный шаг увеличивался от Лх = 0,002 уЕ в начальном участке струи до Ах = 0,1 Уе к основному участку. Следует отметить, что расчет затопленной струи (иЕ=0) мог устойчиво осуществляться только с малыми Ах из-за осцилляций по шагам при вычислении тЕ. Увеличе-

ние скорости внешнего потока до Ue ^0,01 -ит о позволяет существенно увеличить продольный шаг и тем самым сократить общее время расчетов.

На рис. 2 показаны профили полного давления при R = 25 мм, d=8,2 мм, 8 = 30° и нескольких значениях о. На этом и последующих графиках для сравнения приведены также данные экспериментального исследования рассматриваемой пристеночной струи, полученные в ЦАГИ. Совпадение профилей полного давления в расчете и эксперименте достаточно хорошее вплоть до углов поворота струи, соответствующих безотрывному течению (9 = 60°-i-90o в зависимости от а).

На рис. 3 и 4 представлены поперечные профили статического давления р/ра. Разница между расчетом и экспериментом в этом случае заметно больше. В начальном сечении струи (см. рис. 3) это вызвано тем, что средняя кривизна линий тока меньше соответствующей радиусу цилиндра из-за конусности используемого в эксперименте сопла и конструктивного сдвига среза сопла относительно поверхности цилиндра. Последнее вызывает обратную кривизну линий тока и соответствующее повышение статического давления в пристеночном слое. Резкое (относительно расчетного) увеличение статического давления для 0=1,8 при 0 = 30° (см. рис. 4) связано с началом отрыва. Отметим, что профили статического давления наиболее чувствительны к отрыву.

Профили относительной скорости (отнесенной к максимальной в начальном сечении), пересчитанные по измеренным значениям р0 и р и полученные в расчетах, показаны на рис. 5 для сг=1,5 и нескольких значений угла 0. Соответствие данных удовлетворительное, если принять во внимание, что температура струи в эксперименте не измерялась и при пересчете считалось Т0—const. В слое смешения, особенно на его внешней границе, точность пересчета скорости невелика из-за того, что полное и статическое давления близки к атмосферному и становится существенной даже малая погрешность их измерения.

В целом сравнение с экспериментом показывает применимость предложенной расчетной методики при безотрывном обтекании струей

поверхности. Условия отрыва в расчетах (тю = 0) получаются только ( др\

при увеличении в случае уменьшения кривизны поверхности

> oj или при достаточной величине ^ в градиентном

спутном потоке. При обтекании затопленной струей цилиндра расчетное давление на поверхности монотонно возрастает, трение падает, однако отрыва при больших углах, наблюдаемого в эксперименте, не происходит. Это объясняется, по-видимому, неучетом изменения кривизны линий тока в уравнении (6) при исключении конвективных членов с поперечной скоростью. Использование последних соответствует следующему приближению по ei и требует существенного усложнения метода расчета.

5 «Ученые записки ЦАГИ» № 3

65

Рис. 4

У MM

10

5

О 0,2 О,U 0,6 0,3 1,0 1,2 и/и.,

пах

Рис. 5

Автор выражает признательность И. Н. Соколовой за проведение экспериментальных исследований.

1. Акатнов Н. И. Распространение плоской ламинарной струи несжимаемой жидкости вдоль твердой стенки. — Труды ЛПИ, техническая гидромеханика, 1953, № 5.

2. Акатнов Н. И., Сюй Мян ь-Ф ы н. Плоская полуограниченная струя по криволинейной поверхности. — 1962, ПМТФ, № 6.

3. Коробко В. И. Теория неавтомодельных струй вязкой жидкости. Изд-во Саратовского гос. университета, Саратов, 1977.

4. Kuwabara Kohseh. The laminar jet along curved wall. — J. of the Japan Soc. for aeron. and space sciences., 1980, vol. 28, N 323.

5. Петров А. В., Ш ел о м о в с к а я В. В. Метод расчета коэффициента импульса струи, потребного для ликвидации отрыва потока на профиле крыла. — Труды ЦАГИ, 1979, вып. 1977.

6. Irwin Н. P. А. Н., Smith P. Arnot. Prediction of the effect of streamline curvature on turbulence. — Physics of Fluids, 1975, vol. 18, N 6.

7. L j u b о j a M., R о d i W. Calculation of turbulent wall jets with an algebraic Reynolds stress model. — ASM.E J. of Fluid Eng., 1980, vol. 102.

8. G i b s о n M. M., R о d i W. A. Reynolds—stress closure model of turbulence applied to the calculation of a highly curved mixing layer. — J. of Fluid Mech., 1981, vol. 103.

9. Артюх JI. Ю., К а ш к a p о в В. П., Локтионова И. В. О численном моделировании турбулентной полуограниченной струи.—Тепло-массоперенос в жидкостях и газах, Алма-Ата, 1982.

10. Gibson М. М., Yonnist В. A. Modeling the curved turbulent wall jet. —AIAA J, 1982, vol. 20, N 12.

11. S harm a R. N., Patankar S. V. Numerical computations of wall-jet flow. —Jnt. J. Heat Mass Transfer, 1982, vol. 25, N 11.

12. Пантакар С., Сполдинг Д. Тепло- и массообмен в пограничных слоях:—М.: Энергия, 1971.

ЛИТЕРАТУРА

13. В а н-Д а й к М. Теория сжимаемого пограничного слоя во втором приближении. «Исследование гиперзвуковых течений». — М.: Мир, 1964.

14. Denny V., Landis R. An improved transformation of the Patankar-Spalding type for numerical solution of two-dimentional boundary layers flow. — Int. J. Heat Mass Transfer, 1971, vol. 14.

15. Cebeci Т., Smith A. M. O. Analisis of turbulent boundary layers. — Acad. Press, 1974, V.

16. Абрамович Г. H., Крашенинников С. Ю., Секун-д о в А. Н. Турбулентные течения при воздействии объемных сил и не-автомодельности. — М.: Машиностроение, 1975.

17. С а м а р с к и й А. А. Теория разностных схем. — М.: Наука,

1977.

18. Сыч В. М. Оптимальное распределение сеточных узлов в задачах струйной газодинамики.—Труды ЦАГИ, 1982, вып. 2169.

19. Завьялов Ю. С., Квасов Б. И., Мирошниченко В. Л. Методы сплайн-функций. — М.: Наука, 1980.

Рукопись поступила 9/XI 1983 г.