CC BY
51
УЧЕНЫЕ ЗАПИСКИ Ц А Г И Т о м XII 19 8 1
№ 3
ЧИСЛЕННЫЙ РАСЧЕТ ВНУТРЕННИХ ТЕЧЕНИЙ ВЯЗКОГО ГАЗА С ИСПОЛЬЗОВАНИЕМ УРАВНЕНИЙ
НАВЬЕ—СТОКСА
Т. Д. Асланов, А. П. Биркин, В. В. [Ценников
Приводятся результаты расчетов течений вязкого сжимаемого газа в плоском канале постоянного сечения и в сопле Лаваля при числах Яе ~ 100 с использованием уравнений Навье—Стокса. В качестве метода численного решения использовался метод установления с применением схемы расщепления исходных уравнений по физическим процессам и пространственным переменным [1]. Расчеты проведены для различных перепадов давления, в том числе для перепадов давления, когда на выходе реализуется сверхзвуковой режим течения.
При расчете внутренних течений сжимаемого газа при малых и умеренных числах Ие возникает необходимость использования полных уравнений Навье—Стокса. Это обусловлено тем, что в ряде случаев рассматриваемое течение не может быть описано в рамках модели разделения потока на невязкое ядро и пограничный слой или модели „узкого канала“ (это, например, случаи смыкания пограничного слоя, нерасчетного режима течения с ударными волнами, которые могут привести к отрыву течения на стенке).
Использование уравнений Навье—Стокса хотя и усложняет численный расчет и требует увеличения затрат машинного времени, снимает допущения приближенных постановок и позволяет в принципе по единому алгоритму рассчитывать дозвуковые, трансзвуковые и сверхзвуковые области течения, течения при наличии ударных волн и возникающие при этом отрывные течения на стенках сопл.
В настоящей работе приводятся результаты расчетов течения вязкого газа в плоских каналах и соплах, полученные с использованием консервативной разностной схемы применительно к полным уравнениям Навье—Стокса [1].
1. Будем рассматривать задачу об истечении газа из большой емкости через насадок малого поперечного сечения в среду с заданным давлением (рис. 1). Насадок представляет собой либо плоский канал постоянного сечения, либо плоское сопло Лаваля, течение в котором предполагается стационарным.
Для описания течения газа в рассматриваемом насадке используются нестационарные уравнения Навье—Стокса, решение которых будем находить методом установления. Предельное решение
_________|_
Рис. 1
этих уравнений при £-* оо (£— временной параметр) и стационарных граничных условиях будем интерпретировать как стационарное решение поставленной задачи.
Исходная система уравнений Навье—Стокса, записанная в дивергентной форме, имеет вид:
5 (Р“) + £ и)=°>
дх
д (р“) + (р“2) + -57 Ии)
дх
-|г(р *0 + ^(ри*0 + -|гИ2) =
а ,
ж!ры
Ц2 +
")]
+
й2 + и;
где
и=•-—
Ие
дх' ^ =
р = (*— 1) рв, и = Iх (е),
д / ди ду \ д*
. 3 д* \Г дх)' ду У ду] 1
ду
3 дх ду
(!)
У = -
Ке
+
Ие
— —(V —) +— (V—) +-*-{г*а)-±'±(г ^у
3 ду \ ду) дх \ дх) дх \ ду / 3 ду \ дх)
Е—____-__Г— (\
Ие Рг [дх \
('
д ( де
*Тх
+
4 д ( ди \ . д ! ди ,
3 дх) + ду V ду )~}~ 'ду[1>‘и ~дх)
4 д / ду \ . д ( ду \ д
~3~ д^ У дх) дх У ду/ 3 ду
. д ( де
ду V ду
\ , д [ ду\ 2 д / ду
1-1-------1\хи — I — --------------и.к —
3 дх \ ду,
. ~ , ди\ 2 д ( да
~ дх ду) 3~ ду дх
Безразмерные величины связаны с размерными следующими соотношениями:
еа.
УцоХ
р°
Р = -Хгг- , Р1
У
Ут1
и -
V = ■
р1 К)2
О 9
Н
Ие:
Р1 и1 Ууи\ Iх!
где Ь° — время; х°, у° — расстояния, измеряемые вдоль оси и по нормали к плоскости симметрии канала; и°, Vе — продольная и. нормальная составляющие скорости; е°, р°, р°, р.0 — соответственно внутренняя энергия, плотность, давление, коэффициент динамической вязкости; и\, р°, |х°, у°ш1—-значения соответствующих величин в начальном сечении канала (лг, = 0); значение и\, кроме того, отвечает начальному моменту времени (£ = 0); индекс ча соответствует стенке канала, х — показатель адиабаты, Рг — число Прандтля.
Решение рассматриваемой задачи будем находить для стационарных граничных условий следующего вида:
и-.
и = 0, V — 0.
V = 0, е = е{, р = р1 при х — хи
ди
р = р1 при х = хи у=1,
де
ду
ду
ду
д2 и
дх2
и —0, г> = 0, е = ет при х>*„ У=Уп,
|^г=0 при л = л2, 0 <у<уь Р (хл, У*) = Р2,
■* £=°.
(2).
где х! = 0, = длина канала, р 1=1/(и М?), ех = 1/(*(х — 1)М?) ,
М1 — число М на входе в канал в начальный момент времени.
Выписанные граничные условия соответствуют заданию однородного потока газа во входном сечении, условию симметрии течения относительно оси х, заданной температуре стенки, прилипанию газа на стенке и заданному давлению газа р2 на срезе выходного сечения канала.
Граничные условия для и, V, в при х = х2 при получении численного решения отвечают линейной экстраполяции этих величин в продольном направлении (по внутренним узлам сетки) на границу расчетной области [2—4]. Отметим, что в работе [5] были предложены граничные условия, следующие из асимптотики течения вязкого газа в следе-за решеткой плоских каналов. Однако
можно надеяться, что использование „мягких11 граничных условий экстраполяционного типа при наличии сверхзвукового ядра течения в выходном сечении канала не приведет к заметному искажению действительной картины течения. В то же время это позволяет сократить расчетную область и уменьшить объем вычислений.
Для изложенной выше постановки задачи значение скорости на входе (величина а) и, тем самым, расход газа через поперечное сечение канала будут искомыми. При расчетах на каждом временном слое при этом использовалось условие ди/дх = 0 при х = 0, у = 0, что с учетом выписанных выше граничных условий при х = 0 соответствует однородному потоку газа на плоскости симметрии в окрестности начального сечения канала.
В общем случае криволинейного канала для численного решения системы уравнений (1) с граничными условиями (2) преобразуем область течения с криволинейной границей в прямоугольную область. Для этого введем преобразование координат по формулам
В координатах •£, т) область течения представляет собой единичную полосу, а уравнения (1) запишутся в виде
(Ли р) 4~ (Ли РИ) + (Уи> Р®) — О»
— — (УтР) + -Щ (У и, ЧР) + жа>.и+-|-(у. 1»»1')+^(у.р»о)--4| + 1/,> | (3)
±[м(е+^)]+±[ми{е + ^]] +
+% [у’>(е+71Н]+1 V)+~ (у. т - ви
где
V — Уш Щ 11} =-----------------------------------
Уги
, т 1(4 д / ди\ 4 д [ . ди\ 4 0 ( , ди \ ,
1кё {т ~дг ^ йг) ~ Т а?\Ут^ з~ ~дч ® ^ Ж)
д I 1 Г 4
да\ д / -дх>\ 2 д ( дь
дп ) + д!) 3
1 д (Уш *Л\
3 дг, ^ <?-Г] /1 ’
штрих означает дифференцирование по продольной координате; операторы 1/], Е1 имеют вид, аналогичный их.
2. При численном решении систем (1) или (3) с соответствующими граничными условиями используется схема расщепления уравнений по физическим процессам и пространственным перемен-
ным [1], включающая при расчете каждого временного шага (я+1) следующие четыре этапа.
/ этап-.
1
1
дЦ
Утю
]п 0 у
П + ~Г
ди . ди
+ ип
дч
3 р" Ие дЧ
IX»
ди
1 , 1 _ п + — —П + —
ду . ду
+ ип-
1
дt
дЧ
РлИе дгс,
Ж \ Г
____п 4- -
де
1
ип
гп+ 4
де
дЧ
1
Р Уи>
[ЬПе-
дЧ
<й"+^ д! )
рл Ие Рг дЧ
? Ут
Р" Уто
—П + -Т
де
51
и2 + У2\п Тп
1*-ип Ы — &1%
(4)
Дифференциальные операторы Ьр, Ьи, Ьу и Ье определяются по данным на предыдущем временном слое п и имеют вид:
и = -щ- (уя ри) + (уи рм),
1* = -^ (У* Р“2) + (У» Р®«) + (У™ Р) — ^ (У« V) —
II этап\
уш ри (е +
д_
д;
ц2 1)2 V
■
+
дт)
Уш рву (? -+
и2 + и2\'
}
+
+ -Ж (У* иР) + -щ(Ут»™Р) — Е 1 •
По аналогии с дифференциальными уравнениями расщепляются и граничные условия задачи.
При установлении ((-» ос) данная схема расщепления уравнений в отличие от схемы расщепления [6] является консервативной.
п+ 2
Нахождение вектора / 4 {р, и, V, е) на каждом этапе (г==1,
2, 3, 4) в каждой точке расчетной области производится методом конечных разностей. Область интегрирования покрывается сеткой с узлами
%1 = 1 А?, ^ = /Дт) (I = О, I, 2, ... , Ь;
/ = 0, 1, 2, . . Л/),
где Д|, Дт]—шаги сетки в продольном и поперечном направлениях;
Ь, N—граничные номера узлов сетки.
Первые производные относительно величин /с чертой, входящие в выражения (4) —(7), аппроксимируем с первым порядком точности с учетом знака компонентов скорости; дифференциальные операторы /,р, Ьи, Ье аппроксимируем с первым или вторым порядком точности по несимметричной схеме [1]. При установлении точность получаемого решения определяется точностью аппроксимации операторов Ьр, Ьи, Ье.
Получаемые на каждом этапе системы алгебраических уравнений для определения в узлах сетки вектор-функции / решаются методом скалярной прогонки. На IV этапе при этом целесообразно
дип + 1 . ' п
использовать соотношение —^-------1- у-ш у\ —^— = 0, вытекающее из
второго и третьего уравнений системы (7).
4—„Ученые записки ЦАГИ“ № 3
49
После определения в каждой точке расчетной области вектор-функции /п+1 значение вектор-функции /п +1 определяется по формуле
/п + 1 = /" 4- /п+1.
Задавая тем или иным способом начальные данные при t = 0, продолжаем расчет до установления течения. Критерием уста-
д/
новления считается выполнение условия шах
дt
< е, где з — малая величина.
Прежде чем переходить к результатам конкретных численных расчетов, отметим, что с использованием описанного метода был проведен расчет тестовой задачи, соответствующей автомодельному течению газа в плоском расширяющемся канале с прямолинейными стенками с числом М=1 и числом 1?е=25 [7]. При этом имело место хорошее соответствие результатов численного расчета и точного решения.
3. Описанный выше метод был использован для получения решения поставленной в п. 1 задачи о течении вязкого газа в насадке, представляющем собой канал постоянного сечения.
Газ принимался двухатомным (х—1, 4), Рг = 0,71. Зависимость коэффициента вязкости от температуры полагалась степенной
Тп, л = 0,76).
Значение плотности на стенке рассчитывалось не из нестационарного уравнения неразрывности, а по давлению на стенке, найденному экстраполяцией по внутренним узлам сетки со вторым порядком точности.
Начальные данные при t=■0 задавались таким образом, что во всей области течения, кроме стенки, продольная составляющая скорости и температура полагались постоянными; поперечная составляющая скорости всюду полагалась равной нулю; р=р(х) — линейная функция. Температура стенки при расчетах принималась равной температуре торможения потока газа при £ = 0.
Были выполнены численные расчеты течений в канале постоянного поперечного сечения при 1 = 8, Ие —50, М: = 0,4236 и четырех различных значениях параметра противодавления к = р2/р1 = 0,7; 0,55; 0,4 и 0,3. Расчеты всех вариантов проводились с первым порядком точности с параметрами сетки /, = 50, А/ = 17 и шагом по времени, изменяющимся от 0,01 до 0,05. Число временных шагов при этом не превышало 400 (е = 10~3 -ь 10-4).
Отметим, что при увеличении А/ —числа узлов сетки в поперечном направлении — следует ожидать уменьшения скорости газа на входе в канал (расхода газа), поскольку при заданном перепаде давления потери на трение будут больше.
При проведении этих расчетов применялось преобразование координат, приводящее к сгущению сетки в пристеночной области. Так, вместо у] использовалась переменная С, определяемая формулой
.. . ( £==1_ 1п|1 -|-а(1 —•<;)]
1п (1 а) ’
где а = 10.
На рис. 1 приведены полученные из расчетов распределения по х продольной составляющей скорости и, отношения давлений р/ри и числа М, отвечающих плоскости симметрии течения. Из приведенных на рисунке данных видно, что при £<0,7 число М
в ядре потока на срезе канала больше единицы, т. е. скорость газа в плоскости симметрии на срезе канала превышает звуковую. При уменьшении к с 0,4 до 0,3 расход газа через поперечное сечение канала практически не меняется, что свидетельствует о близости режима течения при к = 0,3 к режиму „запирания11. Последнее проявляется также в том, что отличие в поведении всех величин по х для указанных случаев заметным образом сказывается лишь в окрестности выходного сечения канала.
Вблизи входного сечения наблюдается немонотонность в поведении кривых р/р1(х), что обусловлено, по-видимому, заданием однородного потока газа на входе в канал и резкой перестройкой течения. При этом на всей длине, за исключением окрестностей входного и выходного сечений, давление практически постоянно по сечению канала.
В качестве иллюстрации картины течения на рис. 2 приведены профили продольной составляющей скорости и(п) в сечениях канала по х, отмеченных значением I (х — 1Ах), отвечающие случаю Л = 0,3.
Поперечная скорость вблизи входного сечения направлена от стенки к плоскости симметрии канала, а при подходе к выходному сечению меняет знак.
Кроме рассмотренных выше случаев, проведен расчет течения газа в канале с длиной /=16 при Ие = 100, к —0,3. Начальные данные, температурные условия на стенке и параметры сетки задавались аналогичным образом.
Представление результатов расчета этого случая и случая, •отвечающего 1 = 8 при Р!е = 50, к — 0,3, по переменной подобия х/Ке [8] показало, что они хорошо коррелируются одной кривой.
Рис. 3
4. Проведены расчеты течения вязкого газа в сопле Лаваля,, контур которого задавался следующим образом:
длина сопла 1 = 4. Площади входного и выходного сечений рассматриваемого сопла одинаковы и равны утроенной площади минимального сечения. Теплофизические постоянные газа, начальные данные при 1 = 0 и температурные условия на стенке принимались такими же, как при расчетах течения в канале постоянного сечения. Параметры сетки I — 50, Ы=\7, сетка в переменных Е, -ц полагалась равномерной.
Расчеты соответствуют значениям Ие = 50 и 250, £ = 0,3 и 0,1. При £ = 0,3 численное интегрирование исходных уравнений проводилось с первым порядком точности, при £ = 0,1 — как с первым,, так и со вторым порядком точности.
На рис. 3 приведены полученные из расчетов распределения по х величин и, р\рх и М в плоскости симметрии сопла. Кривые, отмеченные крестиками, получены с использованием схемы второго порядка точности.
Приведенные на этом рисунке данные указывают на то, что при £ = 0,1 для обоих значений числа Ие в сопле реализуется
сверхзвуковой режим течения с монотонным изменением газодинамических величин по длине сопла (на причину немонотонного характера изменения величины р!р^ по л в окрестности входа ука-
Рис. 5
зывалось в разделе 3). Число М в центре выходного сечения сопла равно 2,47 при Re —250 и 1,93 при Re = 50. Для этих случаев значения числа Rе* == , определенные по парамет-
рам потока в горле сопла, равны соответственно—128 и ~23,5. Сечения, где число М на оси равно единице, сдвинуты вправо от минимального сечения сопла.
Сравнение результатов расчета, полученных с первым и. вторым порядком точности, показывает, что наибольшее их отличие имеет место в окрестности выходного сечения и не превышает 2%. Отметим, что расход газа по х выдерживался постоянным с точностью до 1;5% при расчетах с первым порядком точности и с точностью до долей процента при расчетах со вторым порядком точности.
Для случая Re = 250, & = 0,1 на рис. 4 представлены йрофили u(ri), где n = y/yw, в различных сечениях сопла. Изменение давления в поперечном сечении сопла достигает 30—70%. Указанное обстоятельство необходимо иметь в виду при использовании приближенной постановки задачи о течении вязкого газа в каналах, когда величина давления принимается зависящей только от продольной координаты [8]. v
Расчеты течения в сопле при k = 0,3 для обоих значений числа Re (см. рис. 3) свидетельствуют о немонотонном изменении газодинамических величин по длине за горлом сопла. При этом имеются области сверхзвукового течения, вниз по течению от которых давление возрастает и поток тормозится. Последнее приводит при Re = 250 к отрыву течения от стенки в области выходного сечения сопла.
На рис. 5 приведены профили и(т[) для случая Re = 250, k = 0,3, в котором реализуется отрывное течение.
В заключение отметим, что в рассмотренных случаях течения газа в сопле Лаваля значение скорости на входе в канал при заданном числе Re практически оставалось постоянным; при изменении параметра противодавления k.
ЛИТЕРАТУРА
3. Ковеня В. М., Яненко Н. Н. Разностная схема на подвижных сетках для решения уравнений вязкого газа. „Журн. вычисл. математ. и математ. физики*, т. 19, № 1, 1979.
2. Cleinfe М. С. Computation of two-dimensional, viscous nozzle flow. „А1АА J.“, vol. 14, N 3. 1976.
3. К у з н e ц о в а Л. В., Павлов Б. М. О расчете течения вязкого газа в соплах Лаваля. Сб. „бычислительные методы и программирование”, т. XXVII, ВЦ МГУ, 1977.
4. Пирогов В. Б., СевериновЛ. И. О расчете внутренних течений вязкого сжимаемого теплопроводного газа. ВИНИТИ,
№ 3359—77 Деп.
5. Быркин А. П., Щ е н н и к о в В. В. Расчет дозвуковых течений вязкого газа в плоских каналах и в следах. „Журн. вычисл. матем. и матем, физ.“, т, 19, № 1, 1979.
6. Березин Ю. А., Ковеня В. М., Яненко Н. Н. Об одной неявной схеме расчета течения вязкого теплопроводного газа. Сб. „Численные методы механики сплошной среды", т. 3, № 4. Новосибирск, „Наука", 1972.
7. Быркин А. П. Об одном точном решении уравнений Навье— Стокса для сжимаемого газа. ПММ, т. 33, № 1, 1969.
8. Быркин А. П„ Межиров И. И. О расчете течения вязкого газа в канале. „Изв. АН СССР, МЖГ“, 1967, № 6.
Рукопись поступила 25jXII 1979 г.
