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

УДК 519.6

Г. И. Щепановская

МАТЕМАТИЧЕСКОЕ И ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ТЕЧЕНИЙ ВЯЗКОГО ТЕПЛОПРОВОДНОГО ГАЗА*

Предлагается алгоритм численного решения модифицированных уравнений Навье-Стокса для одномерного движения вязкого теплопроводного газа. Проведены тестовые расчеты. Реализована задача о распространении теплового импульса в газе. Апробированная компьютерная модель используется для изучения одномерных геодинамических процессов.

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

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

Для аппроксимации полной производной по времени каждого уравнения системы используется метод траекторий, который заключается в аппроксимации субстанциональной производной с помощью разностной производной назад по времени вдоль траектории движения частицы. Метод траекторий впервые появился в работах французских и американских ученых для аппроксимации уравнений Навье-Стокса для вязкой несжимаемой жидкости с первым порядком аппроксимации. Наибольшее теоретическое и практическое развитие метод получил для одного уравнения переноса массы [2]. Для решения систем алгебраических уравнений используется многосеточный метод с внешними итерациями по нелинейности.

Модификация уравнений Навье-Стокса обеспечивает повышение точности приближенного решения. Как следует из тестовых расчетов, абсолютная погрешность приближенного решения уменьшается в разы по сравнению с аналогичной погрешностью для немодифицированных уравнений, при этом разностная схема остается первого порядка. Применение комбинации метода траекторий и метода конечных элементов позволяет построить экономичный алгоритм с вычислительной точки зрения.

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

ё р ди

—+р— = 0, ё дх

ёи дР + дтхх

Л дх дх

ёе п ди д<7

р—+Р— = Qt ——+Ф,

& дх дх

Р = Р(р, е), ц = ц(р, е),

(1)

(2)

(3)

(4)

где р - плотность; и - проекция вектора скорости на ось х; Р - давление; ц - динамический коэффициент вязкости; е - внутренняя энергия; Q - внешний поток тепла от внешних источников; тензор напряжений тхх, проекция теплового потока дх и диссипативная функция Ф выражаются следующим образом:

4 1 ди

Тхх = Ц - ’ Чх =

3 Яе дх

У

де ■Ц—, РгЯе дх

2

(5)

где Яе - число Рейнольдса, Рг - число Прандтля, У = 1,4.

Редукция уравнений. Преобразуем уравнения (1) и (3) к новому виду. Для этого, учитывая неотрицательность плотности и внутренней энергии, введем следующие функции:

(6) (7)

Подставим (6) в уравнение неразрывности (1) и полу-

р = ст , е = £2.

ё^ + ст’ * = 0. (8)

М дх

Дифференцируя по t, имеем

„ ё ст 2 ди

2ст + ст2— = 0.

ёt дх

Далее, сокращая последнее уравнение на 2ст, получим

ё ст 1 ди

— +—ст— = 0. (9)

ёt 2 дх

* Работа выполнена при финансовой поддержке РФФИ (грант № 11-01-00224), проект № 89 СО РАН.

чим

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

ё (е2) + ддх Л дх

дQ ди

— - Р— + Ф. дґ дх

Преобразуем (10) следующим образом:

, ёе д<7 дQ „ди ^

р2е— + -^ = —- Р— + Ф. & дх дґ дх

(10)

(11)

Далее сократим последнее уравнение на 2е и получим

ё е + 1 ддх 1 дQ Р ди + 1 ф

р ёґ 2е дх 2е дґ 2е дх 2е

(12)

Подставим (7) в выражение для теплового потока Чх из (5) и возьмем производную по х, получим

дЧх

дх

Чх =

2У

2У

де

-це—, РгЯе дх

РгЯе

Ц

