Метод исследования диссипативных свойств разностных схем в эйлеровых координатах

Общепринятыми характеристиками разностных схем являются аппроксимация и устойчивость. Методы исследования этих свойств разностных схем подробно изучены в [1]. Однако в разностных схемах, описывающих сильные разрывы, должен присутствовать механизм диссипации энергии, обеспечивающий возрастание энтропии. В результате его действия может быть создано не только необходимое изменение энтропии, но и излишнее, обусловенное свойствами конкретной разностной схемы. Возникает вопрос, какое изменение энтропии можно считать допустимым. Анализ вклада аппроксимационной вязкости в энтропию для ряда разностных схем приведен в [2]. В качестве критерия диссипативности разностной схемы предложено использовать уравнение производства энтропии на слабой ударной волне [3, 4]

T o AS =

- 12 (d 2 ^p 2) _S ^{A V 3 + O (A V 4} • ^{AS 2 )}.

Таким образом, в зависимости от скорости роста энтропии на ударных волнах, схемы можно разделить на два класса: сильно диссипативные и слабо диссипативные. Физический смысл этого критерия весьма прост: разностная схема является приемлемой, если изменение энтропии из-за погрешностей аппроксимации не превосходит ее изменения в характерных физических процессах, какими являются слабые ударные волны. В сильно диссипативных разностных схемах слабые ударные волны неразличимы на фоне погрешностей. Иными словами, сильно диссипативные разностные схемы имеют низкую ≪разрешающую способность≫, и наоборот, слабо диссипативные РС – высокую [5]. Метод исследования диссипативных свойств разностных схем в лагранжевых координатах [6], заключается в построении уравнения производства энтропии, погрешность аппроксимации которого может быть выражена комбинацией погрешностей аппроксимации уравнений, составляющих разностную схему.

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

1.    Метод исследования

Запишем уравнение скорости изменения энтропии

TdS  de  pdV

dt    dt dt , где T – температура, S – энтропия, ε – удельная внутренняя энергия, P – давление, V – удельный объем. Поскольку вдоль траектории каждой частицы сохраняются еҷ лагранжевы координаты, то полная производная какой-либо величины совпадает с еҷ частной производной в этих координатах. В соответствии с [6] погрешность аппроксимации уравнения производства энтропии имеет вид

∂S

^!. =■ ■

