Диссертация (Консервативные дискретизации одномерных квазигазо- и квазигидродинамических систем уравнений), страница 5

1.9: Течение над холмом для тестов (а)–(в): поведение средней полной энергии E (tot)во времени2622.1Глава 2. Пространственная дискретизация одномерной квазигазодинамической системы уравнений с общими уравнениями состояния и уравнение балансаэнтропииОдномерная квазигазодинамическая система уравнений c общими уравнениями состояния и уравнение баланса энтропииКвазигазодинамическая (КГД) система уравнений с одной пространственной переменной в форме [15, 43] состоит из следующих уравнений балансамассы, импульса и полной энергии (в отсутствие массовых сил)∂t ρ + ∂x j = 0,(2.1)∂t (ρu) + ∂x (ju + p) = ∂x Π,(2.2)∂t E + ∂x {(u − w)(E + p)} = −∂x q + ∂x (Πu) + Q.(2.3)Здесь ∂t и ∂x — частные производные по аргументам t > 0 и x ∈ [0, X].Функции ρ > 0, u, E = 0.5ρu2 + ρε — это плотность, скорость и полнаяэнергия газа, а p и ε — это давление и внутренняя энергия.В уравнениях стоят поток массы с плотностью, задаваемой формуламиj = ρ(u − w),ττw = ŵ + u∂x (ρu), ŵ = (ρu∂x u + ∂x p),ρρвязкое напряжение4Π = µ∂x u + ρuŵ + τ u∂x p + ρCs2 ∂x u − (γQ − 1)Q3и тепловой поток, со знаком минус задаваемый формулойp−q = κ∂x θ + τ ρu2 ∂x ε − 2 ∂x ρ − uQ .ρ(2.4)(2.5)(2.6)(2.7)Здесь θ > 0 — абсолютная температура, а µ > 0, κ > 0, τ > 0 — коэффициенты вязкости и теплопроводности и релаксационный параметр, которые в разделах 3.1, 3.2 могут быть произвольными функциями неизвестных.Выражение для Π обобщено на случай общих уравнений состояния согласно[28, 26], причем Cs > 0 — это скорость звука в газе; формулы для Cs > 0 и γQв терминах функций состояния даны ниже.

Величина Q > 0 — это мощностьтепловых источников.27Из этих уравнений следует уравнение баланса внутренней энергии [15, 43,26]∂t (ρε) + ∂x (jε) = −∂x q + Π∂x u − p∂x (u − w) + w∂x p + Q.(2.8)Будем рассматривать общие уравнения состояния газа в формеp = p(ρ, θ), ε = ε(ρ, θ),(2.9)связанные равенством Максвеллаp = θpθ + ρ2 ερ в D0(2.10)и удовлетворяющие условиям термодинамической устойчивости видаpρ > 0, εθ > 0 в D0 ,(2.11)где D0 — область значений пары функций (ρ, θ) в КГД системе (случайD0 = (0, ∞) × (0, ∞) допускается, но не является единственно возможным)и, например, pρ , pθ — частные производные функции p = p(ρ, θ). Если фактически pρ > 0 в D0 , то последние условия можно переписать в известномэквивалентном виде cp > cV ≡ εθ > 0 в D0 , где cp и cV — теплоемкостигаза при постоянном давлении и постоянном объеме [37].

Включение случаяpρ > 0 позволяет охватить одну из простейших моделей двухфазной смесигаз/жидкость.При этом величины Cs > 0 и γQ задаются формуламиCs2 = pρ +θp2θpθ,γ=1+.Qρ2 ε θρεθВведем также энтропию s = S(ρ, ε) посредством формул Гиббсаp1Sρ = − 2 , Sε = .ρθθДля простейшего случая совершенного политропного газа имеемp = (γ − 1)ρε, ε = cV θ, S = S0 − R ln ρ + cV ln ε(2.12)(2.13)(2.14)с постоянными γ > 1, cV > 0 и R = (γ − 1)cV .

