b一维可压ns方程b0814093738.docx

1、b一维可压ns方程b0814093738附录B一维可压缩黏性流动问题的数值解法与计算程序一维可压缩黏性流动是气体动力学中最经典的黏性流动问题，对 它采用迎风型TVD差分算法进行数值求解。同时，为了初学者入门和 练习方便,这里给出了由C语言和Fortran7语言编写的、计算一维可压 缩黏性流动问题的计算程序，供大家学习参考。B-1利用二阶迎风型TVD差分格式求解一维可压缩黏性流动问题 1. 一维可压缩黏性流动问题在两端开口管道中充满了可压缩黏性流体。 当黏性流体以超声速从左向右运动时，一定会在管道中形成一道正激波，如图 B.1所示。1、Ui、Pi和2、U2、P2分别为激波波前和波后的参数。该问题

2、可简化为一 维可压缩黏性流动问题。当数值解达到稳定时，在管道中可求解得到 一道稳定的激波。r1| - k1p2 111 U-1 *2 11II1 耳-13P*-P1 *111111X 0 I2.基本方程组、初始条件和边界条件设流体是黏性流体。一维可压缩黏性流动问题，在数学上可以用TEC T1 2I IPrCpp 2 ,M2ECv 1u2,PrkU L1CpTRe ,cv9 j,a1M2CvM为:f ( B.1)x x其中(B.1b)其中、up和E分别是量纲为一的密度、速度、压力和单位体积总能， s为流体的黏性项。Pr为普朗特数（此处公式中Cp、k、 是有量纲 量），Re为雷诺数，c为比定容热容

3、，Cp为比定压热容，Cp, cv是量纲 为一的量，称为气体绝热指数，a为当地声速。初始条件：在t 0时刻，x,0 1，其他物理量采用线性插值得到。右边界x 1处：x,t0x 1xu 1,t2M211 . . 2 M1T 1,t 2 M211 121 11 M23.二阶精度迎风型TVD差分格式一维N-S方程组中的对流项采用Harten的二阶精度迎风型分格式：n 1 n . nuj uj r h 1J J j 2hnJs12 Xt (B.4)1ji 2 fj fj1 Rj2j2 (B.4a)其中r 。向量在第l个特征方向上分量 为:xi i i.1 gj gj 1 Qj 2la 1j 2l.1 .

4、 1j T j T2 2(B.4b)gj mi nmod 1 1J 2l1 , j 2l l1 1j - j -2 2(B.4c)(B.3)TVD差l lgj 1 gjl.1 j20(B.4d)0,流通量矢量f的非线性Jiacobia系数矩阵A= R ar-1为:非线性Jiacobia系数矩阵A的特征值为:(B.6)非线性Jiacobia系数矩阵A的右特征矢量R、R-1为:一维N-S方程组中的黏性项111R uu a ua(B.7)1 2 u2h ua hua11 u21u12a22 a2 aR1u1 u21 u112a2( B.8)2a4a22a22au1 u21 u112a4a22a22a

5、2a2t采用二阶精度中心差分格式4.计算结果分析采用用C语言和Fortran77语言对一维可压缩黏性流动问题编制了计算程序, 并对雷诺数Re 50的流动进行了计算，计算结果如图 B.2和图B.3所示。图B.2 Fortran7语言程序得到 图B.3 C语言程序得到图B.2和图B.3是计算达到稳定后激波间断位置和密度（）、速度（u）、压的结果 的结果力（P）和单位质量内能（e）的分布。由上述计算结果中可以看出，采用二阶精度迎 风型TVD差分格式计算一维可压缩黏性流动问题得到的数值解和经典文献中的 结果是完全一致的。计算结果表明，迎风型TVD差分格式能够精确地捕捉激波间 断，计算效果较好。由于本问

6、题中黏性较大（Re 50），所以计算得到的激波比较光滑，有一定的宽度。一维可压缩黏性流动问题的解是连续、光滑的B-2 一维可压缩黏性流动问题的数值计算源程序1. C语言源程序/Upwi ndTVD_1D.c/ /二阶迎风型TVD差分格式求解一维可压缩黏性流动问题（C语言版本）/#includestdio.h#includemath.h#definegama1.4#defineTt5.0#defineim201/ 网格数/ 全局变量： doubleQ3im,Qold3im;/Q:rou,rou*u,E doublerouim,uim,pim,Tim,Eim,aim; doublePr,Re,cv

