一维完全气体的Riemann问题求解

已于 2023-04-18 09:30:00 修改
阅读量2.2k 收藏 15
点赞数 3
分类专栏： 传统数值方法 文章标签： 算法
于 2021-05-04 16:13:22 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/forrestguang/article/details/116401291
版权
传统数值方法 专栏收录该内容
35 篇文章 18 订阅 ¥199.90 ¥299.90
订阅专栏
本文详细介绍了利用LF、MacCormack和Roe格式对一维完全气体欧拉方程进行数值求解的过程。通过解决两个具体的Riemann问题，展示了不同格式的计算步骤和数值实验结果。在LF和MacCormack格式中，通过CFL条件调整时间步长确保稳定性，而Roe格式涉及矩阵运算，可能导致数值爆炸。代码分析部分给出了问题1和问题2的Python实现。
摘要由CSDN通过智能技术生成

一维完全气体欧拉方程的引入

在这里我们将利用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​ρ=∂x​u=∂x​P=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γ−3​u22γ−2​u3−γ−1a2u​​1(3−γ)u23−2γ​u2+γ−1a2​​0γ−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,tm​ax]成立，其中 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)=21​u(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+21​n​−F^j−21​n​).
在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+21​n​=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γ−3​u22γ−2​u3−γ−1a2u​​1(3−γ)u23−2γ​u2+γ−1a2​​0γ−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−ua​1u21​u2​1u+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,1​1−u/a−2a3,1​a3,1​​a3,2​−1/a1/a−2u3,2​a3,2​​1/(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))
Galerkin码农选手
关注
  • 3
    点赞
  • 踩
  • 15
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 5
    评论
