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

Актуальной проблемой при численном моделировании сверхзвуковых течений газа является корректное воспроизведение скачков уплотнения. В настоящее время для расчета таких течений широко используются методы, построенные на так называемых TVD схемах [1, 2], использующих подход Эйлера к рассмотрению движения сжимаемой сплошной среды. Одним из недостатков данных методов является большая численная вязкость, приводящая к размазыванию скачков уплотнения (фронтов ударных волн). Для решения этой проблемы Годунов [3] и др. предлагали использовать подвижные сетки, что существенно усложняет задачу, особенно при переходе к рассмотрению двумерных и трехмерных областей течения. В данной работе, как и в [4], предлагается метод, использующий подход Лагранжа к рассмотрению движения сжимаемой сплошной среды и фиксированную однородную сетку. Как будет показано ниже на тестовых задачах, у предлагаемого метода отсутствует численная диффузия на ударных волнах.

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

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

∂ρ ∂ρ ∂u

+u +ρ =0 ∂t∂x∂x,

∂ u ∂ u1 ∂ p

+u+=0 ∂t∂xρ∂x,

∂ε∂εp∂u

+u+=0 ∂t∂xρ∂x.

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

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

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

С.К. Годуновым [3] был предложен метод решения данных уравнений, заключающийся в подстановке в явную по времени центральную конечно-разностную схему решения задачи о распаде разрыва. В данной работе предлагается метод, использующий подход Лагранжа для рассмотрения движения жидкости. Сначала для каждой пары ячеек решается задача о распаде разрыва предложенным в [3] методом, в результате получаем значения “больших” переменных: R – плотность, U – скорость, E – внутренняя удельная энергия, P – давление, D – скорость распространения ударной волны или волны разрежения, – справа и слева от границы ячеек. Величины справа и слева от границы разрыва будем обозначать индексами R и L. При этом процессы конвекции и “акустики” рассматриваются отдельно, т.е. производится расщепление по физическим процессам, так как местная “акустическая” скорость может в разы отличаться от конвекционной скорости течения. Если шаг по времени для расчета процесса конвекции оказывается меньше шага по времени для процесса “акустики”, то конвекция не рассчитывается до первого расчёта этапа “акустики”, а за- тем выполняется столько расчетов этапа конвекции, сколько было пропущено.

2. СХЕМА МЕТОДА ДЛЯ ЭТАПА “АКУСТИКИ”

то для всех ячеек сетки, для которых выполняется условие x 1 <  X j <  x₂ , решение в следующий момент времени k будет тривиальным

P

k Rj

= 2(P , ^k - ¹

+ P

k - 1

Ri + 1

),

На этапе “акустики” в численном методе аналогично работе [4] рассматриваются две волны, бегущие влево и вправо, причём левая перено-

U_r ? = ¹ (U_r ? ^-1 + U _r? J?),

Rj 2 ^х Ri            Ri + 1 ⁷ ,

сит величины w 1 ={P L ,U L ,R L ,E L } с местной “акустической” скоростью C L = D L - U L , а пра-

R ^k = — (R ? ^{- 1} + R ^{k - 1} ) Rj 2 ^v Ri           Ri+1 * ,

вая - w 2 ={P,U,R,E} RRRR

со скоростью

CR = Dr - UR .

RRR

Предварительно для каждой пары ячеек i и i+1 с набором данных { p i ,u i , p i , s i } и { p_{i + 1},u_{i + 1}, p _{i + 1}, s i _{+ 1} } решается задача о распаде разрыва предложенным Годуновым в [3] методом,

в результате получаются значения больших пере-

менных {P, U,R,E,D,R ,E,D} , которые LLL RR R записываются в ячейки сетки следующим образом

^Ч = ^{P ,} U L" = ^{U ,} R^ = ^RL ^,

E l? ₊i = E l ^, C L?;¹ = D l - U ^, (3) они составляют вектор w 1 волны, бегущей в отрицательном направлении оси x (влево); и

p ^k-1 = P и ^{k -1} = U R ^k-¹ = R PRi ^{P , U}Ri ^{U ,} RRi RR ^,

F ^{k - 1} - F  C k ^{- 1} - D - U

