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

Том XXXIII

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

№ 3—4

УДК 519.63:517.96

629.7.015.3.036:533.697.2

РАЗРАБОТКА МЕТОДА ПОСТРОЕНИЯ АДАПТИВНОЙ МНОГОБЛОЧНОЙ СЕТКИ ДЛЯ РАСЧЕТА СЛОЖНЫХ ДВУМЕРНЫХ ТЕЧЕНИЙ НА БАЗЕ ПОЛНОЙ СИСТЕМЫ УРАВНЕНИЙ ЭЙЛЕРА

А. В. ЛЫСЕНКОВ

Рассмотрена явная схема решения уравнений для метода построения адаптивных сеток Брекбилла — Зальцмена. Предложен оригинальный способ расчета коэффициентов метода.

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

В настоящее время существует много методов, посвященных решению задач, связанных с расчетом сложных внешних и внутренних течений газа. Одним из распространенных является метод, основанный на технологии многоблочных структурированных сеток. Этот подход является оптимальным для описания течений со сложной топологией (в том числе — неодносвязных расчетных областей). Для повышения надежности расчета необходимо, чтобы сетка была адаптирована не только к геометрии расчетной области, но и к структуре течения. Так, для более детального описания областей с большими градиентами параметров желательно, чтобы в этих областях ячейки сетки были более мелкими. Такие сетки называются адаптивными.

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

Для численного решения дифференциальных уравнений при адаптации сетки использована явная численная схема. Во всех работах по адаптации сеток, известных автору, а также в большинстве работ по численному интегрированию уравнений Навье — Стокса для решения параболических уравнений используются неявные численные методы. Главным аргументом в пользу неявных схем остается их абсолютная устойчивость (при правильном выборе схемы) и возможность расчета с любым шагом по времени. Однако большие шаги по времени порождают большие численные погрешности; при этом неявные схемы требуют на порядок больших компьютерных ресурсов, чем явные. При использовании неявных схем все усилия вычислителя сосредоточены на решении громоздких систем линейных уравнений для всех ячеек расчетной области.

Явный же подход позволяет сосредоточить внимание не на решении алгоритмических проблем, а на физических и математических процессах, протекающих в каждой отдельной

ячейке; явные схемы предоставляют гораздо большие возможности по управлению процессом расчета на всех стадиях его развития. Опыт работы последних лет группы вычислительной аэродинамики НИО-1 ЦАГИ показал, что явные схемы позволяют решать уравнения Навье — Стокса с той же скоростью, что и неявные, при гораздо меньших затратах оперативной памяти компьютера. В настоящей работе показано, что явный подход может быть столь же эффективно использован и в задачах адаптации расчетных сеток.

Здесь открываются большие возможности. Так, исследовалась адаптация сетки к структуре численного решения, уже полученного на квазиравномерной сетке; наиболее оптимальным, с точки зрения точности получаемого результата, является адаптация сетки на каждом шаге расчета, когда сетка движется вместе со структурой численного решения. Для того чтобы процесс пересчета сетки на каждом шаге не занимал много времени, нужны эффективные, быстрые, дешевые алгоритмы адаптации сетки; здесь явные схемы имеют преимущество перед неявными.

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

Описание метода Брэкбилла — Зальцмена. На рис. 1 изображена одна ячейка расчетной сетки. Пусть (х, у) — координаты в физической плоскости, а (^, ц) — криволинейные

координаты, связанные с сеточными линиями. Масштаб вдоль осей £, и ц выбирается таким образом, чтобы в плоскости (^, ц) каждая ячейка являлась квадратом со стороной 1. (На

практике (^, ц) - индексы, которыми нумеруются ячейки сетки.) Тогда координаты вектора АВ (рис. 1) в физической плоскости приближенно равны (х^; у^ = (ц^; -цх3), где

3 = = х,1’Т) - хТ)1’, — якобиан преобразования координат; вектор

СО « (хТ); _)’Т] ) = (-^,,7; с,х-1^; площадь ячейки Я « | АВ х С!) ^ |./ |.

Рис. 1. Ячейка расчетной сетки и системы координат

Основной задачей является построение сетки, которая будет отвечать следующим требованиям:

1. Сетка должна сгущаться в областях больших градиентов параметров потока. Пусть У/ (х, у) — градиент какого-нибудь характерного газодинамического параметра. Потребуем, чтобы мелкие ячейки сетки были сосредоточены в тех областях, где значение У/ велико. Для идеальной сетки произведение | У/ | -Б « | У/ | -3 должно быть одинаковым для всех ячеек. Это достигается, если минимизировать интеграл 1У = Д^ | У/12 32й^ . Вместо этого в методе Брэкбилла — Зальцмена предлагается минимизировать интеграл

