计算流体力学简介(八)——激波管问题

最新推荐文章于 2024-04-26 16:55:17 发布
最新推荐文章于 2024-04-26 16:55:17 发布
阅读量1w 收藏 57
点赞数 8
分类专栏： 偏微分方程数值解法 # 计算流体力学 C/C++
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/artorias123/article/details/104063044
版权
问题描述

求解如下方程
u = [ ρ ρ u E ] , F = [ ρ u ρ u 2 + p u ( E + p ) ] , E = p γ − 1 + 1 2 ρ u 2   ∂ u ∂ t + ∂ F ∂ x = 0 u=\left[\begin{matrix} \rho \\ \rho u\\ E \end{matrix}\right], F=\left[\begin{matrix} \rho u\\ \rho u^2+p\\ u(E+p) \end{matrix}\right], E=\frac{p}{\gamma-1}+\frac{1}{2}\rho u^2\\ \ \\ \frac{\partial u}{\partial t}+\frac{\partial F}{\partial x}=0 u=⎣⎡​ρρuE​⎦⎤​,F=⎣⎡​ρuρu2+pu(E+p)​⎦⎤​,E=γ−1p​+21​ρu2 ∂t∂u​+∂x∂F​=0
的间断初值问题称为激波管问题也叫黎曼问题。
该问题描述的是无粘理想气体在中一根管道中由压力或者速度驱动在突变的初始条件下的变化过程。
通过求解通量的雅各比矩阵 A = ∂ F ⃗ ∂ u ⃗ A=\frac{\partial {\vec F}}{\partial {\vec u}} A=∂u ∂F ​可以将前面的控制方程写成
∂ u ∂ t + A ∂ u ∂ x = 0 \frac{\partial u}{\partial t}+A\frac{\partial u}{\partial x}=0 ∂t∂u​+A∂x∂u​=0
其中
u = [ ρ ρ u E ] , A = [ 0 1 0 − 3 − γ 2 u 2 ( 3 − γ ) u γ − 1 − u ( 1 γ − 1 c 2 + 2 − γ 2 u 2 ) 1 γ − 1 c 2 + 3 − 2 γ 2 u 2 γ u ] u=\left[\begin{matrix} \rho \\ \rho u\\ E \end{matrix}\right], A=\left[\begin{matrix} 0&1&0 \\ -\frac{3-\gamma}{2}u^2&(3-\gamma)u&\gamma-1\\ -u(\frac{1}{\gamma -1}c^2+\frac{2-\gamma}{2}u^2)&\frac{1}{\gamma -1}c^2+\frac{3-2\gamma}{2}u^2&\gamma u \end{matrix}\right] u=⎣⎡​ρρuE​⎦⎤​,A=⎣⎡​0−23−γ​u2−u(γ−11​c2+22−γ​u2)​1(3−γ)uγ−11​c2+23−2γ​u2​0γ−1γu​⎦⎤​
其中 c = γ p ρ c=\sqrt{\gamma\frac{p}{\rho}} c=γρp​ ​是声速
对于该方程存在一个特殊关系
以守恒变量表示的守恒形式和非守恒形式的方程之间满足雅各比矩阵 A A A有齐次性
即 d F = A d u dF=Adu dF=Adu且 F = A u F=Au F=Au

任何与 A A A相似的矩阵 B B B以及相应的可逆变换 P P P满足 A = P − 1 B P A=P^{-1}BP A=P−1BP，都可以得到一个新的等价形式
∂ ( P u ) ∂ t + B ∂ ( P u ) ∂ x = 0 \frac{\partial (Pu)}{\partial t}+B\frac{\partial (Pu)}{\partial x}=0 ∂t∂(Pu)​+B∂x∂(Pu)​=0
于是根据原变量 ρ , u , p \rho,u,p ρ,u,p可以得到相应的等价形式
∂ u ∂ t + A ∂ u ∂ x = 0 u = [ ρ u p ] , A = [ u ρ 0 0 u 1 ρ 0 ρ c 2 u ] \frac{\partial u}{\partial t}+A\frac{\partial u}{\partial x}=0\\ u=\left[\begin{matrix} \rho \\ u\\ p \end{matrix}\right], A=\left[\begin{matrix} u&\rho&0 \\ 0&u&\frac{1}{\rho}\\ 0&\rho c^2&u \end{matrix}\right] ∂t∂u​+A∂x∂u​=0u=⎣⎡​ρup​⎦⎤​,A=⎣⎡​u00​ρuρc2​0ρ1​u​⎦⎤​
通过原变量形式的方程可以很容易计算出矩阵 A A A的特征值，从而得到与 A A A相似的对角阵，将方程解耦，得到三个方程。分别求出三个方程的特征线并对应三个变量，沿三条特征线三个变量的值不发生改变，这三个变量称为黎曼不变量。