E Ri      E R ^{, C}Ri D R ^{U ,}

E k ¹ /т^-1 k —1     t^-1 k—1\

Rj = 2^(E Ri ⁺ E Ri + 1 ) .            ⁽⁹⁾

Здесь e 1 = 10 ^{- 5} — допуск, связанный с погрешностью вычисления “больших” величин в итерационном методе Годунова.

• 2)    Если выполняется условие

• 2.1)    если

( ^CRi + 1 ^{- C}Ri ) <  ^e1 ,           ⁽¹⁰⁾

тогда

(X 2 - X 1 ) >  e₂ ^{,              (11)}

где e₂ = 10 - 8 - допуск, связанный погрешностью машинной арифметики, то имеем “слабую” ударную волну, которая не обгоняет решение из впереди стоящей ячейки, поэтому решение определится как

которые, в свою очередь, составляют вектор w 2 волны бегущей в положительном направлении оси x (вправо). В случае волны разрежения переменой D присваивается скорость более медленной характеристики волны разрежения. Например, для правой волны разрежения

D R = U + c_w - "\' ^: (u^ - U) ,   (5)

где

p k _ p k -1

^Г Rj    ^rRi ,

U_R ^k = U^{k - 1}, R_R ^k = R ^k-¹ Rj Ri ^, Rj Ri

E Rk = E r? - ¹ , если X 1 <  X j <  (X 1 + X 2 ) / 2 ,

p k _ p k -1

^r Rj    ^rRi + 1 ,

TT k —TT k-1  p k — p k -1

^URj    ^URi + 1 , R Rj    ^RRi + 1 ,

Рассмотрим волну, бегущую в положительном направлении оси x (вправо). Ячейки рассматриваются попарно. Координаты ячеек после перемещения на этапе акустики для “правой” волны определятся как x1 = xi + CRk-1 дtc ,

^X2 = ^Xi + 1 ⁺ О^ д t c .

Решения задачи на этапе “акустики” ищется

в следующем виде:

1) Если выполняются условия

, k-l_p ? -11

Ri + 1 P Ri | < ^e 1 ,

|^CRi + 1 ^CRi | <  e 1 , ⁽⁸⁾

E Rk = Erm , если ( ^X1 ^{+ X}2^)/2 < ^Xj <  ^X2 ; (12) 2.2) если

(X2 - X1) < e2,(13)

то решение определится как

P r‘ = P r‘ - 1 , U r‘ = U r‘ - ¹ , R r‘ = R r‘ - 1 ,

E r‘ = E Rk - ¹ , если ( U Rk - ¹ - U Rk - 1 ) < 0 и

(X1 -h/2) < Xj < (X1 + h/2),(14)

так как имеем ударную волну, и

P r? = P Rk + 1¹ , U r? = U Rk - ¹ , R Rk = R Rk;¹ ,

E r? = E r?;¹ , если ( U_R ?;1 - U r? ^- 1 ) > 0

и (x2 -h/2) < Xj < (x2 + h/2) ,(15)

так как имеем волну разрежения.

• 3)    Если ( ^CR? ₊₁ 1 ^{- C}r? ¹ ) >  ^e1 , то /тт ?-1 TT ?-11

• a) если I U_{Ri + 1} - U_Ri I <  e 1 , то имеем ударную волну

P_R k = P^ - 1 , U_R k = U^- 1 , R_R k = R_R k - 1 Rj Ri ^, Rj Ri ^, Rj Ri

ERk = Er^-1, если (Xi -h/2) < Xj < (x + h/2), p k _ p k-1 tj k _ и k1 R k - R k1

^Г Rj ^rRi + 1 , ^URj U RI + 1 , ^lvRj ^lvRi + 1

E Rk = Erm , если (X i + h/2) <  X j <  X 2 . (16)

Данный вариант возможен в силу следующего отношения “акустической” скорости ударной волны к скорости звука невозмущенного течения перед её фронтом, которое будет иметь вид

Cr“ = Dr - U = P,, (K + 1)P + K-1

C r" c ,, RrN 2 k P ii 4 k ⁽¹⁷⁾

или

7     —

K-1

⁺ IT , (18)

