Научная статья на тему 'Сравнение схем с расщеплением потока для численного решения уравнений Эйлера сжимаемого газа'

Сравнение схем с расщеплением потока для численного решения уравнений Эйлера сжимаемого газа Текст научной статьи по специальности «Математика»

CC BY
298
69
i Надоели баннеры? Вы всегда можете отключить рекламу.
Ключевые слова
УРАВНЕНИЯ ЭЙЛЕРА / СХЕМЫ С РАСЩЕПЛЕНИЕМ ПОТОКА / СХЕМА STEGERWARMING / СХЕМА VAN LEER / СХЕМА AUSM+-UP / СХЕМА TORO-V´AZQUEZ / СХЕМА HLLC / ЗАДАЧА ДВОЙНОГО МАХОВСКОГО ОТРАЖЕНИЯ

Аннотация научной статьи по математике, автор научной работы — Смирнова Н.С.

Статья посвящена сравнению схем с расщеплением потока для уравнений Эйлера сжимаемого газа. Рассматриваются четыре схемы с расщеплением потока: StegerWarming, van Leer, AUSM+-up и Toro-V´azquez. Проводится серия вычислений на ряде одномерных тестов, а также задаче двойного Маховского отражения. Результаты расчетов позволяют сделать вывод о преимуществе методов van Leer и AUSM+-up над другими подходами c расщеплением потока.

i Надоели баннеры? Вы всегда можете отключить рекламу.
iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

Текст научной работы на тему «Сравнение схем с расщеплением потока для численного решения уравнений Эйлера сжимаемого газа»

УДК 519.63

Н. С. Смирнова

Московский физико-технический институт (государственный университет)

Сравнение схем с расщеплением потока для численного решения уравнений Эйлера сжимаемого газа

Статья посвящена сравнению схем с расщеплением потока для уравнений Эйлера сжимаемого газа. Рассматриваются четыре схемы с расщеплением потока: Steger-Warming, van Leer, AUSM+-up и Toro-Vázquez. Проводится серия вычислений на ряде одномерных тестов, а также задаче двойного Маховского отражения. Результаты расчетов позволяют сделать вывод о преимуществе методов van Leer и AUSM+-up над другими подходами с расщеплением потока.

Ключевые слова: уравнения Эйлера, схемы с расщеплением потока, схема Steger-Warming, схема. Vein Leer, схема AUSM+-up, схема Toro-Vázquez, схема HLLC, задача двойного Маховского отражения.

N. S. Smirnova Moscow Institute of Physics and Technology (State University)

Comparison of flux splitting schemes for numerical solution of the compressible Euler equations

In this paper, flux splitting schemes for the compressible Euler equations are studied. We consider four flux splitting schemes, viz. Steger-Warming, van Leer, AUSM+-up and Toro-Vázquez. A comparison of the schemes is conducted using one-dimensional tests and the double Mach reflection problem. The numerical results show the advantage of Vein Leer and AUSM+-up schemes over other methods with flux vector splitting.

Key words: Euler equations, flux vector splitting, Steger-Warming scheme, van Leer scheme, AUSM+-up scheme, Toro-Vázquez scheme, HLLC scheme, double Mach reflection.

1. Введение

В газовой динамике широко используется модель идеального газа [1], движение которого описывается с помощью уравнений Эйлера. Для получения численного решения уравнений Эйлера сжимаемого газа рассматриваются противопоточные схемы. Дискретизация уравнений на сетке осуществляется с учетом направления распространения возмущения, что позволяет учитывать физические свойства течения. Существуют два основных подхода для выявления направления потока: методы типа Годунова [2], [3] и методы с расщеплением потока [4].

Первый подход основывается на модели, в которой поток на границе между ячейками аппроксимируется с помощью точного или приближённого решения задачи о распаде произвольного разрыва. К схемам типа Годунова относятся HLL (Harten-Lax-van Leer) [5], HLLC (Harten-Lax-van Leer-Contact) [6] и др.

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

© Смирнова Н. С., 2018

(с) Федеральное государственное автономное образовательное учреждение высшего образования «Московский физико-технический институт (государственный университет)», 2018

частиц на движущиеся «назад» и «вперед» применяют декомпозицию вектора невязкого потока на несколько частей. В настоящее время этот подход активно применяется в аэродинамике, так как имеет ряд преимуществ. Он достаточно быстро определяет направление распространения возмущения и имеет сравнительно простое расширение для расчетов реальных газов. Кроме этого, методы расщепления потока хорошо применимы для неявных схем, которые при решении задач аэродинамики являются предпочтительнее явных методов. Но уменьшенная сложность таких схем по сравнению со схемами типа Годунова приводит к ухудшению разрешения разрывов, в особенности контактных волн. К схемам с расщеплением потока относятся схемы: Steger-Warming [7], van Leer f8], AUSM+-up (Advection-Upstream-Splitting-Method) [9], Toro-Vázquez [10], [11] и др.

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