д ( де + е—I ц—

дх і дх

(13)

С учетом (13) и выражения для диссипативной функции Ф из (5) уравнение (12) примет вид

ё е у

ёґ РгЯе

ц(де е і дх

д ( де +--------1 Ц—

дх і дх

Р ди 2 ц(ди

2

(14)

ё ст 1 ди

----+—ст— = 0,

ёґ 2 дх

ёи дР + дт

ёґ дх дх

ё е у

ёґ РгЯе

д ( де +--------1 Ц—

дх і дх

ц( де е ідх

Р ди 2 ц(ди

2е дх 3Яе е і дх

(16)

(17)

(18)

Замыкают систему уравнений (16)-(18) алгебраические соотношения для давления и динамического коэффициента вязкости совершенного газа

Р = Р(ст,е), ц = ц(ст,е). (19)

Дискретизация. В качестве расчетной области возьмем единичный отрезок О = [0,1] с границей Г, состоящей из концов отрезка. Введем равномерную сетку х, = ¡И, I = 0,1,п с шагом И = 1 / п (п > 2).

Обозначим множество узлов области О:

Пи = № = (X), I = 0,1,п},

введем множество внутренних узлов

Ои = {5 = (X), = 1,-, п -1}

и два граничных узла

Ги = {5,- = (х,.), I = 0, п}.

В результате расчетная область О разбилась на п интервалов.

Для каждого узла 5 е ОИ введем базисную функцию ф, (х):

/-ч К1 -їх-х1), хє ^-^xi+l], Фг (х)=<р

І0, х є (xг-l, х+l),

(20)

которая равна единице в 5 и равна нулю во всех остальных узлах О .

Будем искать приближенное решение в виде

2е ^‘ 2е дх 3Яе е ^ дх,

Замечание. Обратим внимание на появившиеся множители ц/е и Р/е в уравнении (14), которые, как показывают расчеты, являются естественными «регуляризаторами». Например, для совершенного газа уравнение состояния (4) запишется как

Р = Р(у-1)е или Р =ст2(у-1)е2. (15)

Поскольку внутренняя энергия всегда положительна и больше единицы по отношению к ее величине на бесконечности, то множитель 1 /е «гасит» растущее как е2 давление. Аналогичны рассуждения относительно ц/е . Для совершенного газа, как следует из формулы Сазерленда, динамический коэффициент вязкости является степенной функцией от внутренней энергии.

Будем решать систему уравнений, преобразованную к следующему виду:

ст

!(ґ, х) = £стг (ґ )фг (х),

і=0

и (ґ, х) = ^ иі (ґ )фі (х),

і=0 п

еп (ґ, х) = £ег (ґ )Фі (х).

(21)

(22)

(23)

і=0

После дискретизации по пространству уравнения (16) получаем дискретный аналог уравнения неразрывности:

ёст

1

-± ^стг (и1+1 - щ-1) = °

ёґ 4п

I = 1,п -1.

(24)

Дискретный аналог уравнения количества движения (17) принимает следующий вид:

ёи 2 1

р' + 7ТГ~Б~((иі - и-1 )(Ц-1 + Ц) -

ёґ 3П Яе

1

2П

- (иі+1 - и )(Ц/ +Ц+1)) =

(Р/-1 - Р/+1), ' =1,..., п -1.

(25)

По аналогии с предыдущим уравнением выпишем дискретный аналог уравнения энергии (18) в следующем виде:

- «-+1-е-)2+<6'-е'-1)2)+

1 У

</+4х1) = ст^+\ ик+‘(х1) = и^1 и ек+Чх1) = е*+1 выпишем каждую систему.

Для уравнения переноса массы система алгебраи-

~ к+1

ческих уравнений относительно неизвестных ст

имеет вид

И I 1 /„к „,к \ к+1 _ И к, \гк\

+ (и/+1 и/-1) |ст/ = ст (х1),

т 4 ) т (34)

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

Для уравнения количества движения система алгебраических уравнений относительно неизвестных в узлах разбиения ик+1 имеет вид

_к+1 .,к+1

к+1

„к+1 /

к+1

((еі -еі-1)(Ц'-1 +Ц') - (е'+1 -еі )(Ц' +Ц'+1)) =

— (иі+1 - иі-1) + “77"Г ((иі+1 - иі) + (иі - иі-1) ),

2П2 РгЯе 1 Р

4П е,

3П Яе еі і = 1,..., п -1.

(26)

Метод траекторий. Введем равномерную сетку по времени с шагом т = / т :

йт= {tk : tk = кт, к = 0,..., т}.

Производную по времени (субстанциональная) в уравнении (24) аппроксимируем с помощью разностной производной назад по времени вдоль траектории движения частицы следующим образом [2]:

(27)

+4х1) -стк (Хк) ёt т

где Х<к - решение при € = кт уравнения

ёХ

— = иЦ *, X), X ((к + 1)т) = хг. (28)

ёt

Для численного решения уравнения (28) применим метод Эйлера, вычисление производится с (к + 1)-го шага назад по времени. Полагая, что частица в промежутке времени (?к, tk+!) движется равномерно, получаем

