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

Вестник Челябинского государственного университета. 2013. № 25 (316). Физика. Вып. 18. С. 57-63.

И. В. Березанский, В. С. Суров

УЗЛОВОЙ МЕТОД ХАРАКТЕРИСТИК ДЛЯ РАСЧЕТА ЗАДАЧ МНОГОЖИДКОСТНОЙ ГИРОДИНАМИКИ

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

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

1. Введение

В литературе описан ряд подходов, с использованием которых в переменных Эйлера решается задача определения положений контактных границ между различными идеальными сжимаемыми жидкостями. В одних межфазные границы моделируются невесомыми маркерами, перемещающимися в пространстве [1—2], в других для фиксации контактных границ используются методы VOF [3], level set [4], концентраций [5-6] . Для локализации контактных границ между идеальными газами применяется 7-модель [7] . В у-модели возможен расчет совместного движения не только идеальных газов, но также сред, описываемых с использованием более сложных уравнений состояния (см . [7]) . В [8] у-модель использовалась для сред с уравнениями состояния вида Ван-дер-Ваальса, в [9] — с уравнениями состояния Ми — Грюнайзена. В последнее время для локализации контактных границ часто используется а-модель [10-11] . Заметим, что применение а-модели ограничено в отличие от используемой в настоящей работе обобщенно-равновесной или ОР-модели, предложенной в [12], двумя сжимаемыми средами .

В случае применения ОР-модели задача ставится следующим образом [13]. Необходимо рассчитать совместное течение п различных идеальных сжимаемых жидкостей (газов), разделенных друг от друга контактными границами. При этом движение каждой из сред, которые занимают некоторые конечные не обязательно связные объемы V / = 1, ... , п), описывается уравнениями Эйлера . Для решения поставленной задачи необходимо, чтобы в начальный момент времени (? = 0) объемные доли фракций в смеси удовлетворяли равенствам а.. = 5.. (Ухе V), где а.. — объемная доля .-й фракции в/-м объеме; 5. — дельта Кронекера. Следует отметить, что в ОР-модели концентрации компонентов в смеси должны быть всегда строго больше нуля, поэтому в расчетах нулевые начальные концентрации компонентов заменяются очень малыми значениями порядка 10-8, что не снижает точности вычислений .

2. Обобщенно-равновесная модель смеси

Уравнения, описывающие совместное движение сред, следующие [14]:

где

—р + (u • V)p + рdiv u = 0, —^ + (u V)u + ^gradp = 0, —p + (u • V)p + pc2div u = 0;

—p0 —tt

-p- + (u • V)p0 + pG div u = 0, -—1 + (u-V)a,. + a,. (1 - Gt )div u = 0, i = 1, ..., n -1, (1)

G"=P0

( " -1

V^P? У

л

p-Pc2 v P? P ^ ,

Обозначения к (1)-(2) приведены в конце статьи.

Для одномерных течений систему (1) перепишем в векторно-матричной форме

где

А =

и =

и

Р

а1

Р0

а„

и-1 У

и -р 0 0 . .. 0 0"

0 и - 1 0 . . 0 0

р

0 -р с2 и 0 . . 0 0

0 -р№ 0 и . . 0 0

0 «,(6,-1) 0 0 . . 0 0

0 -Р 1&-! 0 0 . . и 0

0 ап-Л^п-1 _1) 0 0 . . 0

(3)

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

1 п

(4)

1 ^ 0 е = Р,е,'.

Р /=1

Как показано в [14], система уравнений (3) является гиперболической с корнями характеристического уравнения

X, = и - с, X = ... = X = и, X , = и + с.

1 ’2 2п ’ 2п+1

Если при описании поведения компонентов смеси использовать уравнения состояния вида

-,<Р.Р,")-- Р-С,<й° ~Р--) = ^ р( < -■-1) Р

где в. = 1/(т,- - ^ d i = с2.,в Ь. = <р.,- ^ р* с.. —

константы, индивидуализирующие -ю фракцию), то выражение для удельной внутренней энергии смеси (4) примет вид

й-1

X а. (рв* - +к )+рвп+к

1

