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

Вестник Челябинского государственного университета. 2010. № 24 (205) Физика. Вып. 8. С. 15-22.

В. С. Суров, Е. Н. Степаненко

СЕТОЧНЫЙ МЕТОД ХАРАКТЕРИСТИК ДЛЯ РАСЧЁТА ТЕЧЕНИЙ ОДНОСКОРОСТНОЙ МНОГОКОМПОНЕНТНОЙ ТЕПЛОПРОВОДНОЙ СРЕДЫ

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

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

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

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

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

С. К. Годунова, Мак-Кормака, Лакса — Венд-роффа и других подобных им невозможно.

Гиперболическая модель гетерогенной теплопроводной среды. Рассмотрим «-компонентную теплопроводную смесь с первыми т сжимаемыми фракциями. Система уравнений, описывающая течение односкоростной гетерогенной среды, в которой т-й компонент является несущим (все остальные фракции «вкраплены» в него), а также учтены силы межфракционного взаимодействия, имеет вид [5] ф

дл

+ &у(р и) = 0;

+ (и-У)и

дЛ

- grad р = Е;

д

дл

і 1 I |2 рі є+2 и1

+div

дw

[ 1 |2 р ^

р є + — и + — и + w

1 2 р V

= Е • и;

дл да, Р° дл

ті-------+ (и • 1 + х grad Т + w = 0;

П (кФі)

(а,- р°и )=Е '7-к; (1)

к=1

дє,

р| "дг"+(и •У)є, 1 +

+

а,р

р,

п (кфі)

п (кфі)

Ек -[1+(" 'у)р-

к=1

Е(к+& )-[є, - 2 и|2 ^;

к=1

п (к ф)

к=1

і = 1, ... , т - 1;

Р

да

~дГ

Й1У (а 1и ) = Ео X -}к,

р./ к=1 1 = т + 1, ... , «,

где Е — плотность массовой силы; } — интенсивность превращения массы из /-й фракции в 1-ю на единицу объема смеси; р — давление; Q/j — тепловыделение в единицу времени на единицу объема смеси, возникающее вследствие превращения /-й фракции в 1-ю; Щ — количество тепла в единицу времени на единицу объема смеси, которое поступает /-й фракции от 1-й вследствие излучения; Т — осредненная температура; и — вектор скорости; w — вектор плотности теплового потока; а/ — объемная доля; е, — удельная внутренняя энергия; т — коэффициент тепловой релаксации смеси; р —

О ■ -

плотность смеси; р, — истинная плотность /-й

фракции; . = а, р0 — приведенная плотность

/-го компонента; % — коэффициент теплопроводности смеси.

Поведение сжимаемых фракций описывается калорическими и термическими уравнениями состояния, которые для /-й фракции в общем случае имеют вид в, = в,(р, р0) и Т1 = Т1 (р, р0), поэтому выражения для удельной внутренней энергии смеси

Система определяющих уравнений (1) с учетом соотношений (4)-(5) при отсутствии массовых сил, фазовых и химических превращений, а также переноса тепла за счет излучения для одномерных плоских течений приводится к виду

др др ды ды ды 1 д р

— + ы—- +р— = 0; — + ы—+-------------— = О;

дt дх дх дt дх р дх

др др 2 ды ттд w

—- + ы— + рс2— + Н-----------= 0;

дt д х д х д х

дw дw , др , д р

----+ ы------------+ кЕ —+ к^ — +

дt дх р дх дх

т—1

+

Хка^ + к., ^р/

о Л

/=1

д х х

дt

да

дt

(6)

М. + „ М. +р°0! = 0;

дt дх / дх

да/ -ы^ + а/(1 — а,)^ = 0, / = 1, ... , т - 1;

д х д х

ды „ ,

а.— = 0, 1 = т + 1, ... , «,

1 д х

где

«

в=- / ,р- в,, р-,=1

(2)

к =1^, кр =^дТ, ка =-*5Т Е х Зр х др

х За1

а также средней температуры, определяемой формулой

П

т = > а, Т,, /=1

могут быть записаны соответственно как в = в (р, р, а1, р.,...,

ат—1,рт—1, ат+!,•••,ап)

Т = Т (р, р, а1, р!,-, а т—1,

(3)

(4)

к- 0 = 1-^... ка дТ

р1 х ’ ат—1 х дат—1

хдр,

кр0 =

гт-1

X зт

х Зр0т—1

к =Х^^ к =х дТ

™ А.а

х За

т+1

х Зап

Соответствующие выражения для Н, а/ и адиабатической скорости звука с имеют вид

/ л —1 / Лп—1

де + ^ де, др до

Н =

,=1

1

а=

р,

Г

чдР.у

' де, ^ —1

1др. V

^ а, де & е

р, да, др, Л

Р. — .с 2

р, др

рт—1,а т+1,-,а п ).

(5)

и

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

р - с,

є, =

рв, + (,

где

в = У(У, -1), Л, = 4в,, ь, = 4р*,:

то выражение (2) запишется как

є = -

рВт + Ьт + т—1

+ Еа, (т + ('т — Л,тр,0 )

,=1

п

Е

у =т+1

0

а урує у

— Лт .

Здесь

В,т = В, — Вт , (,т = (, — (т , Л,т = Л, — Лт .

Соответствующие выражения для Н, О, и с примут вид

Н = р

В„

т—1

Е

,=1

а, ((тВ' — Ь'Вт )

+ рВ,

О = Рс В, — р .

' ь + рв, ;

В простейшем случае бинарной смеси идеального газа с показателем адиабаты у и несжимаемой второй фракцией газодинамическая часть системы уравнений (6) имеет вид

др др ди ди ди 1 д р

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

дл дх дх дл дх р дх

др д р 2 ди ТТдw п

— + и— + рс2— + Н----------= 0; (7)

дл д х д х д х

да да ,дип

— + и------(1 — а)— = 0;

дл дх д х

где с = — адиабатическая скорость звука

У ар

в смеси, Н = ——1)р, а — объемная доля газа в а

смеси.

Выражение для закона Фурье с тепловой релаксацией, учитывая соотношение Т = Т (р, р, а), которое является следствием (5), пере-

пишем как

дw дw , др , др , да 1 „

+ и-------+ кр----+ кр----+ ка-----+— w = 0, (8)

дл дх р дх р дх а дх т

где

к X дТ к X дТ % дТ

кр Л , кр Л , ка Л .

^ я~ ^ т да

х д^ р х др

Характеристическое уравнение системы (7)-(8) имеет пять действительных корней: ы ± с1, ы, ы ± с2 [5], где

с =

2]с2 + крН + . Іс4 + Н

кр (2с2 + крН) + 4

1

Р

1

2 V + kpH —. Ic4 + H

kp (2c2 + kpH) + 4 Отметим, что C1 определяет скорость рас-

f — ка (1 — а) Л P „

пространения газодинамических возмущений, а с2 — тепловых. Характеристические соотно-

шения на характеристиках ^ = u ± C1 и

^ = u ± C2 могут быть получены из уравнения

§— u — P О О

О §—u — V P О

0 — pc2 § — u О

О 1 — а О § —u

— kP О — kP — ка

d p du

— u —— p— dt dt

du і dp

dt p dt

2 du dp „ dw

— pc2-------------u — — H—

dt dt dt

. du d а

(1 — а)---------------u------

dt dt

,^w — k dP — k dP . kp k p

т p dt p dt

d а dt

= 0.

(9)

После вычисления определителя из (9) получим следующие выражения, справедливые вдоль характеристических направлений ^ = ы ± с1 и £ = ы ± с2:

(u ± C1)

(c2 — c2)

ка d а + к d p ±

т

±c,

kpPc + kpP — ка(1 — а)

kpc12 + kp —

ка (1 — а)

(C2 — C )

±c.

11 11 wc2

ка d а + kp d p±—

т

kpc2 + kp —

kpPC 2 + kpP — ка (1 — а)

ка (1 — а)

те сжимаемости жидкости. Выражение для средней температуры (З) дает

T =-

а2 p

[p —p0 (l —а)] R

+ (1 — а)?0, (13)

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

[ ± c1)dp + Hdw] = 0. (10)

kP=-

а 2х P

tR[P — P2(1 — а)]

kp =

2

а Х

tR[P — P2(1 — а)]

ка =-

а p(2p + аp2) R[P — P00(1 — а)]

2 — T0

В частности, для водно-воздушной смеси при нормальных условиях и объемной доле

[(u ±c2)dp + Hdw] = 0. (11) газовой составляющей а = 0,1 (p2 = 1000 кг/м ),

Вдоль траекторной характеристики справедливо выражение

d p d а

— +-------= 0.

p 1 — а

(12)

При рассмотрении пузырьковой жидкости без потери точности можно считать, что температура несжимаемой фракции постоянна, т. е. Т2 = const, поскольку ее доля в смеси значительна. Это предположение снимается при уче-

значения скоростей ^ , 02 и с равны 38,51, 8,48 и 39,44 м/с соответственно. В расчетах коэффициент теплопроводности смеси определяется из выражения

Х = -((1 +Р2Х 2 ) р' ’ где х = 2,58 10-2 кг-м/(с3-К) — для воздуха, X = 60,2 10-2 кг-м /(с3-К) — для воды, а коэффициент х = 100 с.

и

Метод характеристик для гетерогенной среды. Опишем применение метода характеристик для фиксированной сетки. Для решения задачи достаточно рассмотреть способ определения значений искомых параметров в узле (хк, tn+\) по известным данным в узлах конечноразностной сетки на п-м временном слое. Решение поставленной задачи получим с использованием следующей итерационной процедуры. На «нулевой» итерации (5 = 0) полагаем, что значения искомых функций в точке (хк, tn+\) совпадают с данными из (хк, Лп), поэтому характеристические направления ёх/^ = и ± С1,

dx|dt = и , ёх/ё = и ± С2, отмеченные на рис. 1 прямыми линиями, исходящими из узла (хк, tn+\), аппроксимируются выражениями

хк - х^ = (и(^ + с,*))Дґ;

(*) , „(*)>

хк - х^ = (и(^ + с,*))Дґ;

хк - хС) = и(^Дґ;

х^ - х« = (и(*) - с,*))Дґ;

хк - х« = (и(^ - с,*))Дґ,

ГДе Д = ґп+1 - Ґп .

х[5) = хк - (и(5) + с^))Лt,

хС5) = хк - и(5)Лt; (14)

х« = хк - (и(5) - с1 5, х£ = хк - (и(5) - с25))Лt..

Параметры р, и, р и а в найденных точках хь , хь^, хс, хК1, хк^ находятся интерполяцией по их известным значениям в узлах хк-1, хк и хк+1. Соотношения вдоль характеристических направлений из (10)-(12) перепишем в конечно-разностном виде как

А11р( )(хк, п+1) + А12и( )(хк, ^-п+1) +

+ А13р(5+1) (хк, 1:п+1) + Л4а(5+1) (хк, ^+0 +

+Аи М ) (хк, tn+l) = ,

А21р(5 1)( хк, 1:п+1) + А22и (5 1)( хк, 1:п+1) +

+А23 р (5+1}( хк, ^) + А24а(5+1)( хк, ^+1) +

+А25м( )(хк, tn+1) = ^Я1,

А51р(5+1)(хк, tn+l) + А^а^Чхк, tn+l) = ^ (15)

А31р( )(хк, ^п+1) + А32и ( )(хк, ^-п+1) +

+ А33р(5+1} (хк, 1:п+1) + А34а(5+1} (хк, 1:п+1) +

+ А35 М )(хк, ^п+1) = «V

Рис. 1. Схема пересчета параметров с п-го временного слоя на (п+1)-й с использованием сеточного метода характеристик

Из последних равенств определяем положения точек пересечения характеристик с прямой ґ = ґп (см. рис. 1):

х^) = хк - (и(^ + с,* ))Дґ;,

А4іР(* 1)( хк, Ґп+1) + А42и ( * 1)( хк, ґп+1) +

+ А43 Р (*+1)( хк, ґп+і) + А44а (*+1)( хк, ґп+1) +

+ А45^( *+1)( хк, ґп+і) = V,

-------где

А11 =[кр (с12 - с2)]

А12 = |С1

А13 =

2 "!)(*)

крРС + крР- ка (1 -а),

кА + кр-

( *)

А14 =[ка (с1 - С2)

(*) Ь :

А15 = '

Н

: + с,

к с, + к -

і^р^ і і д.р

ка (1 -а)

(*)

Ь

(с,2 - с2)с1м>

(*)

Дґ + А,,рЬ ) +

+ А12иЬ*) + А13 РЬ*) + А14аЬ*1) + А15^ь*); А21 =[кр (с2 - с2)]^