T 2 (d ² P\ (dV V

- 12 [Vv^) _s\-at) ⁺ -

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

∂V ∂U     ∂U ∂P ∂t ∂m ω 1 , ∂t ∂m ω 2 , d~ (e + 0, 5U 2 ) +     (PU ) — “ 3 . ∂t           ∂m

Тогда погрешность аппроксимации уравнения производства энтропии ω _s имеет вид

∂S Tdf — “s — P“1 + “3 — U “2,                      (3) где U – скорость.

Запишем систему уравнений в эйлеровых переменных

∂ρ ∂ρU dt + dx — “1-                                (4) dpU + «eUU +dp—“2,                      (5) ∂t     ∂x    ∂x dpE dpEU dPU

∂t ∂x ∂x где E = £ + 0, 5U2 — удельная полная энегия.

Покажем, что для системы (4) – (6) выражение (3) имеет тот же вид. Для уравнения неразрывность проведем преобразования, доказывающие возможность перехода между двумя формами. Начнем с его записи в эйлеровых координатах dp dpU

∂t ∂x

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

% + u^ + ^ = 0.

∂t ∂x ∂x

Запишем первые два слагаемые как полную производную и сделаем замену р = V dp   dU    d V   1 dU п

^{+ р}я^— = ⁰,   ⁺     = ⁰

dt    ∂x      dt   V ∂x

После этого остается привести уравнение к требуемому виду путем несложных действий

1 dV £ dU =0

V ² dt V ∂x ^.

С учетом m = px, dm = рдх, уравнение (7) примет вид dV - ∂U dt ∂m

Окончательно получим вид уравнения в лагранжевых координатах

∂V _- ∂U

∂t ∂m

Более кратко изложим аналогичные преобразования уравнения движения и энергии. Рассмотрим уравнение движения в эйлеровых координатах dpU + dpUU + dP = 0 ∂t ∂x ∂x

После преобразований получим

, + U^ + u^ + pU⁹/ + ^a-P = 0. ∂t ∂t ∂x ∂x ∂x

Слагаемые 2 и 3 образуют уравнение неразрывности, домноженное на U. Разделим уравнение на ρ для выделения субстанциональной производной dU UdU 1 dP =0 ∂t ∂x ρ ∂x .

Окончательно получим

dU dP = 0

∂t ∂m .

Закончим рассмотрение эквивалентности уравнением энергии. Проведем аналогичные действия. Запишем его в эйлеровых координатах dpE + dpEU + dPU _ 0 ∂t ∂x ∂x pf + E% + E^U + Puf + dPU _0, ∂t ∂t ∂x ∂x ∂x dE dPU _ 0

∂t ∂m .

Таким образом, подстановка выражений (4) – (6) в (3) приводит к получению уравнения (1). Это показывает справедливость нахождения погрешности аппроксимации уравнения производства энтропии через комбинацию погрешностей исходных уравнений разностной схемы для методов в эйлеровых координатах.

2.    Исследование одной разностной схемы в эйлеровых координатах

Реализуем данный подход применительно к методу крупных частиц [7]. Для упрощения ограничимся рассмотрением одномерного случая, а также примем направления потоков слева направо с их вычислением с использованием формул первого порядка точности. Рассмотрим разностные уравнения неразрывности, движения и энергии pn+1 - p?   p? (Ui1 + U+1) - p^iUi- + Un)„ т +               2h’

?^п +1т pn+i u?    ?^птт^п(тт^п frⁿ ?ⁿ nⁿ (fiⁿ    nⁿ} Pⁿ Pⁿ

Pi Ui - Pi Ui  , Pi ui (Ui + ui+i) — pi-iUi-i(Ui-i + ui )  , Pi+i - Pi-i т                             2h2h

P?+1E"+1 — p? E?   p, (Ui1 + U;+i)Ei — pi-1 (U?1 + U? )E?-i т +                2h

+ (P i ¹ + P ?+i )(U ? + U ?+i ) — (P i-i + P^UH i + U n ) _₀

+                         4h'

Идея МКЧ подразумевает использование предварительных значений скорости U и энергии E .

U i ⁿ _ U i ⁿ

-

(P ?, 2p ? h ^{i +1}

— P i-1 ),

E i ⁿ

_ E i ⁿ

-

(P n +1 / 2 U i +1 / 2

-

nn τ

^P i-1/2 ^U i-1/2 ) ?_nh = ρ _i

_ E i — (₄(P ? + P i+X U i ¹ + U n+1 ) — ₄(P -1 + P^U n-i + U ? )) p ? h.

Отметим, что в настоящей работе рассматривается метод крупных частиц без использования дополнительной искусственной вязкости, поэтому значения P i+ 1 и U n _{+ 1} вычисляются как ^{2 2}

n

^P i+ 2

U + 2

1,

_2(P ⁿ + P i+1 ), _₂(U n + U n +1 ).

Теперь приступим к разложению разностных уравнений в ряды Тейлора в точке (i, n) и выделению погрешностей аппроксимации. Для упрощения данного процесса, начнем с уравнений предварительной скорости и энергии. Ограничимся слагаемыми первого порядка малости по времени и пространству иn = Un --7T + O(t2h2), i i ρin ∂x

En = Etn — /(4Pn^h + 4UП"Ph + O(t2, h2)- i i    4hpny i dx i dx J v 7

Вместе с этим потребуются первые и вторые производные предварительных величин по пространству. В результате получим

^ = ^ - I^P +      + O(T 2 ,^)

dx dx p dx2   p2 dx dx d2U d2U т d3P т dpd2P т d2pdP т dpd2P 2т dP dp dp o( 2 h2) dx2 dx2 p dx3 + p2 dx dx2 p2 dx2 dx p2 dx dx2 p3 dx dx dx T , ’

EE_ dE - 2 PPPU - Pd^U  PdUdp - Ud ² P  UdPdp _{O( 2} h _{2 )}

dx dx p dx dxT p dx2 T p2 dx dxT p dx2 T p2 dx dxT ’’ д2E = d2E  3 д2PdU   3 dP d2U 4 dpdUdP   Pd3UUd dx2   dx2   p dx2 dxT  p dx dx2 т ^ p2 dx dx dxT   p dx3 тp dx

Pd2 U   U d2 P   Pdpd2 U   Udpd2 P   P dU д2 p   UdPd2p p2 dx2    p2 dx2    p2 dx dx2    p2 dx dx2    p2 dx dx2    p2 dx dx2

PdUdpdp _ UdPdpdp 2 h2) p4 dx dx dxT p4 dx dx dxTT

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

p n+_ -p n + T	tX^U+X-pXXXi+U^ 2h	dp , 1 dP , 1 / nmfrn, ' dt+ 2 m T + 2h (p i (2U i +	∂U h + ∂x
1 d 2 U + 2 dx2	h 2 ) - (p n - dph + 1 ^h2)^’ - dx    2 dx 2	^w                         ^w dUh + 1 ?Uh 2 )) + O(T 3 ,h 3 ) dx    2 dx 2	=
=	9nn dUhM nrp      dU dp    nd d2p h^MCK^ h*\ 2h (2p д   ' 2 x' dxdxh   U i dx 2 h) + "