£ = — р

1=1

а соотношения для О. и скорости звука (2) перепишутся как

Рс в - Р ^ + РВ1

с =

ж-1 а. (Ь В — Ь В )

О |Х ¡У п I I п '

В» + х—I---------п---

= Ь + РВ1

\Ьп + Р

1 + Вп —X

1а. (Ьш + РВп )

¡=1 Ь1+РВ1

где В = В - В , d = d - d , Ь = Ь - Ь .

т I п т I п т I п

Характеристические соотношения вдоль характеристических направлений dxldt = и ± с могут быть получены из уравнения

£-И о. 1 0 0 . . 0 .. с1р <1и -и р — Л Л

0 \—и 1 Р 0 . . 0 .. ¿и 1 с1р -и Л р Л

0 -рс2 %-м 0 . . 0 .. 2 йи с1р -рс и — Л А

0 -р& 0 %-м . . 0 .. ^и ^Р\ — р, Сг, и А А

0 04(3-1) 0 0 . . 0 .. „ ч А* с/вц -а, (1 - С,) и—- 1 1 А А

о : -Р1&-1 0 0 .. . £,-и .. р\в "_1 Я_1 А А

0 -1) 0 0 . . 0 .. с1и (1ап -а ,(1-и„ ,) м—- П~1К п~1’ Л Л

где ^ = dxldt.

После вычисления определителя имеем выражения

Ср ± рсСи = 0, (5)

справедливые вдоль характеристических направлений Сх/С = и ± с .

Вдоль траекторной характеристики Сх/С = и выполняются равенства

йр - е1 йр = 0, йа, -

Р/

а,. (1 - О,)

йр —‘-йр = 0, / = 1, ..., п -Р

1,

(6)

которые следуют из системы (3) .

При интегрировании системы уравнений (3) применялся узловой метод характеристик (УМХ), который позволяет с минимальной погрешностью рассчитывать задачи многожидкостной гидродинамики. Заметим, что при использовании других известных методов возможно «искажение» численного решения, что будет показано на примере 3, где дополнительно проводились расчеты с помощью метода Куранта — Изаксона — Риса (КИР) [15].

р() -р(хт) +рх)с«)

3. Узловой метод характеристик для смеси

УМХ относится к бессеточному типу. Для его описания достаточно рассмотреть способ определения значений искомых величин в узле (х *■ ) по известным данным в узлах на т-м временном слое . Для регулярных узлов решение поставленной задачи получим с использованием следующей итерационной процедуры [16-17]. На «нулевой» итерации (V = 0) полагаем, что значения искомых переменных в точке (хк, 1т+1) совпадают с данными из (хк, 1^, поэтому характеристические направления dxldt = и, dx/dt = и ± с аппроксимируются выражениями:

хк - хС = и'Ы, хк - х^ = (и + су)Л^ хк - х^ = иуМ, где Лt = tm+1 - tm. Из последних равенств определяем положения точек пересечения характеристик с прямой t = tm (рис . 1а):

х^ = хк - (и + cv)Лt, хС = хк - иуЛ^

х; = х - (и - су)м.

(7)

Параметры (р, и,р, р0, а1, ..., р0и-1, аи-1)(0) в найденных точках (х хс, х;)(0) находятся интерполя-

и (х]

,т+1 к

цией по их известным значениям в узлах х^ , хк и

= о,

|'+1 і'

) - и(хт)

/ т+1 \ / т\ / т\ / т\

