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

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

Том XXXVI

200 5

№ 3 — 4

УДК 533.6.011.3/.5 + 519.6

НОВЫЙ МЕТОД ИСПОЛЬЗОВАНИЯ ПРИНЦИПА МИНИМАЛЬНЫХ ПРИРАЩЕНИЙ В ЧИСЛЕННЫХ СХЕМАХ ВТОРОГО ПОРЯДКА АППРОКСИМАЦИИ

С. В. МАТЯШ

Предлагается новый метод расчета приращений значений функций в численной схеме Годунова — Колгана — Родионова, основанный на принципе минимальных приращений (ПМП). При вычислении приращений компонент скорости принцип применяется только к длине вектора, а на направление не распространяется. Данный подход приводит к уменьшению энтропийных слоев, уменьшению чувствительности сильных разрывов к изломам сетки, повышению устойчивости схемы и, соответственно, расширению возможности использования методов ПМП типа Ван Лира.

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

В одномерных задачах ПМП реализуется следующим образом. Обозначим пространственную координату буквой х. Пусть центр расчетной ячейки имеет индекс 1, центр соседней ячейки справа — (/ +1), слева — (i -1). Обозначим левую границу ячейки индексом 1 -12, правую — индексом i +1/2. Тогда значение приращения параметра/ для границы i +1/2 вычисляется по формуле:

где введено обозначение Аt+y2а = at+1 - a, а limiter (Aj, A2 ) — функция, ограничивающая

значение приращения. Первый вариант такой функции был предложен в работе Колгана [1].

Несколько модифицированный вариант функции Колгана был детально исследован в работе [2] и получил название «minmod»:

Как показано, например, в работе [3], minmod (А1, А2 ) — наиболее диссипативная из всех функций-ограничителей, при использовании которых выполняется ТУБ-свойство схемы. Среди всех таких функций minmod(А1, А2) порождает максимальное сглаживание неоднородностей

(1)

(2)

численного решения. Существуют и менее диссипативные функции-ограничители. Один из наиболее удачных классов функций-ограничителей предложен Ван Лиром [4]:

limiter(А1, А2 ) = VL(А1, А2 ) =

sign(А1 )• min PjA^, р|А2

|а 2h

А1 ^А2 < 0,

А1 •А 2 > 0.

(3)

Здесь коэффициент Ре[1; 2]. В оригинальной работе [4] Р = 2. Применение данной функции

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

Обобщение ПМП на многомерные задачи оказывается проблемой, не имеющей однозначного решения [5]. В данной работе для сопоставления результатов используется следующий способ трехмерного обобщения ПМП (см. работу [6]).

Будем предполагать, что параметр / распределен линейно вдоль индексных осей, которые проходят через центры ячеек и центры граней ячеек. Тогда получаем следующую формулу для приращения параметра / вдоль индексной оси I:

fi+1/2 fijk

AR ?

• limiter

Ai-12f

Af-1?

A L?

Ai+1/2f

\

A R ?

aL+1?

(4)

A Lr

= \?,jk - Г-1/2

A R ?

= r

i+1/2 ijk •

где введены обозначения: ? = (x; y; z),

Существуют и другие способы трехмерного обобщения ПМП. Один из известных методов был предложен Тилляевой [7]. Этот способ менее чувствителен к изломам индексных осей, чем вышеизложенный. Однако метод (4) гораздо более экономичен и, в отличие от метода Тилляевой, всегда гарантирует, что f+^2 min (f, f+1); max (f, f+1 )J. Достаточно полный

обзор рассматриваемых методов приведен в работе [5].

2. Новая реализация ПМП для векторных параметров. Недостаток известных трехмерных обобщений ПМП заключается в том, что ко всем величинам полей, в том числе и компонентам вектора скорости V = (u; v; w), эти методы применяются как к независимым

величинам. Очевидно, что вектор скорости является единым физическим объектом и это необходимо учитывать в алгоритме.

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

однако, этого не происходит. Значения коэффициента давления ер в некоторых точках могут

отличаться на величину 0.06 при числе Маха набегающего потока 0.1. На рис. 1 приведено поле разности функции cp в пространстве, для основной и повернутой расчетных сеток около

профиля Жуковского. Профиль рассчитывался при а = 3°, у = 15°, М = 0.1 с использованием метода (4) с функцией-ограничителем (2), примененного к скорости покомпонентно.

На рисунке видно, что зона наибольшего расхождения сосредоточена в районе поверхности профиля. Это сказывается в конечном итоге на расхождении значений аэродинамических коэффициентов cx и cy. В первом случае cx = 0.0102, cy = 0.9457, а в случае повернутой сетки