问题求解

在这里插入图片描述
一根两端封闭的管道内左右两侧有不同的气体，中间以一物体将两侧气体分隔开，左侧气体满足 ρ L = 1 , u L = 0 , p L = 1 \rho_L=1,u_L=0,p_L=1 ρL​=1,uL​=0,pL​=1右侧气体满足 ρ L = 2 , u L = 0 , p L = 2 \rho_L=2,u_L=0,p_L=2 ρL​=2,uL​=0,pL​=2现在突然撤去中间分隔物，两侧气体将在压力作用下相互混合，最终达到平衡状态。
这里以二阶迎风的通量差分格式计算。
由原变量表达的方程中空间导数项的系数矩阵 A A A的特征值为 u + c , u − c , u u+c,u-c,u u+c,u−c,u分别表示三个黎曼不变量的传播速度，其中本问题中显然一定有亚音速的区域，因此必然有 u + c > 0 , u − c < 0 u+c>0,u-c<0 u+c>0,u−c<0说明方程本身既有向左传播的信息又有向右传播的信息。这时就需要根据信息传播方向决定使用哪一侧的迎风格式。
通量分裂的方法有很多，这里使用Lax-Friedrichs分裂
前面提到了雅各比矩阵的特征值是 u , u + c , u − c u,u+c,u-c u,u+c,u−c那么根据当地的 ρ , u , p \rho,u,p ρ,u,p可以很容易计算出三个特征值的值，找到三者中模最大的特征值 ∣ λ ∣ m a x |\lambda|_{max} ∣λ∣max​然后构造新的雅各比矩阵 A + = A + ∣ λ ∣ m a x I 2 , A − = A − ∣ λ ∣ m a x I 2 A^+=\frac{A+|\lambda|_{max}I}{2},A^-=\frac{A-|\lambda|_{max}I}{2} A+=2A+∣λ∣max​I​,A−=2A−∣λ∣max​I​这样必然有 A = A + + A − A=A^++A^- A=A++A−且 A + A^+ A+的全部特征值大于0， A − A^- A−的全部特征值小于0，利用守恒型变量雅各比矩阵的齐次性有
F + = A + u = F + ∣ λ ∣ m a x u 2 , F − = A − u = F − ∣ λ ∣ m a x u 2 F^+=A^+u=\frac{F+|\lambda|_{max}u}{2},F^-=A^-u=\frac{F-|\lambda|_{max}u}{2} F+=A+u=2F+∣λ∣max​u​,F−=A−u=2F−∣λ∣max​u​
于是得到分裂的守恒方程
∂ u ∂ t + ∂ F + ∂ x + ∂ F − ∂ x = 0 \frac{\partial u}{\partial t}+\frac{\partial F^+}{\partial x}+\frac{\partial F^-}{\partial x}=0 ∂t∂u​+∂x∂F+​+∂x∂F−​=0
式中第二项表示向右传播的信息，第三项表示向左传播的信息，分别使用 3 I − 4 E − 1 + E − 2 2 Δ x \frac{3I-4E^{-1}+E^{-2}}{2\Delta x} 2Δx3I−4E−1+E−2​和 − E 2 − 4 E 1 + 3 I 2 Δ x -\frac{E^2-4E^1+3I}{2\Delta x} −2ΔxE2−4E1+3I​进行离散。
时间使用单步推进。
Δ t = 0.005 , Δ x = 0.1 \Delta t=0.005,\Delta x=0.1 Δt=0.005,Δx=0.1
具体代码如下