Работа организована следующим образом: в п. 2 рассматриваются уравнения Эйлера и их дискретизация, в п. 3 описываются четыре схемы с расщеплением потока (Steger-Warming, van Leer, AUSM+-up, Toro-Vázquez), n. 4 посвящен результатам расчетов на ряде тестовых задач и сравнению данных схем.

2. Уравнения Эйлера

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

д и + д f + д g + д ^

dt дх ду дх

0.

(1)

Вектор и является вектором консервативных переменных, а Е, С и Н — физические потоки в X, У и Z направлениях соответственно, которые представляют собой функции зависимости от консервативных переменных и задаются таким образом:

U

( Р \ pu pv pw \ E J

F

/ pu \

pu2 + p puv puw V u(E + p) /

G

pv puv pv2 + p pvw V v(E + p) J

(

H

\

pw puw pvw pw2 + p V w(E + p) )

где р — плотность газа, и, V, и — скорости газа в направлениях X, У и Z, р — давление газа, аЕ — его полная энергия:

V 2

Е = р ■ (— + е), V2 = и2 + V2 + и)2,

где е — внутренняя энергия газа.

Для замыкания уравнений необходимо уравнение состояния газа, которое задает внутреннюю энергию. Для совершенного газа уравнение состояния имеет следующий вид:

e(p,p) =

p

p(Y -1):

здесь y — показатель адиабаты.

Для дискретизации уравнений Эйлера сжимаемого газа используется метод конечных объемов [4]. Этот метод требует расчета численного потока по нормали к рассматриваемой ячейке, элемента расчетной сетки. В данной работе проводились вычисления на структурированной декартовой сетке, поэтому далее будем рассматривать лишь поток F в Х-направлении. Для трехмерной декартовой сетки конечно-объемная схема с потоком в Х-направлении записывается в виде консервативной схемы:

Un+1 = Un - Д(F, +1 - F, 1), (2)

где численные потоки F,, i, F- i через границы ячейки i могут быть найдены с помощью

г+ 2 г 2 " схем типа Годунова или схем с расщеплением потока.

В представленной работе были получены численные решения уравнений Эйлера в одномерной и двумерной постановке. Из трехмерного случая можно легко перейти к двумерному и одномерному: в первом случае приняв скорость w равной нулю, а во втором случае обнулив скорости v и w в векторе пот ока F.

3. Схемы с расщеплением потока

В данном разделе рассмотрим четыре метода разбиения потока для расчета численных потоков в консервативной схеме (2): метод Steger-Warming [7], метод van Leer [8], метод AUSM+-up [9] и метод Toro-Vázquez [10], [11]. Описание первых двух схем также приведено в книге Е.Ф. Topo [4].

3.1. Метод Steger-Warming

Метод Steger-Warming [7] был развит как обобщение противопоточных методов на дозвуковые области течений газа. Данная схема расщепления основана на декомпозиции вектора потока на два вектора F+ и F v которых матрицы Якоби имеют неотрицательные и неположительные собственные значения соответственно.

Рассмотрим диагональную матрицу А, составленную из собственных чисел X, матрицы Якоби A = Щи вектора потока. Можно представить матрицу Л в виде суммы матрицы

Л+

у которой собственные числа X+ неотрицательные, и матрицы Л - с неположительными

Xi

Л = Л+ + Л ".

Стегер и Уорминг предложили взять собственные числа X± равными

X± = 1(X, ±\Хг\).

A

A = КЛК -1,

К — матрица из соответствующих собственных правых векторов. Далее, используя свойство однородности вектора потока, получаем расщепление газодинамического потока:

F = AU = КЛК - 1U = К(Л+ + Л - )К - 1U = (A+ + A-)U = F+ + F -.

Собственные числа матрицы Якоби вектора потока в системе уравнений Эйлера (1) таковы:

X1 = u — a, X2 = X3 = X4 = u, X5 = u + a.

Составляющие газодинамического потока определяются следующим образом:

/ X± + 2(y — 1)X± + X± \

(u — a)X± + 2(y — 1)uX± + (u + a)X± F± = P- vX± + 2(y — 1)vX± + vX±

27 wX± + 2(7 — 1)wX± + wX±

\(H — ua)X± + (7 — 1) V2X±± + (H + ua)X±)

где H — энтальпия, определяемая следующим выражением:

E + p u2 + v2 + w2 a2 H =-=----+

P 2 7 - 1

Запишем численный поток через грань i + 2 для схемы (2), применив идею построения противопоточных схем:

F+ 2 = F+(U?) + F"(Un+1). (3)

Численный поток F+i на границе между двумя ячейками i к i + 1 интерпретируется как

композиция компоненты потока F+, движущейся вперед из левой ячейки i, и компоненты потока F—i, идущей назад из правой ячейки i + 1.

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

3.2. Метод van Leer

Ван Лир [8] предложил сделать расщепление потока для уравнений Эйлера, как и в методе Steger-Warming, на два таких слагаемых, что матрица Якоби первого слагаемого обладает неотрицательными собственными числами, а матрица Якоби второго - неположительными собственными числами. Ван Лир сформулировал условия, которым должно удовлетворять данное расщепление:

1) Компоненты f + и f- потоков F± повторяют симметрию f относительно числа Маха, т.е f+(M) = ±f-(-M), если f(M) = ±f(-M).

2) Матрицы Якоби

i + dF+ Л 9F"

Л + — _ Л - — _

дИ ' дИ

должны быть непрерывными.

3) При дозвуковом течении расщепленные потоки вырождены, т.е А+ и А_ имеют нулевое собственное число.

4) Потоки Е±должны быть непрерывны:

Е+ = Е, если М > 1; Е" = Е, если М < 1.

5) Составляющие Е± — полиномы от числа Маха наименьшей степени.

Условия (1) и (2) были введены для устранения особенностей в звуковой точке и для получения более гладкого численного решения. Условие (4) гарантирует, что данное расщепление в сверхзвуковых областях сводится к стандартному противопоточному методу. Наконец, условие (5) позволяет однозначно найти расщепленние газодинамического потока.

Исходя из вышеизложенных требований, ван Лир построил расщепление для системы уравнений Эйлера (1):

F± = ± 1 pa(1 ± M)2

1

2a ( Y-Y ( 2

2a (y-1 m ± 1)

v w

Vy2-t (^ M ± 1)2 + 2 (v2 + w2))

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

здесь число Маха M =

Как и в схеме Steger-Warming, расчет численного потока на границе осуществляется с помощью схемы (3).

Ван Лир оценивает устойчивость данной схемы с помощью неравенства

C (, ■ + a) < 2Y + \M\(3 — Y) Ccfl = ДЖ(u + a) <-7+3-.

Заметим, что Ccfi = Ccfi(M). Возьмем показатель адиабаты для воздуха y = 1.4:

cm* = 1, \м\ = 1,

Cf и 0.636..., \м\ = 0.

M

для данного метода, чем для схемы Годунова или Steger-Warming, v которых Ccfi < 1. Следовательно, для получения устойчивого численного решения с помощью схемы van Leer будет требоваться больше времени.

3.3. Метод Аи8М+-ир

Основная идея семейства схем АШМ (Аёуесйоп-ирб^еат-БрШ^-МеШоё) заключается в представлении конвективной волны и акустической волны как два физически различных процесса. Невязкий поток в схеме АиБМ+-ир [9] расщепляется на конвективный поток и поток с давлением:

р = рс + р = т ф + Р,

где

т = ри, Ф = (1 и V ш Н )Т, Р = ( 0 р 0 0 0 )Т.

В зависимости от процесса первое или второе слагаемое играет ключевую роль. Например, при нагреве стенки в пограничном слое поток массы оказывает большее влияние, а поперек разрыва - давление.

Конвективная часть представляет собой произведение скалярной величины-массы на вектор Ф. Числовой поток на границе между двумя ячейками можно записать следующим образом:

^ т + \тI + т — \т\ т

1 =-— Ф+ +-— Ф- + р]^,

г+2 2 2