K + 1- ( k- 1) R R

( k + 1) Eil - ( k- 1)

D r - U ₌_P ii_ c_II      R_R

( k+ 1) K

ческой” скоростью не обгоняло (не затирало) решение впереди идущей ударной волны, если таковая существует.

Для волны, бегущей в отрицательном направлении оси x (влево) и переносящей величины w1 ={PL, UL,RL, EL} , условия и выражения, записанные выше, получаются аналогичным образом.

После перемещения левой и правой волн в каждой ячейке сетки мы имеем наборы величин kkkk    kkkk w1={PL, UL,RL, EL} и w2={PR, UR,RR, ER}, для которых снова решается задача о распаде разрыва предложенным в [3] методом. Причём в начальных данных задачи значения w2 считаются расположенными слева от границы разрыва, а значения w1 - справа. В результате снова получаем значения “больших” переменных {Pk*, Uk*,RLk*, ELk*,RRk*, ERk*}.

Так как давление и конвективная скорость не испытывают скачка в распаде разрыва, то их значения после подшага “акустики” определятся как где величины с индексом II относятся зоне перед фронтом ударной волны.

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

ше единицы, а потом неограниченно возрастать.

В частности, тестовый случай ударной вол ны, рассмотренный Рое [1], относится к вариан ту            /тт k-1 TT к

b) Если ( U_{Ri + 1} - U_Ri I ^ e 1 , то имеем вол ну разрежения и решение для X 1 <  X j <  X₂ оп

-

-

-

-

ределяется из условия сохранения инвариантов Римана

k

U Rj

= U к-1 U Ri

^{k - 1}_ТТ ^{к - 1}

Ri + 1     ^u Ri /

(xj

^X2 ^{- X}1

^{- X}1 )

^k E Rj

TE RT + ^EE - ^E^ (X j - X , ) ^X 2 ^{- X} 1                  J ’

R R ^k j

p ^{k - 1} R Ri

k

ERj к-1

E

^Ri J

+ Rr“

k

ERj к-1

^ERi + 1 _v

k- 1

/2

p ^k = P ik* , u ^k = U ^k* . (20)

Значения для плотности и энергии выбираются исходя из физического смысла решения. Если P i^k* >  p k ^{- 1} , то выбирается решение для ударной волны. В случае если имеем две ударные волны или пограничный случай, когда ударные волны отсутствуют, то решение находится следующим образом

k k* k k* k* k* k-1

Pi =RRi , £i =ERi , если RL >RR >pi k* k* или RR >pi >RL , k k* k k* k* k* k-1

Pi =RLi , £i =ELi , если RR >RL >pi или R,kk*>pk-1>RRk*.                    (21)

Если P^* < pk-1, то выбирается решение для волны разрежения. В случае если имеем две волны разряжения или пограничный случай, когда волны разряжения отсутствуют, то решение находится следующим образом k k* k k* k-1 k* k*

^p i =R Ri , ^£ i =E Ri , ^{если p} i >R R > ^RL k* k*

или R L > p i >R_r ,

k* k*                    k*

P l = ^RLi , £i =E ti , ^если P l >R l >R

k* k*

или R R > p i >R L

.

k*

Р *^к = ( k- 1)R R^k E R^k .

При выполнении условий 1) – 3), если оказывается, что в ячейку сетки на данном подшаге уже было записано решение, то оно заменяется “новым”, когда давление “нового” решения больше давления решения уже имеющегося в ячейке сетки. Кроме того, при записи решения в ячейки сетки также проверяется условие, чтобы распространение данных решений с местной “акусти-

3. СХЕМА МЕТОДА ДЛЯ ЭТАПА КОНВЕКЦИИ

В предлагаемом методе на этапе конвекции ячейки также рассматриваются попарно. Координаты ячеек после перемещения на этапе конвекции определяются как

^X1 = ^Xi ⁺ u i" ¹ Д t u , ^X2 = ^Xi + 1 ⁺ u i +i ¹ Д t u . (23)

Решения задачи на этапе конвекции отыскиваются в следующем виде:

1) если выполняются условия

I u i + 1 - u il < e 1 . lp i + 1 - pj < e 1 . l^S i + 1 -S il < e 1 .(24)

