拉瓦爾噴管簡(jiǎn)介

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

控制方程

拉瓦爾噴管由于管道形狀發(fā)生變化因此需要修改方程，這里重新推導(dǎo)連續(xù)性方程。
首先提一下原本一維流動(dòng)的連續(xù)性方程，在管道中任意位置取一控制體。控制體左側(cè)的變量為${\rho }_0,u_0,A_0$，右側(cè)為${\rho }_1,u_1,A_1$。則控制體內(nèi)的總質(zhì)量變化應(yīng)為
$\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)$
動(dòng)量方程和能量方程方法基本類似，其中能量方程和連續(xù)方程分析完全相同，動(dòng)量方程略為復(fù)雜一些，在安德森的計(jì)算流體力學(xué)中有詳細(xì)推導(dǎo)，具體方程如下
$$\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}$$
這里假設(shè)拉瓦爾噴管的外形是雙曲線滿足，$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}$$

數(shù)值方法

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

邊界條件

這里剛好趁著這個(gè)問題再說(shuō)一下無(wú)反射邊界。激波管問題中提到對(duì)于無(wú)粘可壓縮流動(dòng)，流動(dòng)中有三個(gè)黎曼不變量，分別以$u,u+c,u?c$傳播，對(duì)于亞音速流動(dòng)而言必有$u+c>0,u?c<0$也就是說(shuō)流動(dòng)總有從左到右和從右到左兩個(gè)方向的信息，在邊界上就是必然有從邊界外到內(nèi)場(chǎng)的信息，和內(nèi)場(chǎng)傳到邊界外的信息。如果這些信息在邊界上與內(nèi)場(chǎng)不匹配就會(huì)產(chǎn)生反射，從而在內(nèi)場(chǎng)出現(xiàn)雜波。

亞音速入口

對(duì)于亞音速入口而言，特征速度為$u,u+c$的黎曼不變量是由邊界外進(jìn)入內(nèi)場(chǎng)的，而$u?c$是由內(nèi)場(chǎng)傳到入口上游去的。$u,u+c,u?c$三個(gè)特征速度對(duì)應(yīng)的黎曼不變量分別為$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$這三個(gè)參數(shù)和三個(gè)黎曼不變量一一對(duì)應(yīng)。前面提到邊界的三個(gè)黎曼變量?jī)蓚€(gè)變量由邊界直接認(rèn)為給出，一個(gè)由內(nèi)場(chǎng)決定。換言之就是入口的$\rho ,u,p$中只能固定兩個(gè)量，另一個(gè)量由內(nèi)場(chǎng)的值和邊界給出的兩個(gè)值計(jì)算出來(lái)。我這里就固定$\rho ,p$兩個(gè)變量，固定$\rho =1,p=2$，速度利用黎曼不變量得出。記內(nèi)場(chǎng)第一個(gè)點(diǎn)的${\rho }_0,u_0,p_0$計(jì)算出的$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}}$，由邊界和內(nèi)場(chǎng)的$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$兩個(gè)特征速度的黎曼不變量由內(nèi)場(chǎng)推出，而$u?c$由邊界給定，換言之就是$\rho ,u,p$中兩個(gè)變量由內(nèi)場(chǎng)計(jì)算，一個(gè)變量給定，這里給定$p=1$，于是由
內(nèi)場(chǎng)最后一個(gè)點(diǎn)上$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$

超音速出口

超音速出口相比于亞音速出口就容易多了，因?yàn)槌羲贄l件下$u,u+c,u?c$都為正，因此信息只出不進(jìn)，于是直接有超音速出口邊界為${\rho }_b={\rho }_n,u_b=u_n,p_b=p_n$

本文使用的邊界

本文的入口是亞音速入口，而出口由于初場(chǎng)給定的是全場(chǎng)速度為0，密度壓力均為1，因此開始流體受壓力作用從0開始加速，開始時(shí)出口是亞音速的，但是隨著流動(dòng)不斷發(fā)展速度不斷增大，最終達(dá)到超音速，于是在出口處做一次判斷，內(nèi)場(chǎng)最后一個(gè)點(diǎn)的$u_n<c_n$時(shí)使用亞音速出口，否則使用超音速出口。

數(shù)值計(jì)算

計(jì)算代碼如下

#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,//空間點(diǎn)數(shù) NS=15000, SKIP_STEP=1500;//時(shí)間步數(shù) const double rb=-5,l=10,//計(jì)算域左邊界，計(jì)算域長(zhǎng)度 dt=0.001,//時(shí)間步長(zhǎng) 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'; }