N = ( 0 nx ny nz 0

\T

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

Поток массы

Поток массы через границу между ячейками рассчитывается по противопоточному правилу:

Л/Г Г рь, м1/2 > 0,

т = М1/2Й1/Н рл, м^ > 0.

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

М1/2 = М + + М- + Мр,

M + и M выражаются аналогично формулам в методе van Leer с помощью полиномиальной функции от величины M:

M± = ( 2(M±\M|), \M| ^ 1,

M s ±4(M ± 1)2 ± 1 (M2 - 1)2, \M\ < 1.

Vn u • n

М =

а1/2 а1/2

Третье слагаемое Мр — это диффузионный член, который вводится для улучшения численных расчетов при низких и трансзвуковых значений числа Маха:

Мр = - К шах(1 - аМ2, 0) ,

/а Ра1/2

М 2 = МЬ + МЯ - = Рк + РЯ 2 ' Р 2 '

С помощью функции шах диффузионный член исчезает при больших числах Маха и влияет на расчет числа Маха на границе при М2 < 1. Область значений чисел Маха, в которой этот член появляется, регулируется с помощью параметра а, который принимает положительные значения меньше единицы. Коэффициент Кр можно варьировать в промежутке от нуля до единицы, функция /а определяется таким образом:

Д(М0) = Мо(2 - Мо),

где

М2 = шт(1, шах(М2, М^).

Мсо — «сг^-ой» параметр, обычно берется такого же порядка, как число Маха на бесконечности М^ (в невозмущенном пространстве). Мсо не должно обращаться в ноль, иначе это будет приводить к численной ошибке.

Скорость звука на границе между ячейками выражается так:

Й1/2 = шт(йь, ая),

а*2 а*2

шах(а*,У+У шах(а*, —V-)'

а*

числа Маха равняется единице:

а*2 = 2(7 — 1)

7 + 1

Потоковая часть с давлением

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

Р = Р +|а • Рь + Р~ |а • РЯ + Ри, здесь Р ± — полином пятой степени от числа Маха:

H.

P ± = ^ 1

1 (1 ± sign(M)), \M\ ^ 1,

4(M ± 1)2(2 ^ M) ± aM(M2 - 1)2, \M\ < 1,

з , . _ е2л г 3 3

а = !б<—4 + Г>/Д € — 4; 16].

Ри

Ри = —Ки • Р + • Р- • (рь + Ря) • (/аа1/2) • (V- — ¥+).

+

ными следующим величинам:

а = 1, Мсо = М^, Кр = 0.25, Ки = 0.75.

3.4. Метод Toro-Vázquez

Topo и Васкес-Сендон [10], [11] предложили метод расщепления, в котором поток F(U) представляется в виде суммы конвективной части A(U) и части с давлением P(U):

F(U)

( pu \

pu2 + p puv puw

\u( 1 pv2 + pe + p)J

A(U) + P(U),

где A(U) и P(U) задаются столбцами:

A(U) = u

p

pu pv pw

V2 pv V

uK(U), P(U)

0

p

0 0

\u(pe + p)J

Обратим внимание, что конвективная часть не содержит переменной давления. Все пере-р

Система уравнений Эйлера разбивается на две системы:

д^ + дх 0

—U — P

dt dx

0.

Рассмотрим последовательно систему с потоковой частью Р(и) и систему с потоковой частью А(и).

Система с потоковой частью Р(и)

Запишем данную систему в терминах физических переменных

д д д V + B(V)d V

dt

dx

0,

где

V = B

1

p 0 0 0 0 0

u 0 0 0 0 i p

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

v , B = 0 0 0 0 0

w 0 0 0 0 0

p 0 YP 0 0 u

Ai = ^(u - A) <\2 = Аз = A4 = 0 < A5 = ^(u + A),

где

A = Л ¡u2 + 4YP.

p

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

Далее рассмотрим локальную задачу Римана:

шУ + B(V) iV = 0,

' Vl = vn, Ж < 0, Уд = vn+1, x> 0.

V(x, 0) =

Собственным числам Л2,з,4 соответствует контактный разрыв, сопровождающийся скачком плотности со значения рь на ря• Слева и справа от контактной волны идут «нелинейные» волны. Будем искать и* и р* в области между левой и правой волнами. Решение получим путем аппроксимации с помощью инвариантов Римана:

„.* _ cruR-cLuL

ui+i =

- с^ ы - PL),

p¡+1 =

с rPl—С lPr i 1 Cr Cl

Cr-Cl

+ 1 CcRCL (uR - Ul^

где

Cl = Pl(Ul - Al), Cr = рд(Ur + Ar).

Этот способ плохо реализуется при наличии сильных ударных волн для задачи Римана в консервативных переменных и = (р, ри, ру, р'Ш, Е)т. Теперь можно найти численный поток давления

/

P

i+ i

0

\

Р*

i+1 0

0

V Y-1 U*+1 Р*+1

2 " 1 2 /

Конвективная система

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

A(U)

( ри \

2

ри

puv

puw

U puv У

uK(U), K(U)

р

ри pv fw

M pv2/

K

u

Васкес-Сендон предложили 2 алгоритма. Схема TV

Схема основана на представлении конвективного потока:

A(U)

и* i K,

i+1 '

u

i+1

конвективная скорость между узлами г и г + 1, которая находится из задачи

Римана для системы с потоком Р(И). Численный поток рассчитывается по противопоточ-ному правилу:

Ai+1 = и*+1

i

если и *, 1 > 0, i+ 2

Kn+1, если и*+1 ^ 0.

+ i+ 2

Схема TV-AWS

Topo и Васкес-Сендон представили модификацию схемы, предложенной в методе Zha-Bilgen [12] для конвективной системы.

A*+i = A+ + a7+v

где

A" = 1(1 - w)«?K?, A+ = 1(1 + u)u?Kn,

un

w = w«) = . г ==,

( *) Ф + (un)2,

е - малая положительная величина, например, е = 0.1.

Численное решение уравнений Эйлера с помощью метода Toro-Vázquez можно найти, используя консервативную схему (2), где численный газодинамический поток равен следующей сумме:

г+ 2 г+ 2 + г+1

В работе конвективный поток рассчитывался с помощью схемы ТУ.

4. Результаты расчетов

В данном разделе представлены численные решения одномерных и двумерных уравнений Эйлера для сжимаемого газа на примере двух одномерных тестов и задачи двойного маховского отражения.

Для повышения порядка точности результатов в двумерном случае применялась конечно-объемная схема WENO (Weighted Essentially Non-Oscillatory) [13], [14]. Обозначим схемы с потоком X как WTENO-X. Тогда, например, схема с потоком, вычисляемым с помощью метода HLLC, будет обозначаться как WTENO-HLLC.

Схемы WTENO-X используют выпуклую линейную комбинацию шаблонов из трех точек с адаптивными коэффициентами для блока реконструкции. Адаптивный шаблон строится таким образом, чтобы он не пересекал разрыв, иначе возникает эффект Гиббса, т.е. если какой-то шаблон содержит разрыв, то соответствующий весовой коэффициент должен быть близким к нулю. Данная схема WTENO-X имеет 5-й порядок по пространству и 3-й по времени.

4.1. Одномерные задачи

В представленной работе были получены численные решения одномерных уравнений Эйлера с помощью четырех схем с расщеплением потока: Steger-Warming, van Leer, AUSM+-up и Toro-Vázquez. Расчеты данных методов сравниваются с точным решением на примере двух тестовых задач, постановка которых взята из книги Е. Ф. Topo [4].

В качестве начальных данных рассматриваются две области с газом, разделенные перегородкой жо- Слева и справа от мембраны в момент времени t = 0 задаются векторы начальных значений физических переменных Vl = (pl,ul,Pl)t и Vr = (pR,ur,pr)T. В следующий момент времени перегородка хо убирается, и решается задача о распаде разрыва (задача Римана). В зависимости от тестовой задачи с течением времени формируются различные конфигурации разрывов.

Численные и точные решения были найдены в области 0 < х < 1 . Проводилось сравнение полученных результатов при разбиении отрезка [0;1] на M = 100 и M = 1000 ячеек. Число Куранта было взято равным 0.9 (Ccfl = 0.9) для всех методов, кроме van Leer, для него Ccfl = 0.6. Начальные значения физических величин в тестовых задачах приведены в табл. 1.

Т а б .л и ц а 1

Начальные данные для тестовых задач

Тест PL ul PL PR ur PR

1 1.0 0.75 1.0 0.125 0.0 0.1

2 5.99924 19.5975 460.894 5.99242 6.19633 46.0950

Тест 1

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

0.9 0.3 0.7 ,-т 0.6

3 0.5 0.4 0.3 0.2 0.1

О 0.1 0 2 0 3 0.4 0 5 0 6 0 7 0.8 00 1 0 0 1 0.2 0 3 0 4 0.5 0.6 0 7 0 8 0.9 1

Position Position

а) Сетка содержит M = 100 ячеек б) Сетка содержит M = 1000 ячеек

Рис. 1. Решение тестовой задачи 1 с помощью метода Steger-Warming (символ) и точного решения (линия) в момент времени t = 0.2 щи хо =0.3

а) Сетка содержит M = 100 ячеек

б) Сетка содержит M = 1000 ячеек

Рис. 2. Решение тестовой задачи 1 с помощью метода van Leer (символ) и точного решения (линия) в момент времени t = 0.2 щи хо = 0.3

а) Сетка содержит M = 100 ячеек