3) Если

^(ui + 1 - ^ui ) ^ ^e1 , ⁽³⁴⁾ имеем волну разрежения и решение для x₁ <  x j <  x₂ определяется из условия сохранения инвариантов Римана

то для всех ячеек сетки, для которых выполняется условие x₁ <  x j <  x₂ , решение в следующий момент времени k будет тривиальным

u k _j

= ur¹

u

к - 1 i + 1

-

u ^ki ¹

x 2 - x 1

⁽ x j - x 1 )

P j = 2(p + p_{i +} i ) , u j = 2(u + u_{i +} i )

£ к

к ¹ Z^k - 1 „к - lx „к ¹ /„к - 1 „к - 1

^p j = 2 ^{( p} i +^p i + 1 ) , ^s j = ^fo ^{+ s} i + 1 ) . (25)

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

• 2)    Если выполняется условие

⁽ u i + 1 ^- u i ) <  e 1 ,               ⁽²⁶⁾

то имеем ударную волну

• a)    если

(x ₂ - X 1 ) >  C 2 ,             ⁽²⁷⁾

то решение определится как

„ к „к - 1 „к „к - 1 к к - 1

Pj = pi , uj = ui , p j = pi , sk = £к-1, если x1 < xj < (x1 + x2)/2,

„ к „к - 1 „к „к - 1 к к - 1

^Pj = ^pi + 1 ’ u j = u i + 1 ’ ^p j = ^p i + 1

£ к = s k - ¹ , если (x 1 + x 2 )/2 <  x j <  x 2 ; (28)

• b)    если

(x 2 - x) <  e 2 ,             ⁽²⁹⁾

^то            кк- 1      к - 11

• - если p_{i + 1} - p i <  e₁ , сначала определяется граница разрыва

x = x i ■ h/2 + (u ^{k - 1} + и ^к +^- 1¹ ) _д t_c , (30) затем в ячейку j слева от границы разрыва записывается решение

„ к „к - 1 „к „к - 1 „к „к - 1 „к „к - 1

• ^pj = ^pi , u j = u i , ^p j = ^p i , ^S j = ^S i , ⁽³¹⁾ а в ячейку j+1 справа от границы разрыва записываются величины

(          /Тк-¹ /J^-1          ^ 2

TSk^-1+ ^{VSi + 1 Vsi} (x j - x 1 ) ,

( x 2 - x 1               J ,

p j^k

к - 1 p i

V

f _Fк ЛК- 1

^S j

, к - 1

^{+ p} i + 1

^s j

8 ^{к - 1}

s

1 У

К- 1

/2

У

p k = ( К- 1) p k £ к .

При выполнении условий 1)-3), если оказывается, что в ячейку сетки на данном подшаге уже было записано решение, то оно заменяется “новым”, когда давление “нового” решения больше давления решения уже имеющегося в ячейке сетки.

4. ТЕСТИРОВАНИЕ МЕТОДА

Предлагаемый метод тестировался на двух широко известных задачах, отличающихся начальными условиями (НУ).

Первая тестовая задача была в свое время предложена Рое [1] и имеет следующие НУ:

p0 = 1, u0 = 0, pk = 1, если xj < 0, p0 = 0.1, u0 = 0, p0 = 0.125 , если xj > 0. (36)

Область моделирования принималась равной x e [ - 4.5,5.5 ] , сетка содержала 100 ячеек, шаг сетки составлял h = 0.1 . Эти данные аналогичны параметрам моделирования задачи в работе [2].

Шаг по времени для процесса “акустики” выбирался согласно критерию Куранта-Фридри-ха-Леви

h

Д t c ^кс с ,

„ к „к -1 „к „к -1 „к „к -1 „к „к-1 _ZQQX ^pj = ^pi + 1 , u j = u i + 1 , ^p j =^p i + 1 , ^S j = ^S i + 1 , ⁽³²⁾ I к-1 „к -1k„

• - если p_{i + 1} - p_i    >  e₁ , решение определит

ся как к      к -1 к       к -1      к       к -1      к к-1

