计算流体力学简介(九)——拉瓦尔喷管模拟

拉瓦尔喷管简介

计算流体力学简介(九)——拉瓦尔喷管模拟
如图所示拉瓦尔喷管为以收缩-扩张管道,入口速度为亚音速,压缩性较差,在收缩段受管壁收缩挤压作用加速,在最窄的喉部达到音速。随着气体速度增大压缩性逐渐增加,在喉部以后管道扩张使得气体迅速膨胀,密度减小,流速继续增加,达到超音速。

控制方程

拉瓦尔喷管由于管道形状发生变化因此需要修改方程,这里重新推导连续性方程。
首先提一下原本一维流动的连续性方程,在管道中任意位置取一控制体。控制体左侧的变量为ρ0,u0,A0\rho_0,u_0,A_0ρ0​,u0​,A0​,右侧为ρ1,u1,A1\rho_1,u_1,A_1ρ1​,u1​,A1​。则控制体内的总质量变化应为
Aρdxdt=ρ1u1A1ρ0u0A0\frac{A\rho dx}{dt}=\rho_1u_1A_1-\rho_0u_0A_0dtAρdx​=ρ1​u1​A1​−ρ0​u0​A0​
于是有
dρdt=d(ρuA)Adx=d(ρu)dx+AA(ρu)\frac{d\rho}{dt}=\frac{d(\rho uA)}{Adx}=\frac{d(\rho u)}{dx}+\frac{A'}{A}(\rho u)dtdρ​=Adxd(ρuA)​=dxd(ρu)​+AA′​(ρu)
动量方程和能量方程方法基本类似,其中能量方程和连续方程分析完全相同,动量方程略为复杂一些,在安德森的计算流体力学中有详细推导,具体方程如下
ρt=(ρu)x+AA(ρu)(ρu)t=(ρu2+p)x+AAρu2Et=u(E+p)x+AA(u(E+p))\frac{\partial \rho}{\partial t}=\frac{\partial (\rho u)}{\partial x}+\frac{A'}{A}(\rho u)\\ \frac{\partial (\rho u)}{\partial t}=\frac{\partial (\rho u^2+p)}{\partial x}+\frac{A'}{A}\rho u^2\\ \frac{\partial E}{\partial t}=\frac{\partial u(E+p)}{\partial x}+\frac{A'}{A}(u(E+p))∂t∂ρ​=∂x∂(ρu)​+AA′​(ρu)∂t∂(ρu)​=∂x∂(ρu2+p)​+AA′​ρu2∂t∂E​=∂x∂u(E+p)​+AA′​(u(E+p))
这里假设拉瓦尔喷管的外形是双曲线满足,A(x)=x225+1A(x)=\sqrt{\frac{x^2}{25}+1}A(x)=25x2​+1​,A=x25(x2+25){A'}=\frac{x}{\sqrt{25(x^2+25)}}A′=25(x2+25)​x​
u=[ρρuE],F=[ρuρu2+pu(E+p)],G=[ρuρu2u(E+p)],E=pγ1+12ρu2 ut+Fx+AAG=0u=\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], G=\left[\begin{matrix} \rho u \\ \rho u^2\\ 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}+\frac{A'}{A}G=0u=⎣⎡​ρρuE​⎦⎤​,F=⎣⎡​ρuρu2+pu(E+p)​⎦⎤​,G=⎣⎡​ρuρu2u(E+p)​⎦⎤​,E=γ−1p​+21​ρu2 ∂t∂u​+∂x∂F​+AA′​G=0

数值方法

使用前面的激波管问题中带限制器的格式,分裂方式也不需要修改。

边界条件

这里刚好趁着这个问题再说一下无反射边界。激波管问题中提到对于无粘可压缩流动,流动中有三个黎曼不变量,分别以u,u+c,ucu,u+c,u-cu,u+c,u−c传播,对于亚音速流动而言必有u+c>0,uc<0u+c>0,u-c<0u+c>0,u−c<0也就是说流动总有从左到右和从右到左两个方向的信息,在边界上就是必然有从边界外到内场的信息,和内场传到边界外的信息。如果这些信息在边界上与内场不匹配就会产生反射,从而在内场出现杂波。

亚音速入口

对于亚音速入口而言,特征速度为u,u+cu,u+cu,u+c的黎曼不变量是由边界外进入内场的,而ucu-cu−c是由内场传到入口上游去的。u,u+c,ucu,u+c,u-cu,u+c,u−c三个特征速度对应的黎曼不变量分别为R1=cvln(p/ργ),R2=u+2cγ1,R3=u2cγ1R_1=c_vln(p/\rho^\gamma),R_2=u+\frac{2c}{\gamma-1},R_3=u-\frac{2c}{\gamma-1}R1​=cv​ln(p/ργ),R2​=u+γ−12c​,R3​=u−γ−12c​
我们给定入口通常直接给定原变量形式即ρ,u,p\rho,u,pρ,u,p这三个参数和三个黎曼不变量一一对应。前面提到边界的三个黎曼变量两个变量由边界直接认为给出,一个由内场决定。换言之就是入口的ρ,u,p\rho,u,pρ,u,p中只能固定两个量,另一个量由内场的值和边界给出的两个值计算出来。我这里就固定ρ,p\rho,pρ,p两个变量,固定ρ=1,p=2\rho=1,p=2ρ=1,p=2,速度利用黎曼不变量得出。记内场第一个点的ρ0,u0,p0\rho_0,u_0,p_0ρ0​,u0​,p0​计算出的R3=u02c0γ1R_3=u_0-\frac{2c_0}{\gamma-1}R3​=u0​−γ−12c0​​,边界上的R3=ub2cbγ1R_3=u_b-\frac{2c_b}{\gamma-1}R3​=ub​−γ−12cb​​其中c0=γp0ρ0,cb=γpbρbc_0=\sqrt{\gamma\frac{p_0}{\rho_0}},c_b=\sqrt{\gamma\frac{p_b}{\rho_b}}c0​=γρ0​p0​​​,cb​=γρb​pb​​​,由边界和内场的R3R_3R3​相等得到
ub=u02c0γ1+2cbγ1u_b=u_0-\frac{2c_0}{\gamma-1}+\frac{2c_b}{\gamma-1}ub​=u0​−γ−12c0​​+γ−12cb​​这样就得到边界的ρb=1,ub=u02c0γ1+2cbγ1,pb=2\rho_b=1,u_b=u_0-\frac{2c_0}{\gamma-1}+\frac{2c_b}{\gamma-1},p_b=2ρb​=1,ub​=u0​−γ−12c0​​+γ−12cb​​,pb​=2于是就得到了固定压力和密度的亚音速入口的边界条件。

亚音速出口

亚音速出口与亚音速入口相反u,u+cu,u+cu,u+c两个特征速度的黎曼不变量由内场推出,而ucu-cu−c由边界给定,换言之就是ρ,u,p\rho,u,pρ,u,p中两个变量由内场计算,一个变量给定,这里给定p=1p=1p=1,于是由
内场最后一个点上R1=cvln(pn/ρnγ),R2=un+2cnγ1R_1=c_vln(p_n/\rho_n^\gamma),R_2=u_n+\frac{2c_n}{\gamma-1}R1​=cv​ln(pn​/ρnγ​),R2​=un​+γ−12cn​​
边界上R1=cvln(pb/ρbγ),R2=ub+2cbγ1R_1=c_vln(p_b/\rho_b^\gamma),R_2=u_b+\frac{2c_b}{\gamma-1}R1​=cv​ln(pb​/ρbγ​),R2​=ub​+γ−12cb​​
得到固定压力的亚音速出口的边界为ρb=(ρnγ/pn)1γ,ub=un+2cnγ12cbγ1,pb=1\rho_b=(\rho_n^\gamma/p_n)^{\frac{1}{\gamma}},u_b=u_n+\frac{2c_n}{\gamma-1}-\frac{2c_b}{\gamma-1},p_b=1ρb​=(ρnγ​/pn​)γ1​,ub​=un​+γ−12cn​​−γ−12cb​​,pb​=1

超音速出口

超音速出口相比于亚音速出口就容易多了,因为超音速条件下u,u+c,ucu,u+c,u-cu,u+c,u−c都为正,因此信息只出不进,于是直接有超音速出口边界为ρb=ρn,ub=un,pb=pn\rho_b=\rho_n,u_b=u_n,p_b=p_nρb​=ρn​,ub​=un​,pb​=pn​

本文使用的边界

本文的入口是亚音速入口,而出口由于初场给定的是全场速度为0,密度压力均为1,因此开始流体受压力作用从0开始加速,开始时出口是亚音速的,但是随着流动不断发展速度不断增大,最终达到超音速,于是在出口处做一次判断,内场最后一个点的un<cnu_n<c_nun​<cn​时使用亚音速出口,否则使用超音速出口。

数值计算

计算代码如下

#include <iostream>
#include <vector>
#include <cmath>
#define gamma 1.4
#define p_in 2
#define p_out 1
const int NE=100,//空间点数
NS=50000,
SKIP_STEP=500;//时间步数
const double rb=-5,l=10,//计算域左边界,计算域长度
dt=0.001,//时间步长
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;
}
double phi(double r)
{
    if(abs(r)>1) return 1;
    else return abs(r);
}
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++;
    double x=rb-dx;
    for(int i=2;i<w1.size()-2;i++)
    {
        x+=dx;
        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)
                   -dt*dAdx(x)/A(x)*(*(f_p)+*(f_m))*0.5;
        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)
                   -dt*dAdx(x)/A(x)*w2[i]*w2[i]/w1[i];
        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)
                   -dt*dAdx(x)/A(x)*(*(f_p)+*(f_m))*0.5;
        f_p++;
        f_m++;
        
    }
    w1[w1.size()-2]=w1[w1.size()-2]-0.5*(*(f_p)-*(f_p-3))*dt/dx-(*(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-(*(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-(*(f_m+3)-*(f_m))*dt/dx;
    f_p++;
    f_m++;
    int i=w1.size()-2;
    double rho=w1[i],u=w2[i]/w1[i],p=(gamma-1)*(w3[i]-w2[i]*w2[i]/w1[i]*0.5),c=sqrt(gamma*p/w1[i]);
    if(u<c)
    {
        double R1=p/pow(rho,gamma),R23=2/(gamma-1)*c;
        p=p_out;
        rho=pow(p/R1,1/gamma);
        c=sqrt(gamma*p/rho);
        u=u+R23-2*c/(gamma-1);
        w1[w1.size()-1]=rho;
        w2[w2.size()-1]=u*rho;
        w3[w3.size()-1]=p/(gamma-1)+0.5*u*rho*u;
    }
    else
    {
        w1[w1.size()-1]=w1[w1.size()-2];
        w2[w2.size()-1]=w2[w2.size()-2];
        w3[w3.size()-1]=w3[w3.size()-2];
    }
    i=1;
    rho=w1[i],u=w2[i]/w1[i],p=(gamma-1)*(w3[i]-w2[i]*w2[i]/w1[i]*0.5),c=sqrt(gamma*p/rho);
    p=p_in;
    rho=1;
    double R3=u-2*c/(gamma-1);
    c=sqrt(gamma*p/rho);
    w1[0]=1;
    w2[0]=R3+2*c/(gamma-1);
    w3[0]=p/(gamma-1)+0.5*w2[i]*w2[i]/w1[i];
    //w3[0]=p/(gamma-1)+0.5;
}
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=0;
    for(;i<w1.size();i++)
    {
        w1[i]=1;
        w2[i]=0;
        w3[i]=1/(gamma-1);        
    }
    i=1;
    double rho=w1[i],u=w2[i]/w1[i],p=(gamma-1)*(w3[i]-w2[i]*w2[i]/w1[i]*0.5),c=sqrt(gamma*p/w1[i]);
    double R3=u-2*c/(gamma-1);
    p=p_in;
    rho=1;
    c=sqrt(gamma*p/rho);
    u=R3+2*c/(gamma-1);
    w1[0]=1;
    w2[0]=u;
    w3[0]=p/(gamma-1)+0.5*w2[i]*w2[i]/w1[i];
    //w3[0]=p/(gamma-1)+0.5;
}
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';
    }
    //cout<<Q<<'\n';
}

