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

Вычислительные технологии

Том 19, № 4, 2014

К.. U о

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

характеристик

В. С. Суров, И. В. БЕРЕЗАНСКИй Южно-Уральский государственный университет, (Национальный исследовательский университет), Челябинск, Россия e-mail: svs@csu.ru, mynameivanych@gmail.com

Суров В.С., Березанский И.В. К расчёту течений односкоростной вязкой теплопроводной смеси узловым методом характеристик // Вычислительные технологии. 2014. Т. 19, № 4. С. 107-116.

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

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

Surov V.S., Berezansky I.V. The calculation of the flow of single-speed viscous and heat mixture by using the nodular method of characteristics // Computational Technologies. 2014. Vol. 19, No. 4. P. 107-116.

In this work we use single-speed model of viscous and heat heterogeneous media. For describing heat transfer we used Maxwell — Cattaneo law instead of Fourier law. The viscous stresses are described by means of a similar relaxation law. The account of viscosity and thermal conductivity is introduced for the equilibrium state of the whole mixture. Through using the relaxation laws for the source model is still hyperbolic, it allows us to use numerical methods for hyperbolic systems. We use nodular method of characteristics. This method proves to be effective during integration the equations for the single-speed model without considering the heat transfer and the viscosity. In this paper we describe a modification of the classical method of characteristics. We look for the parameters for a new time frame by means an iterative procedure. The calculated results are compared with results of single-speed model when the heat transfer and the viscosity are not considered.

Keywords: single-speed viscous and heat mixture, hyperbolic systems of partial differential equations, nodular method of characteristics, numerical modeling.

Введение

Односкоростные модели многокомпонентных сред используются при моделировании волновых процессов во вспененных жидкостях [1], в полимерах [2], пузырьковых жидкостях [3], для локализации контактных поверхностей в многожидкостной гидроди-

намике [4, 5], при разработке моделей кавитации [6] и в расчётах детонационных явлений [7]. В научной литературе помимо использованной в настоящей работе модели односкоростной смеси [8] представлены и другие гиперболические модели, описывающие течения бинарных смесей [9-13], которые, в отличие от [8], не распространяются на случай с произвольным числом фракций в смеси.

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

В работе [14] при рассмотрении распространения тепла вместо закона Фурье применялся закон Максвелла — Каттанео, учитывающий релаксацию теплового потока и обеспечивающий конечность скоростей перемещения тепловых волн, что в свою очередь связано с принадлежностью системы уравнений к гиперболическому типу.

Вязкие напряжения, как и тепловые эффекты, вводятся в уравнения модели многокомпонентной среды из [8] на уровне смеси в целом. В газовой динамике широко используется упрощённая формула для расчёта вязких напряжений [15], которая для одномерных течений имеет вид

ди

° = ^ (1)

где а — вязкое напряжение, ^ — коэффициент вязкости, и — скорость смеси, х — пространственная переменная. Однако применение этого соотношения приводит к потере гиперболичности исходной системы уравнений. Чтобы остаться в рамках данной системы вместо выражения (1) предлагается использовать соотношение

да да\ ди

т°{ т +идХ) + а = ^дХ' (2)

где Ъ — время, та — время релаксации вязких напряжений смеси. Уравнение (2) есть упрощённый вариант реологического выражения для жидкости Максвелла [16]. С этой точки зрения смесь в целом может рассматриваться как вязкоупругая среда, для которой время релаксации находится из соотношения та = р,/Е, где Е — модуль упругости смеси.

Узловой метод характеристик рассмотрим на примере бинарной смеси, состоящей из идеального газа с показателем адиабаты 7 и несжимаемой фракции с плотностью р3. Для общего случая смеси с произвольным числом компонентов методика расчёта остается той же. Представляемый вариант метода характеристик ранее использовался в расчётах течений, описываемых односкоростной моделью гетерогенной среды как в адиабатическом приближении [17], так и с учётом теплопроводности смеси [18]. Для бинарной смеси с объёмной долей идеального газа а система дифференциальных уравнений модели [8], в которой дополнительно учтены вязкие и теплопроводящие свойства смеси, принимает вид