б) Сетка содержит M = 1000 ячеек

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Рис. 3. Решение тестовой задачи 1 с помощью метода А!18М+-ир (символ) и точного решения (линия) в момент времени £ = 0.2 щи хо =0.3

О 0.1 0 2 0 3 0.4 0 5 0 6 0 7 0.8 00 1 0 0 1 0.2 0 3 0 4 0.5 0.6 0 7 0 8 0.0 1

Position Position

а) Сетка содержит M = 100 ячеек б) Сетка содержит M = 1000 ячеек

Рис. 4. Решение тестовой задачи 1 с помощью метода Toro-Vázquez (символ) и точного решения (линия) в момент времени t = 0.2 щи жо =0.3

На рис. 1, 2, 3 и 4 представлены результаты решения тестовой задачи 1 для четырех схем. На рисунках проиллюстрирована зависимость плотности 1'аза от координаты пространственной оси. Все схемы обеспечивают хорошое разрешение ударной волны и хуже справляются с контактной волной. На рис. 1, 2 и 3 возникает провал в разреженной волне в связи с невыполнением энтропийших) неравенства в звуковой точке. Схема TV (рис. 4) позволяет избежать этой ошибки. При увеличении числа ячеек на координатной оси можно увидеть на 1'рафиках быструю сходимость численных результатов, полученных с помощью четырех методов, к точному решению.

Тест 2

Тест 2 используется для проверки точности и устойчивости численнохх) решения. В результате распада разрыва данной задачи возникают две ударные волны с контактной волной между ними.

а) Сетка содержит M = 100 ячеек

б) Сетка содержит M = 1000 ячеек

Рис. 5. Решение тестовой задачи 2 с помощью метода Stöger-Warming (символ) и точного решения (линия) в момент времени t = 0.035 щи xo = 0.4

а) Сетка содержит M = 100 ячеек

б) Сетка содержит M = 1000 ячеек

Рис. 6. Решение тестовой задачи 2 с помощью метода van Leer (символ) и точного решения (линия) в момент времени t = 0.035 при хо = 0.4

а) Сетка содержит M = 100 ячеек

б) Сетка содержит M = 1000 ячеек

Рис. 7. Решение тестовой задачи 2 с помощью метода А!18М+-ир (символ) и точного решения (линия) в момент времени £ = 0.035 щи хо = 0.4

а) Сетка содержит M = 100 ячеек

б) Сетка содержит M = 1000 ячеек

Рис. 8. Решение тестовой задачи 2 с помощью метода Того-Уагсцюг (символ) и точного решения (линия) в момент времени £ = 0.035 щи хо = 0.4