#include <iostream>
#include <vector>
#include <cmath>
#define gamma 1.4
const int NE=100,//空间点数
NS=1000,
SKIP_STEP=10;//时间步数
const double rb=-5,l=10,//计算域左边界，计算域长度
dt=0.005,//时间步长
dx=l/NE;
using namespace std;
void F(vector<double> &_F,double w1,double w2,double w3)
{
    double u=w2/w1,t=u*w2;
    _F[0]=w2;
    _F[1]=(3-gamma)*t/2+(gamma-1)*w3;
    _F[2]=(1-gamma)/2*u*t+gamma*u*w3;
}
void F_div(vector<double>::iterator &f,const vector<double> &F,double w1,double w2,double w3)
{
    *f=F[0];
    f++;
    *f=F[1];
    f++;
    *f=F[2];
    f++;
}
double max(double x1,double x2)
{
    if(x1>x2) return x1;
    else return x2;
}
void advance(vector<double>& w1,vector<double>& w2,vector<double>& w3,vector<double>& F_p,vector<double>& F_m)
{
    vector<double> tF(3,0);
    double l=0;
    vector<double>::iterator f_p=F_p.begin(),f_m=F_m.begin();
    for(int i=0;i<w1.size();i++)
    {
        F(tF,w1[i],w2[i],w3[i]);
        double u=w2[i]/w1[i],p=(gamma-1)*(w3[i]-w2[i]*w2[i]/w1[i]/2),c=sqrt(gamma*p/w1[i]);
        l=max(max(abs(u+c),abs(u-c)),l);
        F_div(f_p,tF,w1[i],w2[i],w3[i]);
        F_div(f_m,tF,w1[i],w2[i],w3[i]);
    }
    for(int i=0;i<w1.size();i++)
    {
        F_p[3*i]+=l*w1[i];
        F_m[3*i]-=l*w1[i]; 
        F_p[3*i+1]+=l*w2[i];
        F_m[3*i+1]-=l*w2[i];  
        F_p[3*i+2]+=l*w3[i];
        F_m[3*i+2]-=l*w3[i]; 
    }
    f_p=F_p.begin()+3,f_m=F_m.begin()+3;
    w1[1]=w1[1]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx;
    f_p++;
    f_m++;
    w2[1]=w2[1]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx;
    f_p++;
    f_m++;
    w3[1]=w3[1]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx;
    f_p++;
    f_m++;
    for(int i=2;i<w1.size()-2;i++)
    {
        
        w1[i]=w1[i]-0.25*(3*(*f_p)-4*(*(f_p-3))+*(f_p-6))*dt/dx
                   +0.25*(*(f_m+6)-4*(*(f_m+3))+3*(*(f_m)))*dt/dx;
        f_p++;
        f_m++;
        w2[i]=w2[i]-0.25*(3*(*f_p)-4*(*(f_p-3))+*(f_p-6))*dt/dx
                   +0.25*(*(f_m+6)-4*(*(f_m+3))+3*(*(f_m)))*dt/dx;
        f_p++;
        f_m++;
        w3[i]=w3[i]-0.25*(3*(*f_p)-4*(*(f_p-3))+*(f_p-6))*dt/dx
                   +0.25*(*(f_m+6)-4*(*(f_m+3))+3*(*(f_m)))*dt/dx;
        f_p++;
        f_m++;
    }
    w1[w1.size()-2]=w1[w1.size()-2]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx;
    f_p++;
    f_m++;
    w2[w2.size()-2]=w2[w2.size()-2]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx;
    f_p++;
    f_m++;
    w3[w3.size()-2]=w3[w3.size()-2]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx;
    f_p++;
    f_m++;
    w1[0]=w1[1];
    w2[0]=-w2[1];
    w3[0]=w3[1];
    w1[w1.size()-1]=w1[w1.size()-2];
    w2[w2.size()-1]=-w2[w2.size()-2];
    w3[w3.size()-1]=w3[w3.size()-2];

}
struct val
{
    const vector<double> &w1,&w2,&w3;
    val(const vector<double>& _w1,const vector<double>& _w2,const vector<double>& _w3):w1(_w1),w2(_w2),w3(_w3){};
};
void init(vector<double> &w1,vector<double> &w2,vector<double> &w3)
{
    int i=1;
    for(;i<w1.size()/2;i++)
    {
        w1[i]=1;
        w2[i]=0;
        w3[i]=2/(gamma-1);
    }
    for(;i<w1.size()-1;i++)
    {
        w1[i]=1;
        w2[i]=0;
        w3[i]=1/(gamma-1);
    }
    w1[0]=w1[1];
    w2[0]=-w2[1];
    w3[0]=w3[1];
    w1[w1.size()-1]=w1[w1.size()-2];
    w2[w2.size()-1]=-w2[w2.size()-2];
    w3[w3.size()-1]=w3[w3.size()-2];
}
ostream& operator<<(ostream& out,const val& Q)
{
    for(int i=1;i<Q.w1.size()-1;i++)
    {
        double rho=Q.w1[i],u=Q.w2[i]/Q.w1[i],p=(gamma-1)*(Q.w3[i]-Q.w2[i]*Q.w2[i]/Q.w1[i]/2);
        out<<i*dx+rb-dx<<'\t'<<rho<<'\t'<<u<<'\t'<<p<<'\n';
    }
    return out;
}
int main()
{
    vector<double> w1(NE+3),w2(NE+3),w3(NE+3),F_p(3*NE+9),F_m(3*NE+9);
    val Q(w1,w2,w3);
    init(w1,w2,w3);
    cout<<w1.size()-2<<'\t'<<NS/SKIP_STEP<<'\t'<<rb<<'\t'<<l<<'\n';
    cout<<Q<<'\n';
    for(int i=0;i<NS;i++)
    {
        advance(w1,w2,w3,F_p,F_m);
        if(i%SKIP_STEP==0)cout<<Q<<'\n';
    }
}