专栏目录
订阅专栏
常微分方程算法之编程示例一（欧拉法）
欢迎各位大佬光临“寒舍”！请多多指点，共同进步！
06-21 815
主要对本专栏介绍的常微分方程差分格式进行算例补充。
c++:一维黎曼问题多种格式求解
weixin_44575323的博客
03-13 4475
c++:多种格式求解一维黎曼问题）一、一维黎曼问题简介新的改变功能快捷键合理的创建标题，有助于目录的生成如何改变文本的样式插入链接与图片如何插入一段漂亮的代码片生成一个适合你的列表创建一个表格设定内容居中、居左、居右SmartyPants创建一个自定义列表如何创建一个注脚注释也是必不可少的KaTeX数学公式新的甘特图功能，丰富你的文章UML 图表FLowchart流程图导出与导入导出导入 一、一维黎曼问题简介 一无限长管道内部充满理想气体（无粘），中间有一不计厚度的隔膜，初始时刻 t_0，左侧气体的速度、密
一维黎曼问题
05-29
应用c 与fortran 语言解决一维黎曼问题
Ri_solve:黎曼问题的求解器-matlab开发
06-01
该函数求解关于理想多方气体的一维不可压缩欧拉方程的初始数据不连续性的黎曼问题。 计算也可以矢量化
MacCormack差分格式的全局误差分析
LSEC小陆的博客
04-30 2392
MacCormack差分格式的全局误差分析 问题 请分析MacCormack差分格式的全局误差： un+12=un+hun′u_{n+\frac{1}{2}} = u_n + h u&#x27;_nun+21​​=un​+hun′​ un+1=12(un+un+12+hun+12′)u_{n+1} =\frac{1}{2} \left( u_n + u_{n+\frac{1}{2}} + ...
一维理想气体Euler方程 有限体积求解 附python代码
Tankrun1997的博客
09-07 709
一维有限体积方法求解sod激波管问题 python 代码
1D Riemann problem
Leslie Lee's blog includes programming in math, physics, and engineering.
06-27 1333
1D Riemann problem 出自 Riemann Solvers and Numerical Methods for Fluid Dynamics, Pg 44, Example 2.5 空气动力学中的一维黎曼问题是一个特殊的初值问题，由欧拉方程与分段初值条件组成 初值条件，t=0t=0t=0时： x≤0时,ρ=ρL,u=uLx \leq 0时, \rho=\rho_L, u=u_Lx≤0时,ρ=ρL​,u=uL​ x>0时,ρ=ρR,u=uRx > 0时, \rho=\rho_R,
Riemann_accurate_激波管精确解_激波管_一维激波管问题精确解_
10-01
这个问题涉及到在一维管道中，当一个理想的气体经历突然的压力或速度变化时，如何通过数值方法准确地模拟出内部流动的复杂动态过程。这里的"Riemann_accurate"指的是使用黎曼问题的精确解来处理这种情况，这是一种...
一维激波管问题的数值求解代码.rar.rar
06-11
《一维激波管问题的数值求解代码》是一个针对流体力学中经典问题的编程实践，主要涉及一维流动中的激波理论及其数值模拟。一维激波管问题在工程和物理研究中具有重要地位，它可以帮助我们理解和预测气体在高压、高速...
附录A 一维黎曼问题(A).doc
08-10
一维 问题，即激波管问题，是一个典型的一维可压缩无黏气体动力学问题，并有 解析解。对它采用二阶精度 两步差分格式进行数值求解。同时，为了初学者入门和练习方便,这里给出了用 语言和 编写的计算一维 问题的计算...
sod_激波管精确解_激波管_精确解_一维python_一维激波管_
09-30
这个主题涉及到一维流动问题的数值求解，通常用于测试和验证数值方法的准确性和稳定性。在这个压缩包中，包含两个Python脚本——sod_roe.py和sod_exact.py，它们可能是用于实现激波管问题的精确解和Roe平均线性化的...
一维黎曼问题的精确解
03-21
主要采用fortran编译了一维黎曼问题的精确解，可用于计算空气动力学计算
Riemann_Solver(SWE)_Riemann求解器_浅水方程matlab_一维浅水波方程_
10-02
一维浅水波方程的Riemann求解器算法，新版matlab程序
可压缩准一维管道无粘流动-欧拉方程-MacCormack.rar_channel flow _supersonic_亚网格_可压
07-15
本题利用一维欧拉方程求解可压缩一维无粘管道流动，并针对出口不同的条件，出口亚音速和出口超音速两种不同条件，将流道进行网格划分，利用MacCormack进行差分求解，得到管道内的总压、马赫数、总焓、内能的分布，并给出计算过程中残差收敛的过程。
maccormack.zip_maccormack格式前台阶扰流
07-15
二维经典前台阶扰流的有黏二阶精度maccormack格式的代码
一维Sod激波管问题Riemann间断解算例及python代码
03-30
对一维Sod激波管问题的分析和数值求解，内附python代码及注释。本算例中初值已经给定，首先通过理论推导和代入数值分析确定了算例中属于左膨胀波右侧激波情形，而后根据已知信息进行编程，用二分法求解未知压强，进而计算所有参数并绘制图像。
采用有限差分法求解一维欧拉方程：三阶R-K时间推进——三阶迎风差分（Steger-Warming）
最新发布
ocean1702的博客
01-05 1411
采用有限差分法求解一维欧拉方程，并求解Sod激波管问题，采用三阶R-K时间推进——三阶迎风差分（Steger-Warming）来离散控制方程。
欧拉方程推导
要把所有那些负面的信息当做对自身行为的评价，而不是对自身价值的评判。
08-21 642
的方程叫做欧拉微分方程。
Riemann问题精确解及程序实现
热门推荐
meiguanhua的博客
02-24 1万+
一维无粘流动初始间断的演化问题，激波管问题（Sod激波管问题），由于该问题属于典型的间断问题，且有精确解存在，故广泛用于对比验证CFD中离散格式、数值方法的准确性，意义嘛还是蛮重要的。 1 问题描述 2 求解方法 3 解法流程 4 程序实现 5 计算结果
一维riemann 求解cfd初值
05-10
一维Riemann求解CFD初值是指在计算流体力学计算前，利用Riemann问题的方法计算出初始状态的数值，最终得到数值模拟结果。Riemann问题即初始状态下区域内不同介质之间的界面处，存在一个按照所处介质的初始状态时由高压到低压扩散的波。利用这个波，可以得到在这个波之前和之后所有介质的气体动力学参数，如压力、密度和速度等，从而得到计算所需的初值条件。 具体的求解过程分为以下步骤： 1. 定义初始状态：初始状态需要定义各点上气体的压力、密度和速度等参数。这些参数定义需要准确，因为初始状态决定了计算过程中的计算结果。 2. 定义Riemann问题：Riemann问题是在介质之间产生的波，是一个初始状态下介质间的界面处由高到低压扩散的波。通过定义初始状态下气体在不同介质中的初始状态，可以求出Riemann问题的解。 3. 求解Riemann问题：求解Riemann问题需要使用波的数学表达式，根据波的传播速度来计算波的位置。通过定义初始状态下气体在不同介质中的状态，可以计算出波的传播速度和波的位置。波的位置决定了不同介质之间的物理参数。 4. 计算求解后的初值：根据求解Riemann问题得到的解，可以计算不同介质之间的气体动力学参数，如压力、密度和速度等，从而得到初值条件。 5. 进行CFD计算：得到初值条件后，可以进行CFD计算。CFD计算过程中需要用到各个计算点上空气的气体动力学参数，如压力、密度和速度等，这些参数是前面通过Riemann问题求解得到的。CFD计算是求解一组偏微分方程，通过迭代得到最终的数值解。 总之，一维Riemann求解CFD初值是计算流体力学的重要步骤之一。通过求解Riemann问题，可以得到初值条件，确保CFD计算的准确性和可靠性。
评论 5
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

Galerkin码农选手

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付：¥1
获取中
扫码支付

您的余额不足，请更换扫码支付或充值

打赏作者

实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值