Результаты расчетов четырех схем, представленные на рис. 5, 6, 7 и 8, сравниваются с точным решением. На полученных графиках зависимости плотности газа от координаты сильная правая ударная волна хорошо прорабатывается во всех методах. Левая ударная волна лучше разрешена для схем van Leer (рис. 6) и AUSM+-up (рис. 7) в сравнении с методами Того-Vázquez (рис. 8) и Steger-Warming (рис. 5). Данные алгоритмы плохо выделяют контактные волны при заданном разбиении [0;1] на M = 100 ячеек. При уменьшении шага по координатной оси X численные результаты, полученные с помощью четырех способов расщепления потока, сходятся к точному решению.

4.2. Задача двойного маховского отражения

Задача двойного маховского отражения [14], [15] тестирует численные схемы на способность обнаружить и разрешить неустойчивые контактные разрывы, а также на склонность к проявлению эффекта «карбункула». Задача не имеет точного аналитического решения, лишь экспериментальные подтверждения возникновения сложной ударно-волновой структуры [16].

В данной работе задача формулируется следующим образом. Вычислительная область представляет собой прямоугольник размерами [0; 4] х [0; 1]. В начальный момент времени t = 0 сильная ударная толпа с числом Маха M = 10 движется направо под углом 60° к горизонтальной оси X из положения x = 6. Перед ударным фронтом давление и плотность газа задаются равными: p = 1 и р = 1.4. Расчет ведется с числом Куранта CFL = 0.45.

Программа, получающая численное решение задачи двойного маховского отражения, взята из статьи [14]. Программа была модифицирована добавлением численных расчетов потоков van Leer, WENO-AUSM+-up и WENO-TV.

На рис. 9 20 представлены результаты, полученные с помощью методов WENO-HLLC, WENO-VanLeer, WENO-AUSM+-up и WENO-TV в момент времени t = 0.2. Расчеты были произведены на последовательности сеток с размерами 1920 х 480, 3840 х 960 и 5760 х 1440. Схема Steger-Warming здесь не рассматривалась, так как в одномерном случае она показала результаты на тестовых задачах хуже остальных.

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

ком уплотнения, который становится слабее и исчезает при приближении к контактному разрыву.

Наибольшую трудность у численного метода вызывает точность разрешения треугольной области с контактными разрывами. Схемы ^УЕМО-НЬЬС, \УЕРТО-УапЬеег, А118М+-ир п \¥ЕМО-ТУ дают схожие результаты на сетке 1920 х 480 (рис. 9-12). Сравнивая полученные результаты на более подробной сетке 3840 х 960, можно отметить, что схема ^¥Е]ЧО-НЬЬС (рис. 13) наиболее точно прорабатывает контактные разрывы, а схема \¥ЕМ"0-ТУ (рис. 16) не справляется с неустойчивостью. При дальнейшем измельчении сетки схемы ^УЕМ^-УапЬеег (рис. 18) и \¥Е]ЧО-Аи8М+-ир (рис. 19) лучше разрешают треугольную область, чем ^УЕШ-НЬЬС (рис. 17) и \VENO-TV (рис. 20).

Данная задача в связи с наличием сильных ударных волн вызывает проявление в схемах «карбунку.л»-неустойчивости. При решении теста на сетке 3840 х 960 методом \VENO-TV (рис. 16) в ножке Маха явно виден нефизический излом, который и является проявлением эффекта «карбункула», что приводит к значительному искажению течения. При использовании схемы \¥Е]ЧТО-НЬЬС (рис. 13) не возникает ярко выраженного эфекта «карбункула», но наблюдается изгиб в ножке Маха. Схемы \¥Е]ЧО-УапЬеег (рис. 14) и \¥Е]ЧО-Аи8М+-ир (рис. 15) более устойчивы к данному эффекту. На подробной сетке 5760 х 1440 схема \УЕ]ЧО-НЬЬС (рис. 17) дает значительный изгиб ножки Маха, который ухудшает разрешение треугольной области с контактными разрывами. Схема '\УЕ]чЮ-А118М+-ир (рис. 19) начинает проявлять склонность к «карбункулу».

Стоит подчеркнуть, что контактные разрывы - физические неустойчивые структуры. В реальности такие завихрения возникают из-за наличия в среде вязкости и теплопроводности, т.е такие структуры можно получить при численном решении уравнений Навье-Стокса. К уравнениям Эйлера, обладающим нулевой физической вязкостью, можно перейти при устремлении числа Рейнольдса к бесконечности в уравнениях Навье-Стокса. При численном решении уравнений Эйлера в качестве вязкости выступает схемная вязкость. Методы расчета уравнений Навь^Стокса также обладают схемной вязкостью, но при малых числах Рейнольдса физическая вязкость подавляет схемную. При увеличении же числа Рейнольдса схемная вязкость начинает играть большую роль. Задача двойного маховсокого отражения позволяет оценить схемную вязкость. Методы с большой схемной вязкостью проявляют сильные диссипативные свойства и размывают неустойчивые структуры. На основе представленных результатов можно заключить, что подход ^АТ^ЖО-НЬЬС обладает меньшей схемной вязкостью, чем методы \¥Е]ЧО-УапЬеег, \¥Е]ЧО-Аи8М+-ир и \¥ЕШ-ТУ.