计算结果如下
蓝线是压力，黑线是密度，红线是速度
在这里插入图片描述
由于我也没有激波管解析解的结果。所以也没办法判断具体结果差多少，不过从网上看到的一些激波管问题的解来看基本规律应该没有太大的问题，不过由于二阶迎风格式的耗散较小，因此间断处有一些震荡。
要想消除震荡可以利用限制器在间断附近使用耗散更大的低阶通量或者使用ENO和WENO格式来抑制震荡。

现在先测试一下一阶格式，使用如下advance函数

void advance(vector<double>& w1,vector<double>& w2,vector<double>& w3,vector<double>& F_p,vector<double>& F_m)
{
    vector<double> tF(3,0);
    double l=0;
    vector<double>::iterator f_p=F_p.begin(),f_m=F_m.begin();
    for(int i=0;i<w1.size();i++)
    {
        F(tF,w1[i],w2[i],w3[i]);
        double u=w2[i]/w1[i],p=(gamma-1)*(w3[i]-w2[i]*w2[i]/w1[i]/2),c=sqrt(gamma*p/w1[i]);
        l=max(max(abs(u+c),abs(u-c)),l);
        F_div(f_p,tF,w1[i],w2[i],w3[i]);
        F_div(f_m,tF,w1[i],w2[i],w3[i]);
    }
    for(int i=0;i<w1.size();i++)
    {
        F_p[3*i]+=l*w1[i];
        F_m[3*i]-=l*w1[i]; 
        F_p[3*i+1]+=l*w2[i];
        F_m[3*i+1]-=l*w2[i];  
        F_p[3*i+2]+=l*w3[i];
        F_m[3*i+2]-=l*w3[i]; 
    }
    f_p=F_p.begin()+3,f_m=F_m.begin()+3;
    w1[1]=w1[1]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx;
    f_p++;
    f_m++;
    w2[1]=w2[1]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx;
    f_p++;
    f_m++;
    w3[1]=w3[1]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx;
    f_p++;
    f_m++;
    for(int i=2;i<w1.size()-2;i++)
    {
        w1[i]=w1[i]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx;
        f_p++;
        f_m++;
        w2[i]=w2[i]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx;
        f_p++;
        f_m++;
        w3[i]=w3[i]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx;
        f_p++;
        f_m++;
        
    }
    w1[w1.size()-2]=w1[w1.size()-2]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx;
    f_p++;
    f_m++;
    w2[w2.size()-2]=w2[w2.size()-2]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx;
    f_p++;
    f_m++;
    w3[w3.size()-2]=w3[w3.size()-2]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx;
    f_p++;
    f_m++;
    w1[0]=w1[1];
    w2[0]=-w2[1];
    w3[0]=w3[1];
    w1[w1.size()-1]=w1[w1.size()-2];
    w2[w2.size()-1]=-w2[w2.size()-2];
    w3[w3.size()-1]=w3[w3.size()-2];
}