Ясно также, что p = Rρθ, аCs2 = γ(γ − 1)ε и γQ = γ.Из уравнений баланса массы (3.1) и внутренней энергии (2.8) можно вывести уравнение баланса энтропии [26] q κ(∂ θ)2 4 µ(∂ u)2 ρŵ2xx∂t (ρs) + ∂x (js) = ∂x − ++++2θθ3θτθ2τ pρτρεθpQQτQθθ∂x u + u∂x θ −+1−. (2.15)+{∂x (ρu)}2 + 2ρθθρεθ2ρεθθ4ρθεθ28Сумма всех слагаемых правой части, кроме дивергентного первого ∂x (−q/θ),представляет собой производство энтропии. При этом сумма второго и третьего слагаемых плюс слагаемое Q/θ представляет собой производство энтропии по Навье-Стоксу, а остальные слагаемые содержат множитель τ иявляются релаксационными. Ясно, что первые пять слагаемых производстваэнтропии всегда неотрицательны, а последнее неотрицательно при выполнении условия τ Q 6 4ρθεθ .

Отметим, что для производства энтропии возможныи другие формы записи [28, 26].Указанное свойство неотрицательности производства энтропии сохраняетсилу и при µ > 0, κ > 0, τ > 0 (для справедливости последнего слагаемоеρŵ2 /(τ θ) надо переписать в виде τ (ρu∂x u + ∂x p)2 /(ρθ)).2.2Применение к решению одномерных уравнений Эйлера реального газаВ отсутствие массовых сил и источников тепла и для λ = 0 КГД системаупрощается и может быть записана в виде∂t ρ + ∂x j = 0,(2.16)∂t (ρu) + ∂x (ju + p) = ∂x Π,j(E + p)∂t E + ∂x= −∂x q + ∂x (Πu)ρ(2.17)(2.18)с величинами E = 12 ρu2 + ρε иj := ρ(u − w), w =τ∂x (ρu2 + p),ρ4Π = µ∂x u + τ ρu2 ∂x u + 2u∂x p + ρCs2 ∂x u ,3p2−q = κ∂x θ + τ ρu ∂x ε − 2 ∂x ρ ,ρ(2.19)(2.20)(2.21)где j — плотность потока массы.Удобно использовать уравнения состояния газа также в альтернативном(эквивалентном) видеp = P (ρ, ε), θ = Θ(ρ, ε).Заметим, что тогда формула для скорости звука принимает видpCs2 = Pρ + 2 Pε .ρ29(2.22)(2.23)Система уравнений Эйлера для невязкого нетеплопроводного газа не содержит диссипативных слагаемых, т.е.

слагаемых с коэффициентами τ , µ иκ. При дискретизации диссипативные слагаемые рассматриваются как искусственные регуляризаторы. Следуя [15], с. 64 и [43], с. 345, свяжем их формуламиµcpαSµ = αS τ ρpρ (ρ, θ), κ ==τ ρcv Cs2 .(2.24)αPαPЗдесь αS и αP — числа Шмидта и Прандтля (положительные постоянные).Параметр релаксации τ выбирается ниже в зависимости от шага пространственной сетки и Cs .Построим явную двухслойную по времени и симметричную по пространству разностную схему. Введем равномерную сетку по пространству с узламиxi = (i − 1)h, 1 6 i 6 N и шагом h = X/(N − 1), где определены основные искомые функции ρ, u и ε, а также p и θ.

Введем также вспомогательную сеткус узлами xi−1/2 = (i − 0.5)h, 1 6 i 6 N − 1, где определены вспомогательные величины j, w, Π и q. Введем симметричные усреднения и центральныеконечные разности (как для основной, так и вспомогательной сеток):va,i−1/2 =wi+1/2 − wi−1/2vi − vi−1vi−1 + vi, δvi−1/2 =, δwi =.2hhЗначения ρ, u, ε и p на вспомогательной сетке вычисляются с использованиемусреднения (·)a .Введем неравномерную сетку по времени 0 = t0 < · · · < tM = tfin с шагами∆tm = tm+1 − tm .Следуя в основном [20], аппроксимируем КГД систему (2.16)–(2.18) двухслойной явной симметричной по пространству разностной схемой с использованием разностей вперед по t и центральных разностей по x:ρ̂ − ρ+ δj = 0,∆tρcu − ρu+ δ(jua + pa ) = δΠ,∆tÊ − EE a + pa= δ(−q + Πua ).+δ j∆tρaВеличины, помеченные “крышкой” b·, относятся к верхнему слою по времени.После вычисления значений ρ̂, ρcu и Ê на верхнем слое по времени полагаемû =ρcuÊ 1, ε̂ = − û2ρ̂ρ̂230в предположении ρ̂ > 0.

