Применение модифицированного полулагранжевого метода типа Годунова для моделирования разбега одномерных волн разрежения в сверхзвуковом диапазоне

Моделирование сжимаемых течений, и в частности, разбега одномерных волн разрежения в сверхзвуковом диапазоне, является сложной задачей из-за изменения скорости потока в широких пределах. При этом одна граница волны разрежения движется с очень большой скоростью, а вторая – с малой скоростью. Было предложено множество методов типа Эйлера с высоким разрешением для точного определения параметров плавного потока, например: метод искусственной вязкости [1, 2], метод уменьшения полной вариации (TVD) [3], существенно бесколебательная схема (ENO) [4] и взвешенная существенно бесколебательная схема (WENO) [5]. Одним из недостатков методов [15] является неточное определение профиля скорости при моделировании разбега волн разрежения при больших скоростях.

В работе [6] автором данной статьи был предложен метод, использующий подход Лагранжа к рассмотрению движения сжимаемой сплошной среды и фиксированную однородную сетку. Позже в [7] была описана реализация граничных условий для твердой стенки. В этих работах было показано, что у данного метода отсутствует численная вязкость на ударных волнах. В статье автора [8] предлагалась модификация этого метода, заключавшаяся в изменении формулы расчёта плотности сплошной среды для волн разрежения. Было установлено, что данная формула позволяет получать хорошие результаты для всех видов тестовых задач, решенных в [6, 7]. В то время, как в работе [7] для получения успешных результатов использо вались две разные формулы: линейный закон рас-

пределения параметров потока в волне разрежения для задачи Шу и Ошер об ударной трубе [9] и нелинейный закон из [6] для остальных рассмотренных задач. В работах [10, 11] описывается модификация метода, предложенного в [6], [7]. Нововведение заключается в учете дополнительного условия о направлении ударной волны на этапе расчёта подшага акустики. При этом данное новшество, внесенное в схему метода, не влияет на результаты уже рассмотренных в [7, 8] тестовых задач, а проявляется только при численном решении второй тестовой задачи с более сильной ударной волной, приведенной в работах [10, 11].

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

1. ОСНОВНЫЕ ПОЛОЖЕНИЯ И ИДЕЯ МЕТОДА

Уравнения Эйлера одномерной задачи движения сжимаемой сплошной среды в размерных переменных имеют вид

∂ρ   ∂ρ∂

+u +ρ=

∂t    ∂x∂

∂u∂u1∂p +u+=0,(1)

∂ t ∂ x ρ∂ x

∂ε ∂ε p ∂ u

+ u += 0 .

∂ t ∂ x ρ∂ x

Здесь u – скорость, ρ – плотность, ε – внутренняя удельная энергия, p – давление, x – координата, t – время.

Система уравнений (1) замыкается с помощью уравнения состояния для идеального газа p=(κ-1)ρε,

где κ – показатель адиабаты.

Как и в предыдущих статьях [6-8], [10, 11], в данной работе предлагается следующий метод: этапы конвекции и акустики рассматриваются отдельно, разбиение на физические процессы производится потому, что локальная акустическая скорость (скорость звука)

с = к— \ р может сильно отличаться от конвективной скорости потока. Таким образом, на акустическом этапе решается система

∂ ρ   ∂ u

+ ρ = 0 , ∂ t ∂ x

∂ u 1 ∂ p

+= 0

∂ t ρ ∂ x

∂ ε p ∂ u

+= 0

∂t  ρ∂x и уравнение (2). Система (3) получается из (1) отбрасыванием конвективных членов.

Если переменные перед частными производными заменить их постоянными значениями, мы можем получить из (3) характеристическую форму

^{∂ w 1} - c ^{∂ w 1} = 0 ,

∂ t

∂ x

^∂ w² = 0 , ∂ t

∂ w ∂ w

3+c 3=0, ∂t∂x где c * к -1

W1 =--; Р + и--ге, кр* c * кр*

w 2 = ^{р -} -2Г ^е,

c *

к - 1 е .

c *

w₃ = кр*

Здесь звездочка указывает на постоянные значения. Уравнения (4) описывают движение двух волн со скоростями - c и c и их соответствующими значениями переноса w 1 и w₃ (5).

Систему уравнений (4) можно решить, используя предложенный лагранжев подход. Однако в этой статье мы используем точный решатель Годунова [12] для решения полной нелинейной системы (1), (2), а затем вычитаем конвективную скорость из решения, чтобы получить решение уравнения (3). Более подробное описание метода приводится, в частности, в работе [11].

Шаг по времени _Д t_c на этапе акустики выбирался согласно критерию Куранта:

h

Д t c = ^k c —,