Е =

Л Еуй-йц, Еу = м>32, (1)

где положительная весовая функция м(х, у) ~| V/12. Как будет показано ниже, функцию w(х, у) можно подобрать так, чтобы ограничить отношение максимальной площади ячейки к минимальной.

2. Ячейки сетки должны быть как можно ближе к квадратным. На идеальной сетке с

---- о --------- о 2 2 2 2

„ | АВ | +1 СО | х(+х*+у(+у*

квадратными ячейками отношение ----------------« —-----!--2----1 одинаково для всех ячеек.

Поэтому нужно минимизировать следующий интеграл:

2 2 2 2

I, = Ц^дец, ^ х + Хц +у + Уц. (2)

Интеграл (2) можно представить в виде

I, =Ц[(V-)2 + ^ц)2]йхйу .

Поэтому в работе [2] условие минимизации интеграла (2) интерпретируется как условие гладкости сетки, которое требует наименьшей амплитуды колебаний -и ц при изменении аргументов х и у, т. е. предельно ровного «рельефа» сетки.

3. Сетка должна быть как можно ближе к ортогональной: углы между сеточными линиями должны быть близки к 90°. Мерой ортогональности может служить величина

АВ ■ СО ~ х=хТ) + >у’Т) (для ортогональной сетки она равна нулю). Тогда сетка будет максимально

близка к ортогональной при минимальном интеграле

10 = ЦЕ0йц = (х-хц + у-Уц )2 . (3)

Перечисленные требования к сетке, противоречат друг другу и не могут быть выполнены полностью. Для достижения компромисса можно минимизировать сумму интегралов (1) — (3), взятых с некоторыми весами:

1 = 1 о + ^г1г + К1., • (4)

Уравнения сетки и их решения. Для минимизации интеграла (4) используется известный в вариационном исчислении метод Эйлера [3]. Интеграл (4) представляется в виде

I =

Е й-йц,

где Е = Ео + . Тогда система уравнений метода Эйлера выглядит следующим образом:

дЕ д д

= 0,

дх д<- дх- дц дхц

= 0.

ц

дЕ д дЕ д дЕ

ду д- ду- дц дуц

(5)

Соотношения (5) должны выполняться в каждой ячейке расчетной сетки.

Подробный вывод уравнений для рассматриваемой вариационной задачи дется в работе[2]. Конечный результат записывается в виде:

+ (2Ж = &, (6)

где

Б =

(41 ( х 1

, V =

V ^2 V У ,

Ь —

аЪ да дЪ ’

а Л\, Л2, Л3 — матрицы 2 х 2, зависящие от первых производных (Х;, х^, у;, у^). Вектор К —

это вектор координат узлов сетки. В вектор Б входят производные от | V/ |2 по х и у.

Для решения уравнения применяется псевдонестационарный метод. Для этого к уравнению

(6) добавляется искусственный нестационарный член:

( + А1( + (2(4л + АЭЬЛЛ = Б ’

'4л

ЛЛ"

(7)

т д где Ь, = -

Решение уравнения (7) будет совпадать с решением уравнения (6) при , ^ да. В качестве начального приближения выбирается квазиравномерная сетка, не адаптированная к решению.

В отличие от [2], в настоящей работе предлагается использовать для решения уравнения (7) явную схему с центрально-разностной аппроксимацией вторых производных по пространству:

(

Уп+1 = Vй. +т

1,1 1,1

УПЛ - 2УП + УП

‘ -1,1

а;

‘, 1 + У ‘+1,1 , А ' ‘,1 -1

2 + А3

V" 1 - 2УП + V

‘,1 +У ‘,1+1

Ал2

+Ап

Уп - Уп Уп - Уп

У‘+1,1+1 У‘+1,1-1 У‘-1,1+1 У‘-1,1-1

2 а;

2 а;

(8)

2Ал

Решение считается установившимся, когда величины, входящие в уравнение (8), перестанут меняться со временем. Для устойчивого счета необходимо, чтобы шаг по времени (т) удовлетворял следующему ограничению:

1

т =

1 1

— + — Тк Т„

Т% <

к

4

к

4тах (Л11, Л12 )

т <-------------------- ,

л 4тах (Л21, Л22 ) ’

где Л п и Л12 — собственные числа матрицы Л], а Л 21 и Л 22 — матрицы Л2.

Методы расчета коэффициентов. В методе Брэкбилла — Зальцмена есть несколько управляющих параметров: коэффициенты Хи Хж и весовая функция w(х, у) .

Возможны различные способы расчета коэффициентов. В работе [2] для этого используется оценка порядков величин интегралов в уравнении (4). Можно показать, что

I ~1,I ~ wh4,Iо ~ к4 .

(9)

Для того чтобы все три условия (гладкость, ортогональность и сгущение) выполнялись в равной мере, необходимо уравнять порядки слагаемых в уравнении (4). Тогда

X „ = wk4 , Xо = w .

(10)

В настоящей работе вместо оценок (12) предлагается другой способ, основанный на непосредственном вычислении интегралов в уравнении (4). Тогда коэффициенты можно вычислить следующим способом:

X, =Р^, Х0 =р0/° Ро,Р5 ~0(1).

V V

2

Как показывает опыт, если взять коэффициенты р0 и равными единице, то адаптация будет очень слабая. Для расчетов, которые будут описаны ниже, были выбраны значения Ро = 0,5, р5 = 0,05 .

В настоящей работе используется функция w(х, у), которая зависит от градиента плотности газа, поскольку эта величина хорошо отслеживает все детали течения: волны сжатия и разрежения, контактные разрывы, энтропийные слои и т.п.

^, 1 = Е + | Vpj,у |2 . (12)

Положительная величина Е меняет отношение максимального значения функции ^ к минимальному и таким образом управляет силой сгущения (если бы адаптация сетки

w

проводилась путем минимизации только интеграла (1), то отношение —тахк было бы равно

w

тт

отношению площадей наибольшей и наименьшей ячеек). Неотрицательные коэффициенты с^ позволяют корректировать функцию ^ для сгущения сетки возле границ и изломов геометрии, а также в окрестности тех особенностей течения, которые представляют особый интерес.

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

Первый подход состоит в том, что граничные узлы в ходе установления решения уравнения

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

Результаты. Для тестирования метода адаптации сетки в настоящей работе используется задача о двумерном сверхзвуковом течении в канале сжатия (рис. 2). Верхняя поверхность канала плоская, нижняя — клин, переходящий в плоскую стенку. На входе — сверхзвуковой поток (М = 2,88). Относительная высота клина к^ё =0,25; относительные длины клина и канала вдоль

оси х Ьёё = 1,5 и = 5 (все размеры отнесены к ширине канала). Течение считается невязким. Данная задача моделирует течение в воздухозаборнике ПВРД. Расчеты проводились при помощи комплекса программ «80ЬУБЯ3» [4]. В этом комплексе для численного решения полной системы уравнений Эйлера используется явная монотонная схема Годунова — Колга-на — Родионова [5]—[7], имеющая 2-й порядок аппроксимации по всем переменным.

Рис. 2. Расчет двумерного течения в канале (поле числа Маха)

Расчеты проводились в несколько этапов: строилась квазиравномерная сетка, на ней проводился расчет течения, затем сетка адаптировалась и расчет течения повторялся еще раз. Возможно и большее количество итераций.

В рассматриваемой тестовой задаче расчет на квазиравномерной сетке, не адаптированной к структуре течения, дает решение, которое имеет ряд недостатков. Имеются два наиболее существенных недостатка. Первый состоит в том, что скачки уплотнения и другие детали структуры течения слишком сильно «размазываются» из-за численной диссипации схемы. Второй — это наличие энтропийного слоя. Это нефизичный слой с пониженным полным давлением, который возникает у твердых поверхностей за скачками уплотнения и волнами разрежения. Адаптация сетки нужна для уменьшения ошибок, связанных с этими явлениями.

Расчеты показывают, что уменьшение ширины «размазывания» скачка происходит при любом способе адаптации сетки и зависит от параметра Е в функции w. В настоящей работе параметр Е вычислялся таким образом, чтобы отношение максимального значения весовой функции w к минимальному было равно 100.

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

4).

5

ё

С?

у/Н

а)

—г~ М

у/Н

б)

0,6 0,8

Рис. 3. Сравнение расчетов с разными способами распределения граничных узлов (поля полного давления и распределение полного давления поперек скачка уплотнения): а — неподвижные граничные узлы; б — координатные линии перпендикулярны границам; в — граничные узлы распределены

подобно первому уровню узлов —————- — адаптивная сетка; —- — — — неадаптивная сетка

Рис. 4. Влияние сгущения сетки возле изломов геометрии (поля полного давления): а — расчет без сгущения возле изломов геометрии; б — расчет со сгущением возле изломов

геометрии

Представляет интерес сравнение расчетов на неадаптивной сетке и на сетке, адаптированной один раз. Адаптация проводилась наилучшим методом, т.е. узлы по границе распределялись подобно первому внутреннему слою, и использовалось сгущение возле изломов геометрии.

На рис. 5 сравниваются полученные в этих расчетах поля чисел Маха и полного давления. Видно, что расчеты на адаптивной сетке гораздо лучше, чем на неадаптивной. Во-первых, скачки уплотнения менее «размазаны» (см. распределение чисел Маха поперек скачка на рис. 6): ширина «размазывания» уменьшилась приблизительно в 4 раза. Во-вторых, нефизичные потери полного давления в энтропийном слое уменьшились примерно на 40% (см. отличие пунктирной кривой от сплошной в левой части графика на рис. 6).

Рис. 5. Сравнение расчетов на адаптивной и неадаптивной сетках: а— неадаптивная сетка; б — поле полного давления на неадаптивной сетке; в — поле числа Маха на неадаптивной сетке; г— адаптивная сетка; д — поле полного давления на адаптивной сетке; е — поле числа Маха на адаптивной сетке

-0,2 О 0,2 0;4 0,6

У / н

Рис. 6. Сравнение расчетов на адаптивной и неадаптивной сетках (распределения числа Маха в сечении, перпендикулярном скачку уплотнения):

---------------адаптивная сетка;--------равномерная сетка

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

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

При втором подходе многоблочная сетка адаптируется так, что грани могут свободно перемещаться, но некоторые углы блоков фиксируются. Для каждого блока граничной считалась первая внутренняя сеточная линия соседнего блока, если такой есть. То есть границы блоков, которые не являются границами области, участвовали в адаптации как внутренние сеточные линии.

При использовании первого метода (рис. 7, а) оказалось, что нужно настраивать алгоритм для каждого блока отдельно, так как для разных блоков необходимы различные коэффициенты метода адаптации. Кроме того, фиксированные грани блоков ограничивают «свободу» узлов, что приводит к ухудшению характеристик адаптации.

Рис. 7. Два подхода к адаптации многоблочной сетки: а — адаптация многоблочной сетки по блокам; б — адаптация многоблочной сетки, как одноблочной

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

В результате изложенного можно сделать следующие выводы:

1. Предложена явная схема для решения уравнений метода построения адаптивных сеток Брэкбилла — Зальцмена.

2. Предложен оригинальный способ расчета коэффициентов метода.

3. Проведено сравнение различных способов распределения граничных узлов. Наилучшие результаты получены в случае, когда граничные узлы распределялись подобно первому внутреннему слою и узлы дополнительно сгущались к изломам геометрии.

4. Сравнение результатов тестовых расчетов на адаптивной и неадаптивной сетках показало, что ширина «размазывания» скачков уплотнения уменьшилась в 4 раза, а нефизичные потери полного давления в энтропийном слое уменьшились на 40%.

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

Автор выражает благодарность своему научному руководителю, кандидату физикоматематических наук Власенко В. В. за помощь в подготовке настоящей работы.

ЛИТЕРАТУРА.

1. Brackbill J. U. and Saltzman J. S. Adaptive zoning for singular problems in two dimensions//J. of Computational Phisics.— 1982. Vol. 49.

2. Крейс Р. И., Теймз Ф. К., Хасан Х. А. Построение адаптирующихся сеток с помощью вариационного метода Брэкбилла-Зальцмена//Аэрокосмическая техника.—1987,

№ 1.

3. Федорюк М. В., Обыкновенные дифференциальные уравнения.— М.: Наука.—

1980.

4. V l a s e n k o V. V., S a b e l n i k o v V. A. Numerical simulation of inviscid flow with hydrogen combustion after shock waves and in detonation waves//AIAA-94-3177.— 1994.

5. Годунов С. К., Забродин А. В., Иванов М. Я., Крайко А. Н., Прокопов Г. П. Численное решение многомерных задач газовой динамики.— М.:

Наука.— 1976.

6. Колган В. П. Применение принципа минимальных значений производной к построению конечно-разностной схемы для расчета разрывных решений газовой динамики//Ученые записки ЦАГИ.— 1972. Т. III, № 6.

7. Родионов А. В. Монотонная схема 2-го порядка аппроксимации для сквозного расчета неравновесных течений//ЖВМ и МФ.— 1987. Т. 27, № 4.

Рукопись поступила 10/V 2001 г.