CC BY
102
Вычислительные технологии
Том 1, № 2, 1996
МОДЕЛИРОВАНИЕ ГАЗОДИНАМИЧЕСКИХ ПРОЦЕССОВ В КАМЕРАХ СГОРАНИЯ НА МНОГОПРОЦЕССОРНОЙ МАШИНЕ**
Н. А. Ваганова, О. О. Коврижных, О. Б. ХАйРУЛЛИНА Институт математики и механики УрО РАН Екатеринбург, Россия
Комплекс программ СОКОЛ-М для многопроцессорной вычислительной системы МВС позволяет проводить серии расчетов течений газа при решении квазистацио-нарной задачи течений газа в реальных каналах, а также при различных параметрах входного потока. Параллельный расчет позволяет сократить время расчета задач в 1.6Ж раз при счете N вариантов на 3Ж процессорах с одновременной визуализацией всех вариантов в разных окнах экрана.
При проектировании камер сгорания двигателей, представляющих собой каналы сложных конфигураций, важным моментом является знание параметров и структуры течений газа в этих камерах. Математическое моделирование процессов движения газов в каналах сложных геометрий связано с преодолением многих трудностей: большой диапазон изменения скоростей потока — при учете соплового объема от дозвукового до сверхзвукового; большие числа Рейнольдса; сильно меняющаяся по времени форма поверхности горения; образование в дозвуковой области многих вихревых зон с замкнутыми линиями тока, вызванное взаимодействием встречных потоков. Успех расчетов, позволяющих передать качественные и количественные особенности описанных явлений, связан как с выбором соответствующих надежных математических моделей, так и с созданием эффективных алгоритмов и программ, позволяющих исследовать большое число достаточно разнообразных конструкций.
Задача о математическом моделировании процессов движения газов в камерах весьма сложных геометрий рассматривается в двумерной стационарной осесимметричной постановке. При этом обычно требуется рассчитать серию вариантов течения для различных фиксированных геометрий поверхности горения. С учетом небольших скоростей горения допущение о стационарности течения в определенные моменты времени оправдано и используется в большинстве работ.
Газовый поток в камерах сгорания характеризуется большими числами Рейнольдса, поэтому для расчета течения газа в камере сгорания вводятся некоторые упрощения. Следуя [1]—[3], считаем: поверхности каналов (рис. 1) состоят из участков пористой поверхности (Г1), через которые вдувается газ, твердых стенок (Го) и участков выхода газа (Г2); размеры пограничных слоев, нарастающих вдоль стенок, малы по сравнению с поперечными
*© Н. А. Ваганова, О. О. Коврижных, О. Б. Хайруллина, 1996.
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №94-01-00361-а.
размерами каналов; погранслои не взаимодействуют между собой; вдуваемый газ однороден, но неизоэнергетичен. После таких допущений для численного моделирования газодинамических процессов в каналах можно использовать модель идеального газа, течение которого удовлетворяет уравнениям Эйлера.
Использование модели идеального газа неочевидно. Однако, если на входы канала подается завихренный (изо- или неизоэнергетический) поток газа, то и в рамках модели идеального газа внутренние течения газа могут иметь замкнутые вихревые зоны. В работе [4] для плоского случая приведен класс точных вихревых решений с замкнутыми линиями тока, соответствующий такой ситуации. Проведенные расчеты по различным вязким и невязким моделям среды показали удовлетворительное согласование результатов, в том числе и в вихревых зонах [5].
При расчете течений газа используется комбинированный метод, при этом область во всем газовом тракте разбивается на три подобласти — дозвуковую (x < Xd), трансзвуковую (xt < xd < X < xts < xs) и сверхзвуковую (x > xts) с перекрытием подобластей [6].
Для расчета дозвукового течения на входных участках рассматриваются два типа входных граничных условий, когда задаются:
плотность р0 и скорость Vo в направлении нормали n к границе Гі, поток массы р0V0 в направлении нормали n к границе Г1 и энтальпия H.
На выходе задается расход газа pVn|г2 = v2(s) > 0, на твердых стенках — условие непротекания Vn|Го = 0, на оси симметрии — условия симметрии
Кроме этого, задается давление торможения Р0 или давление Р* в критическом сечении. Граничные условия должны удовлетворять соотношению баланса массы
По мере выгорания топлива геометрии каналов сильно меняются, хотя параметры вдуваемого газа остаются неизменными. При моделировании процессов, протекающих в камере, возникают две задачи: отработка параметров исходных данных и отработка геометрии канала. Поэтому для определения структуры течения газа в реальной камере необходимо рассчитать серию вариантов либо для одной и той же геометрии при разных входных
г
Рис. 1. Схема расчетной области канала.
dv/ = dp = dp = 0
dr dr dr
V2(s), 0,
S Є Г2, s Є Го.
данных, либо для различных геометрий поверхности горения, т. е. для фиксированных моментов времени, при одних и тех же параметрах вдуваемого газа. Расчет первой серии вариантов позволяет правильно скорректировать давление торможения, заданное из эксперимента с некоторой погрешностью. При расчете второй серии вариантов можно установить, на каком этапе процесса горения образуются замкнутые вихревые потоки, которые при натурных испытаниях приводят к разрушению конструкций. Для расчета каждого варианта из-за сложности конфигураций каналов требуется большое количество точек расчетной сетки, а следовательно, и значительное время счета.
Для исследования структуры течения газа в камере сгорания поставленная задача, как ив [7], решается в переменных функция тока ^ — вихрь ш. Исходная система уравнений газовой динамики расщепляется на две подсистемы: гиперболического типа для функции тока ^ (при известных вихревой функции и энтальпии) и эллиптического типа для вихревой функции ш и энтальпии Н (при известной функции тока).
Конфигурации каналов достаточно сложны, и расчет течения в таких областях на прямоугольных сетках вызывает определенные трудности [7]. Для преодоления их в предлагаемом методе использованы оптимальные криволинейные сетки [8]—[10], обладающие хорошими аппроксимационными свойствами, позволившие снять ограничения на класс рассматриваемых конфигураций каналов.
Для численного моделирования уравнения Эйлера записываются в интегральной форме в криволинейной системе координат (91,92) [11]
с
1/ 1 |_ргА V 1 - 3 51 - РТд
А ^ - Аз ^ I ¿,2
д91 д,2
где
ш / ,
- — ¿,1 + — ¿92
г V д<91 д,2
с
Ша
р /дУ2 , дУ2 ,
о —¿91 + -д— “92
2 V д,1 д,2
с
с
22
Н ( дт"“^ + |^“92 д,1 д,2
(1)
(2)
(3)
А = х1 + гь А = х2 + г2, Аз = Ж1Ж2 + Г1Г2,
А = х1г2 — х2г1, xi = дх/д^, г\ = дг/д, (г = 1, 2).
При этом компоненты вектора скорости У1,У2, функция тока <^, вихревая функция ш, энтальпия Н, давление Р и плотность р удовлетворяют соотношениям
У _ 1 (\ У _ 1 (
У = — ртд (д^Х2 — а92х7'= — РТд (д,!г2 — д,2г1
ш
д91
Р = Ро —
Х2
дУ1 дУ2 дУ2
т;—Х1 + т;— г2 —— д,2 д,1 д,2
р д V2 ш
г1
1
А:
Ь(М0,М)
“91—
0
1
/
2 дд2 г дд2
(5)
Ь(М0,М)
(6)
Р = Ро + / гК'92)р(£'92) [У1(£'92)г1&92) — у2&92)х1Й'92)] +
210
31
+ / г(^1,п)р(91,п) [У1 (91,п)г2(91 ,п) — ^2(91, пМ^ъ п)] ¿П,
— произвольная подобласть, ограниченная контуром С, из заданной расчетной области П; Ь — кривая, соединяющая точку М0 с точкой М из П, 7 — показатель адиабаты.
Для расчета установившихся дозвуковых течений газа использовался конечно-разностный итерационный последовательный алгоритм для поля скоростей. При дискретизации
(1)-(6) использовались специальная аппроксимация, учитывающая направление потока и особенности криволинейности сеток [11], согласованная аппроксимация для расчета давления [7], позволившая избежать паразитных флуктуаций, а также прямой экономичный способ обращения матриц для решения регуляризованной системы линейных разностных уравнений. На каждой итерации последовательно рассчитывались функция тока р (1), компоненты вектора скорости V(У1, У2) (4), давление Р (5), вихревая функция ш (2) и энтальпия Н (3), плотность р (6). Этот алгоритм [11], реализованный в комплексе программ СОКОЛ (на основе комплекса программ РДТ ИММ [12]) для нескольких типов граничных условий, позволил эффективно рассчитывать как сжимаемые, так и несжимаемые потоки с присутствием многочисленных вихревых зон.
На решение систем разностных уравнений для функции тока и вихревой функции уходит основная доля процессорного времени (около 99 %). Оказалось, что при определенных условиях можно рассчитывать функции тока, вихря и энтальпию одновременно на разных процессорах. При распараллеливании алгоритма по типам уравнений для решения каждой из подсистем выделяется по процессору, и итерационный процесс строится по схеме, представленной на рис. 2 (стрелки обозначают пути передачи данных).
На первом процессоре рассчитывается функция тока р (1). На втором — последовательно вычисляются координаты вектора скорости У, У2 (4), вихревая функция ш и энтальпия Н (2)-(3), давление Р (5), плотность р (6). При вычислении какой-либо из этих функций значения остальных берутся с предыдущей итерации с соответствующего процессора. Еще один процессор выделяется для графической обработки результатов каждой итерации.
Работа многопроцессорной машины сначала была имитирована на одном процессоре. В соответствующие модули комплекса программ СОКОЛ были внесены необходимые изменения. Полученный комплекс программ СОКОЛ-М (многопроцессорный) поставлен на многопроцессорную вычислительную систему МВС. Во время работы параллельной программы были проведены замеры времени, которые показали, что при использовании для расчета течений газа двух процессоров коэффициент ускорения равен 1.5 - 1.7 и зависит от количества точек расчетной сетки. Чем больше точек, тем выше коэффициент. По сравнению с теоретическим (1.8 - 1.95) коэффициент оказался ниже, и связано это обстоятельство с потерей времени при обмене данными между процессорами. Для расчета
I процессор II процессор
ф° 0 т 0 0 т 0 0 о У>і, со , Н , Р р
ф1 у}, УІ1^\н\р\р1
1 1
ф2 у/\ф уЬ К22, оД#2,Р2,Р2
фк-1 X
к-ї'Х сох / к-1 \Ф
фк т К К К К К К У\, У2 , со , Н, Р , р
Рис. 2. Схема распараллеливания.
серии вариантов можно задействовать процессоров, где N — число вариантов. Каждый
третий процессор проводит быструю графическую обработку расчета варианта на каждой итерации и выводит результат обработки в соответствующее окно экрана. Параллельный расчет задач позволяет сократить время расчета примерно в 1.6N раз при счете N вариантов на 3N процессорах с одновременной визуализацией всех вариантов в разных окнах экрана.
Комплекс программ СОКОЛ-М позволяет проводить различные серии расчетов:
а) решение квазистационарной задачи течений газа в камерах сгорания;
б) корректировку экспериментально заданного давления торможения;
в) методические расчеты.
Для определения структуры и параметров газодинамических течений в камере сгорания двигателя (в квазистационарной постановке) необходимо рассчитать течение газа в серии областей различных геометрий, соответствующих разным моментам времени горения. На рис. 3, а представлены две конфигурации области на моменты времени ¿1, ¿2; на рис. 3, б, в — каpтины линий тока газа, рассчитанных на эти моменты времени. Видно, что в момент времени ¿1 течение газа было безвихревым, а в момент времени в головной части камеры сгорания двигателя образовался замкнутый вихревой поток.
При расчете газодинамического процесса в камере сгорания с сопловой частью оказалось, что если течение рассчитывать по заданным граничным условиям только в дозвуковой зоне без учета параметров в трансзвуковой, то могут возникнуть следующие ситуации: произойдет запирание потока в камере или в критическом сечении число Маха не достигнет единицы. При некотором изменении давления торможения, взятого из эксперимента, эти явления исчезнут. На основе методики А. Д. Рычкова [6] и его программы расчета течения в трансзвуковой зоне разработан алгоритм корректировки давления торможения, который заключается в поочередном расчете течений в дозвуковой и трансзвуковой зонах и корректировке безразмерной плотности (а в конечном счете — давления торможения) во входных данных для расчета дозвукового течения из условия равенства расхода газа его приходу. При последовательном расчете эта процедура занимает довольно большое время счета — от одного до нескольких часов.
Рис. 3. Расчет квазистационарной задачи: а — конфигурации расчетных областей; б, в — линии тока газа на моменты времени ¿1,^2-
При параллельном расчете N вариантов для одной и той же геометрии с одними и теми же входными условиями, но с разным давлением торможения (заданным с некоторым шагом) этот процесс сокращается по времени практически до расчета одного варианта. На рис. 4 приведен пример расчета в модельной области течения газа при следующих параметрах: 7 = 1.4,Я = 300, = 1, газ вдувается по поверхности АВ с постоянной
скоростью У0 = 1, в торце канала на ОБ скорость газа распределена по косинусоидальному закону [11], так что скорость газа на оси симметрии в точке Б больше скорости в точке О в 5.5 раз, выход газа — в сопловой части ЕЕ. Рассчитывалось 11 вариантов течения газа при заданном давлении торможения Р0 от 100000 до 200000 с шагом 10000.
Из рис. 4 видно, что картина линий тока в камере в диапазоне изменения давления торможения Р0 от 100000 до 150000 неустойчива. Расчеты показали, что давление торможения Р0 = 150000 обеспечило в критическом сечении число Маха, равное 1, и при увеличении Р0 структура течения не меняется.
Рис. 4. Расчет течений с целью корректировки давления торможения при Ро = 100000 (а), 120000 (б), 130000 (в), 150000 (г).
С целью исследования влияния сжимаемости жидкости на структуру течения в ряде каналов проводились расчеты при одних и тех же граничных условиях для несжимаемой и сжимаемой жидкостей. В некоторых областях в случае несжимаемой жидкости течение в результате расчета получалось безвихревым, в случае сжимаемой — образовывались циркуляционные зоны. На рис. 5, б представлен расчет течения несжимаемой жидкости,
Рис. 5. Методические расчеты: а — расчетная сетка, б — несжимаемая жидкость, в — сжимаемая жидкость.
а на рис. 5, в — с учетом сжимаемости. Сравнение подобных картин показывает, что в случаях большого диапазона изменения числа Маха, а здесь оно изменялось от 0.005 до 0.8, необходимо учитывать сжимаемость жидкости.
Аналогичные расчеты течений газа проводились в одних и тех же каналах для разных типов входных граничных условий: по заданным приходу массы газа, энтальпии и по заданным плотности и скорости вдува газа. Расчеты показали, что хотя плотность и скорость при первом наборе данных на участках вдува мало отличались от постоянных, заданных при втором наборе данных, в первом случае, как правило, рассчитывалась неизоэнергети-ческое течение, а во втором — изоэнергетическое. В результате в ряде областей по этим входным данным рассчитывались разные структуры течений — при расчете течений газа по входным данным второго типа условий в некоторых областях образовывались циркуляционные зоны. Сравнение аналогичных расчетов показали, что при построении картин реальных течений важно выбирать правильную постановку граничных условий.
Параллельный расчет задач позволяет сократить не только затраты машинного времени при счете, но и затраты на подготовку исходных и обработку полученных данных.
Сравнительно малое число процессоров, а также существенно меньшие затраты времени, используемые для расчета одного варианта задачи, позволяют рассчитывать на многопроцессорной вычислительной системе одновременно серию задач, связанных между собой.
Список литературы
[1] Райзберг Б. А., Ерохин Б. Т., САмсонов К. П. Основы теории рабочих процессов в реактивных системах на твердом топливе. Машиностроение, М., 1972.
[2] ЛипАнов А. М., Бобрышев В. П., Алиев А. В. и др. Численный эксперимент в теории РДТТ. Наука, Екатеринбург, 1994.
[3] Ахмадеев В.Ф., Гусева Г. Н., Козлов Л.Н., Сухинин С. В., Хайруллина О. Б. Гидродинамические источники акустических колебаний в камерах сгорания. ЦНИИНТИКПК, М., 1990.
[4] КАпцов О. В. Новые решения двумерных стационарных уравнений Эйлера. Прикладная математика и механика, 54, вып. 3, 1990, 409-415.
[5] Ахмадеев В. Ф., Сидоров А. Ф., Спиридонов Ф. Ф., Хайруллина О. Б. О трех методах численного моделирования дозвуковых течений в осесимметричных каналах сложной формы. Моделирование в механике, 4(21), №5, 1990, 15-25.
[6] Рычков А. Д. Математическое моделирование газодинамических процессов в каналах и соплах. Наука, Новосибирск, 1988.
[7] ХАкимзянов Г. С., Яушев И. К. Итерационный метод расчета двумерных дозвуковых установившихся внутренних течений идеальной сжимаемой жидкости. Препринт №4-87, СО АН СССР, ИТПМ, Новосибирск, 1987.
[8] Хайруллина О. Б. Метод расчета блочных оптимальных сеток в двумерных многосвязных областях. Вопросы атомной науки и техники, Сер. Мат. моделирование физических процессов, 1, 1992, 62-66.
[9] Хайруллина О. Б. Построение блочно-регулярных оптимальных сеток. Там же, 1, 1994, 19-25.
[10] Хайруллина О. Б. Ускорение сходимости итерационного процесса при построении блочно-регулярных сеток. Там же, 1, 1995, 54-59.
[11] Хайруллина О. Б. Расчет стационарных дозвуковых вихревых потоков идеального газа в осесимметричных каналах сложных геометрий. Там же, 3, 1990, 32-39.
[12] Хайруллина О. Б. РДТ ИММ - комплекс программ расчета стационарных дозвуковых течений в каналах сложных конфигураций. В “Вычисл. технологии”, 1, №2, ИВТ СО РАН, Новосибирск, 1992, 327-333.
Поступила в редакцию 12 марта 1996 г.
