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

拉瓦尔喷管简介

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

控制方程

拉瓦尔喷管由于管道形状发生变化因此需要修改方程，这里重新推导连续性方程。
首先提一下原本一维流动的连续性方程，在管道中任意位置取一控制体。控制体左侧的变量为${\rho }_0,u_0,A_0$，右侧为${\rho }_1,u_1,A_1$。则控制体内的总质量变化应为
$\frac{A\rho dx}{dt}={\rho }_1{u}_1{A}_1-{\rho }_0{u}_0{A}_0$
于是有
$\frac{d\rho }{dt}=\frac{d(\rho uA)}{Adx}=\frac{d(\rho u)}{dx}+\frac{A^{\prime }}{A}(\rho u)$
动量方程和能量方程方法基本类似，其中能量方程和连续方程分析完全相同，动量方程略为复杂一些，在安德森的计算流体力学中有详细推导，具体方程如下
$$\begin{gathered} \frac{\partial \rho }{\partial t}+\frac{\partial (\rho u)}{\partial x}+\frac{A^{\prime }}{A}(\rho u)=0 \\ \frac{\partial (\rho u)}{\partial t}+\frac{\partial (\rho u^2+p)}{\partial x}+\frac{A^{\prime }}{A}\rho u^2=0 \\ \frac{\partial E}{\partial t}+\frac{\partial u(E+p)}{\partial x}+\frac{A^{\prime }}{A}(u(E+p))=0 \end{gathered}$$
这里假设拉瓦尔喷管的外形是双曲线满足，$A(x)=\sqrt{\frac{x^2}{25}+1}$，$A^{\prime }=\frac{x}{\sqrt{25(x^2+25)}}$
$$\begin{gathered} u=\left[\begin{array}{c}\rho \\ \rho u\\ E\end{array}\right],F=\left[\begin{array}{c}\rho u\\ \rho u^2+p\\ u(E+p)\end{array}\right],G=\left[\begin{array}{c}\rho u\\ \rho u^2\\ u(E+p)\end{array}\right],E=\frac{p}{\gamma -1}+\frac{1}{2}\rho u^2 \\ \frac{\partial u}{\partial t}+\frac{\partial F}{\partial x}+\frac{A^{\prime }}{A}G=0 \end{gathered}$$

数值方法

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

边界条件

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

亚音速入口

对于亚音速入口而言，特征速度为$u,u+c$的黎曼不变量是由边界外进入内场的，而$u-c$是由内场传到入口上游去的。$u,u+c,u-c$三个特征速度对应的黎曼不变量分别为$R_1=c_{v}ln(p/{\rho }^{\gamma }),R_2=u+\frac{2c}{\gamma -1},R_3=u-\frac{2c}{\gamma -1}$
我们给定入口通常直接给定原变量形式即$\rho ,u,p$这三个参数和三个黎曼不变量一一对应。前面提到边界的三个黎曼变量两个变量由边界直接认为给出，一个由内场决定。换言之就是入口的$\rho ,u,p$中只能固定两个量，另一个量由内场的值和边界给出的两个值计算出来。我这里就固定$\rho ,p$两个变量，固定$\rho =1,p=2$，速度利用黎曼不变量得出。记内场第一个点的${\rho }_0,u_0,p_0$计算出的$R_3=u_0-\frac{2c_0}{\gamma -1}$，边界上的$R_3=u_b-\frac{2c_b}{\gamma -1}$其中$c_0=\sqrt{\gamma \frac{p_0}{{\rho }_0}},c_b=\sqrt{\gamma \frac{p_b}{{\rho }_b}}$，由边界和内场的$R_3$相等得到
$u_b=u_0-\frac{2c_0}{\gamma -1}+\frac{2c_b}{\gamma -1}$这样就得到边界的${\rho }_b=1,u_b=u_0-\frac{2c_0}{\gamma -1}+\frac{2c_b}{\gamma -1},p_b=2$于是就得到了固定压力和密度的亚音速入口的边界条件。