7、,cp,Ma,dx,dt;/ voidinitial()doublexl,xr,x;doubleul,Tl,ur,Tr;/ 进出口的 u，T 值 inti;dx=1.0/(im-1);dt=1.0e-6;Pr=0.72;Re=50.0;Ma=2.0;cv=1.0/(gama*(gama-1.0)*Ma*Ma);cp=gama*cv;xl=0.0;xr=1.0;ul=1.0;Tl=1.0;ur=(2.0/(gama-1.0)+Ma*Ma)/(gama+1.0)/(gama-1.0)*Ma*Ma);Tr=2.0*gama/(gama+1.0)*Ma*Ma-(gama-1.0)/(gama+1.0)

8、;Tr=Tr*(gama-1.0)/(gama+1.0)+2.0/(gama+1.0)/Ma/Ma);for(i=0;i=im-1;i+)x=i*dx;roui=1.0;ui=(x-xl)/(xr-xl)*(ur-ul)+ul;Ti=(x-xl)/(xr-xl)*(Tr-Tl)+Tl;pi=roui*Ti/(gama*Ma*Ma);Ei=roui*(cv*Ti+0.5*ui*ui);for(i=0;i0.1)result=fabs(x);elseresult=(x*x/0.1+0.1)/2.0;returnresult;/ doubleminmod(doublex,doubley)double

