CFD simulation in Laval nozzleSIAE 090441313AbstractWe aim to simulate the quasi one dimension flow in the Laval nozzle based on CFD computation in this paper .We consider the change of the temperature ,the pressure ,the density and the speed of the flow to study the flow.The analytic solution of the flow in the Laval nozzle is provided when the input velocity is supersonic.We use the Mac-Cormack Explicit Difference Scheme to slove the question.Key words :Laval nozzle ,CFD,throat narrow.Contents Abstract .1Introduction .2Simulation of one-dimensional steady flow.3 Basis equations .3Dimensionless .10Mac -Cormack Explicit Difference Scheme.11Boundary conditions .13Reference .13Annex .14IntroductionLaval nozzle is the most commonly used components of rocket engines and aero-engine, constituted by two tapered tube, one shrink tube, another expansion tube.Laval nozzle is an important part of the thrust chamber. The first half of the nozzle from large to small contraction to a narrow throat to the middle. Narrow throat and then expand outwards from small to big to the end. The gas in the rocket body by the front half of the high pressure into the nozzle, through the narrow throat to escape by the rear half. This architecture allows the speed of the air flow changes due to changes in the jet cross-sectional area, the airflow from subsonic to the speed of sound, until accelerated to transonic. So, people flared nozzle called transonic nozzle. Since it was invented by the Swedish Laval, also known as Laval nozzle. Analysis of the principle of the Laval nozzle. The rocket engines of the gas flow in the combustion chamber under pressure, after the backward movement of the nozzle into the nozzle . At this stage, the gas movement follow the principle of the fluid moves in the tube , the small cross-section at the flow rate large sectional large flow velocity, thus accelerating airflow.Laval nozzle When you reach the narrow throat, the flow rate has exceeded the speed of sound. Transonic fluid movement they no longer follow the principle of cross-section at small velocity, at a flow rate of small cross-section large, but on the contrary the larger cross-sectional flow faster. The gas flow speed is further accelerated to 2-3 km / sec,equivalent to 7-8 times the speed of sound, thus creating a great thrust. The Laval nozzle fact played the role of a flow rate Enlargement Device. In fact, not just rocket engines, missile nozzle is this horn shape, so the Laval nozzle weapons has a very wide range of applications.Simulation of one-dimensional steady flow1.Basis equationsAs we know,Laval nozzle is a zooming nozzle flow channel to narrow further expansion.Allows the airflow to further accelerate to reach the speed of sound at the throat into a supersonic flow.Now,we want to simulate the quasi one-dimension flowing.Firstly,we will analysis on theory.The flow is isentropic,so we can apply the following equations.(1)Continuity equation:In the flow, we need to consider the following physical quantities.The pression ,the temperature ,the speed of the fluid and the cross-section .They are respectively represented by P,T,u,A. We apply the conservation of the mass.we will obtain this equation.And then we get(2)Equation of momentum(in the direction of the axis)According to the theory of momentum:The simplification of this equation is (3)Energy equationIdeal gas equation of stateR is ideal gas constant,R=8.314J/g/K.M is the masse per mole.(4)The equation of ThermodynamicsBecause the flow is isentropic,so dS=0And we use the equation of momentum,we haveCombine with others equations,we result withWe called u the speed of sound,we noted a. We apply the continuity equationWe defined the Mach number If we have the relation as We have the figure So M1,supersonicIf dA0.If dA0,we have du0.M1,subsonicIf dA0.If dA0,we have du0.This is the reason why this architecture allows the speed of the air flow changes due to changes in the jet cross-sectional area, the airflow from subsonic to the speed of sound, until accelerated to transonic.We have the consequence as followsThen we replace P and T in this equation.The consequence will becomeTo simplifyIn this equation,the variable is the much number,as the speed of the flow is from subsonic to supersonic ,so we can suppose that there exist a critical section where M equal to 1.ThenFigure This section is called narrow throat.The same method,we can obtainFigure We know the section in narrow throat is minimum.we can judge that the function attains the maximum or not 2 Dimensionless Combining CFD with one-dimension flow theory,we make the variables dimensionless.According to the condition initial which is given .We note We use the variable dimensionless to represent the equations.And the equations have the following changes (1)Continuity equation(2)Equation of momentum(3)Energy equation3.Mac-Cormack Explicit Difference SchemeThen we use the Mac-Cormack Explicit Difference Scheme,the principal of this theory is using the surrounding points to present differential parts of a point and we consider the question with one dimension.The distance between two points is h.So we can use two points adjacent to present the differential parts.Using this method,we make an estimation and correct the error.Estimationcorrect the error Intermediate value Then the equation has the following changeAt the moment t ,we will know the value in the whole plan .And we define 4. Boundary conditions Hyperbola equation has two characteristics lines.When one of the characteristics lines enter the flow zone .We admit a parameter to be fixed ,otherwise when one of the characteristics go out the flow zone ,we admit a parameter to be a variable depends the time.Applying this theory ,we can determine the boundary conditions.Reference ：1章利特，高鐵瑜，夏慶鋒，徐廷相.拉瓦爾噴管內(nèi)的準(zhǔn)一維定常流動.中國科技論文在線。2王平，李昌平，陳柏松.基于CFD數(shù)值模擬的拉瓦爾噴管流場分析.航空計算技術(shù)2012年7月第42卷第4期。3王如根，趙瑞賢，李全通.基于實際發(fā)動機(jī)拉瓦爾噴管的流場分析.99學(xué)術(shù)會議空軍工程學(xué)院飛機(jī)推進(jìn)系統(tǒng)實驗室。4周文祥，黃金泉，周人治.拉瓦爾噴管計算模型的改進(jìn)及其整機(jī)仿真驗證.航空動力學(xué)報2009年11月第24卷第11期。5王克印，韓星星，張曉濤，劉耀鵬，陳吉潮.縮擴(kuò)型超音速噴管的設(shè)計與仿真.中國工程機(jī)械學(xué)報2011年9月第9卷第3期。Annex 1:Figure 1程序：x=0:0.1:5;a=1.398+0.374*tanh(0.8*x-4);plot(x,a)Figure 2程序：gama=1.33;M=0:0.01:2;A=(1./M).*(1+(gama-1).*M.2./2)./(gama+1)./2).(gama+1)./(2.*(gama-1);Xlabel(variable M);ylabel(variable A);Plot(M,A)Figure 3程序：gama=1.33;M=0:0.01:2;T=(gama+1)./2)./(1+(gama-1.).*M.2./2);P=(gama+1)./2)./(1+(gama-1).*M.2./2).(1./(gama-1);rho=(P./T).*(gama-1)./2)./(1+(gama-1).*M.2./2).(gama./(gama-1);x=T;P;rho;y=M;M;plot(M,T,M,P,M,rho);subplot(221);plot(M,T);xlabel(variable M);ylabel(variable T)；subplot(222);plot(M,P);xlabel(variable M);ylabel(variable P)；subplot(2,2,3:4);plot(M,rho);xlabel(variable M);ylabel(variable rho)；模擬程序（未完成）：M1=1.5;%input(來流馬赫數(shù)M1=);P1=47892.4;%input(來流氣體壓強(qiáng)P1=);rho1=1.222;%input(來流氣體密度rho1=);gama=1.4;%input(比熱比gama=);R=8.314;%input(氣體常量R=);C=1.5;%input(科朗數(shù)C=);T1=293;%input(來流氣體溫度T1=);a1=sqrt(gama*R*T1);V1=M1*a1;L=10; %input(噴管長度L=);I=300;%input(等分步數(shù)I=);N=1000;%input(時間步數(shù)N=);d_t=0;e=0;delta_x=1/(I-1);A1=1.398-0.347*tanh(4);A=zeros(I,1);V=zeros(I,N);rho=zeros(I,N);T=zeros(I,N);ahead_V=zeros(I,N);ahead_rho=zeros(I,N);ahead_T=zeros(I,N);ahead_der_V=zeros(I,N);ahead_der_rho=zeros(I,N);ahead_der_T=zeros(I,N);der_V=zeros(I,N);der_rho=zeros(I,N);der_T=zeros(I,N);ave_der_V=zeros(I,N);ave_der_rho=zeros(I,N);ave_der_T=zeros(I,N);delta_t=zeros(I,N);a=zeros(I,N);e=zeros(I,N); for n=1:N ahead_V(1,n)=M1; ahead_rho(1,n)=1; ahead_T(1,n)=1; V(1,n)=M1; rho(1,n)=1; T(1,n)=1;endfor i=1:I V(i,1)=(1.5+1.09*i*delta_x)*sqrt(T(i,1); rho(i,1)=1-0.3146*i*delta_x; T(i,1)=1-0.2314*i*delta_x; A(i,1)=(1.398+0.347*tanh(0.8*(i-1)*delta_x*L-4);endfor n=1:N for i=2:(I-1) der_V(i,n)=-V(i,n)*(V(i+1,n)-V(i,n)/delta_x-1/gama*(T(i+1,n)-T(i,n)/delta_x+T(i,n)/rho(i,n)*(rho(i+1,n)-rho(i,n)/delta_x); der_rho(i,n)=-rho(i,n)*(V(i+1,n)-V(i,n)/delta_x-V(i,n)*(rho(i+1,n)-rho(i,n)/delta_x-rho(i,n)*V(i,n)*(log(A(i+1,1)-log(A(i,1)/delta_x; der_T(i,n)=-V(i,n)*(T(i+1,n)-T(i,n)/delta_x-(gama-1)*T(i,n)*(V(i+1,n)-V(i,n)/delta_x+V(i,n)*(log(A(i+1,1)-log(A(i,1)/delta_x); a(i,n)=sqrt(gama*R*T(i,n)/a1; delta_t(i,n)=C*delta_x/(V(i,n)+a(i,n); end d_t=min(delta_t(i,n); for i=2:1:(I-1) ahead_V(i,n+1)=V(i,n)+der_V(i,n)*d_t; ahead_rho(i,n+1)=rho(i,n)+der_rho(i,n)*d_t; ahead_T(i,n+1)=T(i,n)+der_T(i,n)*d_t; ahead_der_V(i,n)=-ahead_V(i,n+1)*(ahead_V(i,n+1)-ahead_V(i-1,n+1)/delta_x-1/gama*(ahead_T(i,n+1)-a