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

УЧЕНЫЕ ЗАПИСКИ Ц А Г И Том XIII 1982

№ 4

УДК 532.526.011.55

МЕТОД РАСЧЕТА СВЕРХЗВУКОВЫХ СОПЛ ПРИ СИЛЬНОМ ВЛИЯНИИ ВЯЗКОСТИ

О. В. Денисенко

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

При достаточно больших числах Ие, когда толщина пограничного слоя мала по сравнению с поперечным размером невязкого ядра, задача о профилированном сопле сводится к расчету поправки на толщину вытеснения пограничного слоя. Более сложным является случай, когда толщина пограничного слоя сравнима или даже больше невязкого ядра потока. Здесь уже согласно работе [1] существенное влияние на течение оказывают различные эффекты второго порядка малости: поперечная и продольная кривизна, скорость скольжения и скачок температуры на стенке и т. д. Кроме того, в уравнении для поперечной составляющей импульса следует учитывать центробежную силу. В результате градиент давления поперек пограничного слоя отличен от нуля, и при расчете маршевым методом необходимо применять специальные меры, чтобы построить решение, соответствующее условиям отсутствия возмущений на выходе сопла [2—4].

Предлагаемый в данной работе метод является по существу модификацией метода работы [1]. Вместо „естественных координат“ были использованы переменные Мизеса (функция тока — продольная координата), что позволило в расчетной области получить фиксированную сетку и добиться автоматического выполнения условия равенства расходов газа в заданном для невязкого и в профилированном для вязкого газа соплах. Применение для расчета пограничного слоя неявного метода второго порядка аппроксимации по поперечной координате дало возможность значительно увеличить шаг по продольной координате (по сравнению с работой [1]), т. е. сократить время расчета. В настоящее время

существуют различные методы регуляризации упрощенных уравнений вязкого газа для устранения в дозвуковой области некорректности задачи Коши при интегрировании системы маршевым методом [3, 4]. Применительно к уравнениям пограничного слоя второго порядка точности в данной работе предложен способ аппроксимации продольного градиента давления поперек пограничного слоя, который по существу является одним из возможных методов регуляризации интегрируемой системы уравнений. Используемая методика позволила получить не только градиентное сращивание численных решений в вязкой и невязкой областях, но и добиться совпадения решений уравнений пограничного слоя с решением невязких уравнений в области малого влияния вязкости.