計(jì)算結(jié)果如下

由于我沒有拉瓦爾噴管的解析解，這里僅僅驗(yàn)證了一下拉瓦爾噴管的關(guā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}$$
將上式左右兩端相減繪制曲線如下，最終結(jié)果基本滿足該關(guān)系，因此應(yīng)該問題不大。

拉瓦爾噴管工作狀態(tài)與正激波

拉瓦爾噴管根據(jù)噴管道出口流體壓力($p_e$)和背景($p_a$)壓力關(guān)系，可以分為多種工作狀態(tài):
1.$p_a=p_e$理想工作狀態(tài)，這時(shí)氣體恰好完全膨脹，管道出口為超音速流動(dòng)，并且離開管道后也不發(fā)生膨脹和壓縮，外部擾動(dòng)無(wú)法影響管內(nèi)流動(dòng)狀態(tài)。
2.$p_a>p_e$欠膨脹狀態(tài)，氣體沒有完全膨脹，管道出口為超音速流動(dòng)，氣體離開管道后會(huì)繼續(xù)膨脹一段，直至壓力達(dá)到背景壓力，這時(shí)同樣外部擾動(dòng)無(wú)法傳入內(nèi)部
3.$p_a<p_e$過膨脹狀態(tài)，氣體在出口處過膨脹，當(dāng)$p_a,p_e$差別不大時(shí)，管口仍為超音速流動(dòng)，但是由于出口處氣體壓力小于背景壓力，因此會(huì)在管口處形成一系列斜激波，最終使得流出氣體壓力和背景壓力相等；隨著$p_a,p_e$差值的增大，在管道出口處的斜激波會(huì)逐漸變?yōu)檎げ?#xff0c;繼續(xù)增加$p_e?p_a$的值，正激波會(huì)向管內(nèi)移動(dòng)使得速度由超音速變?yōu)閬喴羲?#xff0c;這時(shí)外部擾動(dòng)會(huì)影響管內(nèi)流動(dòng)。隨著出口背景壓力逐漸增大正激波位置逐漸前移，當(dāng)正激波移至喉部時(shí)達(dá)到臨界狀態(tài)，繼續(xù)增加出口背景壓力管內(nèi)流動(dòng)將全部保持在亞音速條件下。

根據(jù)前面的分析，理想狀態(tài)和欠膨脹狀態(tài)出口都是超音速流動(dòng)，因此管道內(nèi)部計(jì)算不受外部狀態(tài)影響，至于背景壓力如何管道內(nèi)流動(dòng)都是不變的。唯一有區(qū)別的是管道外是否發(fā)生膨脹，管內(nèi)流動(dòng)始終相同。

過膨脹狀態(tài)則略有不同，由于管內(nèi)可能會(huì)形成一道正激波，正激波后流動(dòng)全部是亞音速流動(dòng)，亞音速出口外部擾動(dòng)將影響內(nèi)部流動(dòng)，這是外部壓力將使得管內(nèi)流動(dòng)狀態(tài)發(fā)生改變。

在計(jì)算上的區(qū)別就是，欠膨脹和理想狀態(tài)最終都是超音速出口，而過膨脹出口總是亞音速的。

這里我試著算了一下過膨脹狀態(tài)

這是入口壓力1.1，出口壓力1時(shí)的狀態(tài)，可以看到一開始管內(nèi)出現(xiàn)了激波，說(shuō)明這時(shí)是過膨脹狀態(tài)，在喉部后側(cè)存在一道正激波，使得速度變?yōu)閬喴羲佟５俏以囁懔艘幌虏徽撛趺礈p小入口壓力只要出口壓力是1，入口壓力大于1管道內(nèi)就不可能出現(xiàn)全部亞音速的情況，似乎正激波最終會(huì)固定到$x=2$的位置處。具體原因暫時(shí)不太清楚，目前估計(jì)是入口邊界的原因，似乎固定壓力和密度的亞音速入口決定了喉部的速度。我使用真實(shí)的氣體參數(shù)作為參數(shù)進(jìn)行計(jì)算，入口設(shè)置壓力為1.001倍的標(biāo)準(zhǔn)大氣壓，密度由氣體狀態(tài)方程利用壓力和溫度計(jì)算而來(lái)，溫度設(shè)為25°C(298.15K)，出口是標(biāo)準(zhǔn)大氣壓，這是計(jì)算收斂到全場(chǎng)亞音速的結(jié)果，我懶得做無(wú)量綱了這里就不展示了。

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。