cx = 0.0098, cy = 0.9463.

Рис. 1

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

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

Рассмотрим предлагаемый алгоритм на примере одного из сеточных направлений — I. Пусть центр ячейки имеет индекс і, центр соседней ячейки справа — (і +1), слева — (і _1). Индекс центра правой грани рассматриваемой ячейки і +1/2, левой і _1/2. Следует отметить, что

в рассматриваемом алгоритме координаты центра каждой ячейки равны средним арифметическим координат ее вершин, а центр грани — среднее арифметическое вершин грани. Поэтому центр ячейки находится на равном расстоянии от центров каждой пары противолежащих гра-

ней

Л)" Г

Л* Г

Пусть вектора скорости в центрах ячеек равны \і_1, \і, \і+1 соответственно.

На первом этапе определим направления скоростей на боковых гранях ячейки і следующим образом:

Х:-У2=Ч-Л% Х:+1/2=Ч+Л%

где У 1+у2, У1_у2 — скорости на гранях,

АУ/ = ^(аУь+АУк),

ау; =

лг-1/2У

Аг_1?

Агь?

Ак?

АУр =

лг+1/2У

Агк?

лгЬ+1?

Ак?

Единичные векторы вдоль скорости на гранях имеют вид

- _ У-1/2 -

ег-1/2 “ -, 5 ег+1/2

У'+1/2

В случае, когда

V'

У-1/2

= 0, полагается е = 0.

У'+1/2

На втором этапе расчета приращений компонент вектора скорости применим ПМП к длинам векторов:

У

= Ншкег^А Ук ; А Уь ),

(6)

где

Ут

А-1/2 V

Агк_1? + Дгь?

Ак?

У»

^г+1/2 V

Агк? + д*1?

Ак?

В качестве функции-ограничителя в (6) можно использовать, например, функцию тттос1(А1, А2) (2) или уь(а15а2) (3).

На третьем этапе, имея направления на гранях, получим векторы скорости на гранях:

Уг-1/2 _ ®г—1/2 ‘ ( У <1 1 У

У+1/2 = ®г+1/2 ‘ ( У <1 + У

(7)

Найдем новые приращения вектора скорости вправо и влево от центра рассматриваемой ячейки

АУь=Уг-Уг_1/2,

ЛУк=Уг+1/2-Уг,

На четвертом этапе симметризуем полученные приращения:

Аиг

АУг =\ Аг,

(8)

АУЬ+АУК У/+1/2 ~ Уг-1/2

(9)

Это и есть искомое приращение компонент скорости вдоль ОДНОГО из сеточных направлений. Для скалярных величинр, р, Т и т. д. применяется формула (4). Для анализа данный

алгоритм сложен. Однако в одномерном случае (V = и) можно показать, что схема имеет второй

порядок аппроксимации. Будем полагать шаг сетки постоянным. Тогда А,и' =

Аг-1/2и + Аг+1/2и

одномерное приращение, соответствующее величине АV| в формуле (5).

Далее А,|и| = ^Ншкег(А,-^2 |и|; А,+^2 |и|) — одномерное приращение модуля скорости,

соответствующее величине а|V в формуле (6).

Тогда, подставляя полученные выражения в формулы (7) — (9), получим следующее выражение для приращения скорости:

Аиь + Аик _ и,+12 - и, -1/2 _ 1

Аи = -

, Лил + А,|и| , Лил -А,|и

(иг + А,и') ^----------^-(и, - Аги'У г 1 г 1

и + а и

и, - Аи'\

\и1 + А и и, - А и +а и' и, +А и и, - А и

и, + А и' и, - а и: и, +А и] и - А и]

(10)

Полагая |А ,и'|□ |и,|, |А,|и|| 0 |и ,|, находим

( л ( л

1 Л 1 1 и и, + А и 1 + V и, у \и,\ 1+ V и1 у

и, +А и] | | ы 1+а и' и \и,\ (1+а“'1 V и, у

=1+ам _Аи:+о

и I и,

( _ _ 2 Л

А, и'

_ и, _

V у

(

здесь учтено, что для малых приращений

, А,и' А,и'

1+—— = 1 + ——

и и, у

. Аналогично,

и, - А \и\ л А, и А и' ^ I 'I *1 1= 1 -_ил + _,_ + о

( _ _ 2 Л

А, и:

_ и, _

V у

Подставив эти выражения в (10), получим:

Аи = и,

Аи-АУ_+о

\и I и,

А и'

2 Л’

-А и'

1 + 0

2

А и'

(11)

1 ди

Теперь представим А,и' = — И-----------------+ О(И2 ). Если, кроме того, приращения А,-^^ \и\ и

2 дх ' ' '

и и А,+12 |и|

I I 1 д и / 2\

имеют одинаковые знаки, то А, и = — И^—1 + О(И ), где И — шаг сетки.

2 дх \ '

2 дх

_ 1 д I и| 1 ди

В тех областях решения, где знак и не меняется, имеем — - = —— и, следовательно,

\и\ дх и дх

О (И2 ). Выражение (11) преобразуется к следующему виду:

4

(И2)+1И * П + О(И2 )] = 1И^+о(И1)

V / 2 дх 1_ V /_ 2 дх V /

Аи = О. .

дх

Это и означает, что схема имеет второй порядок аппроксимации.

Если же в окрестности данной ячейки скорость меняет знак, то нарушается наше исходное предположение |А,и: □ |и,-1, |А, |и|| □ |и, |. Тогда необходимы иные оценки вместо (11). В этом

/л / \ к| + А,|и| |и,| -А,|и|

случае имеем и, = О (А,и) = О (И) и в наихудшей ситуации ^—1--------------||—--------,= О (1).

| и,- + А,и | | и,- — А,и |

Поэтому в районах смены знака скорости порядок аппроксимации схемы снижается до первого. Кроме того, он снижается до первого в тех областях, где знак и не меняется, но А,\и\ и

Аг+1/2 |и|

имеют разные знаки, и А, |и| = 1ишкег(а,-^ |и|; А,+^2 |и|) = 0. Заметим, однако, что

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

3. Результаты тестов. Рассмотрим теперь примеры применения нового алгоритма расчета приращений скорости на конкретных задачах. Данная методика проверялась на разнообразных задачах — от расчета обтекания профилей, клиньев, расчета взаимодействующих скачков до изучения сложных компоновок при числах Маха от 0.1 до 5 и выше.

Как было показано в предыдущем разделе, новая реализация ПМП отличается от классических реализаций (1) — (3) даже в одномерном случае. Поэтому вначале были рассмотрены стандартные тесты для одномерных течений невязкого газа — всевозможные варианты задачи Римана о распаде произвольного разрыва. Результаты, полученные по предлагаемой методике, совпали с высокой точностью с результатами, полученными по методу (2), примененному к скорости покомпонентно.

Далее был рассмотрен двумерный аналог задачи Римана — задача о взаимодействии двух косых скачков уплотнения. Эта задача представляет интерес, так как скачки уплотнения пересекают линии расчетной сетки. Расчетная область имеет форму прямоугольника размером 1х 0.5 и покрыта равномерной сеткой с квадратными ячейками (200 х100 ячеек). Сквозь левую

границу втекает горизонтальный равномерный поток газа с числом М1 = 2.5. Сквозь верхнюю границу под углом -10° втекает поток с параметрами за косым скачком уплотнения, в котором поток поворачивается на угол 10°. Сквозь нижнюю границу под углом +20° втекает поток с параметрами за косым скачком уплотнения, в котором поток поворачивается на угол 20°.

Расчеты были проведены с использованием метода (4) с функцией-ограничителем (2), примененного к скорости покомпонентно, и с использованием нового метода. Результаты расчетов представлены на рис. 2. На рис. 2, а, б изображены поля полного давления (отнесенного к параметрам набегающего потока) как поля, на которых особенно сильно проявляются ошибки расчета. Видно, что из левого нижнего и левого верхнего углов расчетной области, где в зоне поворота потока помещается лишь одна ячейка расчетной сетки, идут нефизичные «висячие» энтропийные слои. Также хорошо видны малые пространственные колебания полного давления, порождаемые переходом косого скачка с одного слоя ячеек на другой. Для сравнения методов рассмотрим распределения плотности (рис. 2, в) и полного давления (рис. 2, г) в сечении х = 0.4. На рисунке видно, что новый метод дает более быстрый рост параметров в скачках уплотнения и область размазывания скачков численной вязкостью имеет более резкие границы (значение плотности быстрее выходит на «полку»). При этом, однако, новый метод дает слабую немонотонность распределения плотности в окрестности скачков. Распределение полного давления в окрестности скачков и контактного разрыва оказывается немонотонным в обоих методах, но при этом новый метод порождает меньший «заброс» полного давления на контактном разрыве и больший — на верхнем скачке уплотнения. Нефизичные «висячие» энтропийные слои при использовании нового метода оказываются более узкими и менее

I I I I I I I I I I у М

о 01 ог о» о« о.5

Рис. 2. Поля полного давления р0/р01:

а — метод (4) + (2); б — данный метод (1 — рассматриваемое сечение); в — плотность р/рь г — полное давление р0/р01; ------------------------точное решение;-— метод (4) + (2),-— новый метод

интенсивными (см. рис. 2, г). В целом, точность методов оказывается близкой, но новый метод менее диссипативен.

Наиболее характерным примером уменьшения энтропийного слоя при использовании предложенного алгоритма является обтекание профиля КЛСЛ0012 на режиме М = 0.15, а = 22°. Расчетная сетка в данном примере имела 400 ячеек на профиле, 200 — до границы. Расстояние до границы расчетной области было равно 50 хордам профиля. Этот режим интересен тем, что при полностью дозвуковом обтекании профиля и низком числе Маха набегающего потока на носике профиля местное число Маха превышает 0.94. Теоретически сопротивление такого профиля при обтекании невязким сжимаемым потоком должно быть равно нулю.

Метод (4) с функцией-ограничителем (2), примененный к скорости покомпонентно, дает следующие результаты для аэродинамических коэффициентов:

сх = 0.0137, су = 2.5582.

Методика, изложенная выше, дает такие значения:

сх = 0.0074, с у = 2.5983.

При использовании в формуле (6) функции-ограничителя (3) значение сх еще сильнее приближается к нулю:

сх = 0.0045, су = 2.6177.

Подъемная сила при этом растет, что соответствует лучшему разрешению течения на носике профиля.

На рис. 3 приведены поля полного давления для рассмотренных случаев. Рис. 3, а соответствует случаю использования метода (4) с функцией-ограничителем (2), примененного к скорости покомпонентно; рис. 3, б представляет поле, полученное по новому методу с функцией-

ограничителем (2), рис. 3, в — поле, полученное по новому методу с функцией-ограничителем (3). Очевидно, что приграничный слой с ошибочным значением р0 существенно меньше в случае учета векторной природы скорости.

Для иллюстрации проблемы взаимодействия решения с неоднородностями расчетной сетки рассмотрим обтекание уступа с клином 5 = 26.6° при М = 3. На рис. 4, а, б, в приведены поля полного давления р0 для тех же вариантов ПМП, что и на рис. 3. Расчетная сетка имела ярко выраженный излом сеточных линий над точкой перехода поверхности клина в горизонтальное направление, так как горизонтальные сеточные линии были параллельны поверхности клина. Число ячеек на поверхности клина было равно 50. Ячейки имели стороны равной длины.

Рис. 3

Рис. 4

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

Кроме приведенных примеров, в ряде случаев, например при расчете обтекания клина 5 = 45° при М = 5 с отошедшим скачком, предлагаемый метод обеспечил устойчивость счета и глубокую сходимость решения, тогда как стандартные методы не позволили получить стационарное решение.

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

уменьшение сопротивления, связанного с потерями полного давления в пристеночных слоях из-за ошибок схемы;

уменьшение чувствительности сильных разрывов к изломам расчетной сетки; уменьшение размазывания скачков уплотнения;

улучшение сходимости решения и увеличение устойчивости счета на критических для схемы режимах;

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

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

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

ЛИТЕРАТУРА

1. Колган В. П. Применение принципа минимальных значений производной к построению конечно-разностной схемы для расчета разрывных решений газовой динамики //

Ученые записки ЦАгИ. — 1972. Т. III, № б.

2. Harten А. High resolution schemes for hyperbolic conservation laws // J. of Computational Physics. — 1983. Vol. 49, N 3.

3. S we by P. K. High resolution schemes using flux limiters for hyperbolic conservation laws // SIAM J. of Numerical Analysis. — 1984. Vol. 21, N 5.

4. Leer van B. Towards the ultimate conservative difference scheme. Part V: A second-order sequel to Godunov’s method // J. of Computational Physics. — 1979. Vol. 32, N 1.

5. Куликовский А. Г., Погорелов Н. В., Семенов А. Ю. Математические вопросы численного решения гиперболических систем уравнений. — М.: Физматлит. —

2001.

6. Bosniakov S., Fonov S., Jitenev V., Shenkin A., Vlasenko V., Yatskevich N. Method for noise suppressing nozzle calculation and first pesults of its implementation // J. of Propulsion and Power. — 1998. Vol. 14, N 1.

7. Тилляева Н. И. Обобщение модифицированной схемы С. К. Годунова на произвольные нерегулярные сетки // Ученые записки ЦАГИ. — 198б. Т. XVII, № 2.

Рукопись поступила 9/IV 2004 г.