一维CFD模拟仿真.doc

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章利特，高铁瑜，夏庆锋，徐廷相.拉瓦尔喷管内的准一维定常流动.中国科技论文在线。2王平，李昌平，陈柏松.基于CFD数值模拟的拉瓦尔喷管流场分析.航空计算技术2012年7月第42卷第4期。3王如根，赵瑞贤，李全通.基于实际发动机拉瓦尔喷管的流场分析.99学术会议空军工程学院飞机推进系统实验室。4周文祥，黄金泉，周人治.拉瓦尔喷管计算模型的改进及其整机仿真验证.航空动力学报2009年11月第24卷第11期。5王克印，韩星星，张晓涛，刘耀鹏，陈吉潮.缩扩型超音速喷管的设计与仿真.中国工程机械学报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(来流马赫数M1=);P1=47892.4;%input(来流气体压强P1=);rho1=1.222;%input(来流气体密度rho1=);gama=1.4;%input(比热比gama=);R=8.314;%input(气体常量R=);C=1.5;%input(科朗数C=);T1=293;%input(来流气体温度T1=);a1=sqrt(gama*R*T1);V1=M1*a1;L=10; %input(喷管长度L=);I=300;%input(等分步数I=);N=1000;%input(时间步数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