亚音速出口

亚音速出口与亚音速入口相反$u,u+c$两个特征速度的黎曼不变量由内场推出，而$u-c$由边界给定，换言之就是$\rho ,u,p$中两个变量由内场计算，一个变量给定，这里给定$p=1$，于是由
内场最后一个点上$R_1=c_{v}ln(p_n/{\rho }_n^{\gamma }),R_2=u_n+\frac{2c_n}{\gamma -1}$
边界上$R_1=c_{v}ln(p_b/{\rho }_b^{\gamma }),R_2=u_b+\frac{2c_b}{\gamma -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$

超音速出口

超音速出口相比于亚音速出口就容易多了，因为超音速条件下$u,u+c,u-c$都为正，因此信息只出不进，于是直接有超音速出口边界为${\rho }_b={\rho }_n,u_b=u_n,p_b=p_n$

本文使用的边界

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

数值计算

计算代码如下

#include <iostream>
#include <vector>
#include <cmath>
#define gamma 1.4
#define pa_in 6
#define p_in 2
#define u_in 1
#define p_out 1
const int NE=100,//空间点数
NS=15000,
SKIP_STEP=1500;//时间步数
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);
}
double A(double x)
{
    return sqrt(x*x/25+1);
}
double dAdx(double x)
{
    return x/sqrt(625+25*x*x);
}
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]*0.5),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]; 
    }
    double x=rb;
    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++)
    {
        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-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++;
    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<-1e-7)
    {
        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);
    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]=rho;
    w2[0]=u*rho;
    w3[0]=p/(gamma-1)+0.5*rho*u*u;
    
}
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]=p_out;
        w2[i]=0;
        w3[i]=p_out/(gamma-1); 
    }
    i=1;
    w1[0]=1;
    w2[0]=0;
    w3[0]=p_in/(gamma-1)+0.5*w2[0]*w2[0]/w1[0];
}
ostream& operator<<(ostream& out,const val& Q)
{
    double x=rb-dx*2;
    for(int i=0;i<Q.w1.size();i++)
    {
        x+=dx;
        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);
		double c=sqrt(gamma*p/rho);
        double Ma=u/c;
        out<<i*dx+rb-dx<<'\t'<<rho<<'\t'<<u<<'\t'<<p<<'\t'<<c<<'\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()<<'\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';
}

计算结果如下

由于我没有拉瓦尔喷管的解析解，这里仅仅验证了一下拉瓦尔喷管的关系式
$$\begin{gathered} (Ma-1)u^{\prime }/u=A^{\prime }/A \\ (\frac{A}{A^{\ast }}{)}^2=\frac{1}{M{a}^2}[\frac{2}{\gamma +1}(1+\frac{\gamma -1}{2}M{a}^2){]}^{\frac{\gamma +1}{\gamma -1}} \end{gathered}$$
将上式左右两端相减绘制曲线如下，最终结果基本满足该关系，因此应该问题不大。

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

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

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

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

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

这里我试着算了一下过膨胀状态

这是入口压力1.1，出口压力1时的状态，可以看到一开始管内出现了激波，说明这时是过膨胀状态，在喉部后侧存在一道正激波，使得速度变为亚音速。但是我试算了一下不论怎么减小入口压力只要出口压力是1，入口压力大于1管道内就不可能出现全部亚音速的情况，似乎正激波最终会固定到$x=2$的位置处。具体原因暂时不太清楚，目前估计是入口边界的原因，似乎固定压力和密度的亚音速入口决定了喉部的速度。我使用真实的气体参数作为参数进行计算，入口设置压力为1.001倍的标准大气压，密度由气体状态方程利用压力和温度计算而来，温度设为25°C(298.15K)，出口是标准大气压，这是计算收敛到全场亚音速的结果，我懒得做无量纲了这里就不展示了。