^ 1

Теперь заменим производные предварительной скорости по пространству ∂∂Ux на ранее полученное выражение (7). В результате получим

A+1 - Я    + йЫ-еЫй^^ = т                     2h dp , ndU , rmdp , 1 дР     d2P   1 dUdpu   1Tnd2pu      2 h2}

= 77 ' P^P ^{+ U}i^+ ⁺ ^^^{T - T}^2 ^- 77777^{h -} U^Ui ^{h + O(T} , ^h )

∂t i ∂x i ∂x

2 dt2      dx2    2 dx dx 2 i dx2          ’

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

„ 1 = - ¹ ^ т + т^ + ¹ ^h + ¹ U * &  + О(т ². h ² ). 2 dt ² dx ²   2 dx dx 2 dx ²

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

ρ n+1 _{U i} n+1 _-

τ

p * u *   P * U * (U * + U ) - P —Ur-A U — + U * )  P ^

+                2h

-

2h

n

= 0.

Разложение в ряд Тейлора в точке (i, n) будет иметь вид n∂U n∂ρ

^{P i} dt ⁺ U i dt

-

2u *

ⁱ ∂x ²

₊ 1 uⁿ—

+ 2 i dt ² ∂U ∂P

-

∂x ∂x τ

. 4 1 * d^ dUdp       dU dp т + 2P" a*т + •... + 2' ^ dx + u"U dx

-

1 * dUdU,   1 *_TT* d ² U,

2 ^{P i} dx dx ^h 2 ^Pi U i dx²

-

l _и * dUdp.

2 i dx dx

-

- ¹ U * U * 1 2 4 h + ^ + О(т ², h ² ) = 0. 2      dx ²     dx

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

„ 2 = -

1 Un ^d 2 ^p _T

2 U i dt² ‘

-

¹ P n

2 i^г

∂ ² U   ∂U ∂ρ

∂t ^{2 τ}

-

∂t ∂t

* d²P   dUdP 1 * dUdU,

^{T + 2U}" ~А~2^{т +} ТЕ "а^{-т + P}Pi ЕЕ ТЕ ' ' dx ²     dx dx 2 dx dx

+1 Pi U *S h + 3U * ^^h + 1U * U * ^Ph + О(тh²)

2 i dx²     2 i dx dx 2 i i dx²

Остается лишь провести вышеприведенные действия в отношении уравнения энергии. Рассмотрим его разностный вид