cs где cs – максимальная акустическая скорость разрежения, а kc – целое число. Если, например, kc = 3 , тогда волна распространяется на три ячейки сетки на подшаге акустики.

Модификация предлагаемого в данной статье метода по сравнению с работами [10, 11] заключается в ведении условия, чтобы распространение решений w ₁ и w ₃ с местной «акустической» скоростью не обгоняло (не затирало) не только решение впереди идущей ударной волны, но и границу раздела разбегающихся волн разрежения, если таковые существуют. При этом для соседних ячеек i и i + 1 проверяется условие

( и <- 10 ^{- 4}) и ( и_м >  10 ^{- 4}) .

На этапе конвекции решаем систему только с конвективными членами д Р   д р

—+ и — = 0 ,

д t

д x

д и    д и

+ и = 0 ,

д и

д t

д x

де

де   де

— + и— = 0 ,

д t     д x

также как систему уравнений (4).

Шаг по времени на стадии конвекции выбирался аналогичным (6) образом: ^h

Д ^t u k u      ,

us где us – максимальная конвективная скорость волны разрежения, а ku – тоже целое число. Например, если ки = 2, тогда волна распространяется на две ячейки сетки на подшаге конвекции. Если

_Д t_c ^ _Д t_u , тогда маршевый шаг _Д t выбирается, как

= п л

Д cc "с Д ,

Д tu nu Д t , где nc и nu – целые числа. Более подробно численная схема метода приводится в работах [6-8], [10, 11].

• 2.    ВЫЧИСЛИТЕЛЬНЫЙ ЭКСПЕРИМЕНТ

В этом разделе мы применяем предложенный метод для решения задачи разбега одномерных волн разрежения с числом Маха M = 2.2 , которая имеет следующие начальные условия:

p0 = 101325.0, и0 =-748.64, р = 1.225 ,если x < 0, j                        j                                                    j pj = 101325.0, uj = 748.64, р* = 1.225, если xj > 0,                   (9)

Область моделирования принималась равной x g [ - 50;50 ] , сетка содержала 1000 ячеек, шаг сетки составлял h = 0.1 . Шаг по времени для процесса «акустики» выбирался согласно критерию Куранта-Фридриха-Леви (6), где k c = 11, т.е. скоростная часть волны разрежения на этапе «акустики» перемещалась со скоростью c s = 340.3 на 11 ячеек сетки, а самая медленная часть - со скоростью c = 190.56 на 6 ячеек сетки.

Шаг по времени для процесса конвекции выбирался по аналогичному правилу (8), где k u = 24 , т.е. скоростная часть волны разрежения на этапе конвекции перемещалась на 24 ячейки сетки со скоростью u s = 748.64. Маршевый шаг равнялся _д t = 0.0032. Этапы «акустики» и конвекции выполнялись на каждом шаге по времени, т.е. _д t_c = _д t_u = _д t .

Результаты расчётов для трёх методов: схемы TVD Хартена [3], метода предыдущей работы автора [11] и предложенного в данной работе, - в момент времени t = 0.032 показаны на рис. 1-5, Предложенная модификация полулагранжевого численного метода с условием (7) позволяет получить более точный метод в то время, как в методе, описанном в [11], более медленная граница волны разрежения распространяется с большей скоростью, чем точное решение. С помощью схемы TVD не удаётся получить центральную зону разрежения с постоянными минимальными параметрами течения (см. рис. 3).

На основании полученных результатов можно сделать вывод, что предлагаемая схема метода имеет преимущество по учёту качественных особенностей течения и точности решения по сравнению со схемой TVD Хартена [3] и методом предыдущей работы автора.

Рис. 2. Распреде л ение скоро с ти в задаче о разбеге о дномерны х волн разр е жения (9):

- точн о е решение, --- • ~ - решение предыдущ е й работы а в тора [11],

“ “ • - п редлагаем ы й метод,--- ' • - р е шение схе м ы TVD [3]

Рис. 1. Распредел е ние плотн о сти в задаче о разбеге одномерн ы х волн раз р ежения (9): - точн о е решение, --- • ~ - решение предыдущ е й работы а в тора [11],

“ “ • - п редлагаем ы й метод,--- ' • - р е шение схе м ы TVD [3]

Рис. 3. Распре д еление ск о рости (уве л иченный м а сштаб)

в задаче о р а збеге одномерных вол н разрежен и я (9):

- точн о е решение, ---- • - - решение предыдущ е й работы а в тора [11],

- пре д лагаемый м етод,---- • • - реш е ние схемы T VD [3]

Рис. 4. Распределение давле н ия в задаче о разбеге о дномерны х волн разр е жения (9): - точн о е решение, --- • ~ - решение предыдущ е й работы а в тора [11],

“ “ • - предлагаемый метод, — ’ • - решение схемы TVD [3]