Рис. 9. Решение задачи двойного маховского отражения с помощью метода WENO-HLLC. Сетка 1920 х 480. 30 контуров в диапазоне от 2 до 22

1 1

|1п о: 2 6.82759 11.6552 16.4828 21.3103

Рис. 10. Решение задачи двойного маховского отражения с помощью метода "У^ЕМО-УапЬеег. Сетка 1920 х 480. 30 контуров в диапазоне от 2 до 22

Жо: 2 6.82759 11.6552 16.4828 21.3103

Рис. 11. Решение задачи двойного маховского отражения с помощью метода \УЕГЮ-Аи8М+-ир. Сетка 1920 х 480. 30 контуров в диапазоне от 2 до 22

|+ю: 2 6.82759 11.6552 16.4828 21.3103

Рис. 12. Решение задачи двойного маховского отражения с помощью метода "У^ЕМО-ТУ. Сетка 1920 х 480. 30 контуров в диапазоне от 2 до 22

1 1 1 1

1+1 о: 2 6.82759 11.6552 16.4828 21.3103

Рис. 13. Решение задачи двойного маховского отражения с помощью метода "У^ЕМО-НХЬС. Сетка 3840 х 960. 30 контуров в диапазоне от 2 до 22

Мпо: 2 6.82759 11.6552 16.4828 21.3103

Рис. 14. Решение задачи двойного маховского отражения с помощью метода "У^ЕМО-УапЬеег. Сетка 3840 х 960. 30 контуров в диапазоне от 2 до 22

I I I

Жо: 2 6.82759 11.6552 16.4828 21.3103

Рис. 15. Решение задачи двойного маховского отражения с помощью метода \УЕГЮ-Аи8М+-ир. Сетка 3840 х 960. 30 контуров в диапазоне от 2 до 22

I I

|+1о: 2 6.82759 11.6552 6.4828 21.3103

Рис. 16. Решение задачи двойного маховского отражения с помощью метода 'УУЕГЮ-ТУ. Сетка 3840 х 960. 30 контуров в диапазоне от 2 до 22

■ I I ■

гИо: 2 6.82759 11.6552 16.4828 21.3103

Рис. 17. Решение задачи двойного маховского отражения с помощью метода "У^ЕМО-НХЬС. Сетка 5760 х 1440. 30 контуров в диапазоне от 2 до 22

Рис. 18. Решение задачи двойного маховского отражения с помощью метода \УЕ1>ТО-УапЬеег. Сетка 5760 х 1440. 30 контуров в диапазоне от 2 до 22

Рис. 19. Решение задачи двойного маховского отражения с помощью метода \¥Е1МО-Аи8М+-ир. Сетка 5760 х 1440. 30 контуров в диапазоне от 2 до 22

rho: 2 6.82759 11.6552 16.4828 21.3103

Рис. 20. Решение задачи двойного маховского отражения с помощью метода "У^ЕМО-ТУ. Сетка 5760 х 1440. 30 контуров в диапазоне от 2 до 22

5. Заключение

В настоящее время схемы с расщеплением потока активно используются в аэродинамике. В данной работе проведено сравнение схем с расщеплением потока, разработанных в последнее время: Steger-Warming,van Leer, AUSM+-up, Toro-Vázquez, и хорошо себя зарекомендовавшей схемы типа Годунова HLLC, чтобы выявить наиболее эффективный метод расчета. Были получены численные решения уравнений Эйлера для сжимаемого газа на примере одномерных тестов и задачи двойного маховского отражения.

Расчетные исследования позволили сделать следующие выводы. Схема Steger-Warming наименее предпочтительна для использования в задачах газовой динамики. Метод Toro-Vázquez хорошо применим в одномерных задачах, но проявляет «карбункул»-неустойчивость и не разрешает неустойчивые структуры в задаче двойного маховского отражения. В двумерной задаче схема HLLC показывает хорошее разрешение контактных разрывов, но имеет тенденцию к образованию «карбункула». Схемы van Leer и AUSM+-up дают лучшие результаты при моделировании сильных ударных структур.

Литература

1. Эглит М.Э. Лекции по основам механики сплошных сред. М.: Книжный дом «ЛИБРОКОМ», 2010.

2. Годунов С. К. Разностный метод численного расчета разрывных решений уравнений гидродинамики // Математический сборник. 1959. Т. 47, №3. С. 271-306.

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

4. Toro E.F. Riemann Solvers and Numerical Methods for Fluid Dynamics. Springer, 2009.

5. Harten A., Lax P.D., van Leer B. On upstream differencing and Godunov-tvpe schemes for hyperbolic conservation laws // SIAM Review. 1983. V. 25, N 1. P. 35-61.

6. Toro E.F., Spruce M., Speares W. Restoration of the Contact Surface in the Harten-Lax-van Leer Riemann Solver // Journal of Shock Waves. 1994. V. 4, N 1. P. 25-34.