др др ди ди ди 1 д(р — а)

— + и— + р— = 0, — + и— +-----= 0,

дЪ дх дх дЪ дх р дх

др + др + 2 ди + нд№ о {да + д^ + ди дъ и дх р°а дх дх ' П дъ 4 дх) ° ^ дх'

да да /л .ди дШ дШ . др . др .да Ш

"тгт + и---(1 — а)— = 0, + и-- + кр — + + ка — +--

дЪ дх дх дЪ дх дх дх дх тщ

где

'y(p - a) H _ Y - 1 k _ k _ AdT k _ X

, H , kp о , kp о , ka ¿л ,

ap a tw dp tw dp tw da

p — давление, T и W — осреднённые температура и плотность теплового потока, p _ aPg + (1 — a)ps — плотность смеси, х — коэффициент теплопроводности смеси, tw — время тепловой релаксации смеси, нижние индексы g и s означают газ и несжимаемая фракция.

Вывод последнего уравнения системы (3) имеется в [19].

При рассмотрении пузырьковой жидкости без потери точности можно считать, что температура несжимаемой фракции постоянна: Ts _ const, поскольку её доля в смеси значительна. Это предположение снимается при учёте сжимаемости жидкости. Выражение для средней температуры

T _ aTg + (1 — a)Ts

2

a2p

[p — ps(1 — a)]R

+ (1 — a)To,

(4)

где Т0 — начальная температура среды, Я — газовая постоянная. С использованием (4) коэффициенты кр, кр и ка принимают вид

2

k ___a2xp_

Р tw [p — ps(1 — a)]2R,

kp

Xa

Tw [p — ps(1 — a)]R'

ka —

x ( ap(2p + aps)

To

Ш \ [р - - а)]2Я

Вязкость и коэффициент теплопроводности смеси рассчитывались из выражений

V = - [аРаN + (1 - а)р8^з] , X = Р [аРдХд + (1 - а)Р5Хз] . Р р

Аналогично вычисляются коэффициенты релаксации смеси.

Характеристическое соотношение системы (3) определяется из уравнения

5)

£ — u —p 0 0 0 0

0 С — u — 1/p 1/p 0 0

0 —pcl С — u 0 0 —H

0 ^¡Ta 0 С — u 0 0

0 1—a 0 0 С — u 0

— kp 0 — kp 0 ka С — u

или (раскрывая определитель) —

(С — и)2 {(С — u)4 — (ca + Hkp + u2)(u — i)2 + H корни которого u, u, u ± c1, u ± c2, где

u)2kv — kp +

(1 — a)

ci

\

1 ca+w2+Hkp W(ci+w2)2+h

[Hkp + 2(ca + w2)] kp + 4(kp — w2kp —

1 — a

ka)

0

0

a

С2

Л

- ca+w2+Hkp (d+W2)2+h

[Hkp + 2(cl + w2)] kp + 4(kp - w2kp -

1 — a P

ka)

Здесь £ = —, ш2 =-.

т та р

Характеристические соотношения вдоль направлений dx/dt = п ± dx|dt = п ± с2 системы (3) могут быть получены из уравнения (6), в котором вместо последнего столбца в определителе должен быть использован вектор-столбец

/ dр dп \

—п—--р—

-п -р dtl Ма

—п—----г +---—

dt р dt р dt Мр 2 dп ттdW —п— — рс2— — Н-

dt ' а м <а

^а Ма ^ ^ Мп

та dt та dt Ма Мп

—п~АЪ + (1 — а)

W dW Мр Мр Ма

--п—--Кр—--кр—--ка —

МЬ dt dt dt /

\

Раскрывая этот определитель, имеем соотношение

АрМр + АиМп + АрМр + Ла Ма + АаМа + dW + = 0,

где

Ap = twraр£(п — О2 [(u — О2 — (cl + ш2)] к Au = Tw Ta p2C(u — С)3

к 2 к (1 —а) кр + cakp

Р

Ap = —TwTaрС(u — С)2 ^кр + [(u — С)2 — ш2] кр —

(1 — а)

Aa = Ta Tw РС (U — С)'

кр + cakp

Aa = —TwTaРС(П — С)2 [(u — С)2 — (cl + Ш2)] ка,

)\ кр,

ка ,

(1 — а)

Р

ка ,

)] ка,