结果如下
在这里插入图片描述
可以看到使用高耗散的低阶格式间断处的震荡就可以被抑制住。但是耗散非常明显。

这里使用限制器来抑制间断震荡，首先相邻确定差分的比值根据信息的传播方向分成两个方向的比值
r i + 1 2 − = u i + 2 − u i + 1 u i + 1 − u i r i + 1 2 + = u i − u i − 1 u i − 1 − u i − 2 r^-_{i+\frac{1}{2}}=\frac{u_{i+2}-u_{i+1}}{u_{i+1}-u_{i}}\\ r^+_{i+\frac{1}{2}}=\frac{u_{i}-u_{i-1}}{u_{i-1}-u_{i-2}} ri+21​−​=ui+1​−ui​ui+2​−ui+1​​ri+21​+​=ui−1​−ui−2​ui​−ui−1​​
限制器函数 φ ( r ) = m i n ( 1 , ∣ r ∣ ) \varphi(r)=min(1,|r|) φ(r)=min(1,∣r∣)
格式如下
u i n + 1 = u i n − Δ t 2 Δ x ( φ ( r i + 1 2 + ) F i + 1 2 + ( 1 ) − φ ( r i − 1 2 + ) F i − 1 2 + ( 1 ) + ( 1 − φ ( r i + 1 2 − ) ) F i + 1 2 − ( 2 ) − ( 1 − φ ( r i + 1 2 − ) ) F i + 1 2 − ( 2 ) ) u_i^{n+1}=u_i^n-\frac{\Delta t}{2\Delta x} (\varphi(r^+_{i+\frac{1}{2}})F^{+(1)}_{i+\frac{1}{2}}-\varphi(r^+_{i-\frac{1}{2}})F^{+(1)}_{i-\frac{1}{2}}+ (1-\varphi(r^-_{i+\frac{1}{2}}))F^{-(2)}_{i+\frac{1}{2}}-(1-\varphi(r^-_{i+\frac{1}{2}}))F^{-(2)}_{i+\frac{1}{2}}) uin+1​=uin​−2ΔxΔt​(φ(ri+21​+​)Fi+21​+(1)​−φ(ri−21​+​)Fi−21​+(1)​+(1−φ(ri+21​−​))Fi+21​−(2)​−(1−φ(ri+21​−​))Fi+21​−(2)​)
这里做了一个简化，直接令
r i + = u i − u i − 1 u i − 1 − u i − 2 r i − = u i + 2 − u i + 1 u i + 1 − u i r^+_i=\frac{u_{i}-u_{i-1}}{u_{i-1}-u_{i-2}}\\ r^-_i=\frac{u_{i+2}-u_{i+1}}{u_{i+1}-u_{i}} ri+​=ui−1​−ui−2​ui​−ui−1​​ri−​=ui+1​−ui​ui+2​−ui+1​​
用 φ ( r i ) \varphi(r_i) φ(ri​)来决定使用一阶还是二阶格式
即
∂ u ∂ t + φ ( r + ) ∂ F + ( 2 ) ∂ x + φ ( r − ) ∂ F − ( 2 ) ∂ x + ( 1 − φ ( r + ) ) ∂ F + ( 1 ) ∂ x + ( 1 − φ ( r − ) ) ∂ F − ( 1 ) ∂ x = 0 \frac{\partial u}{\partial t}+\varphi(r^+)\frac{\partial F^{+(2)}}{\partial x}+\varphi(r^-)\frac{\partial F^{-(2)}}{\partial x}+(1-\varphi(r^+))\frac{\partial F^{+(1)}}{\partial x}+(1-\varphi(r^-))\frac{\partial F^{-(1)}}{\partial x}=0 ∂t∂u​+φ(r+)∂x∂F+(2)​+φ(r−)∂x∂F−(2)​+(1−φ(r+))∂x∂F+(1)​+(1−φ(r−))∂x∂F−(1)​=0
具体的advance函数如下

