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

Наконец, начальные значения теста 3 имеют вид(−2 при x 6 0.5ρ0 (x) = 1, u0 (x) =, p0 (x) = 82 при x > 0.5и tfin = 0.07. Требуемые значения ε0 (x) выражаются из (2.25).Для всех тестов 1–3 берутся простейшие значения αS = αP = 1 и N =400, α = 0.2 и β = 0.1. Основные искомые функции ρ, u и ε, а такжеp (пространственные графики на момент времени tfin и соответствующиепространственно-временные графики) представлены на рис.

2.1, 2.2, рис. 2.3,2.4 и рис. 2.5, 2.6 для тестов 1, 2 и 3 соответственно. Результаты на рис. 2.1,2.3, 2.5 находятся в хорошем соответствии с полученными в [54] при N = 200.2. Классическая модель газа Ван-дер-Ваальса в форме (2.22) задается следующим образом:P (ρ, ε) =ε + aρ(γ − 1)ρ(ε + aρ)− aρ2 , Θ(ρ, ε) =1 − bρcvс γ = cRv +1 (см., например, [37]), где a, b, R и cv — положительные физическиепостоянные.

Предполагается, что 0 < ρ < b−1 . В соответствии с формулой(2.23) имеемγ(γ − 1)(ε + aρ)− 2aρ.Cs2 =(1 − bρ)2Заметим также, чтоpρ (ρ, θ) =Rθ(γ − 1)(ε + aρ)−2aρ=− 2aρ,(1 − bρ)2(1 − bρ)2и рассматривается только случай термодинамической устойчивости pρ (ρ, θ) >0 (фазовые переходы отсутствуют).Для газа Ван-дер-Ваальса рассмотрим тест 5 из [48].

Значения константтаковы:a = 1684.54, b = 0.001692, R = 461.5, cv = 1401.88.Также X = 10, а начальные данные имеют следующий вид:((250 при x 6 535966778 при x 6 5ρ0 (x) =, u0 (x) = 0, p0 (x) =166.6 при x > 527114795 при x > 5.33Значение времени tfin в [48] отсутствует, поэтому наше значение tfin = 0.005подобрано для получения наиболее близких результатов.Для этого теста снова берутся αS = αP = 1, но теперь N = 700, α = 0.65 иβ = 0.1. Основные искомые функции ρ, u и θ, а также p показаны на рис. 2.7 и2.8. Результаты на рис.

2.7 находятся в хорошем соответствии с полученнымив [48] при N = 500.Полученные результаты подтверждают эффективность построенной разностной схемы.34ρu3120.500110020030040000100200300400300400εp924226201831600141002003004000100200Рис. 2.1: Тест 1 (двучленная модель): пространственные графики для tfin = 0.075Рис. 2.2: Тест 1 (двучленная модель): пространственно-временные графики35ρu23211010020030040000100200300400300400εp401083062041020001002003004000100200Рис. 2.3: Тест 2 (двучленная модель): пространственные графики для tfin = 0.1Рис. 2.4: Тест 2 (двучленная модель): пространственно-временные графики36ρu2110−10.50100200300400−20100200300400300400εp2081961817416201002003004000100200Рис. 2.5: Тест 3 (двучленная модель): пространственные графики для tfin = 0.07Рис. 2.6: Тест 3 (двучленная модель): пространственно-временные графики37ρu26040240302202020010180160035070000700θp7x 103507207103.570036902.501002003004005006007006800350700Рис.

2.7: Тест 4 (модель Ван-дер-Ваальса): пространственные графики для tfin = 0.005Рис. 2.8: Тест 4 (модель Ван-дер-Ваальса): пространственно-временные графики382.3Пространственная дискретизация одномерной квазигазодинамической системы уравнений и дискретное уравнение балансаэнтропии2.1. Рассмотрим теперь более сложную дискретизацию КГД системы (2.1)–(2.3).

Введем произвольную неравномерную сетку ω̄h на [0, X] с узлами 0 =x0 < x1 < · · · < xN = X и шагами hi = xi − xi−1 . Пусть hmax = max16i6N hi .Введем также вспомогательную сетку ωh∗ с узлами xi+1/2 = (xi + xi+1 )/2,0 6 i 6 N − 1 и шагами ĥi = xi+1/2 − xi−1/2 = (hi + hi+1 )/2.На функциях v, заданных на ω̄h , и y, заданных на ωh∗ , введем операторысеточного усреднения, сдвига аргумента и разностные отношения[v]i+1/2 = 0.5(vi + vi+1 ), (v± )i+1/2 = vi±1 , δvi+1/2 =[y]∗i =hi yi−1/2 + hi+1 yi+1/22ĥi, δ ∗ yi =vi+1 − vi,hi+1yi+1/2 − yi−1/2ĥi.Операторы [·], (·)± , δ переводят функции, заданные на ω̄h , в функции, заданные на ωh∗ , а операторы [·]∗ , δ ∗ — функции, заданные на ωh∗ , в функции,заданные на ωh = {xi ; 1 6 i 6 N − 1}.Нам потребуется набор формул, связывающих эти сеточные операторыδ(uv) = δu · [v] + [u]δv,δ∗δ ∗ (y[v]) = δ ∗ y · v + [yδv]∗ ,[u][v] − 0.25h2+ δu · δv = δ ∗ [u] · v + uδ ∗ [v],[y]∗ v = [y[v]]∗ − 0.25δ ∗ (h2+ yδv),(2.27)(2.28)(2.29)(2.30)где функция u задана на ω̄h (см., в частности, [31]).

Здесь и ниже для сокращения количества скобок предполагается, что, например, δu · [v] = (δu)[v](т.е. знак умножения · прекращает действие предыдущих операторов слева).Обратим также внимание на формулу [δv]∗ = δ ∗ [v].2.2. Построим следующую дискретизацию по пространству уравнений баланса массы, импульса и полной энергии (2.1)–(2.3)∂t ρ + δ ∗ j = 0,(2.31)∂t (ρu) + δ ∗ (j[u] + [p]) = δ ∗ Π,(2.32)∂t E + δ ∗ ([u] − w)([E]2 + [p]) − 0.25h2+ δu · δp = δ ∗ (−q + Π[u]) + [Q]∗(2.33)39на ωh , где полная энергия, давление и внутренняя энергия имеют стандартный видE = 0.5ρu2 + ρε, p = p(ρ, θ), ε = ε(ρ, θ).(2.34)Будем использовать следующие дискретизации плотности потока массыj = [ρ]0 ([u] − w),ττ[u]δ(ρu), ŵ =([ρ][u]δu + δp),w = ŵ +[ρ][ρ]вязкого напряжения4Π = µδu + [ρ][u]ŵ + τ [u]δp + [ρ][Cs2 ]4 δu − ([γQ ]5 − 1)Q3и теплового потока (со знаком минус)[p]1−q = κδθ + τ [ρ][u]2 δε − 2 δρ − [u]Q .[ρ](2.35)(2.36)(2.37)(2.38)Основные неизвестные функции ρ, u, E и функции p, ε, θ определены наосновной сетке ω̄h , причем предполагается, что ρ > 0, ε > 0.

Функции j, w, ŵ,Π, q, τ , µ, κ, Q определены на вспомогательной сетке ωh∗ . Такая операторнаязапись дискретизации КГД системы уравнений следует [31].В ней наряду с простейшими усреднениями [ρ], [u], [p] применяются и другие дискретизации величин ρ, p, E, ε, Cs2 , γQ на ωh∗ :[ρ]0 , [p]1 , [E]2 = 0.5[ρ]0 u− u+ + [ρ]0 [ε]3 , [ε]3 , [Cs2 ]4 , [γQ ]5 .Вид дискретизации [E]2 , как и добавок −0.25h2+ δu · δp в уравнении (2.33),следуют [31] и обеспечивают естественную форму дискретного уравнения баланса внутренней энергии (см. ниже).

Остальные аппроксимации пока произвольны, и в ходе дальнейшего анализа станет ясна их роль в записи дискретного уравнения баланса энтропии, содержащего дивергентные и недивергентные сеточные дисбалансные слагаемые. Более того, будет указан иходнозначный выбор, гарантирующий отсутствие недивергентных дисбалансных слагаемых и неотрицательность производства энтропии.Этот выбор окажется достаточно нетривиальным. Чтобы его описать, длятермодинамических функций ψ = ψ(ρ, θ) введем разделенные разности иусреднения по аргументам( ψ(ρ ,θ)−ψ(ρ ,θ)( ψ(ρ,θ )−ψ(ρ,θ )+−+−приρ=6ρпри θ− 6= θ+−+ρ+ −ρ−θ+ −θ−dρ ψ =, dθ ψ =,ψρ (ρ− , θ)при ρ− = ρ+ψθ (ρ, θ− )при θ− = θ+aρ ψ = 0.5 {ψ(ρ− , θ) + ψ(ρ+ , θ)} , aθ ψ = 0.5 {ψ(ρ, θ− ) + ψ(ρ, θ+ )} .40Следующая несложная формула для разностного отношения суперпозициифункций играет ниже важную рольδψ(ρ, θ) = (aθ dρ ψ)δρ + (aρ dθ ψ)δθ.(2.39)Обратим также внимание на интегральные представленияZ 1Z 1(α)dρ ψ = hψρ iρ ≡ψρ ([ρ] , θ) dα, dθ ψ = hψθ iθ ≡ψθ (ρ, [θ](α) ) dα, (2.40)00где, например, [ρ](α) = (1 − α)ρ− + αρ+ (в частности, [ρ](1/2) = [ρ]).2.3.

Вывод уравнения баланса энтропии базируется на уравнениях балансамассы и внутренней энергии. Хотя вывод последнего фактически повторяетданный в [31], кратко воспроизведем его для полноты рассуждений. Сначалаумножим уравнение баланса импульса (2.32) на u. Воспользуемся формулой∂t (ρu) · u = 0.5∂t (ρu2 ) + 0.5∂t ρ · u2 ,применим уравнение баланса массы (2.31) и дважды — формулу разностногодифференцирования произведения (2.28) и получим∂t ρ · u2 = −δ ∗ j · u2 = −δ ∗ (j[u2 ]) + [jδ(u2 )]∗ ,δ ∗ (j[u]) · u = δ ∗ (j[u]2 ) − [j[u]δu]∗ .Далее с учетом элементарных формул [u]2 = 0.5[u2 ] + 0.5u− u+ и 0.5δ(u2 ) =[u]δu выведем уравнение баланса кинетической энергии0.5∂t (ρu2 ) + 0.5δ ∗ (ju− u+ ) + δ ∗ [p] · u = δ ∗ Π · u.(2.41)Вычтем из уравнения баланса полной энергии (2.33) последнее уравнение.По формулам разностного дифференцирования (2.28), (2.29) имеемδ ∗ [u][p] − 0.25h2+ δu · δp = δ ∗ [u] · p + δ ∗ [p] · u,(2.42)δ ∗ (w[p]) = δ ∗ w · p + [wδp]∗ , δ ∗ (Π[u]) = δ ∗ Π · u + [Πδu]∗ .(2.43)В итоге верно следующее дискретное уравнение баланса внутренней энергии∂t (ρε) + δ ∗ (j[ε]3 ) = −δ ∗ q + [Πδu]∗ − pδ ∗ ([u] − w) + [wδp]∗ + [Q]∗ .(2.44)Оно является аналогом дифференциального уравнения (2.8).Приступим к выводу уравнения баланса энтропии.

Согласно формуламГиббса (2.13) имеемp11p∂t (ρs) = ∂t ρ · s + ρ − 2 ∂t ρ + ∂t ε = ∂t ρ · s +− ∂t ρ + ∂t (ρε) − ε∂t ρ .ρθθθρ41Формула δ ∗ (j[s]) = δ ∗ j · s + [jδs]∗ и уравнение баланса массы (2.31) даютp1∂t (ρε) + δ ∗ j + εδ ∗ j + [jδs]∗ .(2.45)∂t (ρs) + δ ∗ (j[s]) =θρНиже полагаем, что s = s(ρ, θ) для единообразия с уравнениями состояния(2.9). По формулам Гиббса (2.13) и равенству Максвелла (2.10) имеемsρ = −pθεθ,s=.θρ2θ(2.46)Согласно формуле (2.39) верно равенствоδs = (aθ dρ s)δρ + (aρ dθ s)δθ.(2.47)Воспользуемся им в (2.45) и преобразуем результат к виду∂t (ρs) + δ ∗ (j[s]) =1{∂t (ρε) + δ ∗ (j[ε]3 ) + pδ ∗ ([u] − w)} + Dh ,θ(2.48)где имеющий чисто сеточное происхождение дисбаланс Dh таковDh :=p ∗pε1δ j − δ ∗ ([u] − w) + δ ∗ j − δ ∗ (j[ε]3 ) + [j {(aθ dρ s)δρ + (aρ dθ s)δθ}]∗ .ρθθθθНиже нам потребуются элементарные формулы 1δv1[v],=,δ =−vv− v+vv− v+(2.49)где v± 6= 0.

Применив четырежды формулу (2.28), перепишем Dh в виде h i hpipε1Dh = δ ∗ j− ([u] − w)+j− [ε]3+ρθθθθ∗ppε1+ −jδ + ([u] − w)δ + j −δ + [ε]3 δ+ j {(aθ dρ s)δρ + (aρ dθ s)δθ} .ρθθθθВ силу формулы (2.27) запишем h ip 1p 1pδ =δ ·+δ .ρθθ ρθ ρКроме того, согласно формуле (2.39) верны равенстваp ppεεεδ = aθ dρ δρ + aρ dθ δθ, δ = aθ dρ δρ + aρ dθ δθ.θθθρθθθ42(2.50)Теперь, применив также формулы (2.35) и (2.49), выводим представлениеDh = δ ∗ B1h + [([u] − w)(Dhρ δρ + Dhθ δθ)]∗(2.51)с сеточными дисбалансными слагаемыми и множителями h ih i ppε1B1h := ([u] − w) [ρ]0−+j− [ε]3,(2.52)ρθθθθ [ρ]0 h p ipp1εaθ dρ ++ aθ dρ − [ρ]0 aθ dρ + [ρ]0 aθ dρ s =Dhρ := −[ρ]0ρθ ρ − ρ+ θθθ ε11 hpi1aθdρ p + [ρ]0+ aθ d ρ s −,(2.53)= 1 − [ρ]0ρθρ− ρ+ θθ 1pε[ε]3Dhθ := −[ρ]0+ 1 aρ dθ − [ρ]0 aρ dθ − [ρ]0+ [ρ]0 aρ dθ s =ρθθθ− θ+ pε1[ε]3= 1 − [ρ]0aρ dθ + [ρ]0 aρ dθ s −−.(2.54)ρθθθ− θ+Отметим, что B1h не обращается pв 0 при простейших ε 1[ρ] 0 = [ρ], [ε]3 = [ε]p(хотя это и происходит при [ρ]0 = ρθ / θ , [ε]3 = θ / θ ).