9、result;if(x*yfabs(y)result=y;elseif(fabs(x)fabs(y)result=x;returnresult;voidQ2U()inti;for(i=0;iim;i+)roui=Q0i;ui=Q1i/roui;Ei=Q2i;Ti=(Ei/roui-0.5*ui*ui)/cv;pi=roui*Ti/(gama*Ma*Ma); ai=sqrt(Ti)/Ma;/ voidUpwindTVD_1D_Solver()inti,l,k;doubledq3im,lamda3im,sigma3im,f3im; doublefl3,fwave3im,fr3,g3im,Gv3im

10、; doublealfa3im,beta3im,Rn33im,R33i m;doubleut,at;/ 临时变量Q2U();for(i=0;iim;i+)f0i=roui*ui;f1i=roui*ui*ui+pi;f2i=ui*(Ei+pi);for(i=0;i=im-2;i+)lamda0i=(ui+ui+1)/2.0;lamda1i=(ui-ai+ui+1-ai+1)/2.0;lamda2i=(ui+ai+ui+1+ai+1)/2.0;for(i=0;i=im-2;i+)for(l=0;l=2;l+)dqli=Qli+1-Qli;sigmali=qz(lamdali)/2.0;/ 计算左右

11、特征向量for(i=0;i=im-2;i+)ut=(ui+ui+1)/2.0;at=(ai+ai+1)/2.0;R00i=1.0;R01i=1.0;R02i=1.0;R10i=ut;R11i=ut-at;R12i=ut+at;R20i=ut*ut/2.0;R21i=at*at/(gama-1.0)+ut*ut/2.0-ut*at;R22i=at*at/(gama-1.0)+ut*ut/2.0+ut*at;Rn00i=1.0-(gama-1.0)*ut*ut/2.0/at/at;Rn01i=(gama-1.0)*ut/at/at;Rn02i=-(gama-1.0)/at/at;Rn10i=ut/

12、2.0/at+(gama-1.0)*ut*ut/4.0/at/at;Rn11i=-(gama-1.0)*ut/2./at/at-1.0/2.0/at;Rn12i=(gama-1.0)/2.0/at/at;Rn20i=-ut/2.0/at+(gama-1.0)*ut*ut/4.0/at/at;Rn21i=(1.0-gama)*ut/2./at/at+1.0/2.0/at;Rn22i=(gama-1.0)/2./at/at;for(i=0;i=im-1;i+)for(l=0;l=2;l+)alfali=0.0;/ 计算 alfafor(i=0;i=im-2;i+)for(l=0;l=2;l+)fo

13、r(k=0;k=2;k+) alfali=alfali+Rnlki*dqki;/ 计算 gfor(i=1;i=im-2;i+)for(l=0;l=2;l+)gli=minmod(sigmali*alfali,sigmali-1*alfali-1);for(i=0;i=im-2;i+)for(l=0;l=2;l+)if(alfali!=0.0)elsebetali=0.0;for(i=0;i=im-2;i+)for(l=0;l=2;l+)for(k=0;k=2;k+)fwaveli=fwaveli+0.5*Rlki*(gki+gki+1-qz(lamdaki+betaki)*alfaki);fo

14、r(i=1;i=im-2;i+)for(l=0;l=2;l+) fll=0.5*(fli-1+fli)+fwaveli-1; frl=0.5*(fli+1+fli)+fwaveli;/ 黏性项的处理for(i=1;i=im-2;i+)Gv0i=0.0;Gv1i=(Ti*(ui+1-2.*ui+ui-1)+(Ti+1-Ti-1 )*(ui+1-ui-1)/4.)*4./3./Re/(dx*dx);Gv2i=(Ti*(ui+1-2.*ui+ui-1)*ui+(Ti+1 -Ti-1)*(ui+1-ui-1)*ui/4.+Ti*(ui+1-ui-1)*(ui+1-ui-1)/4.)*4./3.+(Ti

15、*(Ti+1-2.*Ti+Ti -1)+(Ti+1-Ti-1)*(Ti+1-Ti-1)/4.)*cp/Pr)/Re/(dx*d x); for(i=1;i=im-2;i+) for(l=0;l=2;l+)Qli=Qli+Gvli*dt;/ 边界条件Q0im-1=Q0im-2;Q1im-1=Q0im-1*uim-1;Q2im-1=Q0im-1*(cv*Tim-1+0.5*uim-1*uim-1 );/ doubleerror(doubleQ13im,doubleQ23im)inti,l;doubleerr,ddu;err=1.0e-10;for(i=0;iim;i+)for(l=0;l3;l+)

16、if(Q2li!=0.0)ddu=fabs(Q2li-Q1li)/Q2li);elseddu=fabs(Q2li-Q1li);if(errddu)err=ddu;returnerr;/ voidoutput()doublex,rou1,u1,E1,p1,T1;/ 临时变量inti;FILE*fp;fp=fopen(result.dat,w);fprintf(fp,variables=x,rou,u,p,en);for(i=0;iim;i+)x=i*dx;rou1=Q0i;u1=Q1i/rou1;E1=Q2i;T1=(E1/rou1-0.5*u1*u1)/cv;p1=rou1*T1/(gama*

17、Ma*Ma);fprintf(fp,%15.8e%15.8e%15.8e%15.8e%15.8en,x,rou1,u1,p1,cv*T1);fclose(fp);/ main()doublet,err;intn,i,l;initial();t=0.0;n=0;err=1.0e-7;while(t1.0e-8)t=t+dt;n=n+1;for(i=0;iim;i+)for(l=0;l3;l+)Qoldli=Qli;UpwindTVD_1D_Solver();err=error(Q,Qold);if(n/10000*10000=n)printf(%10dtime:%16.9eerror:%16.9

18、en,n,t,err);output();output();printf(Programend!n);2. Fortran77语言源程序!UpwindTVD_1D.f! !二阶迎风型TVD差分格式求解一维可压缩黏性流动问题 ! （ Fortran77语言版本）! programUPWIND_TVD_1Dimplicitreal*8(a-h,o-z) parameter(mx=201,Tt=5.0)dimensionQ(3,mx),Qold(3,mx)!Q:rou,rou*u,E dimensionrou(mx),u(mx),p(mx),T(mx),E(mx)real*8Macommon/par

19、a_def/Sf,Re,Ma,pr,gama,cp,cv,dx,dt,im im=mxcallInitialize(Q)time=0.0n=01n=n+1time=time+dtdoi=1,imdol=1,3Qold(l,i)=Q(l,i)enddoenddocallUpwindTVD_1D_Solver(Q) write(*,20)n,time,error(Q,Qold)callOutPut(Q)endif20format(1x,i10,time=,e16.9,error=,e16.9,1x) .and.goto1endifcallOutPut(Q)write(*,*)Programend!