p j = p i , u j = u i , p j =p i , ^S j =^S i , если jV   >  p ^k _{+ 1} ¹ и (x 1 - h/2) <  x j <  (x + h/2) ,

_ к      к -1 „к    „к -1     к      к -1     к к-1

• ^pj = ^pi + 1 , u j = u i + 1 ,   ^p j = ^p i + 1 ,   ^S j = ^S i + 1 ,

если p ' <  p k-¹ и (x 2 - h/2) < x j < (x 2 + h/2) . (33)

где C = D — U - “акустическая” скорость распространения ударной волны; к_с = 1 , т.е. ударная волна на этапе “акустики” перемещалась на одну ячейку сетки.

Шаг по времени для процесса конвекции выбирался по аналогичному правилу

h

Д t u ^ки и ,

где U – конвективная скорость распространения ударной волны; k u = 2 , т.е. ударная волна на этапе конвекции перемещалась на две ячейки сетки.

Маршевый шаг равнялся _д t = 0.02020929 .

Этап “акустики” выполнялся через 6 шагов по времени, т.е. _д t_c = 6 _д t , этап конвекции - через 10 шагов _д t_u = 10 _д t .

Результаты расчетов после 110 маршевых шагов по времени показаны на рис. 1-3 в сравне-

Рис. 1. Распределение плотности в задаче о распаде разрыва (36)

– точное решение Годунова [3], – решение Хартена [2], – предлагаемый метод

Рис. 2. Распределение скорости в задаче о распаде разрыва (36):

– точное решение Годунова [3], – решение Хартена [2], – предлагаемый метод

Рис. 3. Распределение давления в задаче о распаде разрыва (36)

– точное решение Годунова [3], – решение Хартена [2], – предлагаемый метод нии с точным решением Годунова [3] и данными работы Хартена [2], полученными им для схемы ULT1C.

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

НУ второй тестовой задачи [2] задавались следующим образом p0 = 3.52773, u0 = 0.69888, pjk = 0.445, если Xj< 0, p0 = 0.571, u0 = 0, p0 = 0.5 , если xj > 0.                            (39)

Область моделирования принималась равной x e [-8,6], сетка содержала 140 ячеек, шаг сетки составлял h = 0.1 . Эти данные были также аналогичны параметрам моделирования задачи в работе [2]. Шаг по времени для процесса “акустики” выбирался согласно критерию (37), где kc = 2 , т.е. ударная волна на этапе “акустики” перемещалась на две ячейки сетки. Шаг по времени для процесса конвекции выбирался по правилу (38), где ku = 3 , т.е. ударная волна на этапе конвекции перемещалась на три ячейки сетки. Маршевый шаг равнялся д t = 0.2 , этапы “акустики” и конвекции выполнялись на каждом шаге по времени, т.е. д tc = д tu = д t.

Результаты расчетов в момент времени t = 2.0 показаны на рисунках 4-6 в сравнении с данными работы Хартена [2], полученными им для схемы ULT1C.

Графики показывают отсутствие численной диффузии на ударных волнах, в отличие от метода Хартена. Как и в первом примере, из-за накопления ошибки округления задний фронт ударной волны запаздывает в показанный момент времени на одну ячейку сетки.

В заключение можно сделать следующие выводы:

• 1)    Преимуществом предлагаемого метода по сравнению с методом Хартена является отсутствие численной вязкости (диффузии) на ударных волнах.

• 2)    С течением времени ударные волны или волны разрежения могут распространяться несколько медленнее или быстрее, чем в точном решении. Это происходит из-за округления положения фронтов волн с точностью до ячейки сетки вследствие применения фиксированной расчетной сетки.

  

  Рис. 4. Распределение плотности в задаче о распаде разрыва (39)

  – точное решение Годунова [3], о – решение Хартена [2], 4- – предлагаемый метод

  
  

  u

  
  
  

  Рис. 5. Распределение скорости в задаче о распаде разрыва (39)

  – точное решение Годунова [3], О – решение Хартена [2], 4- – предлагаемый метод

  
  
  

  Рис. 6. Распределение давления в задаче о распаде разрыва (39)

  – точное решение Годунова [3], о – решение Хартена [2], 4- – предлагаемый метод

  
  

  X