计算结果如下
计算流体力学简介(九)——拉瓦尔喷管模拟
由于我没有拉瓦尔喷管的解析解,这里仅仅验证了一下拉瓦尔喷管的关系式
(Ma1)u/u=A/A(AA)2=1Ma2[2γ+1(1+γ12Ma2)]γ+1γ1(Ma-1)u'/u=A'/A\\ (\frac{A}{A^*})^2=\frac{1}{Ma^2}[\frac{2}{\gamma+1}(1+\frac{\gamma-1}{2}Ma^2)]^{\frac{\gamma+1}{\gamma-1}}(Ma−1)u′/u=A′/A(A∗A​)2=Ma21​[γ+12​(1+2γ−1​Ma2)]γ−1γ+1​
将上式左右两端相减绘制曲线如下,最终结果基本满足该关系,因此应该问题不大。
计算流体力学简介(九)——拉瓦尔喷管模拟

拉瓦尔喷管工作状态与正激波

拉瓦尔喷管根据喷管道出口流体压力(pep_epe​)和背景(pap_apa​)压力关系,可以分为多种工作状态:
1.pa=pep_a=p_epa​=pe​理想工作状态,这时气体恰好完全膨胀,管道出口为超音速流动,并且离开管道后也不发生膨胀和压缩,外部扰动无法影响管内流动状态。
2.pa>pep_a>p_epa​>pe​欠膨胀状态,气体没有完全膨胀,管道出口为超音速流动,气体离开管道后会继续膨胀一段,直至压力达到背景压力,这时同样外部扰动无法传入内部
3.pa<pep_a<p_epa​<pe​过膨胀状态,气体在出口处过膨胀,当pa,pep_a,p_epa​,pe​差别不大时,管口仍为超音速流动,但是由于出口处气体压力小于背景压力,因此会在管口处形成一系列斜激波,最终使得流出气体压力和背景压力相等;随着pa,pep_a,p_epa​,pe​差值的增大,在管道出口处的斜激波会逐渐变为正激波,继续增加pepap_e-p_ape​−pa​的值,正激波会向管内移动使得速度由超音速变为亚音速,这时外部扰动会影响管内流动。随着出口背景压力逐渐增大正激波位置逐渐前移,当正激波移至喉部时达到临界状态,继续增加出口背景压力管内流动将全部保持在亚音速条件下。