ка ) ,

Aw = Tw Ta p(u — С )2S u(u — С ) [(u — С)2 — (cl + ш~ Л —

-H

кр + {(u — С)2 — ш2) kp — к

p 1ьа

Р

кр + C-'^p

2u (1 — а)

At = P^u — С )2 j Ta (u — OW [(u — С)2 — (cl + Ш2)] + Tw a Вдоль траекторной характеристики dx/dt = u справедливы равенства (1 — a)dp + pda = 0, ^dp + Ta pda + padt = 0,

Р

которые следуют из первого, четвёртого и пятого уравнений системы (3)

a

1. Узловой метод характеристик для гетерогенной среды

Для описания узлового метода характеристик достаточно рассмотреть способ определения значений искомых величин в к-м узле (хк, гп+1) по известным данным в узлах на п-м временном слое. Решение поставленной задачи получим с использованием следующей итерационной процедуры. На "нулевой" итерации (V = 0) полагаем, что значения искомых переменных в точке (хк, гп+1) совпадают с данными в точке (хк, гп), поэтому характеристические направления &х/дк = и, йх/(И = и ± с1 и йх/(И = и ± с2 аппроксимируются выражениями:

хк - хС = иДг, хк - х1г = (и" + с1 )Дг, хк - х12 = (ии + с2)Д1, хк - хиЕ1 = (и" - с1 )Дг, хк - хЪ2 = (и" - с%)Дг,

где Дг = гп+1 - 1п, V — номер итерации. Из последних равенств определяем положения точек пересечения характеристик с прямой I = 1п (рис. 1):

х^ = хк - (и + с)Дг, х12 = хк - (и + с2)Дг, хС = хк - иДг, хг1 = хк - (ии - с1 )Дг, хг2 = хк - (ии - с2)Дг. (10)

Параметры (р, и, р, а, а, Ш)(0) в найденных точках (хы1, хы2, хс, хг1, хг2)(0) находятся интерполяцией по их известным значениям в узлах хк-1, хк и хк+1. Соотношения (8)-(9) перепишем в конечноразностном виде как

АРы Ри+1 + + АрЬ1 р"+1 + Л1Ь1 а-+1 + Л^ аи+1 + Л^ Ш"+1 = Х^,

АрЬ2 Ри+1 + Л^ и"+1 + ЛрЬ2 р"+1 + ЛаЬ2 а-+1 + ЛаЬ2 аи+1 + Л^2 Ш"+1 = Х[2, ЛРС Ри+1 + и+1 + Л;ср"+1 + ЛиаС а»+1 + Л^ а^+1 + Л^с Ш "+1 = Х£,

Рис. 1. Расчётная схема УМХ для регулярных узлов

BVc PV+1 + BVc uv+1 + BVcpv+1 + BVc aV+ + BVc aV+1 + Bwc WV+1

YV

-V+1 _l AV „,V+1 _l ÂV ,rV+1 _l лг^1 _l ЛV a "+1 + WU+1 = X"

A;Rl pv+1 + A"uRi u "+1 + ä;Ri pV+1 + AvaRi aV+1 + A..

V

aRi1

Ri'

ä;r2 pv+1 + AUr2 uV+1 + ävpr2 pv+ + AIr2 aV+1 + ÄaR2 aV+1 + A^ WV+1 = XV

где коэффициенты, например, для характеристического направления Мх/М^ = п + С1 имеют вид

АРь1 = — Штарс1(п + С1) [с? — (са + и2)] К

P\Li

ÄVaLi = -TWТар c1(u + С1)

kp + cakp

(1 - а) (

P '

Li

ÄpL1 = -TWТаpcj(u + С1) ^kp + (c? - и2) kp - (1 р а ka^j

Li

ÄVaL1 = Та TW pcj(u + d)

k 2 k (1 - a) k kp + cakp k0

р

Li

ÄUaL1 = -TwTapc\(u + c1) [c? - (cl + UJ2)] kj"

Li

ÄWl1 = - Tw Ta fx?A uc1 [c? - (cl + и2)] - H kp + (c? - и2) kp - ——— k0

(Ät)Ll = fcl(u + c1 К Ta c1 W (cl + и2 - cl) + Tw a

2, (1 - a) 7

k -I- r2k — -_-k

р

Li

Li

Аналогичные формулы выписываются и для других характеристических направлений.

Решая систему (11) (при v = 0) относительно переменных (р, u, p, a, a, W)(1), найдем уточнённые значения искомых функций в точке (xk, tn+1 ). Затем по этим данным из выражений (10) вычисляются новые координаты (xli , Xl2 , xc, Xri , Xr2 )(1), которые в свою очередь используются для определения (р, u, p, a, a, W)(2) из (11), где необходимо принять v = 1. Описанный итерационный процесс продолжается вплоть до сходимости.

Расчёты по адиабатическому варианту модели смеси проводились с использованием численного метода Куранта — Изаксона — Риса (Courant—Isacson—Rees) [20]:

Uk+1 — Uk Uk _Uk

U^AIUk + (^^-<1/2 U±A-U +

Uk Uk

ui ui-1

Ax

um+1/2 = 1 (um - um+1) + 2 i^-1 [sisn (л)] ^im (um - um+o, m=г,г -1.

2

Здесь

U

р

u p

a

A

u 0 0

р u

р^а

00 1/р о

u0

л

±

2(л ±|Л|),

у 0 а - 1 0 u /

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

и

а

и

и

v

0

2. Результаты численного моделирования

С использованием описанных выше методов рассчитаны две задачи о распаде произвольного разрыва в пузырьковой жидкости. В первом случае параметры смеси до распада были следующими: слева от диафрагмы (х < 0) — р0Ы = 0.15 МПа, и0Ы = 0, а0ы = 0.12, рёоь = 1.19 кг/м3; справа от диафрагмы (х > 0) — р0я = 0.1 МПа, и0к = 0, а0к = 0.12, рёок = 1.19 кг/м3. Плотность жидкой фракции р8 = 1000 кг/м3,

■40

-0,00015 к

-20

20 АО х

\

а

0,136 0,128 0,120 0,112 0,104

-40 -20 0 20 40 х

40

-20

- \ г

0 20 40 х

40 х

Рис. 2. Зависимости р/р0я(х), и(х), Т(х), а(х), а(х), Ш(х) для пузырьковой жидкости к моменту времени Ь = 0.6 с, полученные при учёте неравновесных процессов (сплошные линии) и в адиабатическом приближении (штриховые линии)

вязкость = 10-3 кг/(м • с), коэффициент теплопроводности Хб = 60.2 • 10-2 кг^м/(с3^К), времена релаксации тав = = 10-1 с. Показатель адиабаты газа 7 = 1.4, его вязкость = 1.81 • 10-5 кг/(м • с), коэффициент теплопроводности Xg = 2.58 • 10-2 кг • м/(с3 • К), времена релаксации тag = тwg = 10-1 с. В момент времени г = 0 диафрагма мгновенно удаляется, при этом реализуется режим течения с ударной волной, движущейся вправо, и волной разрежения, перемещающейся со скоростью с1 = 35.5 м/с влево. Скорость

Рис. 3. Зависимости р/ро^(х), и(х), Т(х), а(х), а(х), Ш(х) к моменту времени Ь = 0.19 с, полученные при учёте неравновесных процессов (сплошные и штрихпунктирные линии — соответственно для Хв = 60.2 • 10-2 и 30.1 • 10-2 кг • м/(с3- К)) и в адиабатическом приближении (штриховые линии)

тепловых волн у контактной границы c2 = 8.1 м/с. Нижние индексы означают: 0 — в невозмущённой среде, L и R — для параметров смеси соответственно слева и справа от контактного разрыва.

На рис. 2 представлены результаты расчётов, полученные методом узлового метода характеристик к моменту времени t = О.5 c на сетке из 200 узловых точек, которые сопоставлены с решением рассматриваемой задачи Римана в адиабатическом приближении, рассчитанным на той же сетке, но с использованием конечно-разностной схемы Куранта — Изаксона — Риса. Как видно из рис. 2, д, влияние сил вязкого трения на общую картину течения незначительно. Большее значение оказывают тепловые эффекты.

Во второй задаче — о симметричном разлёте двух одинаковых масс пузырьковой жидкости, значения переменных были следующими: u0l = —1 м/с, u0r = 1 м/с, p0L = p0R = 0.1 МПа, а0ь = a0R = 0.12. Параметры составляющих смесь фракций оставались те же, что и в предыдущей задаче. В этом случае формируются две волны разрежения. На рис. 3 для момента времени t = 0.19 c приведены результаты расчётов, выполненные с помощью узлового метода характеристик на сетке из 200 узловых точек. Штрихпунк-тирной кривой отмечены зависимости, полученные для варианта с уменьшенным в два раза коэффициентом теплопроводности жидкости при неизменных значениях других параметров компонентов смеси. Здесь же представлены результаты расчётов, полученные в рамках адиабатического приближения на той же сетке, но с использованием схемы Куранта — Изаксона — Риса.

Как видно из данных рис. 2 и З, количественные характеристики течения смеси при распаде разрыва существенно зависят от интенсивности теплообменных процессов.

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

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

[1] Суров В.С. Об отражении воздушной ударной волны от слоя пены // Теплофизика высоких температур. 2000. Т. 38, № 1. С. 101-110.

[2] Суров В.С. Расчёт взаимодействия воздушной ударной волны с пористым материалом // Вестник Челябинского гос. ун-та. 1997. № 1. С. 124-134.

[3] Суров В.С. Анализ волновых явлений в газо-жидкостных средах // Теплофизика высоких температур. 1998. Т. Зб, № 4. С. б24-бЗ0.

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

[5] Суров В.С. Узловой метод характеристик в многожидкостной гидродинамике // Там же. 2013. Т. 8б, № 5. С. 1080-108б.

[6] Goncalves E. Numerical study of expansion tube problems: Toward the simulation of cavitation // Comput. & Fluids. 2013. Vol. 72. P. 1-19.

[7] Schoch S., Nikiforakis N., Lee B.J., Saurel R. Multi-phase simulation of ammonium nitrate emulsion detonations // Combustion and Flame. 2013. Vol. 1б0. P. 1883-1899.

[8] Суров В.С. Течение Буземана для односкоростной модели гетерогенной среды // Инженерно-физический журнал. 2007. Т. 80, № 4. С. 45-51.

[9] Wackers J., Koren B. A fully conservative model for compressible two-fluid flow //J. Numer. Meth. Fluids. 2005. Vol. 47. P. 1337-1343.

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

[11] Deledicque V., Papalexandris M.V. A conservative approximation to compressible two-phase flow models in the stiff mechanical relaxation limit // Ibid. 2008. Vol. 227. P. 9241-9270.

[12] Saurel R., Petitpas F. and Berry R.A. Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures // Ibid. 2009. Vol. 228. P. 1678-1712.

[13] Kreeft J.J., Koren B. A new formulation of Kapila's five-equation model for compressible two-fluid flow, and its numerical treatment // Ibid. 2010. Vol. 229. P. 6220-6242.

[14] Cattaneo C. Sur une forme de l'equation de la chaleur elinant le paradoxe d'une propagation instantance // CR. Acad. Sci. 1958. Vol. 247. P. 431-432.

[15] Самарский А.А., Попов Ю.П. Разностные методы решения задач газовой динамики. М.: Наука, 1980.

[16] Лодж А.С. Эластичные жидкости. М.: Наука, 1969.

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

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

[19] Суров В.С. Односкоростная модель многокомпонентной теплопроводной среды // Инженерно-физический журнал. 2010. Т. 83, № 1. С. 132-141.

[20] Куликовский А.Г., Погорелов Н.В., Семёнов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. Изд. 2-е, дополн. и исправл. М.: Физматлит, 2012. 656 с.

Поступила в 'редакцию 12 сентября 2013 г., с доработки — 24 февраля 2014 г.