{-с1

крРс + крР- ка (1 -а)

^( *) }Я ,

А23 =

крс,2 + к-,-

( * )

А25 = '

А24 =[ ка (с2 - с 2) ] ,

Н

- с,

крс2 + кр-

,(*)

5Яі =

(с,2 - с2)с1м>

(*)

Дґ + А21рЯ*) + А22 иЯ*; +

(*)

+ А23 рЯ1) + А24а Я*1) + А25^^,);

А31 =[кр (с2 - с 2) ]

А32 ={с2 [ кр Рс 2 + крР- ка (1 -а) ]}”

А33 =

крс|+кр-

(*)

А34 =[ка (с2 - с2)]

А35 = '

Н

=

: + с-,

(с| - с2)с2 ^

кр4+к,-

(*)

(*)

Дґ + А31рЬ2) + А32иь2) +

+А33 рь2) + А34аЬ*2) + А35^Ь2);

А41 =[кр (с22 - с2)]^

А42 ={-с2 [крРс2 + крР- ка(1 -а)]}^):

А43 =

крс2 + к,-

(*)

А44 =[ка (с2 - с2)]^

А45 = '

Н

- с

(с^ - с2)с2 ^

+ А43 рЯ2) + А44аЯ52) + А45 Мя2 э

А51 = (1 - а)С,