根据前面的分析,理想状态和欠膨胀状态出口都是超音速流动,因此管道内部计算不受外部状态影响,至于背景压力如何管道内流动都是不变的。唯一有区别的是管道外是否发生膨胀,管内流动始终相同。

过膨胀状态则略有不同,由于管内可能会形成一道正激波,正激波后流动全部是亚音速流动,亚音速出口外部扰动将影响内部流动,这是外部压力将使得管内流动状态发生改变。

在计算上的区别就是,欠膨胀和理想状态最终都是超音速出口,而过膨胀出口总是亚音速的。

这里我试着算了一下过膨胀状态
计算流体力学简介(九)——拉瓦尔喷管模拟
这是入口压力1.1,出口压力1时的状态,可以看到一开始管内出现了激波,说明这时是过膨胀状态,在喉部后侧存在一道正激波,使得速度变为亚音速。但是我试算了一下不论怎么减小入口压力只要出口压力是1,入口压力大于1管道内就不可能出现全部亚音速的情况,似乎正激波最终会固定到x=2x=2x=2的位置处。具体原因暂时不太清楚,目前估计是入口边界的原因,似乎固定压力和密度的亚音速入口决定了喉部的速度。

计算流体力学简介(九)——拉瓦尔喷管模拟计算流体力学简介(九)——拉瓦尔喷管模拟 唯我电磁炮永世长存 发布了23 篇原创文章 · 获赞 9 · 访问量 1万+ 私信 关注
上一篇:state processor api


下一篇:libsvm之核函数