_{ρ i}n +1 _{E i}n +1 τ

-

ρ ⁿ E ⁿ

।

p , (и;* + и     e

^w

ww

- P i-1 (U - + U * )E — + 2h

^w

^w

^w

+

p + P + )(U? + u *+i) - ( p — + p *)^—

4h

-       = 0.

Запишем итоговое выражение для погрешности этого уравнения

1 „ d ² E ' "  - 2 ^P " д- . ^т

-

1 _En ^t

2 i dt ²

-

∂ρ∂E n∂P∂U n n∂2U     2∂2P dt dt + Ui dx dx + i Pi dx2 + U dx2 +

, dPdE + E * &P + P * dUdU + U * dPdU + 1 E * dUdP, + U * E * dP dx dx ^г dx^{2 г} dx dx ^г dx dx 2 i dx dx i x dx

+1 P * TUTEh ⁺ 1U E * ^h + 1 P * U *    h +

2 dx dx 2 i i dx²     2 i dx²

+¹ дГдГт + о(т ² h ² )-ρ ∂x ∂x

P ∂ ² P    P ∂ρ ∂P

ρ ∂x ^{2 τ}

-

ρ ² ∂x ∂x

т +

Теперь с помощью уравнения (3) найдем итоговую погрешность ω _s , объединив результаты для отдельных уравнений. Результатом является следующее выражение

X s = Pw i + ^ 3 - Ux 2 =

= _ 1 pn d2 2 i dt ∂ρ ∂E - ∂t ∂t + P n i ∂x ∂x P ∂ 2 P ρ ∂x 2 τ 1 — 2 Р П и

p . Pn d 2 PT+1 PndUdph +1 unpnd 2 ph _ 1 nnd 2 E _ i End 2 p 2   + P i d 2 x ‘ ' 2P i dx dxh 1 2 U i P i dx2     2 p i dt2     2 E i dt 2 n ∂P∂U     n n ∂ 2 U      n n ∂ 2 P   ∂P∂E     n ∂ 2 P + U i dx dx‘ 1 i P i dx 2 ‘      i U i dx 2 ' 1 dx dx ' 1 E i dx 2 + . 1 ™ dU dph. ,- , dph , 1 dU dE 1       d2p 1 d 2 E     (8) 2E dx dxh 1 U i Edxh 1 2 pdx dxh 1 2U i E i dx 2 h 1 2 p i U i dx 2 PdpdP 1 dPdP 1 „ J) 2 p 1    „dUj    dUUdp p 2 dxdx' 1 pdx dx' 1 2 U i U i dt2 ' 1 2 U i p i dt2 ' 1 U i dt dtT n JU^h - 1 p n U n U n ^h - 3 U n U n jUyph - 1 U n U n U n iph. i dx dx     2 i i i dx2     2 i i dx dx     2 i i i dx2

Для возможности сравнения полученной погрешности и допустимого физически обоснованного изменения энтропии на слабой ударной волне, рассмотрим преобразование выражения (2) к виду, содержащему только производные по времени. Для этого используются следующие выражения [4]

∂P ∂P

> = kdV J_s

∂V

∂t ,

d2P fd2P\ /dV\2  /dP\ d2V dt2 = \0V2) S \ dt ) + \dV) s dt2 ’

После ряда действий, рассматриваемое уравнение принимает вид xs =

-

т ² / d ² P\ / dV\3    _

12 VdV ²) Д dt) ⁺ - =

τ ² ∂ ² P ∂ρ

ρ ³ ∂t ² ∂t

-

т ² dP d ² p   т ² dP ( dp \ ²

12p ² dt dt² + 6p 3 dt у dt / +

Сравнение выражений (8) и (9) показывает превышение физически обоснованного изменения энтропии в методе крупных частиц. Это обусловлено наличием в итоговой погрешности слагаемых первого порядка по времени и пространству, в то время как выражение (9) ограничивается членами второго порядка. Следовательно, схема метода крупных частиц является сильно диссипативной на сильных разрывах.

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