拉瓦尔喷管简介
如图所示拉瓦尔喷管为以收缩-扩张管道，入口速度为亚音速，压缩性较差，在收缩段受管壁收缩挤压作用加速，在最窄的喉部达到音速。随着气体速度增大压缩性逐渐增加，在喉部以后管道扩张使得气体迅速膨胀，密度减小，流速继续增加，达到超音速。
控制方程
拉瓦尔喷管由于管道形状发生变化因此需要修改方程，这里重新推导连续性方程。
首先提一下原本一维流动的连续性方程，在管道中任意位置取一控制体。控制体左侧的变量为ρ0,u0,A0，右侧为ρ1,u1,A1。则控制体内的总质量变化应为
dtAρdx=ρ1u1A1−ρ0u0A0
于是有
dtdρ=Adxd(ρuA)=dxd(ρu)+AA′(ρu)
动量方程和能量方程方法基本类似，其中能量方程和连续方程分析完全相同，动量方程略为复杂一些，在安德森的计算流体力学中有详细推导，具体方程如下
∂t∂ρ=∂x∂(ρu)+AA′(ρu)∂t∂(ρu)=∂x∂(ρu2+p)+AA′ρu2∂t∂E=∂x∂u(E+p)+AA′(u(E+p))
这里假设拉瓦尔喷管的外形是双曲线满足，A(x)=25x2+1，A′=25(x2+25)x
u=⎣⎡ρρ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,u−c传播，对于亚音速流动而言必有u+c>0,u−c<0也就是说流动总有从左到右和从右到左两个方向的信息，在边界上就是必然有从边界外到内场的信息，和内场传到边界外的信息。如果这些信息在边界上与内场不匹配就会产生反射，从而在内场出现杂波。
亚音速入口
对于亚音速入口而言，特征速度为u,u+c的黎曼不变量是由边界外进入内场的，而u−c是由内场传到入口上游去的。u,u+c,u−c三个特征速度对应的黎曼不变量分别为R1=cvln(p/ργ),R2=u+γ−12c,R3=u−γ−12c
我们给定入口通常直接给定原变量形式即ρ,u,p这三个参数和三个黎曼不变量一一对应。前面提到边界的三个黎曼变量两个变量由边界直接认为给出，一个由内场决定。换言之就是入口的ρ,u,p中只能固定两个量，另一个量由内场的值和边界给出的两个值计算出来。我这里就固定ρ,p两个变量，固定ρ=1,p=2，速度利用黎曼不变量得出。记内场第一个点的ρ0,u0,p0计算出的R3=u0−γ−12c0，边界上的R3=ub−γ−12cb其中c0=γρ0p0,cb=γρbpb，由边界和内场的R3相等得到
ub=u0−γ−12c0+γ−12cb这样就得到边界的ρb=1,ub=u0−γ−12c0+γ−12cb,pb=2于是就得到了固定压力和密度的亚音速入口的边界条件。
亚音速出口
亚音速出口与亚音速入口相反u,u+c两个特征速度的黎曼不变量由内场推出，而u−c由边界给定，换言之就是ρ,u,p中两个变量由内场计算，一个变量给定，这里给定p=1，于是由
内场最后一个点上R1=cvln(pn/ρnγ),R2=un+γ−12cn
边界上R1=cvln(pb/ρbγ),R2=ub+γ−12cb
得到固定压力的亚音速出口的边界为ρb=(ρnγ/pn)γ1,ub=un+γ−12cn−γ−12cb,pb=1
超音速出口
超音速出口相比于亚音速出口就容易多了，因为超音速条件下u,u+c,u−c都为正，因此信息只出不进，于是直接有超音速出口边界为ρb=ρn,ub=un,pb=pn
本文使用的边界
本文的入口是亚音速入口，而出口由于初场给定的是全场速度为0，密度压力均为1，因此开始流体受压力作用从0开始加速，开始时出口是亚音速的，但是随着流动不断发展速度不断增大，最终达到超音速，于是在出口处做一次判断，内场最后一个点的un<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';
}
计算结果如下
由于我没有拉瓦尔喷管的解析解，这里仅仅验证了一下拉瓦尔喷管的关系式
(Ma−1)u′/u=A′/A(A∗A)2=Ma21[γ+12(1+2γ−1Ma2)]γ−1γ+1
将上式左右两端相减绘制曲线如下，最终结果基本满足该关系，因此应该问题不大。
拉瓦尔喷管工作状态与正激波
拉瓦尔喷管根据喷管道出口流体压力(pe)和背景(pa)压力关系，可以分为多种工作状态:
1.pa=pe理想工作状态，这时气体恰好完全膨胀，管道出口为超音速流动，并且离开管道后也不发生膨胀和压缩，外部扰动无法影响管内流动状态。
2.pa>pe欠膨胀状态，气体没有完全膨胀，管道出口为超音速流动，气体离开管道后会继续膨胀一段，直至压力达到背景压力，这时同样外部扰动无法传入内部
3.pa<pe过膨胀状态，气体在出口处过膨胀，当pa,pe差别不大时，管口仍为超音速流动，但是由于出口处气体压力小于背景压力，因此会在管口处形成一系列斜激波，最终使得流出气体压力和背景压力相等；随着pa,pe差值的增大，在管道出口处的斜激波会逐渐变为正激波，继续增加pe−pa的值，正激波会向管内移动使得速度由超音速变为亚音速，这时外部扰动会影响管内流动。随着出口背景压力逐渐增大正激波位置逐渐前移，当正激波移至喉部时达到临界状态，继续增加出口背景压力管内流动将全部保持在亚音速条件下。
根据前面的分析，理想状态和欠膨胀状态出口都是超音速流动，因此管道内部计算不受外部状态影响，至于背景压力如何管道内流动都是不变的。唯一有区别的是管道外是否发生膨胀，管内流动始终相同。
过膨胀状态则略有不同，由于管内可能会形成一道正激波，正激波后流动全部是亚音速流动，亚音速出口外部扰动将影响内部流动，这是外部压力将使得管内流动状态发生改变。
在计算上的区别就是，欠膨胀和理想状态最终都是超音速出口，而过膨胀出口总是亚音速的。
这里我试着算了一下过膨胀状态
这是入口压力1.1，出口压力1时的状态，可以看到一开始管内出现了激波，说明这时是过膨胀状态，在喉部后侧存在一道正激波，使得速度变为亚音速。但是我试算了一下不论怎么减小入口压力只要出口压力是1，入口压力大于1管道内就不可能出现全部亚音速的情况，似乎正激波最终会固定到x=2的位置处。具体原因暂时不太清楚，目前估计是入口边界的原因，似乎固定压力和密度的亚音速入口决定了喉部的速度。