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

拉瓦尔喷管简介

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

控制方程

拉瓦尔喷管由于管道形状发生变化因此需要修改方程，这里重新推导连续性方程。
首先提一下原本一维流动的连续性方程，在管道中任意位置取一控制体。控制体左侧的变量为ρ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=ρ1u1A1−ρ0u0A0
于是有
dρdt=d(ρuA)Adx=d(ρu)dx+A′A(ρ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+A′A(ρu)=0∂(ρu)∂t+∂(ρu2+p)∂x+A′Aρu2=0∂E∂t+∂u(E+p)∂x+A′A(u(E+p))=0\frac{\partial \rho}{\partial t}+\frac{\partial (\rho u)}{\partial x}+\frac{A'}{A}(\rho u)=0\\ \frac{\partial (\rho u)}{\partial t}+\frac{\partial (\rho u^2+p)}{\partial x}+\frac{A'}{A}\rho u^2=0\\ \frac{\partial E}{\partial t}+\frac{\partial u(E+p)}{\partial x}+\frac{A'}{A}(u(E+p))=0∂t∂ρ+∂x∂(ρu)+AA′(ρu)=0∂t∂(ρu)+∂x∂(ρu2+p)+AA′ρu2=0∂t∂E+∂x∂u(E+p)+AA′(u(E+p))=0
这里假设拉瓦尔喷管的外形是双曲线满足，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∂u∂t+∂F∂x+A′AG=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,u−cu,u+c,u-cu,u+c,u−c传播，对于亚音速流动而言必有u+c>0,u−c<0u+c>0,u-c<0u+c>0,u−c<0也就是说流动总有从左到右和从右到左两个方向的信息，在边界上就是必然有从边界外到内场的信息，和内场传到边界外的信息。如果这些信息在边界上与内场不匹配就会产生反射，从而在内场出现杂波。

亚音速入口

对于亚音速入口而言，特征速度为u,u+cu,u+cu,u+c的黎曼不变量是由边界外进入内场的，而u−cu-cu−c是由内场传到入口上游去的。u,u+c,u−cu,u+c,u-cu,u+c,u−c三个特征速度对应的黎曼不变量分别为R1=cvln(p/ργ),R2=u+2cγ−1,R3=u−2cγ−1R_1=c_vln(p/\rho^\gamma),R_2=u+\frac{2c}{\gamma-1},R_3=u-\frac{2c}{\gamma-1}R1=cvln(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=u0−2c0γ−1R_3=u_0-\frac{2c_0}{\gamma-1}R3=u0−γ−12c0，边界上的R3=ub−2cbγ−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=γρ0p0,cb=γρbpb，由边界和内场的R3R_3R3相等得到
ub=u0−2c0γ−1+2cbγ−1u_b=u_0-\frac{2c_0}{\gamma-1}+\frac{2c_b}{\gamma-1}ub=u0−γ−12c0+γ−12cb这样就得到边界的ρb=1,ub=u0−2c0γ−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两个特征速度的黎曼不变量由内场推出，而u−cu-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=cvln(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=cvln(pb/ρbγ),R2=ub+γ−12cb
得到固定压力的亚音速出口的边界为ρb=(ρnγ/pn)1γ,ub=un+2cnγ−1−2cbγ−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,u−cu,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 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';
}

计算结果如下

由于我没有拉瓦尔喷管的解析解，这里仅仅验证了一下拉瓦尔喷管的关系式
(Ma−1)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γ−1Ma2)]γ−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差值的增大，在管道出口处的斜激波会逐渐变为正激波，继续增加pe−pap_e-p_ape−pa的值，正激波会向管内移动使得速度由超音速变为亚音速，这时外部扰动会影响管内流动。随着出口背景压力逐渐增大正激波位置逐渐前移，当正激波移至喉部时达到临界状态，继续增加出口背景压力管内流动将全部保持在亚音速条件下。

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

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

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

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

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