Х1 = х1 - и (ґ*, хі) т.

(29)

Обозначим и^*, хг) = игк . Значение стк (Хгк) определяется путем линейной интерполяции.

При и > 0

стк (Хк) = стк (хг) + (Хк - хг)ст (х1) -ст (х1 -1). (30)

При и1 < 0

стк (хк) = стк (х1) + (Хк - хі)

хі - хі-1

к „чст (хі+1)-ст (хі)

хі+1 - хі ’

(31)

Производную по времени в уравнениях (25) и (26) аппроксимируем аналогично разностной производной (27), в результате получим

ёи, к+1 ик+1( хі) - ик (хк)

рі^Г *рі Т ’

_ ёеі +1 ек+1(хі) -ек (Хк)

рі^Г ~рі Т .

(32)

(33)

2 1 / к к \

-7ГТТГ- (ці-1 + ц) 3П2 Яе

иы +

к+1 рі 1

2 1

(ц/-1 + 2ц/ + Ц/+1)

т 3П2 Яе

ик+1 +

2 1

3П2 Яе

(35)

ик+1 = иі+1 _

рк+1 1

= 21—ик (Хк ) + — (РД - Рк+!), 1 = 1, ..., п - 1.

т 2И

Система квазилинейных алгебраических уравнений, соответствующая уравнению внутренней энергии, выглядит следующим образом:

1 УЦ//к к \ 1

-(ек-ек-1) -

2п РгЯе ек

к+1 р/

2П РгЯе

^ (цк-1 + цк)

ек+1 + еі-1 +

^ ЬЦ- (2ек-ек+1 -ек-1)+

т 2П2 РгЯе є/

1

^ (2ц/ + Ц/+1 + Ц/-1)

2П РгЯе

1 У ці , к к\ 1

-(ек+1 ~ек/)-

ек+1 +

2П2 РгЯе е/

2П2 РгЯе

(36)

к+1 еі+1

к+1

. рі „к ^кл

-ек (Хк) - 4-\ (ик++11 - игк-+/) + т 4И ек

+ 4 [(и^ - ик+1)2 + (ик+1 - и-)2],

3И2 Яе ек

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

Таким образом, получена консервативная вариационно-разностная схема первого порядка аппроксимации.

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

Тестовые расчеты. Для тестирования полученной математической модели функции и, ст, е задаются следующим явным образом:

-і

Значения ик (Хк) и ек (Хк) вычисляются путем линейной интерполяции, аналогично стк (Хк).

Системы алгебраических уравнений. После аппроксимации производной по времени в дискретных аналогах (24)-(26) получаем системы квазилинейных алгебраических уравнений. С учетом обозначений

и (х, ґ) = х( х -1) ґ,

+

+

+

ст(х, ґ) = х(х -1) ґ +1, е(х, ґ) = х(х -1)2 ґ +1.

(37)

а\\2п = | |(5а)2 ёх | ,

5а = а* - а,-

|(5й)2 ёх = -2 [|| а0 - ао| + |аП - аП\ ] + кг"-агЙ|

і=1

||&

•а|| . п = тах аі - аг-

Для ||$ст||вдп и ||5р||» й можно получить соотноше-

ние

115р1 , п =N1 , п тах1 <+стп

з,П II II»,п '

Норму ||5р||2 п можно выразить через ||5ст^|» п следующим образом:

||8Нк, =[2 0|8р1^ +|8рп

+пЁ118рі II

2

і=1

ІМЦ =[^ ІІНІ,,«|ст0 +ст5

II8

<СТ,|», К +ст.П і +

+ »1|>

і=1

сті

ст,- + ст,-

2 Ї 2

При подстановке функций (37) в исходные уравнения и модифицированные уравнения получаются соответственно правые части /р, /ст, /и, /е и /Е,

которые учитываются в системах уравнений при их численном решении. Нормы погрешности в пространстве Ь2 и Ьт определяются следующим образом:

Нормы погрешностей в табл. 1 и 2 приведены для шага по времени т = 0,000 1 в момент времени t = кт (к = 1000).

Задача о распространении теплового импульса.

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

Яе, число Прандтля Рг, число Маха Мт и ю имеют следующие значения: у = 1,4, Яе = 2 -103, Рг = 0,7, М ю = 4, ю= 0,8 [4-6]. Температура Т определяется из следующего уравнения состояния:

Т =

ум» Р р

(38

В качестве начальных условий задаются условия затухания возмущений в бесконечном удалении от источника. Условия на границе для плотности, скорости и энергии следующие: р| х = р| х=1 = 1,

и\ 0 = и\ 1 = 0, е 0 = е 1 = 1).

1х=0 1х=1 ’ 1х=0 1х=1 у

Графики распределения плотности, температуры, скорости и давления газа при И = 0,002(п = 500), т = 0,0001 представлены на рис. 1-4. Кривые, обозначенные на рисунках как 1, 2, 3, 4, соответствуют расчетным моментам безразмерного времени t = 0,01; 0,05; 0,1; 0,2.

Таблица 1

п 0,05 0,04 0,025 0,02 0,0125 0,01

Н12,п 0,000 044 1 0,000 031 9 0,000 016 1 0,000 011 7 0,000 006 1 0,000 004 6

1|8р| 0,000 088 2 0,000 063 7 0,000 032 1 0,000 023 3 0,000 012 2 0,000 009 2

Іі5^! 2,п 0,000 033 8 0,000 027 8 0,000 019 6 0,000 017 2 0,000 014 2 0,000 013 3

118Н12,п 0,000 067 2 0,000 055 1 0,000 038 8 0,000 034 1 0,000 028 0 0,000 026 4

115Н12,п 0,000 145 0 0,000 106 5 0,000 057 5 0,000 044 1 0,000 027 5 0,000 023 0

Таблица 2

п 0,05 0,04 0,025 0,02 0,0125 0,01

115стЩ 0,000 248 2 0,000 200 6 0,000 128 2 0,000 103 8 0,000 067 0 0,000 054 6

115рЩ 0,000 495 9 0,000 400 7 0,000 256 0 0,000 207 3 0,000 133 8 0,000 109 1

ІМЦ 0,000 211 9 0,000 192 9 0,000 164 1 0,000 154 4 0,000 139 7 0,000 134 8

ІМІ »,п 0,000 420 6 0,000 382 6 0,000 324 9 0,000 305 4 0,000 276 0 0,000 266 2

N1 »,п 0,000 839 1 0,000 692 0 0,000 469 2 0,000 394 5 0,000 281 9 0,000 244 2

Рис. 1

Рис. 2

Рис. 3

Рис. 4

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

Библиографические ссылки

1. Ушакова О. А., Шайдуров В. В, Щепановская Г. И.

уравнений На-

Метод конечных элементов для

вье-Стокса в сферической системе координат // Вестник КрасГУ. 2006. № 4. С. 151-156.

2. Pironneau O. On the Transport-Diffusion Algorithm and Its Applications to the Navier-Stokes Equations // Numerische Mathematik. 1982. № 38. P. 309-332.

3. Самарский А. А., Николаев Е. С. Методы решения сеточных уравнений. М. : Наука, 1978.

4. Шайдуров В. В., Щепановская Г. И., Якубович М. В. Одномерная модель динамики вязкого теплопроводного газа // Материалы XIV Междунар. науч. конф. «Решетневские чтения». Красноярск, 2010. С. 440-441.

5. Шайдуров В. В, Щепановская Г. И. Газодинамическая модель внутреннего строения Земли // Вестник СибГАУ. 2008. № 1. 79-83.

6. Vyatkin A. V., Shaidurov V. V., Shchepanovskaya G. I. Numerical Spherically-Symmetric Simulation of Deep-Seated Geodynamics // J. of Applied and Industrial Mathematics. 2010. Vol. 4. № 2. P. 290-297.

7. Шайдуров В. В., Щепановская Г. И., Якубович М. В. Применение метода траекторий и метода конечных элементов в моделировании движения вязкого теплопроводного газа // Вычислит. методы и программирование. 2011. Т. 12. С. 275-281.

G. I. Shchepanovskaya

MATHEMATICAL AND NUMERICAL MODELLING OF CURRENTS OF VISCOUS THERMALLY CONDUCTIVE GAS

In the work the author offers an algorithm of the numerical decision of the modified equations of Nave-Stoksa for one-dimensional movement of viscous thermally conductive gas. Test calculations are carried out. The task of distribution of a thermal impulse in gas is realized. The approved computer model is used for studying of one-dimensional geodynamic processes.

Keywords: equations of Nave-Stoksa, viscous thermally conductive gas, numerical modeling, a method of trajectories, a method of final elements.

© Щепановская Г. И., 2011