УЧЕНЫЕ ЗАПИСКИ ЦАГИ То м XXIV 199 3
№ 1
УДК 532.525.011.55
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ТЕЧЕНИЯ ГАЗА В СОПЛАХ ЛАВАЛЯ И ТЕЧЕНИЯ ТОРМОЖЕНИЯ В КАНАЛАХ НА ОСНОВЕ ПОЛНЫХ УРАВНЕНИЙ НАВЬЕ—СТОКСА
А. П. Быркин., В. Н. Клюев, А. Ю. Сосунов, А. И. Толстых
В рамках полных уравнений Навье—Стокса выполнено численное моделирование ламинарного течения газа в плоском и двух осесимметричных соплах Лаваля. Осесимметричные сопла спроектированы на получение в рабочей части вакуумной аэродинамической трубы однородного потока воздуха с числом М=6. Проведенное сопоставление расчетных и имеющихся экспериментальных распределений чисел М потока вдоль оси , и в поперечных сечениях сопл показало удовлетворительное их согласие. Апробированный в данных численных экспериментах алгоритм [1| может быть использован для диагностики характеристик ламинарного течения газа в соплах при малых и умеренных числах Ие.
На основе уравнений Навье—Стокса проведено моделирование торможения газа в плоском канале постоянного сечения и в сужающемся канале (Ие= 103-ь 105) при сверхзвуковой скорости на входе (М1>1).
В результате расчетов получена, в частности, картина течения, реализующегося в начальной части канала с наличием скачков уплотнения (исходящих с нижней и верхней передних кромок) и переходящего вниз по потоку в течение с затухающей периодической структурой.
В работе [1] представлен разностный алгоритм высокой точности, основанный на компактных аппроксимациях третьего порядка [2], использовавшийся для расчётов ламинарных и турбулентных течений вязкого газа в каналах, сверхзвуковых и гиперзвуковых соплах аэродинамических труб. В последнем случае применялся метод глобальных итераций для параболизованных уравнений Навье—Стокса.
Основной целью данной работы являются дальнейшие численные исследования внутренних течений, характерных для ряда прикладных задач при наличии факторов, в той или иной степени усложняющих структуру потока. К таким факторам можно отнести, в частности, резкие изменения формы контура, возникновение ударных волн, взаимодействие с пограничным слоем и т. д. В этих случаях целесообразно использовать полные уравнения Навье—Стокса.
Ниже представлены результаты расчетов течения вязкого газа в соплах Лаваля, в том числе предназначенных для вакуумной аэродинамической трубы, и в каналах постоянного и переменного сечения по длине при сверхзвуковой скорости на входе.
Постановка задачи. Метод расчета.
Рассматривается стационарная задача о внутренних ламинарных течениях вязкого газа (канал постоянного или переменного сечения, сопло Лаваля). Для описания течения используются полные нестационарные уравнения Навье—Стокса, решения которых находятся методом установления.
Перейдем от исходной системы координат х, у (декартовой или цилиндрической, см. рис. 1) к системе координат |, Т), где % = Х, Г) = :=У/У^{х), Уи>{х)—контур канала. Сохраняя исходные составляющие скоростей и и V вдоль осей х и у, систему уравнений Навье—Стокса запишем в виде
Ш. Л- 4 <5* дт.
н + 1кг’
(1)
где время, и, /% О, Н, Т — четырехкомпонентные векторы
Р ри
—> ри ( р«2 + р
У=У9УЧ рг> /Л , «2 4- г»2 \ р(х+ 2 -) ! рт >»(*+“’Г)
■У*У
рт
у'ти
рЯУИ - — у\р У XV
Ут
Р_
и2 + V*}
/, , н2 + 1>а\
рда / А Н---------------2—)
0 Л
0 , т= т2
р 9 ш ть
0 Т,
Плоскому случаю соответствует значение параметра v = 0, осесимметричному V = 1; Утю (х) =с1ук(х)/с1х, —у^ци)/уу,\ Т{ = 0, выра-
жения для компонент Т2, Т3, Ть, обусловленных вязкими членами уравнений, из-за громоздкости не выписываем, они приведены, например, в [1].
В уравнениях (1) все величины безразмерные и связаны с размерными (отмеченными индексом °) соотношениями:
пг° II 11 _ ы° ) ^
V0 ^ ’ , л° П- ^0), , р° р —-V, Р1 ■ !А° 1*1
Р! (Щ*
, Ие,
р, р, А, ^ — давление, плотность, энтальпия, динамическая вязкость соответственно; индекс 1 отвечает начальному сечению канала (х = 0), индекс VI) — стенке канала. В качестве характерной скорости задачи №° в случае течения в сопле Лаваля принимается скорость звука газа на входе а\, в случае течения в каналах при сверхзвуковой скорости на входе — значение этой скорости ы,.
Предполагается, что газ является совершенным
■*. -— 1 ,
Р= % Р »
{к — показатель адиабаты) и вязкость ц, является функцией температуры (энтальпии).
Решение рассматриваемой задачи находится для стационарных граничных условий следующего вида:
и = ии V = О, Л = р — рх при \ — 0, 0 < -ц < 1,
^ = ПРИ *>0, *1 = 0,
и — у — 0, Н = И.т(х) при $ > 0, •/]= I,
д^и д2г> д2Н п ее
йр- ” дф ~ йр ПРИ “ ^к’
Р (^Ю 1 ) = Рку
где $к = хк — безразмерная длина канала.
Выписанные граничные условия соответствуют заданию однородного потока газа во входном сечении канала (сопла), условию симметрии течения относительно оси х, заданной температуре стенки, прилипанию газа на стенке и заданному давлению газа рк на срезе канала. В качестве граничных условий при | = |к используются соотношения, отвечающие линейной экстраполяции искомых величин по внутренним узлам сетки на границу расчетной области (в пристеночной области, где скорость газа является дозвуковой, для нахождения плотности привлекается условие заданного изменения статического давления р(у)).
Отметим, что для изложенной постановки задачи в случае сопла Лаваля значение стрости на входе «1 и, следовательно, расход газа через сопло будут искомыми. В случае рассмотрения течения в каналах при «1>й1 (сверхзвуковое значение продольной составляющей скорости на входе) все газодинамические величины р, и, V, И на входе должны быть заданы. Для получения численного решения, кроме того, необходимо задать поля газодинамических величин в начальный момент времени ^ = 0.
В качестве метода численного решения использован численный алгоритм, описанный в работе [1] и основанный на применении компактной схемы третьего порядка аппроксимации конвективных членов (см. [2—3]).
Результаты расчетов
1. На основе алгоритма [1] получено численное решение прямой задачи о течении вязкого газа в плоском (у = 0) сопле Лаваля (сопло 1), контур ую(х) которого представлен на рис. 1. Особенностью данного сопла является наличие протяженного цилиндрического участка в трансзвуковой его части, при этом максимальный угол наклона контура в дозвуковой части составляет 60°, в сверхзвуковой части ~30°. Приведенный на рис. 1 контур сверхзвуковой части сопла характерен для реактивных сопл и сопл химических лазеров [4, 5]. При расчете газ полагался совершенным с показателем х=1,55, число Рг = 0,7, зависимость вязкости от температуры принималась степенной (|х~Тп, п =
п п На\Ут*
= 0,Ь4), число Ке1=---------, определенное по параметрам газа на вхо-
1*1
де и полувысоте минимального сечения равнялось 103. В качестве начальных данных при / = 0 для полей газодинамических величин р, и, V, И. во всей расчетной области, кроме стенки, использовались параметры одномерного течения невязкого газа (на стенке ставились условия прилипания и заданной температуры поверхности; стенка предполагалась охлаждаемой /1Ш//11 = 0,43).
Расчеты проводились на неравномерной по \ и г) сетке с числом узлов 100 x 50, для чего вводились преобразования продольной (с целью сгущения сетки в окрестности критического сечения) и поперечной (для сгущения сетки в пристеночной области) координат. При получении решения шаг по времени т изменялся от 0,001 до 0,5 (максимальное число Куранта составляло К~30), число временных шагов для достижения установления не превышало 300.
На рис. 1 приведены результаты расчета течения в рассматриваемом сопле при значении параметра противодавления на срезе сопла Рк/Р\ = 0,0053. Расчетные данные представлены в виде логарифмических распределений по х статического давления на оси — 1пр0(*), давления на стенке сопла — 1п р-ш{х) и числа М на оси — М0(х). Представленные на рисунке данные показывают, что на оси выходного сечения сопла число М0«5,74. Выше отмечалось, что для рассматриваемой ?адачи о течении газа в сопле Лаваля при заданном перепаде давления между входным и выходным сечениями реализуемый расход газа через сопло является искомой величиной. Заметим, что расход газа по л: выдерживался постоянным с точностью ~0,3%. Число Ие*, определенное по параметрам потока на оси и полувысоте сопла в критическом сечении (М0=1), равно 661.
Сопла1 Сечен11е х = Сечение х=1Б,5 Сечены х=91,В
О 0,5 1,0 1,5 2,0 0 1 Ч ВО 1 ч В
М,и,к,р103
Ш
м
5
Рис. 1
В конце сверхзвуковой части исследуемого сопла давление на стенке превышает давление на оси в несколько раз. Последнее видно на рис. 1, где представлены поперечные профили величин и, h, р, М.
2. Проведено численное моделирование течения газа в двух осесимметричных (v= 1) соплах (сопла 2 и 3), контуры yw(x:) которых рассчитаны с учетом влияния вязкости на получение в рабочей части вакуумной аэродинамической трубы однородного потока воздуха (х=1,4) с числом М = 6 [6] и приведены соответственно на рис. 2 и 3. Они были получены следующим образом. Сначала определялся контур сопла у(х) для невязкого течения газа в результате решения обратной задачи. (Исходный невязкий контур у(х) один и тот же у сопла 2 и 3, он показан на рисунках штриховой линией.) Затем учитывалось влияние вязкости путем добавления к радиусам изоэнтропического контура у толщины вытеснения 6* ламинарного пограничного слоя. Радиус критического сечения обоих сопл равен у* =2,5 мм, параметры торможения газа ро = 5000 Па, Го = 450 К. Расчеты контуров yw{x) выполнены для двух значений температуры стенки сопла (в обоих случаях постоянной по его длине) Tw — 290 К и 80 К, что соответствует охлаждению стенки сопла водой (сопло 2) и жидким азотом (сопло 3). Скольжение и скачок температуры газа на стенке не учитывались (число Кп, рассчитанное по радиусу сопла и условиям на стенке, не превышает 0,03). Число ReL, вычисленное по длине сверхзвуковой части сопла и статическим параметрам газа в его выходном сечении, равно 9300. Столь малое значение Re^ подтверждает, по данным эксперимента [61, наличие на стенке ламинарного пограничного слоя. Из данных, представленных на рис. 2 и 3, следует, что поправка на влияние вязкости соизмерима с радиусом изоэнтропического контура. Радиус выхода сопла 2 равен 41,2 мм, радиус выхода сопла 3 — 28,7 мм, длина сверхзвуковой части обоих сопл 178 мм.
Проведенные в работе [6] испытания указанных сопл в вакуумной аэродинамической трубе показали, что в них реализуются течения, близкие к расчетным. Несмотря на малые размеры невязкого ядра потока в данных соплах, числа М в нем в пределах выходного характеристического ромба практически постоянны и мало отличаются от расчетного значения М = 6.
В настоящей работе при численном моделировании течения воздуха в рассматриваемых соплах, как и при их проектировании, принималось, что х=1,4, число Рг = 0,71, зависимость вязкости от температуры полагалась степенной (л = 0,76). В качестве характерного линейного размера принят радиус минимального сечения соплаyw* =2,5 мм. Число Reb вычисленное по статическим параметрам газа, скорости звука на входе и радиусу yw% для обоих сопл равно 1650. Формирование начальных полей величин р, и, v, h при t = 0 и разностной сетки осуществлялось по аналогии с рассмотренным выше случаем плоского сопла.
На рис. 2 и 3 представлены полученные при расчете продольные распределения давления на оси 1пр0(х) и стенке сопла Inpw(x), числа М на оси — М0(л:); для двух сечений х=180 мм и х = 208 мм (срез сопла) приведены поперечные профили продольной составляющей скорости и, энтальпии h, давления р и числа М.
В случае сопла 2 расчеты выполнены для значения противодавления рк/рі = 0,0006, которое примерно соответствует статическому давлению газа на срезе сопла, отвечающему числу М = 6, а также для противодавлений Рк/Рі = 0,0003 и 0,0012 (в два раза меньших и больших номинального). Несмотря на то, что распределения pw(x) для разных
«о
со
значений рк/р\ отличаются между собой в промежутке по х от 150 до 200 мм (см. рис. 2), распределения по х числа М на оси практически неразличимы. Влияние передачи возмущений вверх по потоку при этом, несмотря на толстый пограничный слой, проявилось в основном на характеристиках течения в пограничном слое. Последнее согласуется с результатами расчета течения в данном сопле с использованием упрощенных уравнений Навье—Стокса [7]. Число Яе* определенное по параметрам потока на оси и радиусу сопла в минимальном сечении, составляет 1060.
В случае сопла 3 расчеты выполнены для значения рк/р1=0,0006, Ие* =1077. На рис. 2 и 3 данные настоящих расчетов в виде распределений чисел М по оси — М0(х) и в поперечном сечении М(у) сопоставляются с данными расчета [7] и эксперимента [6]. Отмечая близкое совпадение между расчетными кривыми (в случае сопла 2), укажем на удовлетворительное согласование результатов расчетов с экспериментом. Отметим, однако, что размер расчетного невязкого ядра потока несколько меньше, чем это получено в эксперименте. При расчетах в невязком ядре течения с большой степенью точности выдерживалось постоянство энтропийной функции Б = р/р;*.
Апробированный в описанных численных экспериментах алгоритм [1] может быть использован при диагностике характеристик течения вязкого газа в соплах Лаваля при числах Не,= 103-М05.
3. Получено численное решение задачи о торможении вязкого газа в плоском ^ = 0) канале постоянного сечения при сверхзвуковой скорости на входе, которая в данном случае являлась масштабом. Исходный вариант расчета соответствует течению в канале с безразмерной длиной хк= 10 при значении числа М.1 = 2, числа 1^1 = 1000 и противодавления рк/р1 = 2. Теплофизические характеристики газа принимались такими же, как и при расчете течений в осесимметричных соплах. Температура стенки полагалась равной температуре торможения газа на входе.
Начальные данные при ^ = 0 задавались таким образом, что во всей области течения, кроме условного пограничного слоя у стенки, все газодинамические характеристики задавались постоянными, причем г>=0 Как и выше, расчетная сетка по I и г] была переменной и содержала то же число узлов (100X50). Для исходного варианта расчета на рис. 4 представлены распределения по л: статического давления на оси (г) = 0)—ро{х), при г)»0,5 — р${х) и стенке — Рю(х), числа М на оси — М.0(х), а также профили величин и, р, Н, М в поперечных сечениях канала. На рисунке, кроме того, приведена картина линий постоянного давления (изобар), откуда видно, что имеет место сложная структура течения газа в канале. Так, в начальной части канала реализуется течение с косыми скачками уплотнения, исходящими с передних кромок нижней и верхней стенок и обусловленными начальным пограничным слоем. Скачки от кромок, взаимодействуя между собой, отражаются от оси и затем размываются, не доходя до стенок (наличие скачков иллюстрируется также поперечными профилями, представленными в средней части рисунка). Далее вниз по потоку видна затухающая периодическая картина изменения газодинамических величин. Судя по распределению М.0(х), число М на оси канала всюду сверхзвуковое, его значение в выходном сечении вследствие торможения уменьшается до 1,52. Отрывные зоны течения у стенок канала при заданном противодавлении отсутствуют.
СО
сл
со
о>
О
Рис. 5
Отметим резкое увеличение давления на стенке канала в окрестности входного сечения, обусловленное повышением давления за формирующимся косым скачком уплотнения.
Проведены расчеты течения в канале при условиях, несколько отличных от исходного варианта. Так, рассматривался случай, когда Pk/Pi=1,5, и случай pu/pi=l,5 и ,>Ск = 5. Сопоставление результатов расчетов для всех трех случаев показало, что картины течения в начальной части канала идентичны, различия в структуре течения проявляются вниз по потоку за размытыми отраженными скачками уплотнения.
Выявленный характер рассмотренного течения в канале вдали от передней кромки качественно согласуется с результатами эксперимента [8], посвященного исследованию течения газа в псевдоскачке (М4 = = 1,62, Re= 106).
4. Ниже представлены результаты расчета течения вязкого газа в сужающемся плоском канале (отношение выходного сечения к входному равно двум). Контур канала yw{x) показан на рис. 5. При расчете задавались число Mi = 6, Число Re± = 105, значение противодавления Рк/Р\ = 1.
Теплофизические характеристики газа, граничные условия на стенке, параметры расчетной сетки соответствуют рассмотренному выше случаю. На рис. 5 приведены полученные распределения Po(x), ps(x), Pw(x), M0(jc), поперечные профили величин и, h, р, М, линии изобар. При течении в канале рассматриваемой конфигурации с принятыми граничными условиями в окрестности сечения х = 8 образуются протяженные зоны повышенного градиента давления, локализованные между осью и стенкой канала. При этом течение у стенки в указанной области имеет отрывной характер. Наличие возвратного течения иллюстрируется показанным на рисунке поперечным профилем продольной составляющей скорости.
В дальнейшем предполагается проведение расчетов течения газа в каналах при сверхзвуковой скорости на входе для ситуаций, при которых скорость на выходе становится дозвуковой.
ЛИТЕРАТУРА
1. Быркин А. П., Толстых А. И. Компактные схемы третьего и четвертого порядков в задачах о внутренних течениях вязкого и невязкого газа//Журнал вычисл. матем. и матем. физики. — 1988. Т. 28, № 8.
2. Толстых А. И. О неявных схемах повышенной точности для систем уравнений//Журнал вычисл. матем. и матем. физики.— 1981.
Т. 21, № 2.
3. Савельев А. Д., Толстых А. И. Метод расчета течений вязкого сжимаемого газа на основе компактных аппроксимаций третьего порядка//Журнал вычисл. матем. и матем. физики.— 1987. Т. 27, № 11.
4. Пирумов У. Г., Росляков Г. С. Течения газа в соплах.—
М.: Изд-во МГУ, 1978.
5. Hyde J. С., Н os ас k G. A., and Osugi G. S. An Investigation of velocity flowfields in chemical laser nozzles//AIAA Paper. — 1973,
N 73-641.
6. Кудрявцева Л. И., Межиров И. И., Пономарев С. П., Якушева В. Л. Экспериментальное исследование осесимметричных профилированных сверхзвуковых сопл при малых числах Ие//Ученые записки ЦАГИ. — 1973. Т. 4, № 3.
7. Быркин А. П., Тимофеева Т. А., Толстых А. И. Применение компактных схем третьего-четвертого порядка для расчета течения газа в соплах с большими сверхзвуковыми числами М на основе упрощенных уравнений Навье—Сгокса//Ученые записки ЦАГИ.— 1988. Т. 19, №6.
8. Т а ш a k i Т., Т о m i t a Y., Y a m a n e R. A study of pseudo-shock (1-st report, Л-type pseud9-shock)//Bylletiri pf JSME, 1970, vol. 13, N 55.
Рукопись поступила 15/IV 1991 e,
7 —«Ученые каннски» № |