CC BY
32
УЧЕНЫЕ ЗАПИСКИ Ц А Г И Том XVIII 198 7
№ 3
УДК 533.6.011.55.011.6
ИСПОЛЬЗОВАНИЕ АДАПТИВНЫХ СЕТОК ПРИ РАСЧЕТЕ ДВУМЕРНЫХ НЕРАВНОВЕСНЫХ ТЕЧЕНИЙ ИДЕАЛЬНОГО ГАЗА
В. Л. Меньшикова
Предложен простой алгоритм построения неравномерных поперек ударного слоя адаптивных сеток, хорошо приспособленный для расчета двумерных задач сверхзвукового обтекания тел невязким воздухом, в котором протекают неравновесные химические реакции. В ходе расчета маршевым методом неравномерная сетка перестраивается несколько раз, в зависимости от перемещения зон высоких градиентов внутри ударного слоя, что позволяет рассчитывать тела большого удлинения без вмешательства программиста в процессе счета. Предложенный метод позволяет получать решения с заданной точностью с вдвое меньшим числом узлов, чем потребовалось бы для получения такого же решения на равномерной сетке.
Одной из характерных особенностей задач гиперзвукового обтекания тел идеальным газом с учетом неравновесных физико-химических процессов является наличие зон больших градиентов параметров внутри ударного слоя. Это обстоятельство заставляет всех авторов, проводящих численные расчеты подобных течений, использовать неравномерные сетки. Однако до сих пор генерация неравномерных сеток была процессом интуитивным, основанным на некотором предварительном знании решения, на аналогии с уже полученными решениями и т. д. Сетка обычно заранее фиксировалась и при необходимости изменялась в процессе счета. Такой подход может давать неплохие результаты при расчете дозвуковых зон течений около гладких затуплений, когда по заданным параметрам невозмущенного потока можно определить заранее, будет ли течение в ударном слое близким к равновесному или замороженному, и таким образом определить, в какой части ударного слоя будет локализована зона больших градиентов, и сгустить узлы расчетной сетки в этой области. Задача становится более сложной при расчете сверхзвуковой области течения. Эти расчеты проводятся, как правило, маршевыми методами. Зоны с большими градиентами в процессе счета перемещаются внутри расчетной области, поэтому предсказать их положение заранее невозможно.
На рис. 1 и 2 в качестве иллюстрации представлены распределения мольной доли окиси азота г]ко поперек ударного слоя на затупленной
пластине с радиусом затупления ^0=19 см в различных сечениях t = L/R0 (L — длина пластины). Пластина обтекается неравновесным воздухом, параметры набегающего потока соответствуют следующему режиму: высота Н — 71 км, скорость 1/<х> = 6,77 км/с, угол атаки а=10° и 35°. Для получения таких решений необходимы неравномерные поперек ударного слоя сетки, но чтобы их строить, надо знать решение заранее. Приведенные результаты получены автором методом прогонки, подробно описанным в работах [1]и[2], в котором используются аналитические формулы перехода к неравномерной сетке. Такой подход позволяет получать результаты с хорошея точностью, но он очень трудоемок. Задачу практически нельзя решить за один запуск в процессор, а приходится решать поэтапно, меняя сетку по ходу решения; в случае же, если решение расходится, что наблюдается при плохом выборе сетки, приходится возвращаться назад к решению задачи, рассчитанному для предыдущих сечений. Практически каждый расчет проводится дважды на каждом отрезке; сначала на некоторой сетке получается промежуточное (как правило, осциллирующее) решение, по которому подбираются более подходящие параметры сгущения сетки, и затем расчет проводится заново на хорошей сетке. При этом выбор параметров сгущения также требует определенного навыка, и нет гарантии, что он делается оптимально. Все сказанное убеждает в необходимости использования в задачах такого рода адаптивных сеток, т. е. сеток, которые генерируются автоматически в процессе счета и при этом настраиваются в зависимости от получаемого решения.
Вопросу построения адаптивных сеток при решении различных газодинамических задач в последние годы уделяется усиленное внимание. Не перечисляя отдельных работ, посвященных этому вопросу, укажем, что достаточно подробный обзор их содержится в работе [3]. Существуют различные подходы к построению адаптивных сеток. Одни из них основываются на принципе минимизации погрешности. Такое решение представлено, например, в работе [4], в которой функция, сгущающая узлы, строится как конечная сумма полиномов Чебышева с коэффициентами, определяемыми из условия минимума в смысле наименьших квадратов меры погрешности. Локальная погрешность аппроксимации при использованной разностной схеме второго порядка оценивалась по третьей производной, определяемой численно. Таким образом, построение функции, обеспечивающей сгущение сетки, сводится к решению задачи нелинейного программирования. Такой подход очень сложен и требует относительно больших затрат машинного времени для построения сетки.
В последнее время чаще используется другой, более простой принцип построения адаптивных сеток. Суть его в одномерном случае можно сформулировать следующим образом: необходимо выбрать узлы сетки так, чтобы некоторая положительная весовая функция W(|) была распределена равномерно на отрезке, т. е. чтобы выполнялось соотношение
s/+i
J* W (£) dX = const
Ч
или в дискретной форме
ht Wt = const, (1)
где hi^Zi+i — %1.
Для выбора функции W(l) в литературе используются разные подходы, один из них исходит из условия равномерности невязки или погрешности решения. В последнем случае в качестве W(|) берется
функция W (£) = /1*||ы*+1||, здесь tif'1 — производная порядка k по I от вектора решения, k — порядок аппроксимации. В работе [5] использовалась функция W ($) = ■/\ +1| и^Ц2. При применении ее в случае одного уравнения условие (1) дает соотношение
У I + ul Sjf = const, которое эквивалентно требованию равенства длин дуг между узловыми точками на кривой решения. В работе [6] предложено в качестве W (£) использовать функции, учитывающие кривизну решения кривой,
UP(5) = (1 + ?\к\)У 1+и?,
J+tiF
где k— (л %.3/2 , или № (Е), определяемую из уравнения
-Wx + °W=*\lfi\ + $\k\,
где а, о, ji — некоторые константы.
При выбранной функции W(|) построение функции |(х), переводящей равномерную сетку по х в неравномерную, осуществляется очень просто. Условие (1) эквивалентно дифференциальному уравнению
\х W (£) = const. (2)
Постоянная в уравнении (2) определяется из условия нормирования. Если при изменении | от 0 до 1 х также изменяется от 0 до 1, то,
используя соотношение X£ = 1/gx, получим
£
Г (I) dl
&
о
Тогда
' ( W (5) dZ
1
W (6) dl
Учитывая, что Ал:=1 /Ы, где N — число интервалов, на которое разбивается отрезок, получаем формулу, определяющую положение узлов неравномерной сетки
1
С «7 (?) dІ
А?,= ------------ • (3)
1 NW(S■t)
Вернемся к поставленной задаче расчета двумерного неравновесного сверхзвукового обтекания. Задача построения адаптивной сетки в этом
случае разбивается на две. Первая задача: по решению, полученному в некотором сечении, построить оптимальное в некотором смысле разбиение отрезка между телом и волной. Вторая задача возникает в связи с тем, что генерировать сетку на каждом шаге в маршевом направлении неэкономично и в этом нет необходимости, так как решение меняется постепенно и сетка, построенная по решению в некотором сечении, может хорошо работать на некотором интервале. Поэтому вто-
Рис. 3
рая задача заключается в том, чтобы определять моменты, когда сетку надо перестраивать. При этом может возникнуть необходимость не только перераспределять узлы, но и увеличивать число узлов, если за счет перераспределения их нельзя достигнуть необходимой точности.
Рассмотрим первую задачу. Прежде всего отметим, что получаемое в задаче решение— семикомпонентный вектор (оь и2, р, Т, т)о, Лк, т^о), где vu и2, р, Т — компоненты скорости, давление и температура; г|о, т)и, ч^о — мольные доли О, Ы, N0; т)о2 и определяются из соотношений материального баланса. Однако для построения адаптивной сетки будем исходить из результатов исследования поведения лишь одной компоненты этого вектора. Такое представляется возможным в связи с тем, что все компоненты вектора решения имеют одни и те же зоны больших градиентов, определяемые физической картиной течения — энтропийными слоями и релаксационными зонами, см. рис. 3. Возможность построения хороших неравномерных сеток по поведению одной компоненты решения подтверждает также опыт работы с неадаптивными неравномерными сетками. В качестве параметра задачи, используемого для настройки сетки, естественно взять наиболее чувствительный к выбору сетки, раньше других начинающий осциллировать параметр. В рассматриваемой задаче им является величина «(|) =т]но (см. рис. 1 и 2).
, Остановимся теперь на выборе весовой функции й^(|) в рассматриваемой задаче. Будем исходить из того, что” необходимо решить задачу с заданной точностью как можно более экономично, т. е. используя возможно меньшее число узлов и как можно реже перестраивая сетку. В расчетах были опробованы три весовые функции И^(£):
\)У I + и\\ 2) А|нц|/А,, где А, = 1/Ы ;
3) (1 +р\к\)У\ + и\;
здесь ы(§)=гц«), «е и — соответственно первая и вторая производные «(!);■£(£)—кривизна кривой и(£).
Оказалось, что первая из перечисленных функций, хорошо зарекомендовавшая себя в ряде задач [5], для рассматриваемой задачи не подходит, так как не способна сгущать точки в областях большой кривизны, характерных для получаемых решений. Функция =
= к\ иц\/к1 использовалась в работе [7] для решения ряда модельных задач, в частности, для решения уравнения Бюргерса с начальными условиями в виде ступеньки. В этой работе показано, что использование такой весовой функции сильно сгущает точки в окрестности скачка, что позволяет увеличить точность расчета или при заданной точности вдвое уменьшить число узлов по сравнению с равномерной сеткой. Но, как показала практика использования этой функции в рассматриваемой задаче, сетка, получаемая с помощью функции Ш (I) = к | |/Л,,
слишком жестко связана с решением и поэтому требует очень частой перестройки. Кроме того, эта сетка очень неравномерна, что приводит к росту погрешности при интерполяции параметров в моменты перестройки сетки.
Наиболее подходящей для рассматриваемой задачи является
функция IV (Е) = (1 + ;а | к |) V 1 + и,\. При этом значение р. должно изменяться в процессе счета в зависимости от значений градиентов и (I) в слое. В сечении, где происходит перестройка сетки, величина |а полагается пропорциональной значению максимального градиента тШах, коэффициент пропорциональности выбирается из условия, чтобы величина [а изменялось от 0 до 0,5 при изменении «5шах от 1 до 10. Величина интервалов Д|г не должна превышать некоторого Д£шах) ТЭК как при больших |А сгущение точек в областях большой кривизны будет очень сильным и в остальной области может не хватить узлов для продолжения счета. Если величина
становится больше Д?твх, то она полагается равной Д£тах, а остальные Д|, соответственно увеличиваются, чтобы выполнилось
N
условие Д^ = 1.
1=1
Описанная процедура позволяет по полученному в некотором сечении решению и{%) строить с помощью функции М^(?) = [1 +
+ ^ («е;шах) | к | ] V 1 + и\ по формуле (3) новое разбиение Д&,. Но знания функции \ (х) только в дискретной форме по ее значениям в узлах хь обычно бывает недостаточно. При решении задачи необходимо знать, как правило, функции | (х) и \х в дробных узлах, а также уметь определять значения обратной функции х (£) в произвольных точках, что требуется, например, при интерполяции параметров с одной сетки на другую. С этой целью по значениям функции ?(хг) строится кубический сплайн по программе, описанной в работе [8], и определяются обратная функция и ее производные в произвольных точках.
Перейдем теперь ко второму вопросу: определение моментов перестройки сетки.
Характерной особенностью функции «(Е) =тшо(Е) является наличие максимума внутри ударного слоя. Максимум ы(|) перемещается внутри интервала [0,1], и его смещение можно принять за критерий, по которому перестраивается сетка. Этот процесс реализуется следующим образом: задается некоторая начальная величина А, и расчет, начиная с некоторого начального сечения, в котором строится соответствующая решению сетка, продолжается до тех пор, пока максимум функции
и{\) не сместится на величину А. Тогда сетка перестраивается по решению в этом сечении.
В процессе счета происходит контроль за качеством получаемого решения: проверяются выполнение законов сохранения и величина осцилляций решения. Если погрешность выполнения интеграла Бернулли превышает некоторую разрешенную величину или осцилляции решения превышают допустимое значение, программой вырабатывается признак, обозначающий, что решение испорчено. В процессе счета происходит запись решения на диск для того, чтобы, в случае необходимости, решение могло быть продолжено с места записи. В случае получения признака порчи решения программа возвращается в точку последней записи и возобновляет счет оттуда с другим, более мелким значением Д. Обычно этого хватает, чтобы продолжить расчет до момента следующей записи решения. Если признак испорченности решения появляется раньше, то программа опять начинает счет с последней записи, но с увеличенным числом узлов. В случае необходимости этот процесс может быть повторен. При этом, естественно, существуют ограничения на максимальное число узлов в рассматриваемом направлении и на минимальное значение А.
Используя адаптивную сетку, построенную по описанной методике,, удалось за один запуск в процессор получить решения, приведенные на рис. 1—3, с той же точностью, которая была при настройке сетки вручную.
Функция сгущения узлов получающаяся в процессе счета и
использующаяся в сечениях, соответствующих рис. 1—3, показана на рис. 4 и 5. Как видно, адаптивная сетка хорошо отражает особенности решения. Сетка в ходе таких расчетов перестраивалась около десяти раз. Использование адаптивной сетки позволяет получить решение с той же точностью, что и на равномерной сетке с вдвое большим числом узлов.
ЛИТЕРАТУРА
1. Бабенко К. И., Косоруков А. Л., Радвогин Ю. Б. Сверхзвуковое стационарное двумерное обтекание воздухом с учетом неравновесных химических. реакций. — Препринт ИПМ АН СССР им. М. В. Келдыша, № 71, 1981.
2. К о с о р у к о в А. Л., Меньшикова В. Л., Р а д в о г и н Ю. Б. Влияние неравновесности обтекания на аэродинамические характеристики затупленной пластины с изломом. — Ученые записки ЦАГИ, 1985, т. 16, № 2.
3. Thompson Т. F. A survey of dynamically—adaptive grid in the solutions of partial differential equations. — AIAA Paper 84—1606.
4. P i e r s о n B. L., Kutler P. Optimal nodal point distribution for
improved accuracy in computational fluid dynamics. — AIAA J., 1980,
vol. 18, N 1.
5. W h i t e, Andrew B. Jr. On the numerical solution of initial boundary — value problems in one space dimension. — SIAM J. of Numer. Analysis, 1982.
6. Ablow С. М., Schechter S. Campylotropic coordinates.— Journal of Computational Physics, vol. 27, 1978.
7. Haase W., Misegades K., N a a r M. Adaptive grid in numerical fluid dynamics. — Internal. J. for Numerical Methods in Fluids, 1985 VI, vol. 5, N 6.
8. Бабиков П. E. Модифицированная версия комплекса программ для вычисления геометрических характеристик таблично заданных поверхностей сложных конфигураций. (АРГОЛА, вып. 13а). — Труды ЦАГИ, 1982, вып. 2134.
Рукопись поступила 6/1 1986
