学习自李新亮的《计算流体力学讲义PPT》第2讲 双曲型方程组及间断解
1 问题描述
一维无粘流动初始间断的演化问题,激波管问题(Sod激波管问题),由于该问题属于典型的间断问题,且有精确解存在,故广泛用于对比验证CFD中离散格式、数值方法的准确性,意义嘛还是蛮重要的,哈。
一无限长管道内部充满理想气体(无粘),中间有一不计厚度的隔膜,初始时刻 t 0 t_0 t0,左侧气体的速度、密度、压力分别为 u 1 u_1 u1、 ρ 1 \rho_1 ρ1、 p 1 p_1 p1,右侧气体的速度、密度、压力分别为 u 2 u_2 u2、 ρ 2 \rho_2 ρ2、 p 2 p_2 p2。不考虑流体粘性,突然将隔膜抽出(也可认为隔膜突然消失),那么到时刻 t t t后,管道内流体的速度、密度、压力将如何分布?
---------------------------------------------------------
u 1 , ρ 1 , p 1 u_1, \rho_1,p_1 u1,ρ1,p1 | u 2 , ρ 2 , p 2 u_2, \rho_2, p_2 u2,ρ2,p2
---------------------------------------------------------
控制方程为一维Euler方程(即无粘可压缩理想流体一维流动的控制方程):
{ ∂ ρ ∂ t + ∂ ( ρ u ) ∂ x = 0 ∂ ( ρ u ) ∂ t + ∂ ( ρ u 2 + p ) ∂ x = 0 ∂ ( ρ E ) ∂ t + ∂ ( ρ E u + p u ) ∂ x = 0 \begin{cases} \displaystyle \frac{\partial{\rho}}{\partial{t}} + \frac{\partial{(\rho u)}}{\partial x}=0 \\ \\ \displaystyle \frac{\partial{(\rho u)}}{\partial{t}} + \frac{\partial{(\rho u^2 + p)}}{\partial x}=0 \\ \\ \displaystyle \frac{\partial{(\rho E)}}{\partial{t}} + \frac{\partial{(\rho E u + pu)}}{\partial x}=0 \end{cases} ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧∂t∂ρ+∂x∂(ρu)=0∂t∂(ρu)+∂x∂(ρu2+p)=0∂t∂(ρE)+∂x∂(ρEu+pu)=0
初始条件:
t = t 0 时 , ( u , ρ , p ) = { u 1 , ρ 1 , p 1 x < 0 u 2 , ρ 2 , p 2 x ≥ 0 t=t_0时,\quad\quad (u,\rho,p)=\begin{cases} u_1,\rho_1,p_1 &x<0\\ u_2,\rho_2,p_2 &x\geq0\\ \end{cases} t=t0时,(u,ρ,p)={ u1,ρ1,p1u2,ρ2,p2x<0x≥0
另外,求解方程需要用到的变量为理想气体比热比 γ = 1.4 \gamma=1.4 γ=1.4,气体守恒能量 ρ E \rho E ρE与气体其他状态参数的关系式
ρ E = p γ − 1 + ρ u 2 2 \displaystyle \rho E=\frac{p}{\gamma-1}+\rho\frac{u^2}{2} ρE=γ−1p+ρ2u2
声速与气体其他状态参数的关系式
c = γ p ρ \displaystyle c=\sqrt{\gamma \frac{p}{\rho}} c=γρp
2 求解方法
流场中可能出现三种波:
间断类型 | 特征 |
---|---|
激波 | 强间断,满足R-H关系式 |
接触间断 | 特殊间断,仅密度突变,而速度和压力不变 |
膨胀波 | 等熵波 |
依据不同的初始条件,可以产生五类不同的情况,分别为
(1)左激波-接触间断-右激波
(2)左膨胀波-接触间断-右激波
(3)左激波-接触间断-右膨胀波
(4)左膨胀波-接触间断-右膨胀波
(5)左膨胀波-右膨胀波
2.1 情况1,左右激波和中间的接触间断
区域1和区域2中物理量与初始状态一致,分别为 u 1 u_1 u1、 ρ 1 \rho_1 ρ1、 p 1 p_1 p1和 u 2 u_2 u2、 ρ 2 \rho_2 ρ2、 p 2 p_2 p2;
待求量为左激波运动速度 Z 1 Z_1 Z1、右激波运动速度 Z 2 Z_2 Z2,注意左激波是向左运动的,右激波是向右运动的,这里把 Z 1 Z_1 Z1和 Z 2 Z_2 Z2的正方向都定义为沿着 x x x轴的正方向。
区域3和区域4中物理量也为待求量,由于接触间断两侧仅密度产生突变,故区域3和区域4中的速度和压力时相同的,将其定义为 u ∗ u^* u∗和 p ∗ p^* p∗,将区域3和区域4的密度分别定义为 ρ ∗ L \rho^{*L} ρ∗L和 ρ ∗ R \rho^{*R} ρ∗R。注意,中间接触间断面的运动速度并不是未知量,因为该间断面两侧速度相同,都为 u ∗ u^* u∗,所以该间断面的运动速度也为 u ∗ u^* u∗。
那么总共有6个未知量,即 Z 1 Z_1 Z1、 Z 2 Z_2 Z2、 u ∗ u^* u∗、 p ∗ p^* p∗、 ρ ∗ L \rho^{*L} ρ∗L和 ρ ∗ R \rho^{*R} ρ∗R。
在激波两侧的物理量满足R-H关系,即质量、动量、能量通量守恒,即
{ ρ 1 ( u 1 − Z ) = ρ 2 ( u 2 − Z ) ρ 1 u 1 ( u 1 − Z ) + p 1 = ρ 2 u 2 ( u 2 − Z ) + p 2 ρ 1 E 1 ( u 1 − Z ) + u 1 p 1 = ρ 2 E 2 ( u 2 − Z ) + u 2 p 2 \begin{cases} \displaystyle \rho_1(u_1-Z)=\rho_2(u_2-Z) \\ \displaystyle \rho_1 u_1 (u_1-Z) + p_1=\rho_2 u_2 (u_2-Z) + p_2 \\ \displaystyle \rho_1 E_1 (u_1-Z) + u_1 p_1=\rho_2 E_2 (u_2-Z) + u_2 p_2 \end{cases} ⎩⎪⎨⎪⎧ρ1(u1−Z)=ρ2(u2−Z)ρ1u1(u1−Z)+p1=ρ2u2(u2−Z)+p2ρ1E1(u1−Z)+u1p1=ρ2E2(u2−Z)+u2p2
那么,对于1-3两区
{ ρ 1 ( u 1 − Z 1 ) = ρ ∗ L ( u ∗ − Z 1 ) ( 1 ) ρ 1 u 1 ( u 1 − Z 1 ) + p 1 = ρ ∗ L u ∗ ( u ∗ − Z 1 ) + p ∗ ( 2 ) ρ 1 E 1 ( u 1 − Z 1 ) + u 1 p 1 = ρ ∗ L E ∗ L ( u ∗ − Z 1 ) + u ∗ p ∗ ( 3 ) \begin{cases} \displaystyle \rho_1(u_1-Z_1)=\rho^{*L}(u^*-Z_1) & (1) \\ \displaystyle \rho_1 u_1 (u_1-Z_1) + p_1=\rho^{*L} u^* (u^*-Z_1) + p^* & (2)\\ \displaystyle \rho_1 E_1 (u_1-Z_1) + u_1 p_1=\rho^{*L} E^{*L} (u^*-Z_1) + u^* p^* & (3) \end{cases} ⎩⎪⎨⎪⎧ρ1(u1−Z1)=ρ∗L(u∗−Z1)ρ1u1(u1−Z1)+p1=ρ∗Lu∗(u∗−Z1)+p∗ρ1E1(u1−Z1)+u1p1=ρ∗LE∗L(u∗−Z1)+u∗p∗(1)(2)(3)
对于2-4两区
{ ρ 2 ( u 2 − Z 2 ) = ρ ∗ R ( u ∗ − Z 2 ) ( 4 ) ρ 2 u 2 ( u 2 − Z 2 ) + p 2 = ρ ∗ R u ∗ ( u ∗ − Z 2 ) + p ∗ ( 5 ) ρ 2 E 2 ( u 2 − Z 2 ) + u 2 p 2 = ρ ∗ R E ∗ R ( u ∗ − Z 2 ) + u ∗ p ∗ ( 6 ) \begin{cases} \displaystyle \rho_2(u_2-Z_2)=\rho^{*R}(u^*-Z_2) & (4)\\ \displaystyle \rho_2 u_2 (u_2-Z_2) + p_2=\rho^{*R} u^* (u^*-Z_2) + p^* & (5)\\ \displaystyle \rho_2 E_2 (u_2-Z_2) + u_2 p_2=\rho^{*R} E^{*R} (u^*-Z_2) + u^* p^* & (6) \end{cases} ⎩⎪⎨⎪⎧ρ2(u2−Z2)=ρ∗R(u∗−Z2)ρ2u2(u2−Z2)+p2=ρ∗Ru∗(u∗−Z2)+p∗ρ2E2(u2−Z2)+u2p2=ρ∗RE∗R(u∗−Z2)+u∗p∗(4)(5)(6)
总共6个方程(1)-(6),6个未知数,可解。
解法如下:
先来看1-3区的关系式(1)-(3),共有4个未知数 ρ ∗ L \rho^{*L} ρ∗L、 u ∗ u^* u∗、 Z 1 Z_1 Z1、 p ∗ p^* p∗,注意 E ∗ L E^{*L} E∗L是其他状态参数的函数,不属于未知量,其与其他量关系为:
ρ ∗ L E ∗ L = p ∗ γ − 1 + ρ ∗ L ( u ∗ ) 2 2 ( 7 ) \displaystyle \rho^{*L} E^{*L}=\frac{p^*}{\gamma-1}+\rho^{*L}\frac{(u^*)^2}{2} \quad (7) ρ∗LE∗L=γ−1p∗+ρ∗L2(u∗)2(7)
接下来要做的是从(1)-(3)中消去 ρ ∗ L \rho^{*L} ρ∗L和 Z 1 Z_1 Z1,找到 u ∗ u^* u∗和 p ∗ p^* p∗的关系式,先由(1)可得
ρ ∗ L = ρ 1 ( u 1 − Z 1 ) ( u ∗ − Z 1 ) ( 8 ) \displaystyle \rho^{*L} = \frac{\rho_1(u_1-Z_1)}{(u^*-Z_1)} \quad (8) ρ∗L=(u∗−Z1)ρ1(u1−Z1)(8)
将上式(8)代入到式(2),消去 ρ ∗ L \rho^{*L} ρ∗L,可得
ρ 1 ( u 1 − Z 1 ) ( u 1 − u ∗ ) = p ∗ − p 1 ( 9 ) \rho_1 (u_1-Z_1) (u_1-u^*) = p^* - p_1 \quad (9) ρ1(u1−Z1)(u1−u∗)=p∗−p1(9)
将式(8)代入式(3),消去 ρ ∗ L \rho^{*L} ρ∗L,并将能量 E 1 E_1 E1和 E ∗ L E^{*L} E∗L用压力、密度和速度来表示,可得
( p 1 γ − 1 + ρ 1 u 1 2 2 ) ( u 1 − Z 1 ) + u 1 p 1 = p ∗ γ − 1 ( u ∗ − Z 1 ) + ρ 1 ( u ∗ ) 2 2 ( u 1 − Z 1 ) + u ∗ p ∗ \displaystyle \left( \frac{p_1}{\gamma - 1}+\rho_1\frac{u_1^2}{2}\right) (u_1-Z_1)+u_1p_1=\frac{p^*}{\gamma - 1}(u^*-Z_1)+\rho_1\frac{(u^*)^2}{2}(u_1-Z_1) +u^*p^* (γ−1p1+ρ12u12)(u1−Z1)+u1p1=γ−1p∗(u∗−Z1)+ρ12(u∗)2(u1−Z1)+u∗p∗
即
[ p 1 γ − 1 + ρ 1 u 1 2 − ( u ∗ ) 2 2