20、end! subroutineInitialize(Q)implicitreal*8(a-h,o-z)real*8Macommon/para_def/Sf,Re,Ma,pr,gama,cp,cv,dx,dt,im dimensionrou(im),u(im),p(im),T(im),E(im),Q(3,im)Re=50.0! 雷诺数Ma=2.0! 马赫数pr=0.72! 普朗特数gama=1.4! 气体常数cv=1.0/(gama*(gama-1.0)*Ma*2)! 比定容热容cp=gama*cv! 比定压热容dt=1.0e-6dx=1.0/(1.0*(im-1)! 空间步长xl=0.0xr=

21、1.0rou0=1.0u0=1.0T0=1.0temp=(gama+1.)/(gama-1.)u1=(2./(gama-1.)+Ma*2)/(temp*Ma*2)T1=(2.*gama/(gama+1)*Ma*2-1./temp)*(1./temp+2./(gama+1.)/Ma*2)doi=1,imrou(i)=rou0xi=(i-1)*dxu(i)=(xi-xl)/(xr-xl)*(u1-u0)+u0T(i)=(xi-xl)/(xr-xl)*(T1-T0)+T0P(i)=rou(i)*T(i)/(gama*Ma*2)E(i)=rou(i)*(cv*T(i)+0.5*u(i)*2)enddo

22、doi=1,imQ(1,i)=rou(i)Q(2,i)=rou(i)*u(i)Q(3,i)=E(i)enddoend! subroutineQ2U(Q,rou,u,p,T,E,a) implicitreal*8(a-h,o-z) real*8Ma common/para_def/Sf,Re,Ma,pr,gama,cp,cv,dx,dt,im dimensionQ(3,im),rou(im),u(im) dimensionp(im),T(im),E(im),a(im) doi=1,im rou(i)=Q(1,i) u(i)=Q(2,i)/Q(1,i)E(i)=Q(3,i) T(i)=(E(i)/

23、rou(i)-0.5*u(i)*2)/cv p(i)=rou(i)*T(i)/(gama*Ma*2) a(i)=dsqrt(T(i)/Ma enddoend! subroutineUpwindTVD_1D_Solver(Q)implicitreal*8(a-h,o-z)real*8Ma,lamda,minmod common/para_def/Sf,Re,Ma,pr,gama,cp,cv,dx,dt,im dimensionQ(3,im),Qold(3,im),dq(3,im),rou(im),u(im),p (im),T(im)dimensionE(im),a(im),f(3,im),fwa

24、ve(3,im),fr(3),fl(3) dimensionGv(3,im),g(3,im),lamda(3,im),sigma(3,im) dimensionalfa(3,im),Rn(3,3,im),R(3,3,im),beta(3,im) callQ2U(Q,rou,u,p,T,E,a)doi=1,imf(1,i)=rou(i)*u(i)f(2,i)=rou(i)*u(i)*2+p(i)f(3,i)=u(i)*(E(i)+p(i)enddodoi=1,im-1lamda(1,i)=(u(i)+u(i+1)/2.0lamda(2,i)=(u(i)-a(i)+u(i+1)-a(i+1)/2.

25、0lamda(3,i)=(u(i)+a(i)+u(i+1)+a(i+1)/2.0enddodoi=1,im-1dol=1,3dq(l,i)=(Q(l,i+1)-Q(l,i) sigma(l,i)=qz(lamda(l,i)/2.0! 定常情况 enddoenddodoi=1,im-1ut=(u(i)+u(i+1)/2.0!ut 临时变量at=(a(i)+a(i+1)/2.0!at 临时变量R(1,1,i)=1.0R(1,2,i)=1.0R(1,3,i)=1.0R(2,1,i)=utR(2,2,i)=ut-atR(2,3,i)=ut+atR(3,1,i)=ut*ut/2.0R(3,2,i)=at*at/(gama-1.0)+ut*ut/2.0-ut*at R(3,3,i)=at*at/(gama-1.0)+ut*ut/2.0+ut*at Rn(1,1,i)=1.0-(gama-1.0)*ut*ut/2.0/at/at Rn(1,2,i)=(gama-1.0)*ut/at/atRn(1,3,i)=-(gama-1.0)/at/atRn(2,1,i)=ut/2.0/at+(gama-1.0)*ut*ut/4.0/at/atRn(2,2,i)=-(gama-1.0)*ut/2./at/at-1.0/2.0/atRn(2,3,i)=(gam