void advance(vector<double>& w1,vector<double>& w2,vector<double>& w3,vector<double>& F_p,vector<double>& F_m)
{
    vector<double> tF(3,0);
    double l=0;
    vector<double>::iterator f_p=F_p.begin(),f_m=F_m.begin();
    for(int i=0;i<w1.size();i++)
    {
        F(tF,w1[i],w2[i],w3[i]);
        double u=w2[i]/w1[i],p=(gamma-1)*(w3[i]-w2[i]*w2[i]/w1[i]/2),c=sqrt(gamma*p/w1[i]);
        l=max(max(abs(u+c),abs(u-c)),l);
        F_div(f_p,tF,w1[i],w2[i],w3[i]);
        F_div(f_m,tF,w1[i],w2[i],w3[i]);
    }
    for(int i=0;i<w1.size();i++)
    {
        F_p[3*i]+=l*w1[i];
        F_m[3*i]-=l*w1[i]; 
        F_p[3*i+1]+=l*w2[i];
        F_m[3*i+1]-=l*w2[i];  
        F_p[3*i+2]+=l*w3[i];
        F_m[3*i+2]-=l*w3[i]; 
    }
    f_p=F_p.begin()+3,f_m=F_m.begin()+3;
    w1[1]=w1[1]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx;
    f_p++;
    f_m++;
    w2[1]=w2[1]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx;
    f_p++;
    f_m++;
    w3[1]=w3[1]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx;
    f_p++;
    f_m++;
    for(int i=2;i<w1.size()-2;i++)
    {
        double r_p=(w1[i]-w1[i-1])/(w1[i-1]-w1[i-2]),
               r_m=(w1[i+2]-w1[i+1])/(w1[i+1]-w1[i]);
               if(w1[i-1]==w1[i-2]) r_p=0;
               if(w1[i+1]==w1[i]) r_m=0;
        w1[i]=w1[i]-phi(r_p)*0.25*(3*(*f_p)-4*(*(f_p-3))+*(f_p-6))*dt/dx
                   -(1-phi(r_p))*0.5*(*(f_p)-*(f_p-3))*dt/dx
                   +phi(r_m)*0.25*(*(f_m+6)-4*(*(f_m+3))+3*(*(f_m)))*dt/dx
                   -(1-phi(r_m))*0.5*(*(f_m+3)-*(f_m))*dt/dx;
        f_p++;
        f_m++;
        r_p=(w2[i]-w2[i-1])/(w2[i-1]-w2[i-2]),
        r_m=(w2[i+2]-w2[i+1])/(w2[i+1]-w2[i]);
        if(w2[i-1]==w2[i-2]) r_p=0;
        if(w2[i+1]==w2[i]) r_m=0;
        w2[i]=w2[i]-phi(r_p)*0.25*(3*(*f_p)-4*(*(f_p-3))+*(f_p-6))*dt/dx
                   -(1-phi(r_p))*0.5*(*(f_p)-*(f_p-3))*dt/dx
                   +phi(r_m)*0.25*(*(f_m+6)-4*(*(f_m+3))+3*(*(f_m)))*dt/dx
                   -(1-phi(r_m))*0.5*(*(f_m+3)-*(f_m))*dt/dx;
        f_p++;
        f_m++;
        r_p=(w3[i]-w3[i-1])/(w3[i-1]-w3[i-2]),
        r_m=(w3[i+2]-w3[i+1])/(w3[i+1]-w3[i]);
        if(w3[i-1]==w3[i-2]) r_p=0;
        if(w3[i+1]==w3[i]) r_m=0;
        w3[i]=w3[i]-phi(r_p)*0.25*(3*(*f_p)-4*(*(f_p-3))+*(f_p-6))*dt/dx
                   -(1-phi(r_p))*0.5*(*(f_p)-*(f_p-3))*dt/dx
                   +phi(r_m)*0.25*(*(f_m+6)-4*(*(f_m+3))+3*(*(f_m)))*dt/dx
                   -(1-phi(r_m))*0.5*(*(f_m+3)-*(f_m))*dt/dx;
        f_p++;
        f_m++;
    }
    w1[w1.size()-2]=w1[w1.size()-2]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx;
    f_p++;
    f_m++;
    w2[w2.size()-2]=w2[w2.size()-2]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx;
    f_p++;
    f_m++;
    w3[w3.size()-2]=w3[w3.size()-2]-0.5*(*(f_p)-*(f_p-3))*dt/dx-0.5*(*(f_m+3)-*(f_m))*dt/dx;
    f_p++;
    f_m++;
    w1[0]=w1[1];
    w2[0]=-w2[1];
    w3[0]=w3[1];
    w1[w1.size()-1]=w1[w1.size()-2];
    w2[w2.size()-1]=-w2[w2.size()-2];
    w3[w3.size()-1]=w3[w3.size()-2];
}