А52 = 0, А53 = 0 А54 = рС, А55 = 0,

^С = А51рС + А54аС. .

Решая систему (15) (при 5 = 0) относительно переменных р(1), и(1), р(1), а(1) и м,(1), найдем уточненные значения искомых функций в точке (хк, tn+l). Затем по этим данным из выражений (14) вычисляются новые координаты хт(1>,

•Ц

хТТ4, хС1), х®, х£>, которые в свою очередь

используются для определения р(2), и(2), р(2), а(2) и м,(2) из (15), где необходимо положить 5 = 1. Описанный итерационный процесс продолжается вплоть до сходимости.

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

(х < 0) - Ро(1) = 0,15 МПа, Ио(1> = 0, ао(1> =

= 0,073 , 70(1 = 3130 К, р00а> = 1,67 кг/м3, у а> =

= 1,4, р00(1> = 1000 кг/м3, хК1> = 2,58 10-2

кг-м /(с3-К), ^2(1> = 60,2 10-2 кг-м/(с3-К),

Т(1> = 100 с; справа от нее (х > 0) - р0(2> = 0,1

МПа, и0СТ = 0, а0^> = 0,1, = 2930 К,

р00(2> = 1,19 кг/м^ У (2> = 1,4 р00(2) = 1000 кг/м3, х1(2> = 2,58 10-2 кг-м/(с3-К), х2(2> = 60,2

(*)

ка (1 -а)

Дґ + А4,р^) + ) +"

(*)

(*)

■/ 4 '

,,(*)•

10 кг-м /(с -К), Т(2) = 100 с. В момент времени ґ = 0 диафрагма мгновенно удаляется, при этом реализуется режим течения с ударной волной, движущейся вправо, и волной разрежения, перемещающейся влево. На рис. 2 представлены результаты расчетов, полученные к моменту времени ґ = 22,5 мс (квадратные точки), которые сопоставлены с решением задачи Римана; полученные на той же сетке, но с использованием конечно-разностной схемы

Р1 Р0<

а)

и

и,

Рис. 2. Зависимости параметров течения при распаде произвольного разрыва к моменту ґ = 22,5 мс, полученные: методом характеристик для теплопроводной смеси (квадратные точки), при отсутствии теплопроводности (круглые точки); методом Куранта — Изаксона — Риса для теплопроводной смеси (штриховые кривые); точное решение для нетеплопроводной смеси

(сплошные линии)

Куранта — Изаксона — Риса [6] (штриховые кривые):

ик+1 - ик +А* Ц+1/2 - Ц/-1/2 =8*

где

Дґ

ик

т+1/2

Д х

4 (и т+и т+і у

2{°[^П (Л)]П-1}‘( (икт - ит+1 ),

т = і, і - 1.

2

Здесь

и =

и

р

а

А =

и Р 0 0 0

0 и 1Р 0 0

0 о а и 0 Н

0 а-1 0 и 0

кР 0 кр ка и у

(

0 ^ 0 0 0

- w| х

Матрицы О, Л, О 1 определяются в соответствии со следующими формулами:

Q =

kp gc2 1 g p k 2ka (c2 - c22) kp gc12 pk 2

1 k - kp g - ka pk g - a kp g - ka pk g 1 a

kp H kp H 0 kp H kp H

Pc1 pc2 pc2 pc1

kp H - kp H 0 - kp H kp H

kpc2 kpc2 2kp(c2 - c2) kpc2 kpc2

kp g - ka g i a g - a g i a kp g - ka

kp (kpH - c\) кр (kpH - c2)

Л =

Q-1 =

V c1 c I

( u - c1 0 0 0 0 >

0 u - c2 0 0 0

0 0 u 0 0 ?

0 0 0 u + c2 0

0 V 0 0 0 u+c J

( pc1 c2 ) k - ■cL c2 ka \ _fL

kp kp H J kp H kp kp

pC2 c2 ^ k - — c22 ka _ fi.

kp kp H J kp H kp kp

0 0 g 0

pC2 ( c 2 k - -2 \ c2 ka c2

kp p H \ J kp H kp kp

pc1 ( c 2' k - ^ \ c12 ka fL

kp p H V / kp H kp kpJ

где g = р/(1 -а). Для сравнения на этом же рисунке представлены результаты расчетов при тех же начальных условиях, но при отсутствии теплопроводности в смеси, также полученные методом характеристик (круглые точки). Здесь же приведено точное решение задачи Римана для нетеплопроводной смеси [7] (сплошные кривые).

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

Выводы: описан метод характеристик,

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

kp (kpH - c2) кр (kpH - c\)

понентной теплопроводной смеси на фиксированных пространственных сетках со «сквозным» расчетом ударных скачков. Рассмотренный подход может быть непосредственно распространен на задачи с несколькими пространственными переменными.

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

1. Зауэр, Р. Течения сжимаемой жидкости / Р. Зауэр. М. : Иностр. лит., 1954.

2. Ляхов, Г. М. Ударные волны в грунтах и в воде вблизи от места взрыва / Г. М. Ляхов, В. Н. Охитин, Ф. Г. Чистов // Приклад. механика и техн. физика. 1972. № 3. С. 131.

3. Хоскин, Н. Э. Метод характеристик для решения уравнений одномерного неустано-вившегося течения / Н. Э. Хоскин // Вычисл. методы в гидродинамике. М. : Мир, 1967.

4. Reverberi, A. P. On the non-linear Maxwell-Cattaneo equation with non-constant diffusivity: Shock and discontinuity waves / A. P. Reverberi [et. al.] // Int. J. Heat and Mass Transfer. 2008. Vol. 51. P. 5327.

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

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

7. Суров, В. С. Задача Римана для односкоростной модели многокомпонентной смеси / В. С. Суров // Теплофизика высоких температур. 2009. Т. 47, № 2. С. 283.