7. Steger J.L., Warming R.F. Flux Vector Splitting of the Inviscid Gasdvnamic Equations with Applications to Finite Difference Methods //Journal of computational physics. 1981. V. 40, N 2. P. 263-293.

8. Van Leer B. Flux-Vector Splitting for the Euler Equations // Eighth international conference on numerical methods in fluid dynamics. 1982. P. 507-512.

9. Liou M.-S. A sequel to AUSM, Part II: AUSM+-up for all speeds // Journal of Computational Physics. 2006. V. 214, N 1. P. 137-170.

10. Toro E.F., Vázquez- Cendón M.E. Flux splitting schemes for the Euler equations // Computers and Fluids. 2012. V. 70. P. 1-12.

11. Toro E.F., Castro С.А., Lee B.J. A novel numerical flux for the 3D Euler equations with general equation of state // Journal of Computational Physics. 2015. V. 303. P. 80-94.

12. Zha G.-C., Bilgen E. Numerical solutions of Euler equations by using a new flux vector splitting scheme.//International Journal for Numerical Methods in Fluids. 1993. V. 17, N 2. P. 115-144.

13. Jiang G.S., Shu C.W. Efficient implementation of weighted ENO schemes // Journal of Computational Physics. 1996. V. 126, N 1. P. 202-228.

14. Titarev V.A., Toro E.F. Finite-volume WENO schemes for three-dimensional conservation laws 11 Journal of Computational Physics. 2004. V. 201, N 1. P. 238-260.

15. Woodward P.R., Colella P. The numerical simulation of two-dimensional fluid flow with strong shocks // Journal of Computational Physics. 1984. V. 54, N 1. P. 115-173.

16. Семенов A.H., Березкина M.K., Красовская И.В. Классификация разновидностей отражения ударной волны от клина. Часть 2. Экспериментальное и численное исследование разновидностей маховского отражения // Журнал технической физики. 2009. Т. 79, №4. С. 52-58.

References

1. Eglith M.E. Lectures on the fundamentals of continuum mechanics. M.: Book house «LIBROCOM», 2010.

2. Godunov S.K. A finite difference method for the computation of discountinuous solutions of the equation of fluid dynamics. Mathematical Sb. 1959. V. 47, N 3. P. 271-306.

3. Kulikovskii A.G., Pogorelov N.V., Semenov A.Yu. Mathematical Aspects of Numerical Solution of Hyperbolic Systems. London: Chapman and Hall/CRC, 2001.

4. Toro E.F. Riemann Solvers and Numerical Methods for Fluid Dynamics. Springer, 2009.

5. Harten A., Lax P.D., van Leer B. On upstream differencing and Godunov-tvpe schemes for hyperbolic conservation laws. SIAM Review. 1983. V. 25, N 1. P. 35-61.

6. Toro E.F., Spruce M., Speares W. Restoration of the Contact Surface in the Harten-Lax-van Leer Riemann Solver. Journal of Shock Waves. 1994. V. 4, N 1. P. 25-34.

7. Steger J.L., Warming R.F. Flux Vector Splitting of the Inviscid Gasdvnamic Equations with Applications to Finite Difference Methods. Journal of computational physics. 1981. V. 40, N 2. P. 263-293.

8. Van Leer B. Flux-Vector Splitting for the Euler Equations. Eighth international conference on numerical methods in fluid dynamics. 1982. P. 507-512.

9. Liou M.-S. A sequel to AUSM, Part II: AUSM+-up for all speeds. Journal of Computational Physics. 2006. V. 214, N 1. P. 137-170.

10. Toro E.F., Vázquez-Cendón M.E. Flux splitting schemes for the Euler equations. Computers and Fluids. 2012. V. 70. P. 1-12.

11. Toro E.F., Castro C.A., Lee B.J. A novel numerical flux for the 3D Euler equations with general equation of state. Journal of Computational Physics. 2015. V. 303. P. 80-94.

12. Zha G.-C., Bilgen E. Numerical solutions of Euler equations by using a new flux vector splitting scheme. International Journal for Numerical Methods in Fluids. 1993. V. 17, N 2. P. 115-144.

13. Jiang G.S., Shu C.W. Efficient implementation of weighted ENO schemes. Journal of Computational Physics. 1996. V. 126, N 1. P. 202-228.

14. Titarev V.A., Toro E.F. Finite-volume WTENO schemes for three-dimensional conservation laws. Journal of Computational Physics. 2004. V. 201, N 1. P. 238-260.

15. Woodward P.R., Colella P. The numerical simulation of two-dimensional fluid flow with strong shocks. Journal of Computational Physics. 1984. V. 54, N 1. P. 115-173.

16. Semenov A.N., Berezkina M.K., Krasovskaya I.V. Classification of shock wave reflections from a wedge. Part 2: Experimental and numerical simulations of different types of Mach reflections. Technical Physics. 2009. V. 54, N 4. P. 497-503.

Llocmynujia e pedaK'n'um 04-12.2017

i Надоели баннеры? Вы всегда можете отключить рекламу.