р(хк ^ - р(хк Л -Р(хк )с(хк )

/ т+1\ / т\

и(хк ^ - и (хк )

р(хГ1) -Р(Хст) -с2(хт)

Р( хГ1 Ґ -Р( Хст )|'

= 0,

р( хт) р( хт )|'

а,. (хт+1)| -а,. (хС )|

_0 / т+1 \1'+^ Л0 / т\|'

Р (хк )\ -Р(хс )

-а, (хт)

1 - О ( хС )

= 0,

Р( хС+1) -р(хС)

= о,

-Р° (хт о (хт)

р(хС+1)Г -р(хС)|'

= о,

(8)

і = 1, ... , п - 1.

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

и

и

х • Соотношения (5)-(6) перепишем в конечноразностном виде как

Решая систему (8) при V = 0 относительно переменных (р, и, р, р01, а1, ..., р0и-1, аи-1)(1), найдем уточненные значения искомых функций в точке (х t ). Затем по этим данным из выражений (7) вычисляются новые координаты (х хС, хк)(1), которые в свою очередь используются для определения (р, и, р, р01, а1, ..., р0и-1, аи-1)(2) из (8), где необходимо положить V =1. Описанный итерационный процесс продолжается вплоть до сходимости.

В алгоритм УМХ нетрудно включить подвижные контактные границы . Будем считать, что в начальный и последующие моменты времени контактная граница располагается в узле хС На «нулевой» итерации (V = 0) полагаем, что значения искомых переменных в точке х^1 совпадают с данными из хС, поэтому характеристические направления dx/dt = и, dx/dt = и ± с аппроксимируются выражениями (рис . 1б)

(х£н)у =*£ +^(и(л£) + и(л£+1)|У)д*, (9)

х+1у=(хт у+[«х)+ф? )]у м,

(хС+1У =(хт У +\и(хт) - фт )]'а^, (10)

которые дают возможность найти (х^1)® и положение точек пересечения крайних характеристик с прямой t = t :

А т

х у=(хт+1у -[«х)+фт )]у х У = (хт+1У -[«(хС) - Фт)]" А/. (11)

Параметры в найденных точках, (хт, х^)(0) находятся интерполяцией по их известным значениям в узлах х хС и хС, х . Вдоль найденных

характеристических направлений выполняются соотношения

которые позволяют получить уточненные значения искомых величин (р, и, р, р1, а1, ..., р^, а )(1) в узле (х”+1)(1). Далее из формул (9) и (11) определяем новые положения контактной границы и точек пересечения крайних характеристик с прямой ґ = ґт. Затем из системы (12), где необходимо положить V = 1, вычисляем значения величин (р, и, р, р0, а1, ., р°п1, аи-1)(2) в точке (х^+1)(2) и т. д . Описанный итерационный процесс продолжается вплоть до сходимости

4. Результаты вычислений

Для иллюстрации применения описанного выше численного метода, используемого для локализации контактных границ, рассмотрено несколько тестовых примеров взятых из [18-19].

Пример 1. Рассмотрена задача Римана при следующих значениях параметров на момент времени ґ = 0 [18]: «слева» (х < 5) от контактной границы, разделяющей различные среды, — (р, и, р, у. р,, с,^ = (50, 0, 7,87, 3, 7,87, 1); и «справа» от нее (х > 5) — (р, и, р, у, р,, с,); = (1, 0, 2, 2, 2, 1) .

На рис . 2 представлены данные численных расчетов течения, полученные с использованием УМХ на равномерной сетке из 400 ячеек в сравнении с точным решением к моменту времени ґ = 1,0 с .

Пример 2. Рассчитана задача Римана при следующих значениях параметров на момент времени ґ = 0 [18]: «слева» (х < 5) от контактной границы, разделяющей различные среды, — (р, и, р, у, р,, с,)) = (1, 0, 7,87, 1,4, 7,87, 1); «справа» от нее (х > 5) — (р, и, р, у, р,, с,); = (10, -2, 8,5, 3, 8,5, 1) .

На рис 3 приведены данные численных расчетов течения, полученные с использованием УМХ на равномерной сетке из 400 ячеек в сравнении с точным решением к моменту времени ґ = 1,0 с .

Р( ХТ 1) - Р(х?) + Р(х? )с(х?)

Р(ХҐ) - Р(х?) -Р(х?)с(х?)

и(х?+1) - и(х?)

и(х?+1)Г - и(х? )Г

= 0, = 0,

і ІУ+1 Ш+\\ Ш\

Р(хс -Р(хс)

= 0,

Р( х?) Р( х?)

а, (х?+1) -а, (х?)

л0 / ?+1 ч ^0/ ш\

Р, (хс П -Р, (хс )

■а, (х?) [1 - С, (х?) ]

Р(х?+1) -Р( х? )

-Р,° (х? )0, (х?)

Р(х?+1) -Р( х? )

= 0,

= 0, і = 1, .

(12)

Рис . 2 . Зависимости параметров течения к моменту времени t = 1,0 с, полученные с помощью УМХ (штриховые кривые); точное решение — сплошные кривые

20

15

10

5

0

0

20

15

10

9 х

0

9 х

Рис . 3. Зависимости параметров течения к моменту времени t = 1. 0 с, полученные с помощью УМХ (штриховые кривые); точное решение — сплошные кривые

-10

-20

0,0 0,3 0,6 0,9 х 0,0 0,3

0,6

0,9

Рис . 4 . Зависимости параметров течения к моменту времени t = 0,3 мс, полученные с помощью УМХ (штриховые кривые); КИР (штрихпунктирные); точное решение — сплошные кривые

р р и

Рис . 5 . Зависимости параметров течения к моменту времени t = 0 . 28 мс, полученные с помощью УМХ (штриховые кривые); точное решение — сплошные кривые

Для следующего тестового примера расчеты проводились также с использованием схемы КИР (см. [15]):

тут+1 ТТт TTm ТТт

U i — U i + Am U i+i/2 ~ U j—1/2 _ о At A Ax _ ’

где

U = (p,u,p, p0, a,, p°2, a2)T,

u^./1=|(u;+u-„)+

+|^-'[Sign(A)]n|'(u;-u;1),

q-i- 1, i.

Здесь Л — диагональная матрица собственных значений матрицы Am, Q — матрица, строками которой являются левые собственные векторы матрицы A .

Пример 3. Рассмотрена задача Римана для трех различных газов при следующих значениях параметров на момент времени t = 0 [19]:

(p, u, р, Y, р*, c*) = (0,8 • 105, 0, 2,5, 1,2, 2,5, 0), х < 0,3,

(p, u, р, y, р*, с*) = (1,0 • 105, 0, 1,5, 1,4, 1,5, 0), 0,3 <х < 0,6,

(p, u, р, y, р*, с*) = (1,2 • 105, 0, 0,5, 1,67, 0,5, 0), х > 0,6 .

На рис . 4 представлены данные численных расчетов течения, полученные с использованием УМХ и метода КИР, к моменту времени t = 0,3 мс на равномерной сетке из 200 ячеек в сравнении с точным решением . Следует отметить при использовании метода КИР положения контактных разрывов и ударных скачков существенно разнятся с точными Не исправляет ситуацию и измельчение расчетной сетки

Пример 4. Рассмотрена задача Римана для трех различных газов при следующих значениях параметров на момент времени t = 0 [19]:

(p, u, р, y, р*, с*) = (1,8 • 105, 0, 1,5, 1,4, 1,5, 0), х < 0,3, (p, u, р, y, р*, с*) = (1,0-105, 0, 1,0, 1,2, 1,0, 0), 0,3<х<0,6, (p, u, р, y, р*, с*) = (0,2 • 105, 0, 0,15, 1,67, 0,15, 0), х>0,6,

На рис 5 представлены результаты численных расчетов течения, полученные с использованием УМХ, на равномерной сетке из 200 ячеек к моменту времени t = 0 28 мс в сравнении с точным решением

5. Выводы

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

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

Обозначения

с — скорость звука в смеси; n — количество фракций в смеси; p — давление; и — скорость смеси; х — пространственная переменная; a — объемная доля; е — удельная внутренняя энергия; р0 — истинная плотность i-й фракции; р — плотность смеси; р = а.р0 — приведенная плотность .

Нижние индексы L и R — для параметров смеси «слева» и «справа» от контактного разрыва (С); k — номер узла; m — номер временного слоя;

Верхние индексы v — номер итерации .

Список литературы

1. Glimm, J . Robust computational algorithms for dynamic interface tracking in three dimensions / J. Glimm, J. Grove, X Li et al . // SIAM J. Sci . Comput 2000. Vol . 21, № 6 . P 2240-2276.

2 . Суров, В . С . Взаимодействие ударных волн с каплями пузырьковой жидкости // Журн . техн. физики. 2001. Т. 71, № 6 . С . 17-22 .

3 . Scardovelli, R. Direct numerical simulation of free-surface and interfacial flow / R. Scardovelli, S . Zaleski // Annu Rev. Fluid Mech. 1999. Vol. 31. P 567-598.

4 . Osher, S . A level set method: An overview and some recent results / S Osher, R P Fedkiw // J Com-put. Phys . 2001. Vol . 169 . P 463-502.

5 . Бахрах, С . М . Расчет газодинамических течений на основе метода концентраций / С. М . Бахрах, Ю . П. Глаголева, М. С. Самигулин и др . // Докл. АН СССР 1981. Т 257, № 3 . С 566-569.

6 Бондаренко, Ю А Расчет термодинамических параметров смешанных ячеек в газовой динамике / Ю . А . Бондаренко, Ю . В . Янилкин // Мат. моделирование . 2002. Т. 15, № 6 . С . 63-81.

7 . Abgrall, R. Computations of compressible multifluids / R. Abgrall, S . Karni // J . Comput . Phys . 2001. Vol . 169 . P 594-623.

8 . Shyue, K. -M . A fluid-mixture type algorithm for compressible multicomponent flow with van der Waals equation of state // J. Comput. Phys . 1999. Vol . 156 . P 43-88 .

9 Shyue K -M A fluid-mixture type algorithm for compressible multicomponent flow with Mie — Gru-

neisen equation of state // J . Comput . Phys . 2001. Vol . 171. P. 678-707.

10 . Allaire, G . A five-equation model for the simulation of interfaces between compressible fluids / G. Allaire, S . Clerc, S . Kokh // J. Comput. Phys . 2002. Vol . 181. P. 577-616.

11. Murrone, A . A five-equation reduced model for compressible two phase flow problems / A Murrone, H Guillard // J. Comput. Phys . 2005. Vol . 202. P 664-698.

12 . Суров В . С . Течение Буземана для односкоростной модели гетерогенной среды // Инженер . -физ . журн. 2007. Т. 80, № 1. С . 75-84.

13 . Суров В . С. О локализации контактных поверхностей в многожидкостной гидродинамике // Инженер . -физ. журн. 2010 . Т. 83, № 3 . С. 518-527.

14 . Суров В . С . Об уравнениях односкоростной гетерогенной среды // Инженер -физ журн 2009 Т. 82, № 4 . С . 45-51.

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

16 Суров В С Об одном варианте метода характеристик для расчета течений односкоростной многокомпонентной смеси // Инженер -физ журн 2010 .Т. 83, № 2 . С . 345-350.

17 . Суров, В . С . Узловой метод характеристик для гиперболических систем / В С Суров, Е. Н . Степаненко, И . В . Березанский // Материалы XVIII междунар . конф . по вычислит, механике и современным прикладным программным системам, г. Алушта, 22-31 мая 2013. М . : Моск. авиац . ин-т. 2013. С . 668-671.

18 . Guoxi, N . A remapping-free, efficient Riemann-solvers based, ALE method for multi-material fluids with general EOS / N . Guoxi, J . Song and W Shuan-ghu // Computers & Fluids . 2013. Vol . 71. P 19-27.

19 . Caia, L . An efficient ghost fluid method for compressible multifluids in Lagrangian coordinate / L . Caia, J . -H . Feng W. -X . Xie // Applied Numerical Mathematics . 2008 . Vol . 58 . P 859-870.

20 Суров, В С Сеточный метод характеристик для расчета течений односкоростной многокомпонентной теплопроводной среды / В . С . Суров, Е . Н . Степаненко // Вестн . Челяб . гос . ун-та. 2010 . № 24 (205) . Физика. Вып. 8 . С . 15-22 .

21 Суров, В С Новые гиперболические модели в механике гетерогенных сред / В С Суров, И В Березанский // Механика композиционных материалов и конструкций : сб . тр . IV Всерос . симпозиума, Москва, 4-6 дек. 2012, Т. 2 . М . : Ин-т приклад, механики РАН, 2012 . С. 252-265 .