Соотношения (2.19)–(2.21) аппроксимируем следующим образом:τaj = ρa (ua − w), w =δ(ρu2 + p),ρa4Π = µa δu + τa ρa u2a δu + 2ua δp + ρa (Cs2 )a δu ,3pa−q = κa δθ + τa ρa ua δε − 2 δρ .ρaПри решении системы уравнений Эйлера величина τ рассматривается какпараметр регуляризации и вычисляется по формулеτ =αh.CsДля выполнения условия устойчивости типа Куранта-Фридрихса-Леви шагпо времени берется в виде∆t = β min16i6Nh.|ui | + Cs,iЗдесь 0 < α < 1 и 0 < β < 1 — числовые параметры, значения которыхопределяются экспериментально для каждого расчета.Применим построенную схему для двух моделей неидеального газа (случайсовершенного политропного газа был подробно протестирован в [20]).

Нижена графиках по оси абсцисс откладывается значение N .1. Сравнительно простая модель с двучленными уравнениями состояниячасто встречается в литературе и описывается следующим образом:ρP (ρ, ε) = (γ − 1)ρε + B−1 ,(2.25)ρ∗где γ > 1, B и ρ∗ — положительные постоянные. Из формулы (2.23) следуетCs2 =γp + BγB (γ − 1)B= γ(γ − 1)ε ++> 0.ρρ∗ρВ отличие от многих других схем, в дополнение к (2.25) нам также понадобится второе уравнение состояния (2.22). Для его получения предположим,что εθ = cv ≡ const > 0, тогда ε(ρ, θ) = cv θ + ε0 (ρ). Следовательно,ρp(ρ, θ) = (γ − 1)ρ[cv θ + ε0 (ρ)] + B−1 ,(2.26)ρ∗31и из равенства Максвелла (2.10) следует обыкновенное дифференциальноеуравнение относительно ε0 (ρ):ρρ2 ε00 (ρ) = (γ − 1)ρε0 (ρ) + B−1 .ρ∗Его общим решением являетсяε0 (ρ) = c0 ργ−1 +BB−γρ (γ − 1)ρ∗с любой постоянной c0 .

Для выполнения свойства ε0 (ρ) > 0 при любом ρ > 0потребуем, чтобы c0 > 0. Тогда точка минимума ε0 (ρ) такова1/γBρmin =γ(γ − 1)c0и таким образомBmin ε0 (ρ) = ε0 (ρmin ) =ρ>0(γ − 1)1ρmin1−ρ∗.Вследствие этого c0 должно удовлетворять неравенствуc0 >B,γ(γ − 1)ργ∗и в итоге полагаем c0 = B/(γ(γ − 1)ργ∗ ).Из формулы (2.26) следует равенствоpρ (ρ, θ) = (γ − 1)cv θ + c0 γ(γ − 1)ργ−1 ,и условия термодинамической устойчивости (2.11) выполнены для всех θ > 0и ρ > 0, если cv > 0 и c0 > 0.Отметим, что поскольку коэффициент κ содержит множитель cv , см.

(2.24),для cv = const, поэтому в здесь и ниже функция cv θ может быть легко исключена из системы уравнений.Для двучленной модели рассмотрим 3 теста из [54] с различными конфигурациями течений газа. Пусть B = ρ∗ = 1, γ = 1.4 и X = 1. В тесте 1начальные значения ρ, u и p таковы:((1при x 6 0.58 при x 6 0.5ρ0 (x) =, u0 (x) = 0, p0 (x) =0.125 при x > 0.50.1 при x > 0.532и tfin = 0.075. В тесте 2 начальные значения таковы:((0.8 при x 6 0.510 при x 6 0.5ρ0 (x) =, u0 (x) = 0, p0 (x) =1 при x > 0.50.1 при x > 0.5и tfin = 0.1.