Δ t = 0.001 , Δ x = 0.1 \Delta t=0.001,\Delta x=0.1 Δt=0.001,Δx=0.1
计算结果如下
在这里插入图片描述
对比前面一阶格式和二阶格式可以看到，使用限制器之后二阶格式的间断震荡被明显抑制住了，而耗散也没有一阶格式那么明显。

遗世独立的理想乡_
关注
  • 8
    点赞
  • 踩
  • 57
    收藏
    觉得还不错? 一键收藏
  • 3
    评论
专栏目录
Riemann问题精确解及程序实现
meiguanhua的博客
02-24 1万+
一维无粘流动初始间断的演化问题，激波管问题（Sod激波管问题），由于该问题属于典型的间断问题，且有精确解存在，故广泛用于对比验证CFD中离散格式、数值方法的准确性，意义嘛还是蛮重要的。 1 问题描述 2 求解方法 3 解法流程 4 程序实现 5 计算结果
一维激波管问题——计算流体力学（Lax-Friedrichs格式）
05-12
本文用python3.6对一维激波管（lax激波管）问题进行实现，用到的格式为Lax-Friedrichs格式，激波管问题的初值采用了lax激波管的初值（具体见代码），推进时间为0.14s，网格数为1000
一维激波管问题——计算流体力学（Roe格式）
03-31
本文用python3.6对一维激波管lax激波管问题进行了计算，使用了Roe一阶迎风格式，激波管问题的初值采用了lax激波管的初值（具体见代码），推进时间为0.14s，网格数为1001，Roe 格式在激波处与解析解符合的较好, 而在接触间断处存在有耗散。
Sod激波管 五阶WENO 流通矢量分裂 三阶RK
最新发布
tokisaygg的博客
04-26 424
以λ+ =（（λ² + epsl²）^0 λ + .5）/2 ，λ- = （（λ² - epsl²）^0 λ + .5）/2 将三个特征值进行分裂，求h j+1/2 在（a，b，c，d，e）下的差分格式，联立三个三阶差分模板（a，b，c），（b，c，d），（c，d，e），求解出权值c1，c2，c3。4.使用特征值以及对应矩阵求出 F+=。1.构造模板 a=f j-2 ，b=f j-1 ，c= f j ，d= f j+1 ， e=f j+2。2.求出特征值 λ1=u ，λ2=u+c ，λ3=u-c。
一维激波管问题——计算流体力学（Lax-Wendroff格式）
05-12
本文用python3.6对一维激波管lax激波管问题进行了计算，使用了Lax-Wendoff格式，激波管问题的初值采用了lax激波管的初值（具体见代码），推进时间为0.14s，网格数为500，相较于Lax-Friedrichs格式在稀疏波与激波处更加符合解析解, 接触间断面处会发生耗散。
matlab-流体Sod激波管问题L-W格式、Roe格式、Van Leer格式（迎风格式）、5阶WENO格式计算其数值解
a_zxswer的博客
07-04 6719
目录Sod 激波管问题1、LW：2、Roe3、Van Leer4、五阶weno源码+文档+演示视频 Sod 激波管问题 针对如下Sod 激波管问题： 分别用L-W格式（搞定）、Roe格式（搞定）、Van Leer格式（迎风格式）、5阶WENO格式（）计算其数值解。 要求：1） 简述求解原则，注明时间推进格式，时间步长； 2） 列出全部计算公式和步骤； 3） 使用三套空间网格数，分别是100,200,400； 4） 画出t=0.14时刻密度、速度及压力的分布，并与精确解进行比较（要求同一套网格四种格式的数
求解一维Sod激波管问题
qq_45885422的博客
11-11 7933
本文是计算流体力学课程的大作业项目，全文基于Matlab代码编写，希望能够给相关专业的同学提供一些帮助。 在查找资料时，发现网上已有很多现成可用的代码，而我在本科时接触的比较多的也是CFD，想特别一些，因此在查论文时，发现了一些可用的方法，相当于复刻论文了。本文基于张涵信院士提出的NND计算格式，采用Steger_Warming矢通量分裂方法。张涵信院士是我国计算流体力学开创者和奠基人之一，长期致力于结合物理分析的计算流体力学理论和应用研究，为我国航空航天事业作出突出贡献。...
matlab-流体Sod激波管问题.zip
07-04
matlab-流体Sod激波管问题L-W格式、Roe格式、Van Leer格式（迎风格式）、5阶WENO格式计算其数值解
一维激波管（Lax shock tube）问题的数值求解
热门推荐
LSEC小陆的博客
05-28 1万+
一维激波管（Lax shock tube）问题的数值求解 文章目录一维激波管（Lax shock tube）问题的数值求解@[toc] 问题描述一般的守恒算法格式精确解常用守恒格式（先介绍标量）Lax-Friedrichs格式Roe格式Lax-Wendroff格式ENO格式简介（先介绍标量）5阶WENO格式空间离散Lax-Friedrichs通量分裂原函数方法重构左值右值空间算子和时间三阶Rung...
Riemann_accurate_激波管精确解_激波管_一维激波管问题精确解_
10-01
一维激波管问题精确解，内设置五种情况，程序内给定的初始条件可改
一维激波管问题的数值求解代码.rar
12-14
一维激波管问题的数值求解的MATLAB程序，包括Lax-Friedrichs, Lax-Wendroff,Roe和ENO等格式。说明文档可参考博客https://blog.csdn.net/lusongno1/article/details/90636623。
一维sod激波管问题weno格式
02-17
一维sod激波管问题weno格式
差分法sod激波管_SOD_sodmatlab_激波管_激波_差分法求解SOD激波管_
10-01
差分法求解SOD激波管matlab程序，可直接使用
北航计算流体力学大作业代码
06-23
北航NACA翼型扰流fortune程序
Lax格式激波管程序
12-10
Lax格式激波管模拟，matlab语言，适合计算流体力学入门。
sod.zip_Sod激波管_matlab_sod激波管问题_流体力学_激波管
09-23
计算sod激波管问题，可以自动判断和自动输出结果，只需要输入激波管问题的初始条件
sob激波管黎曼解matlab程序.rar
04-15
计算流体力学大作业，sob激波管黎曼解matlab程序。
激波管的解法——有限差分方法
qq_43222660的博客
04-14 780
方程及状态参数 进排气管内一维变截面非定常控制方程为: {∂(ρA)∂t+∂(ρuA)∂x=0∂(ρuA)∂t+∂(ρu2A+pA)∂x=−ρGA+p∂A∂x∂(ρe0A)∂t+∂(ρuh0A)∂x=ρq˙A(1) \left\{ \begin{array}{l} \frac{{\partial \left( {\rho A} \right)}}{{\partial t}} + \frac{{\partial \left( {\rho uA} \right)}}{{\partial x}} = 0\\ \f
计算空气动力学sod激波管问题
05-02
空气动力学Sod激波管问题是一个经典的计算流体力学问题。在这个问题中，考虑到在管道中某些初始条件下的气体流动，使用数值方法估算当流动达到稳定状态时，气体流动的压力、密度和速度等参数。这个问题之所以重要，是因为它有助于理解气体流体动力学的基本原理和应用技术。 计算空气动力学Sod激波管问题需要使用数值方法进行近似解，通常使用计算流体力学(CFD) 方法来模拟气体流动。在求解过程中，要确定模型的边界条件和初始条件，选择数值离散化方法和数值计算方案，并设置计算参数和误差范围等。随着计算机技术和数值方法不断发展，对这些参数的选择和优化方法也越来越成熟。 计算空气动力学Sod激波管问题在航空航天、能源工程、材料科学、医学等领域都有很大的应用，例如在发动机喷射器、燃烧室、超声波治疗、血管瘤治疗等方面，都需要使用CFD方法来分析气体流动和热传递。因此，研究计算空气动力学Sod激波管问题具有重要的意义和应用价值。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 3
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

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

余额充值