拉瓦爾噴管簡介

如圖所示拉瓦爾噴管為以收縮-擴張管道，入口速度為亞音速，壓縮性較差，在收縮段受管壁收縮擠壓作用加速，在最窄的喉部達到音速。隨著氣體速度增大壓縮性逐漸增加，在喉部以后管道擴張使得氣體迅速膨脹，密度減小，流速繼續增加，達到超音速。

控制方程

拉瓦爾噴管由于管道形狀發生變化因此需要修改方程，這里重新推導連續性方程。
首先提一下原本一維流動的連續性方程，在管道中任意位置取一控制體。控制體左側的變量為${\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{?\rho }{?t}+\frac{?(\rho u)}{?x}+\frac{A^{\prime }}{A}(\rho u)=0 \\ \frac{?(\rho u)}{?t}+\frac{?(\rho u^2+p)}{?x}+\frac{A^{\prime }}{A}\rho u^2=0 \\ \frac{?E}{?t}+\frac{?u(E+p)}{?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=?\begin{array}{c}\rho \\ \rho u\\ E\end{array}?,F=?\begin{array}{c}\rho u\\ \rho u^2+p\\ u(E+p)\end{array}?,G=?\begin{array}{c}\rho u\\ \rho u^2\\ u(E+p)\end{array}?,E=\frac{p}{\gamma ?1}+\frac{1}{2}\rho u^2 \\ \frac{?u}{?t}+\frac{?F}{?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^{?}}{)}^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)，出口是標準大氣壓，這是計算收斂到全場亞音速的結果，我懶得做無量綱了這里就不展示了。

總結

以上是生活随笔為你收集整理的计算流体力学简介(九)——拉瓦尔喷管模拟的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。