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

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

Похожие темы научных работ по физике , автор научной работы — Матяш С. В.

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

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

УДК 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) порождает максимальное сглаживание неоднородностей

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

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

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

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

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

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

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

где введены обозначения: ? = (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.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

имеют разные знаки, и А, |и| = 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 у М

Рис. 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. Ячейки имели стороны равной длины.

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

Кроме приведенных примеров, в ряде случаев, например при расчете обтекания клина 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. Куликовский А. Г., Погорелов Н. В., Семенов А. Ю. Математические вопросы численного решения гиперболических систем уравнений. — М.: Физматлит. —

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.