一维完全气体欧拉方程的引入
在这里我们将利用LF格式,MacCormack格式和Roe格式对一维完全气体欧拉方程进行数值求解。
{
(
ρ
ρ
u
E
)
t
+
(
ρ
u
ρ
u
2
+
p
u
(
E
+
p
)
)
x
=
0
,
p
=
(
γ
−
1
)
(
E
−
1
2
ρ
u
2
)
,
γ
=
1.4.
\left\{\begin{array}{c} \left(\begin{array}{c} \rho \\ \rho u \\ E \end{array}\right)_{t}+\left(\begin{array}{c} \rho u \\ \rho u^{2}+p \\ u(E+p) \end{array}\right)_{x}=0 ,\\ p=(\gamma-1)\left(E-\frac{1}{2} \rho u^{2}\right), \gamma=1.4. \end{array}\right.
⎩
⎨
⎧
ρρuE
t+
ρuρu2+pu(E+p)
x=0,p=(γ−1)(E−21ρu2),γ=1.4.
这里我们求解Riemann问题,我们将针对以下两个初始条件分别求解方程:
∙
\bullet
∙第一个初始条件参考
U
(
x
,
0
)
=
(
ρ
ρ
u
E
)
(
x
,
0
)
=
{
(
1
,
0
,
2.5
)
T
,
x
<
0.3
,
(
0.125
,
0
,
0.25
)
T
,
x
>
0.3.
\boldsymbol{U}(x, 0)=\left(\begin{array}{c} \rho \\ \rho u \\ E \end{array}\right)(x, 0)=\left\{\begin{array}{ll} (1,0,2.5)^{T}, & x<0.3, \\ (0.125,0,0.25)^{T}, & x>0.3. \end{array}\right.
U(x,0)=
ρρuE
(x,0)={(1,0,2.5)T,(0.125,0,0.25)T,x<0.3,x>0.3.
计算区间为
[
0
,
1
]
[0,1]
[0,1],输出时刻为
t
=
0.2
t = 0.2
t=0.2。
∙
\bullet
∙第二个初始条件参考:
U
(
x
,
0
)
=
(
ρ
u
P
)
(
x
,
0
)
=
{
(
3.857143
,
2.629369
,
10.33333
)
T
,
x
<
−
4
,
(
1
+
0.2
s
i
n
(
5
x
)
,
0
,
1
)
T
,
x
>
−
4.
\boldsymbol{U}(x, 0)=\left(\begin{array}{c} \rho \\ u \\ P \end{array}\right)(x, 0)=\left\{\begin{array}{ll} (3.857143,2.629369,10.33333)^{T}, & x<-4, \\ (1 + 0.2sin(5x),0,1)^{T}, & x>-4. \end{array}\right.
U(x,0)=
ρuP
(x,0)={(3.857143,2.629369,10.33333)T,(1+0.2sin(5x),0,1)T,x<−4,x>−4.
计算区间为
[
−
5
,
5
]
[-5,5]
[−5,5],其中在
x
=
±
5
x=\pm 5
x=±5的边界上有
∂
x
ρ
=
∂
x
u
=
∂
x
P
=
0
\partial_x \rho = \partial_x u = \partial_x P = 0
∂xρ=∂xu=∂xP=0,输出时刻为
t
=
1.8
t = 1.8
t=1.8。
上面的方程是关于守恒量的,为了方便阅读,我们引入记号
φ
=
ρ
u
\varphi = \rho u
φ=ρu,同时代入
p
=
(
γ
−
1
)
(
E
−
1
2
ρ
u
2
)
p=(\gamma-1)(E-\frac{1}{2} \rho u^{2})
p=(γ−1)(E−21ρu2),有
ρ
u
2
+
p
=
(
γ
−
1
)
E
+
3
−
γ
2
φ
2
ρ
.
\rho u^2 + p= (\gamma - 1)E + \frac{3-\gamma}{2}\frac{\varphi^2}{\rho}.
ρu2+p=(γ−1)E+23−γρφ2.
u
(
E
+
p
)
=
γ
E
φ
ρ
+
1
−
γ
2
φ
3
ρ
2
.
u(E + p) = \gamma E \frac{\varphi}{\rho} + \frac{1 - \gamma}{2}\frac{\varphi^3}{\rho^2}.
u(E+p)=γEρφ+21−γρ2φ3.
{
(
ρ
φ
E
)
t
+
(
φ
(
γ
−
1
)
E
+
3
−
γ
2
φ
2
ρ
γ
E
φ
ρ
+
1
−
γ
2
φ
3
ρ
2
)
x
=
0
,
\left\{\begin{array}{c} \left(\begin{array}{c} \rho \\ \varphi \\ E \end{array}\right)_{t}+\left(\begin{array}{c} \varphi \\ (\gamma - 1)E + \frac{3-\gamma}{2}\frac{\varphi^2}{\rho} \\ \gamma E \frac{\varphi}{\rho} + \frac{1 - \gamma}{2}\frac{\varphi^3}{\rho^2} \end{array}\right)_{x}=0 ,\\ \end{array}\right.
⎩
⎨
⎧
ρφE
t+
φ(γ−1)E+23−γρφ2γEρφ+21−γρ2φ3
x=0,
我们给出实验观测量
(
ρ
,
u
,
P
)
T
(\rho,u,P)^T
(ρ,u,P)T的结果,下面我们先把上面两个问题的初值函数图像展示出来。
对于下面即将描述的两个数值格式,特别需要注意的是步长的选取,这里需要根据CFL条件随时调整步长,具体的做法如下,下面我们介绍的格式在空间维度都是固定的步长,只需要考虑时间步长的选取。
对于方程\ref{qiti},我们给出
A
(
U
)
=
∂
F
(
U
)
∂
U
A(U) = \frac{\partial F(U)}{\partial U}
A(U)=∂U∂F(U)表达式\ref{A}和对应的三个特征向量
λ
1
=
u
+
a
,
λ
2
=
u
,
λ
3
=
u
−
a
\lambda_1 = u + a,\lambda_2 = u,\lambda_3 = u - a
λ1=u+a,λ2=u,λ3=u−a,由此我们可以知道矩阵A的谱半径为
∣
u
∣
+
a
|u| + a
∣u∣+a,其中
a
=
(
γ
P
ρ
)
a = \sqrt(\frac{\gamma P}{\rho})
a=(ργP)。根据CFL条件,步长的选取应该满足
τ
h
ρ
(
A
)
<
μ
\frac{\tau}{h}\rho(A) <\mu
hτρ(A)<μ,其中
μ
\mu
μ为CFL数,且小于1。具体的做法就是在计算第
m
+
1
m + 1
m+1层的函数值时,我们需要设定好步长
τ
m
<
μ
h
÷
(
m
a
x
j
∣
u
j
m
∣
+
a
j
m
)
\tau_m < \mu h\div (max_{j}|u_j^m| + a_{j}^m)
τm<μh÷(maxj∣ujm∣+ajm)。
A
(
U
)
=
(
0
1
0
γ
−
3
2
u
2
(
3
−
γ
)
u
γ
−
1
γ
−
2
2
u
3
−
a
2
u
γ
−
1
3
−
2
γ
2
u
2
+
a
2
γ
−
1
γ
u
)
A(U)=\left(\begin{array}{ccc} 0 & 1 & 0 \\ \frac{\gamma-3}{2} u^{2} & (3-\gamma) u & \gamma-1 \\ \frac{\gamma-2}{2} u^{3}-\frac{a^{2} u}{\gamma-1} & \frac{3-2 \gamma}{2} u^{2}+\frac{a^{2}}{\gamma-1} & \gamma u \end{array}\right)
A(U)=
02γ−3u22γ−2u3−γ−1a2u1(3−γ)u23−2γu2+γ−1a20γ−1γu
LF
LF格式参考公式\ref{LF},对于LF格式\ref{LF}而言,我们计算第
m
+
1
m + 1
m+1层的函数值时,根据第
m
m
m层的函数值和CFL条件算出
τ
m
\tau_m
τm,与此同时还需要把边界部分特别拿出来考虑,不妨设空间上一共有N个离散点
x
0
<
…
<
x
N
−
1
x_0 < \ldots < x_{N - 1}
x0<…<xN−1,当我们计算内部的
N
−
2
N - 2
N−2个离散点函数时直接代入LF格式\ref{LF}即可,但是我们计算边界上,比如
j
=
0
,
U
0
m
+
1
j = 0,U_0^{m+1}
j=0,U0m+1时会遇到
U
−
1
m
U_{-1}^{m}
U−1m 的情况,这里就需要做特殊处理。我们事先划定空间区域
[
a
.
b
]
⊇
[
0
,
1
]
[a.b]\supseteq [0,1]
[a.b]⊇[0,1],根据常系数对流方程的知识,我们知道计算数值解的时候,应该需要根据特征线来选取足够大的计算区域
[
a
,
b
]
[a,b]
[a,b]使得
[
a
,
b
]
⊇
[
x
−
a
t
,
x
+
a
t
]
,
∀
x
∈
[
0
,
1
]
,
t
∈
[
0
,
t
m
a
x
]
[a,b] \supseteq[x - at,x + at],\forall x \in [0,1],t \in [0,t_max]
[a,b]⊇[x−at,x+at],∀x∈[0,1],t∈[0,tmax]成立,其中
a
a
a为特征线斜率。在求解上面这个方程的时候同样要有类似的准则,不过我们预先并不知道函数值,所以无法做一个准确的判断,一般的想法是把计算区域
[
a
,
b
]
[a,b]
[a,b]选取尽量大一些。这里我们介绍两种处理方法:
U
j
m
+
1
=
U
j
+
1
m
+
U
j
−
1
m
2
−
τ
m
2
h
(
F
(
U
j
+
1
m
)
−
F
(
u
j
−
1
m
)
)
U_{j}^{m+1}=\frac{U_{j+1}^{m}+U_{j-1}^{m}}{2}-\frac{\tau_m}{2h}\left(F\left(U_{j+1}^{m}\right)-F\left(u_{j-1}^{m}\right)\right)
Ujm+1=2Uj+1m+Uj−1m−2hτm(F(Uj+1m)−F(uj−1m))
第一种处理方法是,针对第一个问题,我们知道函数初值是分片常数,在无穷远处导数为0,因此我们可以把计算区域
[
a
,
b
]
[a,b]
[a,b]选大以后,令
U
x
(
a
,
t
)
=
U
x
(
b
,
t
)
=
0
U_x(a,t) = U_x(b,t) = 0
Ux(a,t)=Ux(b,t)=0,这样的话在利用LF格式的时候就可以直接根据偏导数在边界上为0写出
U
−
1
m
=
U
1
m
,
U
N
m
=
U
N
−
2
m
U_{-1}^m = U_{1}^m,U_{N}^m = U_{N - 2}^m
U−1m=U1m,UNm=UN−2m,问题2正好有偏导数在边界上为0,因此这种处理方法可以直接应用到问题2中。\par
第二种处理方法是,针对第一个问题而言,我们算出第0层,也就是初始条件函数值
U
0
0
,
…
,
U
N
−
1
0
U_{0}^{0},\ldots,U_{N - 1}^{0}
U00,…,UN−10,在计算第一层的时候,我们只计算
U
1
1
,
…
,
U
N
−
2
1
U_{1}^{1},\ldots,U_{N - 2}^{1}
U11,…,UN−21,以此类推,计算第
m
m
m层的时候,我们只计算
U
m
m
,
…
,
U
N
−
m
−
1
m
U_{m}^{m},\ldots,U_{N - m - 1}^{m}
Umm,…,UN−m−1m,这样一直算到
t
m
a
x
t_{max}
tmax,这种做法有几个问题,首先我们知道我们最终需要输出的是
t
m
a
x
=
0.2
t_{max} = 0.2
tmax=0.2时函数在
[
0
,
1
]
[0,1]
[0,1]中的取值,而我们这个做法每做一层运算,区域都在向中间压缩减小,我们需要把计算区域
[
a
,
b
]
[a,b]
[a,b]选取的够大,使得最终算到最后一层的时候那个区域
[
x
m
,
x
N
−
m
−
1
]
⊇
[
0
,
1
]
[x_m,x_{N - m - 1}] \supseteq [0,1]
[xm,xN−m−1]⊇[0,1],这个其实缺少理论保证。
MacCormack格式
MacCormack格式参考
u
(
x
j
,
t
n
+
1
)
=
u
(
x
j
,
t
n
)
+
τ
(
u
t
)
j
n
+
1
2
τ
2
(
u
t
t
)
j
n
+
O
(
τ
3
)
=
1
2
u
(
x
j
,
t
n
)
+
1
2
(
(
u
ˉ
+
τ
u
ˉ
t
)
t
)
j
n
+
O
(
τ
3
)
,
u
ˉ
:
=
u
+
τ
u
t
.
\begin{aligned} u\left(x_{j}, t_{n+1}\right) &=u\left(x_{j}, t_{n}\right)+\tau\left(u_{t}\right)_{j}^{n}+\frac{1}{2} \tau^{2}\left(u_{t t}\right)_{j}^{n}+\mathcal{O}\left(\tau^{3}\right) \\ &=\frac{1}{2} u\left(x_{j}, t_{n}\right)+\frac{1}{2}(\left(\bar{u}+\tau \bar{u}_{t}\right)_{t})_{j}^{n}+\mathcal{O}\left(\tau^{3}\right), \quad \bar{u}:=u+\tau u_{t} . \end{aligned}
u(xj,tn+1)=u(xj,tn)+τ(ut)jn+21τ2(utt)jn+O(τ3)=21u(xj,tn)+21((uˉ+τuˉt)t)jn+O(τ3),uˉ:=u+τut.
对于向量形式,我们就可以对应做修改得到
{
W
j
=
U
j
n
−
τ
h
(
F
(
U
j
+
1
n
)
−
F
(
U
j
n
)
)
,
U
j
n
+
1
=
1
2
(
U
j
n
+
W
j
)
−
τ
2
h
(
F
(
W
j
)
−
F
(
W
j
−
1
)
)
.
\left\{\begin{array}{l} W_{j}=U_{j}^{n}-\frac{\tau}{h}\left(F\left(U_{j+1}^{n}\right)-F\left(U_{j}^{n}\right)\right), \\ U_{j}^{n+1}=\frac{1}{2}\left(U_{j}^{n}+W_{j} \right)-\frac{\tau}{2 h}\left(F\left(W_{j}\right)-F\left(W_{j-1}\right)\right). \end{array}\right.
{Wj=Ujn−hτ(F(Uj+1n)−F(Ujn)),Ujn+1=21(Ujn+Wj)−2hτ(F(Wj)−F(Wj−1)).
同样需要对边界索引做特殊处理,处理的方法和上面的LF格式如出一辙,这里就不赘诉了。
Roe格式
在介绍Roe格式之前,我们先回顾以下拟线性方程差分格式的基本定义,对于拟线性方程组:
U
t
+
F
(
U
)
x
=
0.
U_t + F(U)_x = 0.
Ut+F(U)x=0.
我们基本的差分格式定义如下:
U
j
n
+
1
=
U
j
n
−
τ
h
(
F
^
j
+
1
2
n
−
F
^
j
−
1
2
n
)
.
U_{j}^{n+1} = U_j^n - \frac{\tau}{h}(\hat{F}_{j + \frac{1}{2}}^n - \hat{F}_{j - \frac{1}{2}}^n).
Ujn+1=Ujn−hτ(F^j+21n−F^j−21n).
在Roe格式中,我们定义:
F
^
j
+
1
2
n
=
F
^
(
U
j
n
,
U
j
+
1
n
)
=
F
(
U
j
n
)
+
F
(
U
j
+
1
n
)
2
−
1
2
∣
A
^
j
n
∣
(
U
j
+
1
n
−
U
j
n
)
.
\hat{F}_{j + \frac{1}{2}}^n = \hat{F}(U_j^n,U_{j + 1}^n) = \frac{F(U_j^n) + F(U_{j+1}^n)}{2} - \frac{1}{2}|\hat{A}_{j}^n|(U_{j+1}^n - U_{j}^n).
F^j+21n=F^(Ujn,Uj+1n)=2F(Ujn)+F(Uj+1n)−21∣A^jn∣(Uj+1n−Ujn).
其中
∣
A
^
∣
=
R
∣
D
∣
R
−
1
,
∣
D
∣
=
d
i
a
g
{
∣
λ
1
∣
,
⋅
,
∣
λ
m
∣
}
|\hat{A}| = R|D|R^{-1},|D| = diag\{|\lambda_1|,\cdot,|\lambda_m| \}
∣A^∣=R∣D∣R−1,∣D∣=diag{∣λ1∣,⋅,∣λm∣},其中
R
D
R
−
1
=
A
RDR^{-1}=A
RDR−1=A,整理得到最终的Roe格式为:
U
j
n
+
1
=
U
j
n
−
τ
h
(
F
(
U
j
+
1
n
)
−
F
(
U
j
−
1
n
)
2
+
1
2
(
∣
A
^
j
n
∣
(
U
j
+
1
n
−
U
j
n
)
−
∣
A
^
j
−
1
n
∣
(
U
j
n
−
U
j
−
1
n
)
)
.
U_{j}^{n+1} = U_j^n - \frac{\tau}{h}(\frac{F(U_{j + 1}^n) - F(U_{j-1}^n)}{2} + \frac{1}{2}(|\hat{A}_{j}^n|(U_{j+1}^n - U_{j}^n) - |\hat{A}_{j-1}^n|(U_{j}^n - U_{j-1}^n)).
Ujn+1=Ujn−hτ(2F(Uj+1n)−F(Uj−1n)+21(∣A^jn∣(Uj+1n−Ujn)−∣A^j−1n∣(Ujn−Uj−1n)).
数值实验
下面我们来展示两个问题的数值图像。
关于问题1,我们需要求出时间
t
=
0.2
t = 0.2
t=0.2时,函数在
[
0
,
1
]
[0,1]
[0,1]的取值,这里本人定义的计算区域是
[
−
5
,
5
]
,
h
x
=
0.02
[-5,5],hx = 0.02
[−5,5],hx=0.02。
关于问题2,我们需要求出时间
t
=
1.8
t = 1.8
t=1.8时,函数在
[
−
5
,
5
]
[-5,5]
[−5,5]的取值,这里是固定的边界
[
−
5
,
5
]
[-5,5]
[−5,5],hx = 0.2,这里特别说明,在求解问题2的时候,不知道什么原因,导致Roe格式的CFL数必须要取很小,本人取得是0.3,即使这么处理以后,仍然无法算到1.8秒,格式在计算到中间差不多0.6秒得时候出现了数值爆炸,准确来说是计算声速a的时候出现了对负数开根号的计算,这一点本人实在找不出问题所在,因此本人使用Roe格式的时候只算到了0.5秒。
代码分析
其中CFD-1.py,CFD-2.py分别表示问题1和问题2的数值代码,里面的函数我都写了注释,值得说明的是,本人这两个代码都写了一个matrix函数,这个函数负责把当前的迭代点的取值输入,然后把
∣
A
^
∣
=
R
∣
D
∣
R
−
1
,
∣
D
∣
=
d
i
a
g
{
∣
λ
1
∣
,
⋅
,
∣
λ
m
∣
}
|\hat{A}| = R|D|R^{-1},|D| = diag\{|\lambda_1|,\cdot,|\lambda_m| \}
∣A^∣=R∣D∣R−1,∣D∣=diag{∣λ1∣,⋅,∣λm∣}的
∣
A
^
∣
|\hat{A}|
∣A^∣输出来,根据的原理是:
A
(
U
)
=
(
0
1
0
γ
−
3
2
u
2
(
3
−
γ
)
u
γ
−
1
γ
−
2
2
u
3
−
a
2
u
γ
−
1
3
−
2
γ
2
u
2
+
a
2
γ
−
1
γ
u
)
A(U)=\left(\begin{array}{ccc} 0 & 1 & 0 \\ \frac{\gamma-3}{2} u^{2} & (3-\gamma) u & \gamma-1 \\ \frac{\gamma-2}{2} u^{3}-\frac{a^{2} u}{\gamma-1} & \frac{3-2 \gamma}{2} u^{2}+\frac{a^{2}}{\gamma-1} & \gamma u \end{array}\right)
A(U)=
02γ−3u22γ−2u3−γ−1a2u1(3−γ)u23−2γu2+γ−1a20γ−1γu
λ
1
=
u
−
a
,
λ
2
=
u
,
λ
3
=
u
+
a
,
R
=
(
R
1
,
R
(
2
)
,
R
(
3
)
)
=
(
1
1
1
u
−
a
u
u
+
a
H
−
u
a
1
2
u
2
H
+
u
a
,
)
\begin{array}{c} \lambda_{1}=u-a, \quad \lambda_{2}=u, \quad \lambda_{3}=u+a ,\\ R=\left(R^{1}, R^{(2)}, R^{(3)}\right)=\left(\begin{array}{ccc} 1 & 1 & 1 \\ u-a & u & u+a \\ H-u a & \frac{1}{2} u^{2} & H+u a,\\ \end{array}\right) \end{array}
λ1=u−a,λ2=u,λ3=u+a,R=(R1,R(2),R(3))=
1u−aH−ua1u21u21u+aH+ua,
a
3
,
1
=
u
(
u
a
+
u
2
−
2
H
)
/
(
4
a
H
−
2
a
u
2
)
,
a
3
,
2
=
(
2
H
−
u
2
−
2
a
u
)
/
(
4
a
H
−
2
a
u
2
)
,
R
−
1
=
(
u
/
a
+
a
3
,
1
a
3
,
2
−
1
/
a
1
/
(
2
H
−
u
2
)
1
−
u
/
a
−
2
a
3
,
1
1
/
a
−
2
u
3
,
2
1
/
(
u
2
−
2
H
)
a
3
,
1
a
3
,
2
1
/
(
2
H
−
u
2
)
)
\begin{array}{c} a_{3,1} = u(ua + u^2 - 2H)/(4aH - 2au^2),a_{3,2} = (2H - u^2 - 2au)/(4aH - 2au^2),\\ R^{-1}=\left(\begin{array}{ccc} u/a + a_{3,1} & a_{3,2} - 1/a & 1/(2H - u^2) \\ 1 - u/a - 2a_{3,1} & 1/a - 2u_{3,2} & 1/(u^2 - 2H) \\ a_{3,1} & a_{3,2} & 1/(2H - u^2) \end{array}\right) \end{array}
a3,1=u(ua+u2−2H)/(4aH−2au2),a3,2=(2H−u2−2au)/(4aH−2au2),R−1=
u/a+a3,11−u/a−2a3,1a3,1a3,2−1/a1/a−2u3,2a3,21/(2H−u2)1/(u2−2H)1/(2H−u2)
上面的3个格式,其中LF和Mac格式的实现都相对简单,本人一开始固定时间步长来迭代,结果函数震荡的十分厉害,后来在同学的建议下改成了依据CFL条件每一步修改时间步长的做法,这样才让函数输出稳定了下来,对于LF和Mac格式,本人一开始打算每做一次时间层的迭代,就把空间向中间压缩两个节点的做法,但是不知道为什么这样做得到的解是错的,甚至于是不收敛的。这些代码中,最花时间的是Roe格式,这个格式涉及到矩阵的求逆运算,本来一开始打算定义矩阵A,然后调用Python中自带的矩阵分解函数来求
∣
A
∣
|A|
∣A∣,但是在迭代的过程很难保证这种数值格式产生的
∣
A
∣
|A|
∣A∣是良态的,导致计算特别慢,而且极容易出错,所以本人选择自己手动计算逆矩阵。
CFD-1
import numpy as np
import time
import matplotlib.pyplot as plt
def UB(x):#函数初值
n = x.shape[0]
tmp = np.zeros([n,3])
leq = (x < 0.3).astype('float32')
geq = (x > 0.3).astype('float32')
tmp[:,0] = leq*1 + geq*0.125
tmp[:,1] = leq*0 + geq*0
tmp[:,2] = leq*2.5 + geq*0.25
return tmp
def matrix(U):#计算Roe格式的时候需要求矩阵的分解和绝对值组合
gama = 1.4
eps = 1e-7
n = U.shape[1]
R = np.ones([n,n]);L = np.zeros([n,n])
a = np.sqrt(gama*(gama - 1)*(U[0,2]/(U[0,0] + eps) - 0.5*(U[0,1]**2/(U[0,0]**2 + eps))))
u = U[0,1]/(U[0,0] + eps)
H = (U[0,2] + (gama - 1)*(U[0,2] - 0.5*U[0,1]**2/(U[0,0] + eps)))/(U[0,0] + eps)
R[1,0] = u - a;R[1,1] = u;R[1,2] = u + a
R[2,0] = H - u*a;R[2,1] = 0.5*u**2;R[2,2] = H + u*a
a31 = u*(u*a + u**2 - 2*H)/(4*a*H - 2*a*u**2)
a32 = (2*H - u**2 - 2*a*u)/(4*a*H - 2*a*u**2)
L[0,0] = u/a + a31;L[0,1] = a32 - 1/a;L[0,2] = 1/(2*H - u**2)
L[1,0] = 1 - u/a - 2*a31;L[1,1] = 1/a - 2*a32;L[1,2] = 2/(u**2 - 2*H)
L[2,0] = a31;L[2,1] = a32;L[2,2] = 1/(2*H - u**2)
lam = np.array([u + a,u,u - a])
return R@np.diag(abs(lam))@L
def FF(U):#计算源函数
gama = 1.4
eps = 1e-7
n = U.shape[0]
tmp = np.zeros([n,3])
tmp[:,0] = U[:,1]
tmp[:,1] = (gama - 1)*U[:,2] + 0.5*(3 - gama)*(U[:,1]**2)/(U[:,0] + eps)
tmp[:,2] = gama*U[:,2]*U[:,1]/(U[:,0] + eps) + 0.5*(1 - gama)*(U[:,1]**3)/(U[:,0]**2 + eps)
return tmp
def LAMmax(U):#计算最大特征值
gama = 1.4
eps = 1e-7
a = np.sqrt(gama*(gama - 1)*(U[:,2]/(U[:,0] + eps) - 0.5*(U[:,1]**2/(U[:,0]**2 + eps))))
tmp = abs(U[:,1]/(U[:,0] + eps)) + a
return max(tmp)
def UU(tmp):#把守恒变量替换成原始变量
n = tmp.shape[0]
u = np.zeros([n,3])
gama = 1.4
eps = 1e-7
u[:,0] = tmp[:,0]#rho
u[:,1] = tmp[:,1]/(tmp[:,0] + eps)
u[:,2] = (gama - 1)*(tmp[:,2] - 0.5*tmp[:,1]**2/(tmp[:,0] + eps))
return u
class FD():
def __init__(self,bound,nx):
self.bound = bound#定义空间和时间的边界,比计算区域大
self.hx = (bound[0,1] - bound[0,0])/(nx - 1)#定义步长
self.X = np.linspace(bound[0,0],bound[0,1],nx)#包含计算区域
self.nx = nx
self.gama = 1.4
self.cpu = [0,1]#空间计算区域
self.x = np.arange(self.cpu[0],self.cpu[1] + self.hx,self.hx)#离散空间计算区域
self.xL = np.arange(self.bound[0,0],self.cpu[0] + self.hx,self.hx)#把区域分成3份,这是左
self.xR = np.arange(self.cpu[1],self.bound[0,1] + self.hx,self.hx)#右
def solve(self,mode):
U_old = UB(self.X)
#print(U_old.shape)
U_new = np.zeros_like(U_old)
cfl = 0.9
#---------------------------------------------------
if mode == 'LF':
m = 0
t = self.bound[1,0]
while m < 100:
tau = cfl*self.hx/LAMmax(U_old)
U_new[0,:] = U_old[1,:]
for i in range(1,U_old.shape[0] - 1):
U_new[i:i + 1,:] = 0.5*(U_old[i - 1:i,:] + U_old[i + 1:i + 2,:]) - \
0.5*tau*(FF(U_old[i + 1:i + 2,:]) - FF(U_old[i - 1:i,:]))/self.hx
U_new[-1,:] = U_old[-2,:]
if t <= self.bound[1,1] and t + tau >= self.bound[1,1]:
tau = self.bound[1,1] - t
U_new[0,:] = U_old[1,:]
for i in range(1,U_old.shape[0] - 1):
U_new[i:i + 1,:] = 0.5*(U_old[i - 1:i,:] + U_old[i + 1:i + 2,:]) - \
0.5*tau*(FF(U_old[i + 1:i + 2,:]) - FF(U_old[i - 1:i,:]))/self.hx
U_new[-1,:] = U_old[-2,:]
t += tau
break
else:
t += tau
U_old = U_new
m = m + 1
print('the mode:%s,the epoch:%d,the problem:%d,the time:%.2f'%('LF',m,2,t))
return U_new[self.xL.size - 1:self.xL.size + self.x.size - 1,:]
if mode == 'Mac':
cfl = 0.9
m = 0
t = self.bound[1,0]
while m < 100:
eig = LAMmax(U_old)
tau = self.hx*cfl/eig
W0 = U_old[1:2,:] - tau*(FF(U_old[0:1,:]) - FF(U_old[1:2,:]))/self.hx
W1 = U_old[0:1,:] - tau*(FF(U_old[1:2,:]) - FF(U_old[0:1,:]))/self.hx
U_new[0:1] = 0.5*(U_old[0:1] + W1) - 0.5*tau*(FF(W1) - FF(W0))/self.hx
for i in range(1,U_old.shape[0] - 1):
W0 = U_old[i - 1:i,:] - tau*(FF(U_old[i:i + 1,:]) - FF(U_old[i - 1:i,:]))/self.hx
W1 = U_old[i:i + 1,:] - tau*(FF(U_old[i + 1:i + 2,:]) - FF(U_old[i:i + 1,:]))/self.hx
U_new[i:i + 1,:] = 0.5*(U_old[i:i + 1,:] + W1) - \
0.5*tau*(FF(W1) - FF(W0))/self.hx
W0 = U_old[-2:-1,:] - tau*(FF(U_old[-1:,:]) - FF(U_old[-2:-1,:]))/self.hx
W1 = U_old[-1:,:] - tau*(FF(U_old[-2:-1,:]) - FF(U_old[-1:,:]))/self.hx
U_new[-1:,:] = 0.5*(U_old[-1:,:] + W1) - 0.5*tau*(FF(W1) - FF(W0))/self.hx
if t <= self.bound[1,1] and t + tau >= self.bound[1,1]:
tau = self.bound[1,1] - t
W0 = U_old[1:2,:] - tau*(FF(U_old[0:1,:]) - FF(U_old[1:2,:]))/self.hx
W1 = U_old[0:1,:] - tau*(FF(U_old[1:2,:]) - FF(U_old[0:1,:]))/self.hx
U_new[0:1] = 0.5*(U_old[0:1] + W1) - 0.5*tau*(FF(W1) - FF(W0))/self.hx
for i in range(1,U_old.shape[0] - 1):
W0 = U_old[i - 1:i,:] - tau*(FF(U_old[i:i + 1,:]) - FF(U_old[i - 1:i,:]))/self.hx
W1 = U_old[i:i + 1,:] - tau*(FF(U_old[i + 1:i + 2,:]) - FF(U_old[i:i + 1,:]))/self.hx
U_new[i:i + 1,:] = 0.5*(U_old[i:i + 1,:] + W1) - \
0.5*tau*(FF(W1) - FF(W0))/self.hx
W0 = U_old[-2:-1,:] - tau*(FF(U_old[-1:,:]) - FF(U_old[-2:-1,:]))/self.hx
W1 = U_old[-1:,:] - tau*(FF(U_old[-2:-1,:]) - FF(U_old[-1:,:]))/self.hx
U_new[-1:,:] = 0.5*(U_old[-1:,:] + W1) - 0.5*tau*(FF(W1) - FF(W0))/self.hx
t += tau
break
else:
t += tau
U_old = U_new
m = m + 1
print('the mode:%s,the epoch:%d,the problem:%d,the time:%.2f'%('Mac',m,2,t))
return U_new[self.xL.size - 1:self.xL.size + self.x.size - 1,:]
if mode == 'Roe':
cfl = 0.9
m = 0
t = self.bound[1,0]
U_new = U_old.copy()
while m < 100:
eig = LAMmax(U_old)
tau = self.hx*cfl/eig
mat = matrix(U_old[0:1,:]) + matrix(U_old[1:2,:])
U_new[0:1,:] = U_old[0:1,:] - 0.5*tau*(U_old[0:1,:] - U_old[1:2,:])@mat.T/self.hx
for i in range(1,U_old.shape[0] - 1):
W0 = 0.5*(FF(U_old[i - 1:i,:]) + FF(U_old[i:i + 1,:])) - \
0.5*(U_old[i:i + 1,:] - U_old[i - 1:i,:])@matrix(U_old[i - 1:i,:]).T
W1 = 0.5*(FF(U_old[i:i + 1,:]) + FF(U_old[i + 1:i + 2,:])) + \
0.5*(U_old[i + 1:i + 2,:] - U_old[i:i + 1,:])@matrix(U_old[i:i + 1,:]).T
U_new[i:i + 1,:] = U_old[i:i + 1,:] - tau*(W1 - W0)/self.hx
mat = matrix(U_old[-1:,:]) + matrix(U_old[-2:-1,:])
U_new[-1:,:] = U_old[-1:,:] - 0.5*tau*(U_old[-1:,:] - U_old[-2:-1,:])@mat.T/self.hx
if t <= self.bound[1,1] and t + tau >= self.bound[1,1]:
tau = self.bound[1,1] - t
mat = matrix(U_old[0:1,:]) + matrix(U_old[1:2,:])
U_new[0:1,:] = U_old[0:1,:] - 0.5*tau*(U_old[0:1,:] - U_old[1:2,:])@mat.T/self.hx
for i in range(1,U_old.shape[0] - 1):
W0 = 0.5*(FF(U_old[i - 1:i,:]) + FF(U_old[i:i + 1,:])) - \
0.5*(U_old[i:i + 1,:] - U_old[i - 1:i,:])@matrix(U_old[i - 1:i,:]).T
W1 = 0.5*(FF(U_old[i:i + 1,:]) + FF(U_old[i + 1:i + 2,:])) + \
0.5*(U_old[i + 1:i + 2,:] - U_old[i:i + 1,:])@matrix(U_old[i:i + 1,:]).T
U_new[i:i + 1,:] = U_old[i:i + 1,:] - tau*(W1 - W0)/self.hx
mat = matrix(U_old[-1:,:]) + matrix(U_old[-2:-1,:])
U_new[-1:,:] = U_old[-1:,:] - 0.5*tau*(U_old[-1:,:] - U_old[-2:-1,:])@mat.T/self.hx
t += tau
break
else:
t += tau
m = m + 1
U_new = U_old
print('the mode:%s,the epoch:%d,the problem:%d,the time:%.2f'%('Roe',m,2,t))
return U_new[self.xL.size - 1:self.xL.size + self.x.size - 1,:]
bound = np.array([[-5,5],[0,0.2]])
nx = 501
mode = 'Roe'
fd = FD(bound,nx)
u_pred = UU(fd.solve(mode))
#plt.plot(fd.x,u_pred[:,0],'r*')
plt.plot(fd.x,u_pred[:,0],'r',label = 'rho')
#plt.plot(fd.x,u_pred[:,1],'bo')
plt.plot(fd.x,u_pred[:,1],'b',label = 'u')
#plt.plot(fd.x,u_pred[:,2],'k.')
plt.plot(fd.x,u_pred[:,2],'k',label = 'P')
plt.xlabel('the space:[%.2f,%.2f]'%(fd.cpu[0],fd.cpu[1]))
plt.ylabel('the numerical solutions')
plt.title('the method :%s, at time:%.2f'%(mode,bound[1,1]))
plt.legend(loc = 'upper right')
plt.savefig('%s1.jpg'%(mode))
CFD-2
import numpy as np
import time
import matplotlib.pyplot as plt
def UB(x):#初始条件
n = x.shape[0]
gama = 1.4
tmp = np.zeros([n,3])
leq = (x < -4).astype('float32')
geq = (x >= -4).astype('float32')
tmp[:,0] = leq*3.857143 + geq*(1 + 0.2*np.sin(5*x))
tmp[:,1] = leq*2.629329 + geq*0
tmp[:,2] = leq*10.3333 + geq*1
s = tmp.copy()
s[:,0] = tmp[:,0]
s[:,1] = tmp[:,0]*tmp[:,1]
s[:,2] = tmp[:,2]/(gama - 1.0) + 0.5*tmp[:,0]*tmp[:,1]**2
return s
def matrix(U):#计算Roe格式的时候需要求矩阵的分解和绝对值组合
gama = 1.4
eps = 1e-7
n = U.shape[1]
R = np.ones([n,n]);L = np.zeros([n,n])
a = np.sqrt(gama*(gama - 1)*(U[0,2]/(U[0,0] + eps) - 0.5*(U[0,1]**2/(U[0,0]**2 + eps))))
u = U[0,1]/(U[0,0] + eps)
H = (U[0,2] + (gama - 1)*(U[0,2] - 0.5*U[0,1]**2/(U[0,0] + eps)))/(U[0,0] + eps)
R[1,0] = u - a;R[1,1] = u;R[1,2] = u + a
R[2,0] = H - u*a;R[2,1] = 0.5*u**2;R[2,2] = H + u*a
a31 = u*(u*a + u**2 - 2*H)/(4*a*H - 2*a*u**2)
a32 = (2*H - u**2 - 2*a*u)/(4*a*H - 2*a*u**2)
L[0,0] = u/a + a31;L[0,1] = a32 - 1/a;L[0,2] = 1/(2*H - u**2)
L[1,0] = 1 - u/a - 2*a31;L[1,1] = 1/a - 2*a32;L[1,2] = 2/(u**2 - 2*H)
L[2,0] = a31;L[2,1] = a32;L[2,2] = 1/(2*H - u**2)
lam = np.array([u + a,u,u - a])
return R@np.diag(abs(lam))@L
def FF(U):#计算源函数
gama = 1.4
eps = 1e-7
n = U.shape[0]
tmp = np.zeros([n,3])
tmp[:,0] = U[:,1]
tmp[:,1] = (gama - 1)*U[:,2] + 0.5*(3 - gama)*(U[:,1]**2)/(U[:,0] + eps)
tmp[:,2] = gama*U[:,2]*U[:,1]/(U[:,0] + eps) + 0.5*(1 - gama)*(U[:,1]**3)/(U[:,0]**2 + eps)
return tmp
def LAMmax(U):#计算最大特征值
gama = 1.4
eps = 1e-7
a = np.sqrt(gama*(gama - 1)*(U[:,2]/(U[:,0] + eps) - 0.5*(U[:,1]**2/(U[:,0]**2 + eps))))
tmp = abs(U[:,1]/(U[:,0] + eps)) + a
return max(tmp)
def UU(tmp):#把守恒变量替换成原始变量
n = tmp.shape[0]
u = np.zeros([n,3])
gama = 1.4
eps = 1e-7
u[:,0] = tmp[:,0]#rho
u[:,1] = tmp[:,1]/(tmp[:,0] + eps)
u[:,2] = (gama - 1)*(tmp[:,2] - 0.5*tmp[:,1]**2/(tmp[:,0] + eps))
return u
class FD():
def __init__(self,bound,nx):
self.bound = bound#定义空间和时间的边界,比计算区域大
self.hx = (bound[0,1] - bound[0,0])/(nx - 1)#定义步长
self.X = np.linspace(bound[0,0],bound[0,1],nx)#包含计算区域
self.nx = nx
self.gama = 1.4
def solve(self,mode):
U_old = UB(self.X)
#print(U_old.shape)
U_new = np.zeros_like(U_old)
cfl = 1.0
#---------------------------------------------------
if mode == 'LF':
m = 0
t = self.bound[1,0]
U_new = U_old.copy()
while m < 100:
#print(UU(U_new)[:,0])
tau = cfl*self.hx/LAMmax(U_old)
U_new[0,:] = U_old[1,:]
for i in range(1,U_old.shape[0] - 1):
U_new[i:i + 1,:] = 0.5*(U_old[i - 1:i,:] + U_old[i + 1:i + 2,:]) - \
0.5*tau*(FF(U_old[i + 1:i + 2,:]) - FF(U_old[i - 1:i,:]))/self.hx
U_new[-1,:] = U_old[-2,:]
if t <= self.bound[1,1] and t + tau >= self.bound[1,1]:
tau = self.bound[1,1] - t
U_new[0,:] = U_old[1,:]
for i in range(1,U_old.shape[0] - 1):
U_new[i:i + 1,:] = 0.5*(U_old[i - 1:i,:] + U_old[i + 1:i + 2,:]) - \
0.5*tau*(FF(U_old[i + 1:i + 2,:]) - FF(U_old[i - 1:i,:]))/self.hx
U_new[-1,:] = U_old[-2,:]
t += tau
break
else:
t += tau
U_old = U_new
m = m + 1
print('the mode:%s,the epoch:%d,the problem:%d,the time:%.2f'%('LF',m,2,t))
return U_new
if mode == 'Mac':
cfl = 0.9
m = 0
t = self.bound[1,0]
U_new = U_old.copy()
while m < 100:
eig = LAMmax(U_old)
tau = self.hx*cfl/eig
W0 = U_old[1:2,:] - tau*(FF(U_old[0:1,:]) - FF(U_old[1:2,:]))/self.hx
W1 = U_old[0:1,:] - tau*(FF(U_old[1:2,:]) - FF(U_old[0:1,:]))/self.hx
U_new[0:1] = 0.5*(U_old[0:1] + W1) - 0.5*tau*(FF(W1) - FF(W0))/self.hx
for i in range(1,U_old.shape[0] - 1):
W0 = U_old[i - 1:i,:] - tau*(FF(U_old[i:i + 1,:]) - FF(U_old[i - 1:i,:]))/self.hx
W1 = U_old[i:i + 1,:] - tau*(FF(U_old[i + 1:i + 2,:]) - FF(U_old[i:i + 1,:]))/self.hx
U_new[i:i + 1,:] = 0.5*(U_old[i:i + 1,:] + W1) - \
0.5*tau*(FF(W1) - FF(W0))/self.hx
W0 = U_old[-2:-1,:] - tau*(FF(U_old[-1:0,:]) - FF(U_old[-2:-1,:]))/self.hx
W1 = U_old[-1:0,:] - tau*(FF(U_old[-2:-1,:]) - FF(U_old[-1:0,:]))/self.hx
U_new[-1:0] = 0.5*(U_old[-1:0] + W1) - 0.5*tau*(FF(W1) - FF(W0))/self.hx
if t <= self.bound[1,1] and t + tau >= self.bound[1,1]:
tau = self.bound[1,1] - t
W0 = U_old[1:2,:] - tau*(FF(U_old[0:1,:]) - FF(U_old[1:2,:]))/self.hx
W1 = U_old[0:1,:] - tau*(FF(U_old[1:2,:]) - FF(U_old[0:1,:]))/self.hx
U_new[0:1] = 0.5*(U_old[0:1] + W1) - 0.5*tau*(FF(W1) - FF(W0))/self.hx
for i in range(1,U_old.shape[0] - 1):
W0 = U_old[i - 1:i,:] - tau*(FF(U_old[i:i + 1,:]) - FF(U_old[i - 1:i,:]))/self.hx
W1 = U_old[i:i + 1,:] - tau*(FF(U_old[i + 1:i + 2,:]) - FF(U_old[i:i + 1,:]))/self.hx
U_new[i:i + 1,:] = 0.5*(U_old[i:i + 1,:] + W1) - \
0.5*tau*(FF(W1) - FF(W0))/self.hx
W0 = U_old[-2:-1,:] - tau*(FF(U_old[-1:0,:]) - FF(U_old[-2:-1,:]))/self.hx
W1 = U_old[-1:0,:] - tau*(FF(U_old[-2:-1,:]) - FF(U_old[-1:0,:]))/self.hx
U_new[-1:0] = 0.5*(U_old[-1:0] + W1) - 0.5*tau*(FF(W1) - FF(W0))/self.hx
t += tau
break
else:
t += tau
U_old = U_new
m = m + 1
print('the mode:%s,the epoch:%d,the problem:%d,the time:%.2f'%('Mac',m,2,t))
return U_new
if mode == 'Roe':
cfl = 0.3
m = 0
t = self.bound[1,0]
U_new = U_old.copy()
while m < 150:
eig = LAMmax(U_old)
tau = self.hx*cfl/eig
mat = matrix(U_old[0:1,:]) + matrix(U_old[1:2,:])
U_new[0:1,:] = U_old[0:1,:] - 0.5*tau*(U_old[0:1,:] - U_old[1:2,:])@mat.T/self.hx
for i in range(1,U_old.shape[0] - 1):
W0 = 0.5*(FF(U_old[i - 1:i,:]) + FF(U_old[i:i + 1,:])) - \
0.5*(U_old[i:i + 1,:] - U_old[i - 1:i,:])@matrix(U_old[i - 1:i,:]).T
W1 = 0.5*(FF(U_old[i:i + 1,:]) + FF(U_old[i + 1:i + 2,:])) + \
0.5*(U_old[i + 1:i + 2,:] - U_old[i:i + 1,:])@matrix(U_old[i:i + 1,:]).T
U_new[i:i + 1,:] = U_old[i:i + 1,:] - tau*(W1 - W0)/self.hx
mat = matrix(U_old[-1:,:]) + matrix(U_old[-2:-1,:])
U_new[-1:,:] = U_old[-1:,:] - 0.5*tau*(U_old[-1:,:] - U_old[-2:-1,:])@mat.T/self.hx
if t <= self.bound[1,1] and t + tau >= self.bound[1,1]:
tau = self.bound[1,1] - t
mat = matrix(U_old[0:1,:]) + matrix(U_old[1:2,:])
U_new[0:1,:] = U_old[0:1,:] - 0.5*tau*(U_old[0:1,:] - U_old[1:2,:])@mat.T/self.hx
for i in range(1,U_old.shape[0] - 1):
W0 = 0.5*(FF(U_old[i - 1:i,:]) + FF(U_old[i:i + 1,:])) - \
0.5*(U_old[i:i + 1,:] - U_old[i - 1:i,:])@matrix(U_old[i - 1:i,:]).T
W1 = 0.5*(FF(U_old[i:i + 1,:]) + FF(U_old[i + 1:i + 2,:])) + \
0.5*(U_old[i + 1:i + 2,:] - U_old[i:i + 1,:])@matrix(U_old[i:i + 1,:]).T
U_new[i:i + 1,:] = U_old[i:i + 1,:] - tau*(W1 - W0)/self.hx
mat = matrix(U_old[-1:,:]) + matrix(U_old[-2:-1,:])
U_new[-1:,:] = U_old[-1:,:] - 0.5*tau*(U_old[-1:,:] - U_old[-2:-1,:])@mat.T/self.hx
t += tau
break
else:
t += tau
m = m + 1
U_new = U_old
print('the mode:%s,the epoch:%d,the problem:%d,the time:%.2f'%('Roe',m,2,t))
return U_new
#bound = np.array([[-5,5],[0,0.5]])
bound = np.array([[-5,5],[0,1.8]])
nx = 101
mode = 'Mac'
fd = FD(bound,nx)
u_pred = UU(fd.solve(mode))
#plt.plot(fd.x,u_pred[:,0],'r*')
plt.plot(fd.X,u_pred[:,0],'r',label = 'rho')
#plt.plot(fd.x,u_pred[:,1],'bo')
plt.plot(fd.X,u_pred[:,1],'b',label = 'u')
#plt.plot(fd.x,u_pred[:,2],'k.')
plt.plot(fd.X,u_pred[:,2],'k',label = 'P')
plt.xlabel('the space:[%.2f,%.2f]'%(fd.bound[0,0],fd.bound[0,1]))
plt.ylabel('the numerical solutions')
plt.title('the method :%s, at time:%.2f'%(mode,bound[1,1]))
plt.legend(loc = 'upper right')
plt.savefig('%s2.jpg'%(mode))