1. Сверхзвуковое ламинарное течение вязкого газа в плоском или осесимметричном сопле. Считаем, что сопло, предназначенное для невязкого газа и обеспечивающее необходимое однородное ядро потока, задано, причем контур этого сопла не имеет разрыва кривизны или излома образующей, т. е. является гладким. В предположении, что поток вязкого газа можно подразделить на невязкое ядро и ламинарный пограничный слой, ставится обратная задача: найти форму сопла, которая обеспечивает совпадение течения в невязкой части потока с соответствующей частью течения в сопле известной формы, предназначенном для невязкого газа. Пусть г*, I/*, р*, а*, а*—соответственно радиус критического сечения сопла, скорость, плотность, коэффициент вязкости и скорость звука (Vn. = a.^ в невязком ядре критического сечения.

Введем следующие обозначения: VV№ — полная скорость;

p?*Vl — давление; ррй — плотность; hV\ — энтальпия; ¡j^* — коэффициент вязкости; хг.л, — продольная координата; yr* •—расстояние от оси сопла до линии тока; 6 — угол наклона линии тока к оси сопла; Рг — число Прандтля; f—показатель адиабаты; v = 0 или 1—соответственно плоский или осесимметричный случай; aV* — скорость звука: М= Vja — число М; о г* — толщина пограничного слоя; фр* г4'1 V* — функция тока; Re* = р* V* rj'^ — число Рейнольдса. Предположим, что произведение Re* х достаточно велико. Тогда величину 8¡х можно считать малым параметром и уравнения Навье — Стокса после отбрасывания членов с относительным порядком (3. х)2 преобразуются к следующему виду:

р V cos 6 -f cos б — =

дх дх

р У ^7 ^ р уъ V —j 4- ¿ii + 3i;

+ g -i + ч ;

7 — * и Р = ----------Р h ;

Т

= Re'.2 о У' l/sin 6 ; — = — Rei "ру 'cos ö ; дх dv

Ие12Рг

V вш 6 . , дк

---------\хь V---------

сое 6 д '

_д_

дх

„т т дк \ , ,,, й дк

Iх Р У V — + р У 'V — «• —

(? V , г д .и.

и,--------------V —!-

й : д:

ДЯ г) Йл: с»:

I/2 •

Рг/

При выводе системы уравнений (1) вместо функции тока была введена более удобная для расчета уравнений пограничного слоя нормированная поперечная координата С — Ие*2 [1/(у + 1) — Ф1- Уравнения (1) представляют собой уравнения пограничного слоя второго порядка точности. Группа членов, обозначенных посредством з, и о2, определяет эффект влияния продольной кривизны пограничного слоя на течение, а группа членов, обозначенных и g.2,— прочие эффекты второго порядка. При <з1=а2=§-]=£2=0 и др\д’С,—О уравнения (1) являются обычными уравнениями Прандтля, записанными в переменных Мизеса.

2. Расчет течения невязкого газа. Считаем, что контур сопла невязкого ядра задан. Чтобы исключить особенность в поведении производных от искомых функций по поперечной координате вблизи оси, вместо координаты * введем новую переменную &:

V-}-! _______

= (У+ 1)

(2)

После простых преобразований систему (1) при Ре*= оо представим в виде

¿10 у эш3 0 (М2— 1)соз 6/

дх у(М2_1) (М ? — 1) V \

др_

дх

» ЯШ 6 СОБ 6 Р V2

(м':-- к у

.. д Г др •' - ,>х

?Р_________

дг

р I ' п 0 / у У'1 др

дТ

у у д 0

¿7

Щ - 1

2 1/2 СОЗ 0

/г

Мх = М соз 0 ; р

г!

уз __ 2 ~

7

М.г ■ 7+ 1 2(7-1)

о Н ;

1

д_Ь_

дг

(3)

д

дм

+ 1) Р У‘У СОБ 6.

Граничные условия для системы (3) следующие: на поверхности сопла £ — 1 и 6 =: 6^ (х), где 0^ (х) — заданная функция, а на оси сопла имеют место условия симметрии 1 = 0, 0 = 0 и др/д% — 0. Если в течении Мг>1, то система (3) в полуполосе {0-<£<^1, х>0; является гиперболической и для ее численного решения можно использовать маршевые методы. В данной работе для расчета невязкого ядра применялась явная схема „бегущего счета“ [5]. Схема имеет первый порядок аппроксимации по обеим переменным, и для устойчивого счета должно выполняться следующее необходимое ограничение на соотношение шагов

1)[1 + 1§0(М=-1Г1-'}

р V соэ 0

(4)

V м2 - 1

Здесь х —величина шага по продольной координате, а Д$ — по поперечной. При расчете величина шага тконтролировалась, и при условии выполнения неравенства (4) счет был всегда устойчив. Заметим, что для точного расчета уравнений пограничного слоя (1) приходится использовать значительно меньший шаг по продольной координате, чем шаг, допустимый из критерия устойчивости (4), так как при увеличении числа М величина т растет как ]/" М2—1 и достигает больших значений.

После того, как на новом слое определены угол наклона линий тока 0, давление, скорость, энтальпия и плотность, рассчитывалась функция у (?) (т. е. расстояние от линии тока до оси сопла) путем численного интегрирования методом трапеций интеграла

Как правило, 50— 100 узлов вполне достаточно для расчета невязкого ядра с точностью до 2%, нри этом время расчета невязкого ядра течения на ЭВМ БЭСМ-6 при общем числе шагов около 1500 менее 5 минут. Расчет начинался с критического сечения и предполагалось, что в этом сечении звуковая поверхность является плоской, а значения всех параметров в невязкой части постоянны на этой поверхности. Чтобы не раскрывать особенности в уравнениях (3) при М-> 1, полагалось, что число М в критическом сечении отлично от единицы и равно М2=1 + е,. Остальные величины находились из законов сохранения энергии, массы и уравнения состояния. Величина е, во время расчетов варьировалась, и оказалось, что при г, = 10-2 -ь- 10-5 влияние параметра г, на рассчитываемые величины менее сотой доли процента.

3. Построение профилированного сопла. Рассчитав с помощью схемы „бегущего счета“ невязкое течение, выберем некоторую линию тока Чп невязкого потока. Если при всех д: > 0 на этой линии тока влиянием вязкости можно пренебречь, то, решив уравнения (1) в полуполосе {;6[1, ?„], л: > 0}, мы получим контур профилированного сопла, который обеспечивает, по крайней мере между линией тока 1п и осью сопла, такое же невязкое ядро, как и исходное сопло-Естественно принять эту линию тока %п за внешнюю границу пограничного слоя и потребовать, чтобы на ней скорость V, энтальпия /г, давление р, плотность р, угол наклона 0 и функция у вязкого течения совпадали с соответствующими величинами невязкого потока. Ясно, что при Ие* -> оо величина расхода газа через пограничный слой стремится к нулю и за внешнюю границу пограничного СЛОЯ МОЖНО принять контур ИСХОДНОГО невязкого Сопла (£„ = ^1 = 1)-Но по мере уменьшения числа Ие*, особенно когда величина расхода газа через пограничный слой становится сравнимой с общим расходом газа [Ие* = О (Ю3)|, внешняя граница отодвигается на значительное расстояние от заданного невязкого контура (Ё„<0,5). В качестве критерия выбора линии тока 1п в работе принималось, что отношение вязкого трения на внешней границе пограничного слоя к трению на стенке должно быть меньше некоторого малого параметра е5. При расчете, задавая величину г5 в пределах 10~2 — 10~3, на каждом шаге поставленное условие проверялось, и в случае невыполнения внешняя граница пограничного слоя отодвигалась ближе к оси сопла, а расчет начинался сначала.

(ч+ 1) ыг

р VСОБ в

о

Определение градиента давления при наличии перепада давления поперек пограничного слоя—наиболее ответственная операция численного алгоритма. Дело в том, что по внешнему виду система (1) аналогична уравнениям пограничного слоя на режиме свободного взаимодействия [2], и, следовательно, имеет однопараметрическое семейство решений. Каждое решение из этого семейства можно трактовать как решение, соответствующее некоторым условиям ниже по потоку. Поэтому, если на каждом шаге по л; мы будем итерировать давление (а значит и др/дх), то из-за неизбежных ошибок округления довольно быстро (через несколько шагов по х) получим решение с отрывом или сильным разрежением потока (др!дх=—х) в некоторой точке. Во избежание этого давление на новом слое рассчитывалось до начала итерационного процесса, и все входящие в выражение для др!дС величины, за исключением угла 0, брались с предыдущего слоя. Численное интегрирование велось от внешней границы пограничного слоя с применением метода трапеций для соответствующего уравнения системы (1). Затем на каждой линии тока вязкого течения рассчитывался градиент давления др ';дх пи формуле ( др/дх)\ —(р\ — /?а)/'с, что соответствует аппроксимации первого порядка по продольной координате.

Заметим, что выбранная в качестве внешней границы пограничного слоя линия тока 1п по существу является таковой только в окрестности выходного сечения сопла. Вначале, на разгонном участке, пограничный слой тонок, и большая часть потока в вязкой области является фактически невязкой. С другой стороны, уравнения (1) являются композитными по отношению к уравнениям (3), т. е., кроме членов, входящих в уравнения (3), содержат еще и дополнительные, обусловленные тензором вязких напряжений. Поэтому, если в области пограничного слоя, там, где влияние вязкости мало, давление, градиент давления и угол наклона 6 будут совпадать с соответствующими величинами невязкого потока, то и решение уравнений (1) будет совпадать с решением уравнений (3). Для этого на каждой линии тока вязкого течения С* угол наклона линейно интерполировался по значениям угла 0„ в невязком потоке. Точнее, координаты С* по формуле (2) пересчитывались как бы в невязкие линии тока а затем проводилась линейная интерполяция угла 0Й по рассчитанным из уравнений (3) углам наклона 0„. ГIо найденным значениям 0Йопределялась кривизна дЬ/дх с первым порядком аппроксимации по х: < дЬ/дх ) Ц=( 0* — — 0*)/", и далее рассчитывалось давление на новом слое. Хотя с принятой при выводе системы уравнений (1) точностью угол наклона линий тока 0 в пограничном слое можно считать заданной функцией только продольной координаты (например, в качестве этой функции можно взять угол наклона заданного невязкого контура сопла или угол наклона внешней границы пограничного слоя), вводимая по описанной выше методике зависимость угла наклона © от поперечной координаты позволяет получить не только градиентное сращивание численных решений в окрестности внешней границы пограничного слоя, но и добиться совпадения решения уравнений (1) с решением невязких уравнений в области несущественного влияния вязкости. На рис. 1 приведены профили давления, скорости I' энтальпии в промежуточном сечении л: = 0,54 при расчете сопла на М = 5, Не* = 995. Расчет проводился с числом узлов Л'=51 и Лг1/=81 соответственно для невязкой и вязкой областей.

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

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

системы (3) получаем немонотонное поведение дв/дх и д^/дх2 вдоль линий тока, что при интегрировании уравнений пограничного слоя (1) приводит к колебаниям рассчитываемых характеристик, например, трения на стенке, теплопередачи и тангенса угла наклона профилированного сопла. Чтобы получить гладкое поведение дЬ/дх и д2Ь/дх2, в работе исходный контур невязкого сопла интерполировался кубическими сплайнами с проведением предварительной процедуры сглаживания массива 0а,й.

4. Расчет уравнений пограничного слоя. Для системы уравнений (1) на стенке ставились граничные условия скольжения и скачка температуры. Согласно работам [6 — 8] для упругих сфер и коэффициентов аккомодации и отражения, равных единице, граничные условия скольжения имеют следующий вид:

V-

/ ¡ху'‘У д\/ ' 1}е1/2 \ а д:

/-А / ¡а у‘ V дН Ие12 \ а д"

*„=1,24 у П. ■

1,П

Рг(7 + 1)

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

В случае граничных условий прилипания, чтобы исключить особенность в поведении производных по поперечной координате, делалась замена г) = [2;]1/2. Тогда на поверхности сопла т(=0, 1/=0, Н = В предположении, что вязкость зависит только от температуры, в работе для определения (/г) использовался либо закон Сазерленда, либо интерполяционная формула работы [9].

Для расчета уравнений пограничного слоя (1) применялся неявный метод потоковой прогонки работы [10]. Метод имеет второй порядок аппроксимации по поперечной координате даже на неравномерной сетке. Кроме того, в процессе прогонки автоматически рассчитывается первая производная искомой функции по поперечной координате. Это позволяет простым образом аппроксимировать в выражениях для gl и g■, члены со смешанными производными. Например, на т-]-1 итерации принималось, что

/ГА- \1,т+1 , ,ч1'т(^>т (д/1\1-т ( V1’т - VI)

\ д* ' дх /к 'У К \дк )к \й: )к т

(Г '¿)\т /дУу.т /дУу

1 .ТО* ~

Линеаризация остальных членов в уравнениях (1) проводилась стандартным способом. Для производных по продольной координате х в работе использовалась аппроксимация первого порядка. Итерация начиналась с расчета V, затем 'а и завершалась определением расстояния у от линии тока до оси с подстановкой в выражение для у (С, х) вновь найденных на этой итерации значений скорости и плотности. Функция у находилась интегрированием методом трапеций следующего выражения:

'к

0

У'+1 {--к, *1 )=У'+1 (?«. хг) -

(•' + !) Ие1 2

сг:

V СОБ

Заметим, что величина у (0, х) и представляет собой контур искомого профилированного сопла.

В случае граничных условий скольжения при задании начального профиля в критическом сечении считалось, что толщина пограничного слоя равна нулю и все величины совпадают с соответствующими величинами в невязком потоке. Затем на нескольких шагах по л: (обычно через 20 — 40) поток в пограничном слое постепенно затормаживался до выполнения на стенке условий (5). Для граничных условий прилипания давление в критическом сечении считалось постоянным и равным давлению в невязком газе, а профиль скорости и энтальпии вблизи стенки видоизменялся так, чтобы выполнялись условия прилипания. Плотность и функция у находились из уравнения состояния и закона сохранения массы. С целью проверки влияния вида начального профиля на решение уравнений (1) был проведен ряд расчетов, на основании которых можно утверждать, что при Ие**^ Ю3 течение в пограничном слое не зависит от начальных данных. Точность расчета проверялась дроблением сетки. Оказалось, что для достижения точности в три-четыре значащих цифры во всех рассчитываемых величинах

достаточно 50—80 узлов по поперечной координате в пограничном слое и 40 — 60 узлов в невязком течении при величине шага х не более 0,5. При этом величина ДМ/М, где ДМ—изменение числа М вдоль оси сопла в области характеристического ромба, составляет 1—2% и не уменьшается при дальнейшем дроблении сетки, что объясняется неточностью задания исходного контура bwk и частично сглаживанием. Время расчета варианта порядка 20—60 минут на ЭВМ БЭСМ-6. Заметим, что если бы для расчета системы уравнений (1) применялась, как и в работе [1J, явная схема, то время расчета возросло бы (по порядку) в 0(NV2/Ai) раз, где NV — число узлов в пограничном слое, N — в невязком течении.

5. Результаты расчета. Расчеты проведены для случая f =1,4 и Рг = 0,715. Зависимость вязкости от температуры определялась законом Сазерленда. На рис. 2 приведены контуры сопл с числом М, на выходе сопла равным 8. Число Re было выбрано равным 10°, 105, 2985, 1400, а температурный фактор TJT0 — 0,\2; 0,12; 1 и 0,05 соответственно. Благодаря проведенной стандартной нормировке на числе Re'2 поперечной координаты в пограничном слое метод

эффективно работает во всем рассмотренном диапазоне чисел Re без дробления сетки при увеличении Не*. Более того, как правило, при больших числах Йе* время расчета уменьшалось вследствие некоторого уменьшения числа итераций. При значениях температурного фактора ГГОТ0^0,015 влияние скольжения потока на рассчитанные величины даже при Не* = 1400 мало (менее 1%), что согласуется с результатами работы [1]. Поэтому эти варианты рассчитывались для граничных условий прилипания. Но при Тт/Т0 — \ и Ие* = 2985 скольжение потока оказывает существенное влияние на характеристики сопла (до 10 %). Например, скорость скольжения

на выходе соила достигает значения 0,24 (рис. 3). Случай Не*=1400 и Т,М1Т0 = 0,05 является фактически предельным для рассмативае-мого сопла. Влияние вязкости настолько велико, что поперечный размер невязкого ядра профилированного сопла составляет менее

25% радиуса выходного сечения заданного сопла (профиль скорости приведен также на рис. 3).

На рис. 4 проведено сравнение результатов расчета для контура сопла при М=8, Ие*=2985, Т^Т0= 1 с контуром, полученным

из расчета по методу работы [1J и любезно предоставленным В. В. Михайловым (пунктирная линия). Видно, что совпадение удовлетворительное и расхождение в результатах не превосходит 4%, что согласуется с точностью расчетов по методу работы [1].

В заключение автор считает приятным долгом выразить признательность и благодарность за весьма полезное и ценное обсуждение В. В. Михайлову, В. Н. Гусеву, В. П. Провоторову и П. Е. Бабикову.

ЛИТЕРАТУРА

1. Михаилов В. В. Метод расчета сверхзвуковых сопл с учетом влияния вязкости. „Изв. АН СССР, МЖГ“, 1969, № 1.

2. Ермак Ю. Н. Влияние давления на срезе гиперзвукового сопла на течение внутри сопла. „Ученые записки ЦАГИ“,'т. XIII,

№ 5, 1977.

3. L i n R. Viscous flow over a cone at moderate incidence—1: hypersonic tip region. .Computers and Fluids“, 1973, N 1.

4. Ковеня В. М., Черный С. Г. Решение упрощенных уравнений вязкого газа маршевым методом. Сб. „Численные методы механики сплошной среды“, т. 10, № 1, 1979.

5. Я н е н к о H. Н. Метод дробных шагов решения многомерных задач математической физики. Новосибирск, .Наука“, 1967.

6. Галкин В. С. Об эффектах скольжения при обтекании тел гиперзвуковым слаборазреженным потоком. „Инженерный журнал“, т. 3, вып. 1, 1963.

7. Ко g am М. N. Molecular gas dynamics. In „Annual Review of Fluid Mechanics“, vol. 5, Palo Alto, California, USA, 1973.

8. Рейнольдс М. Измерение профиля скорости в кнудсенов-ском слое для задачи Крамерса. В кн. „Динамика разреженных газов“. М., „Мир“, 1976.

9. Г а л к и н В. С., Н и к о л а е в В. С. О моделировании вязких гиперзвуковых течений в аэродинамических трубах. „Ученые записки ЦАГИ“, т. 1, № 4, 1970.

10. С а м а р с к и й А. А., Николаев В. С. Методы решения сеточных уравнений. М., „Наука“, 1978.

Рукопись поступила 12/1 1981 г.