MATHEMATICAL AND NUMERICAL MODELLING OF CURRENTS OF VISCOUS THERMALLY CONDUCTIVE GAS

Abstract

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 dis- tribution of a thermal impulse in gas is realized. The approved computer model is used for studying of one-dimensional geodynamic processes.

Full Text

Система нестационарных одномерных уравнений d ρ +ρ ∂u = , (1) Навье−Стокса для вязкого теплопроводного газа включает три дифференциальных уравнения в част- 0 dt ∂x ных производных, вытекающих из законов сохране- ρ du ∂P + ∂τxx , (2) ния массы, количества движения и внутренней энер- dt ∂x ∂x гии газа. Предложенная в работе замена искомых ρ de + P ∂u = Q − ∂qx + Φ, (3) функций в уравнениях неразрывности и внутренней dt ∂x t ∂x энергии переводит закон сохранения массы и полной P = P(ρ, e), μ = μ(ρ, e), (4) энергии из нормы пространства L1 в норму гильбер- где ρ – плотность; u – проекция вектора скорости на това пространства L2 . Впоследствии это значительно упрощает обоснование устойчивости и сходимости [1]. Дискретизация по пространству модифицирован- ных уравнений Навье−Стокса осуществляется мето- ось x; P – давление; μ – динамический коэффициент вязкости; e – внутренняя энергия; Q – внешний по- ток тепла от внешних источников; тензор напряжений дом конечных элементов. τxx , проекция теплового потока qx и диссипативная Для аппроксимации полной производной по вре- мени каждого уравнения системы используется метод функция Φ выражаются следующим образом: траекторий, который заключается в аппроксимации τ = 4 1 μ ∂u , q = − γ μ ∂e , субстанциональной производной с помощью разност- xx 3 Re ∂x x Pr Re ∂x (5) ной производной назад по времени вдоль траектории движения частицы. Метод траекторий впервые поя- 4 1 ⎛ ∂u ⎞2 Φ= μ⎜ ⎟ , вился в работах французских и американских ученых 3 Re ⎝ ∂x ⎠ для аппроксимации уравнений Навье−Стокса для вяз- кой несжимаемой жидкости с первым порядком ап- проксимации. Наибольшее теоретическое и практиче- ское развитие метод получил для одного уравнения переноса массы [2]. Для решения систем алгебраиче- ских уравнений используется многосеточный метод с внешними итерациями по нелинейности. Модификация уравнений Навье−Стокса обеспечи- вает повышение точности приближенного решения. Как следует из тестовых расчетов, абсолютная по- грешность приближенного решения уменьшается в разы по сравнению с аналогичной погрешностью для где Re – число Рейнольдса, Pr – число Прандтля, γ = 1, 4. Редукция уравнений. Преобразуем уравнения (1) и (3) к новому виду. Для этого, учитывая неотрица- тельность плотности и внутренней энергии, введем следующие функции: ρ = σ2 , (6) e = ε2 . (7) Подставим (6) в уравнение неразрывности (1) и полу- чим 2 немодифицированных уравнений, при этом разност- d (σ ) + σ2 ∂u = 0. (8) ная схема остается первого порядка. Применение комбинации метода траекторий и метода конечных элементов позволяет построить экономичный алго- ритм с вычислительной точки зрения. Математическая модель. Выпишем дифференци- альные уравнения одномерного вязкого теплопровод- ного газа в виде безразмерных уравнений неразрыв- dt ∂x Дифференцируя по t , имеем 2σ d σ + σ2 ∂u = 0. dt ∂x Далее, сокращая последнее уравнение на 2σ , полу- чим ности, количества движения и уравнения для внут- d σ + 1 σ ∂u = 0. (9) ренней энергии: dt 2 ∂x * Работа выполнена при финансовой поддержке РФФИ (грант № 11-01-00224), проект № 89 СО РАН. 101 Математика, механика, информатика Проделаем аналогичную процедуру для уравнения внутренней энергии (3), т. е. преобразуем уравнение (3) к новому виду с учетом (7). Для этого подставим выражение (7) в (3), в результате получим Замыкают систему уравнений (16)−(18) алгебраи- ческие соотношения для давления и динамического коэффициента вязкости совершенного газа P = P(σ, ε), μ = μ(σ, ε). (19) d (ε2 ) ∂qx ∂Q ∂u ρ + = − P + Φ. (10) Дискретизация. В качестве расчетной области dt ∂x ∂t ∂x возьмем единичный отрезок Ω= [0, 1] с границей Г, Преобразуем (10) следующим образом: ρ2ε d ε + ∂qx = ∂Q − P ∂u + Φ. (11) состоящей из концов отрезка. Введем равномерную сетку xi = ih, i = 0, 1, ..., n с шагом h = 1 / n (n ≥ 2). dt ∂x ∂t ∂x Обозначим множество узлов области Ω : Далее сократим последнее уравнение на 2ε и полу- чим Ωh = {S i = ( xi ), i = 0, 1, ..., n}, ρ d ε + 1 ∂qx = 1 ∂Q − P ∂u + 1 Φ. (12) введем множество внутренних узлов dt 2ε ∂x 2ε ∂t 2ε ∂x 2ε Ωh = {Si = ( xi ), i = 1, ..., n −1} Подставим (7) в выражение для теплового потока qx из (5) и возьмем производную по х, получим и два граничных узла Γ = {S = ( x ), i = 0, n}. q =− 2γ x PrRe με ∂ε , ∂x h i i В результате расчетная область Ω разбилась на n интервалов. ∂q 2γ ⎛ ⎛ ∂ε ⎞2 ∂ ⎛ ∂ε ⎞ ⎞ Для каждого узла S ∈ Ω введем базисную функ- x =− ⎜ μ⎜ ⎟ + ε ⎜ μ ⎟ ⎟. (13) i h ∂x PrRe ⎜ ⎝ ∂x ⎠ ∂x ⎝ ∂x ⎠ ⎟ цию ϕ ( x) : ⎝ ⎠ С учетом (13) и выражения для диссипативной i ⎪⎧(1 − xi − x ), x ∈[ xi −1 , xi +1 ], функции Φ из (5) уравнение (12) примет вид ϕi ( x) = ⎨ ⎪0, x ∈ ( x , x ), (20) d ε γ ⎡ μ ⎛ ∂ε ⎞2 ∂ ⎛ ∂ε ⎞⎤ ⎩ i −1 i +1 ρ − ⎢ ⎜ ⎟ + ⎜ μ ⎟⎥ = которая равна единице в Si и равна нулю во всех ос- dt PrRe ⎢⎣ ε ⎝ ∂x ⎠ ∂x ⎝ ∂x ⎠⎥⎦ (14) тальных узлах Ω . 1 P ∂u 2 μ ⎛ ∂u ⎞2 Будем искать приближенное решение в виде = Qt − + ⎜ ⎟ . 2ε 2ε ∂x 3Re ε ⎝ ∂x ⎠ Замечание. Обратим внимание на появившиеся n σh (t, x) = ∑ σ (t)ϕ ( x), (21) множители μ/ε и P /ε в уравнении (14), которые, как показывают расчеты, являются естественными «регуляризаторами». Например, для совершенного газа уравнение состояния (4) запишется как i =0 n uh (t, x) = ∑ i =0 n ui (t)ϕi ( x), (22) P = ρ(γ −1)e или P = σ2 (γ −1)ε2 . (15) εh (t, x) = ∑ ε (t )ϕ ( x). (23) Поскольку внутренняя энергия всегда положи- тельна и больше единицы по отношению к ее величи- не на бесконечности, то множитель 1 /ε «гасит» рас- i =0 После дискретизации по пространству уравнения (16) получаем дискретный аналог уравнения нераз- тущее как ε2 давление. Аналогичны рассуждения рывности: относительно μ/ε . Для совершенного газа, как сле- дует из формулы Сазерленда, динамический коэффи- d σl + dt 1 σ (u 4h l l +1 − ul −1 ) = 0, (24) циент вязкости является степенной функцией от внутренней энергии. Будем решать систему уравнений, преобразован- ную к следующему виду: l = 1, ..., n −1. Дискретный аналог уравнения количества движения (17) принимает следующий вид: d σ + 1 σ ∂u = 0, (16) ρ dul + 2 1 ((u − u )(μ + μ ) − dt 2 ∂x l dt 3h2 Re l l −1 l −1 l ρ du = − ∂P + ∂τxx , (17) − (ul +1 − ul )(μl + μl +1 )) = (25) dt ∂x ∂x 1 = l −1 − Pl +1 ), l = 1, ..., n −1. d ε γ ⎡ μ ⎛ ∂ε ⎞2 ∂ ⎛ ∂ε ⎞⎤ 2h ρ − ⎢ ⎜ ⎟ + ⎜ μ ⎟⎥ = dt PrRe ⎢⎣ ε ⎝ ∂x ⎠ ∂x ⎝ ∂x ⎠⎥⎦ (18) По аналогии с предыдущим уравнением выпишем 1 P ∂u 2 μ ⎛ ∂u ⎞2 дискретный аналог уравнения энергии (18) в следую- = Qt − + ⎜ ⎟ . 2ε 2ε ∂x 3Re ε ⎝ ∂x ⎠ щем виде: 102 Вестник Сибирского государственного аэрокосмического университета имени академика М. Ф. Решетнева ρ d εl − 1 γ μl ((ε − ε )2 + (ε − ε )2 ) + σk +1 ( x ) = σk +1 , u k +1 ( x ) = u k +1 и εk +1 ( x ) = εk +1 вы- l dt 2h2 PrRe εl l +1 l l l −1 пишем каждую систему. + 1 γ ((ε − ε )(μ + μ ) − (ε − ε )(μ + μ )) = Для уравнения переноса массы система алгебраи- 2h2 PrRe l l −1 l −1 l l +1 l l l +1 ческих уравнений относительно неизвестных имеет вид σk +1 − 1 Pl (u − u ) + 1 1 μl ((u − u )2 + (u − u )2 ), 4h ε l +1 l −1 3h2 Re ε l +1 l l l −1 ⎛ h + 1 u k − u k ⎞ σk +1 = h σk X k , l l ⎜ τ ( l +1 4 l −1 ) ⎟ l τ ( l ) (34) l = 1, ..., n −1. (26) Метод траекторий. Введем равномерную сетку по времени с шагом τ= tfin / m : ⎝ ⎠ l = 1, ..., n −1. Для уравнения количества движения система ал- гебраических уравнений относительно неизвестных в ωτ = {tk : tk = k τ, k = 0, ..., m}. узлах разбиения u k +1 имеет вид Производную по времени (субстанциональная) в ⎡− 2 1 (μk + μk )⎤ u k +1 + уравнении (24) аппроксимируем с помощью разност- ⎢ 3h2 Re l −1 l ⎥ l −1 ной производной назад по времени вдоль траектории ⎣ ⎦ + движения частицы следующим образом [2]: ⎡ ρk 1 2 1 k k k ⎤ k +1 d σl ≈ σ k +1 ( x ) − σk ( X k ) , (27) + ⎢ + ⎣⎢ τ (μ + 2μ + μ ) u + 3h2 Re l −1 l l +1 ⎥⎦ l (35) dt τ + ⎡− 2 1 (μk + μk )⎤ u k +1 = ⎢ 3h2 Re l l +1 ⎥ l +1 где X k − решение при t∗ = k τ уравнения k +1 l ⎣ ⎦ k k 1 k k =dX ∗ = τ u ( X l ) + 2h (Pl −1 − Pl +1 ), l = 1, ..., n −1. = u(t , X ), dt∗ X ((k +1)τ) = xl . (28) Система квазилинейных алгебраических уравне- Для численного решения уравнения (28) применим метод Эйлера, вычисление производится с (k + 1)-го шага назад по времени. Полагая, что частица в про- ний, соответствующая уравнению внутренней энер- гии, выглядит следующим образом: k ⎡ 1 γ μl k k 1 γ k k ⎤ k +1 межутке времени лучаем (tk , tk +1 ) движется равномерно, по- ⎢ ⎢⎣ 2h 2 PrRe (ε − ε ) − εk l l −1 2h (μ + μ ) ε + 2 PrRe l −1 l ⎥⎦ l −1 X k = x − u(t∗ , x ) τ. (29) ⎡ k +1 l l k k k Обозначим l l l u(t∗ , x ) = u k . Значение σk ( X k ) опре- + ⎢ ⎢⎣ τ − 2h2 PrRe (2ε − ε − ε ) + k l деляется путем линейной интерполяции. + 1 γ (2μk + μk + μk )⎤ εk +1 + (36) При u k > 0 2h2 PrRe l l +1 l −1 ⎥ l ⎦ k k ⎡ γ μk k k γ k k ⎤ k + σk ( X k ) = σk ( x ) + ( X k − x ) σ ( xl ) − σ ( xl −1 ) . (30) 1 + ⎢− l (ε −ε ) − 1 (μ +μ )⎥ ε 1 = l l l l xl − xl −1 ⎣⎢ 2h2 PrRe εk + l +1 l 2h2 PrRe l l +1 ⎥⎦ l +1 При u k < 0 ρk 1 k k 1 Pk k +1 k +1 l σk ( X k ) = σk ( x ) + ( X k − x ) σ k ( x ) − σk ( xl ) . (31) = l ε τ k ( X l ) − l 4h εk (ul +1 − ul −1 ) + xl +1 − xl + 1 1 μl [(u k +1 − u k +1 )2 + (u k +1 − u k +1 )2 ], Производную по времени в уравнениях (25) и (26) 3h2 Re εk l +1 l l l −1 аппроксимируем аналогично разностной производной (27), в результате получим l = 1, ..., n −1. Таким образом, получена консервативная вариа- ρ dul ≈ ρk +1 u k +1 ( xl ) − u ( X l ) , (32) ционно-разностная схема первого порядка аппрокси- мации. l dt l τ k +1 k k Для решения систем линейных алгебраических ρ d εl ≈ ρk +1 ε ( xl ) − ε ( X l ) . (33) уравнений с трехдиагональной матрицей использует- l dt l τ ся метод немонотонной прогонки, который отличает- Значения u k ( X k ) и εk ( X k ) вычисляются путем ли- ся высокой вычислительной устойчивостью [3]. Тестовые расчеты. Для тестирования полученной нейной интерполяции, аналогично σk ( X k ) . Системы алгебраических уравнений. После ап- проксимации производной по времени в дискретных аналогах (24)−(26) получаем системы квазилинейных алгебраических уравнений. С учетом обозначений математической модели функции u, σ, ε задаются следующим явным образом: u( x, t ) = x( x −1)2 t, 103 Математика, механика, информатика σ( x, t ) = x( x −1)2 t +1, (37) n −1 ⎛ ∗ 1 2 ⎞ 2 h ⎞ ⎟ ε( x, t ) = x( x −1)2 t + 1. + h∑ ⎜ δσi ∞,h σi + σi ⎟ ⎟ . i =1 ⎝ ⎠ ⎟ ⎠ При подстановке функций (37) в исходные урав- нения и модифицированные уравнения получаются Нормы погрешностей в табл. 1 и 2 приведены для шага по времени τ = 0, 000 1 в момент времени t = k τ соответственно правые части fρ , fσ , fu , fe и fε , (k = 1 000) . которые учитываются в системах уравнений при их численном решении. Нормы погрешности в простран- стве L2 и L∞ определяются следующим образом: 1 δa = ⎛ ∫(δa)2 dx ⎞ 2, Задача о распространении теплового импульса. При решении модельной задачи тепловой импульс задается в окрестности центра расчетной области, газодинамическая постоянная γ , число Рейнольдса 2,h ⎜ ⎟ Re , число Прандтля Pr , число Маха M∞ и ω имеют ∗ h 3 δa = ai − ai , следующие значения: γ = 1, 4, Re = 2 ⋅10 , Pr = 0, 7, 2 h ⎛ ∗ h 2 ∗ h 2 ⎞ n−1 ∗ h M∞ = 4, ω = 0, 8 [4−6]. Температура T определяется ∫(δa) dx = a − a 2 ⎜ 0 0 + an − an ⎟ + h∑ ai − ai , из следующего уравнения состояния: ⎜ ⎟ ⎝ ⎠ i =1 γM2 P δa = max a∗ − ah . T = ∞ . (38 ∞,h i i i ρ Для ние δσ∞,h и δρ∞,h можно получить соотноше- В качестве начальных условий задаются условия затухания возмущений в бесконечном удалении от источника. Условия на границе для плотности, скоро- δρ = δσ max σ∗ + σh . сти и энергии следующие: ρ = ρ = 1 , ∞,h ∞,h i i i x =0 x =1 Норму δρ можно выразить через 2 h δσi ∞,h сле- u = u = 0 , e = e = 1 x =0 x =1 x =0 x =1 дующим образом: Графики распределения плотности, температуры, ⎛ h ⎛ 2 n −1 2 ⎞ 1 2 ⎞ 2 скорости и давления газа при h = 0, 002 (n = 500), δρ = 2,h ⎜ δρ + δρ ∞,h ∞,h ⎟ + h∑ δρi ∞,h ⎟ , τ = 0, 000 1 представлены на рис. 1–4. Кривые, обозна- ⎜ ⎟ ⎜ 2 ⎝ ⎠ i =1 ⎟ ченные на рисунках как 1, 2, 3, 4, соответствуют рас- ⎝ ⎠ δρ = ⎛ =h ⎡⎛ δσ σ∗ + σh ⎞ 2 + ⎛ δσ σ∗ + σh ⎞ 2 + четным моментам безразмерного времени 2,h ⎜ 2 ⎢⎣⎜ 0 ∞,h 0 0 ⎟ ⎜ n ∞,h n n ⎟ ] t = 0, 01; 0, 05; 0,1; 0, 2. ⎝ ⎝ ⎠ ⎝ ⎠ Таблица 1 h0,050,040,0250,020,01250,01 δσ 2,h0,000 044 10,000 031 90,000 016 10,000 011 70,000 006 10,000 004 6 δρ 2,h0,000 088 20,000 063 70,000 032 10,000 023 30,000 012 20,000 009 2 δε 2,h0,000 033 80,000 027 80,000 019 60,000 017 20,000 014 20,000 013 3 δe 2,h0,000 067 20,000 055 10,000 038 80,000 034 10,000 028 00,000 026 4 δu 2,h0,000 145 00,000 106 50,000 057 50,000 044 10,000 027 50,000 023 0 Таблица 2 h0,050,040,0250,020,01250,01 δσ ∞,h0,000 248 20,000 200 60,000 128 20,000 103 80,000 067 00,000 054 6 δρ ∞,h0,000 495 90,000 400 70,000 256 00,000 207 30,000 133 80,000 109 1 δε ∞,h0,000 211 90,000 192 90,000 164 10,000 154 40,000 139 70,000 134 8 δe ∞,h0,000 420 60,000 382 60,000 324 90,000 305 40,000 276 00,000 266 2 δu ∞,h0,000 839 10,000 692 00,000 469 20,000 394 50,000 281 90,000 244 2 104 Вестник Сибирского государственного аэрокосмического университета имени академика М. Ф. Решетнева ρ Рис. 1 Рис. 2 Рис. 3 Рис. 4 В заключение следует отметить, что полученные системы уравнений удовлетворяют законам сохране- ния массы и полной энергии на дискретном уровне, обеспечивая устойчивость дискретного решения по времени. Замена искомых функций в уравнениях не- разрывности и внутренней энергии обеспечивает по- вышение точности приближенного решения и приво- дит к меньшей абсолютной погрешности в норме L2 и L∞ . Применение комбинации метода траекторий и метода конечных элементов не требует согласова- ния триангуляций на соседних временных слоях, что значительно облегчает динамическое разрежение или сгущение триангуляций по времени для оптимизации вычислительной работы или улучшения аппроксима- ции в пограничных слоях и ударных волнах. Для ре- шения систем алгебраических уравнений использует- ся многосеточный метод с внешними итерациями по нелинейности. Совокупность методов позволяет по- строить экономичный алгоритм с вычислительной точки зрения.