附录B 一维可压N-S方程(B).doc

附录B 一维可压缩黏性流动问题的数值解法与计算程序一维可压缩黏性流动是气体动力学中最经典的黏性流动问题，对它采用迎风型差分算法进行数值求解。同时，为了初学者入门和练习方便,这里给出了由语言和语言编写的、计算一维可压缩黏性流动问题的计算程序，供大家学习参考。B-1 利用二阶迎风型差分格式求解一维可压缩黏性流动问题1.一维可压缩黏性流动问题在两端开口管道中充满了可压缩黏性流体。当黏性流体以超声速从左向右运动时，一定会在管道中形成一道正激波，如图B.1所示。和分别为激波波前和波后的参数。该问题可简化为一维可压缩黏性流动问题。当数值解达到稳定时，在管道中可求解得到一道稳定的激波。图B.1一维可压缩黏性流动问题示意图2.基本方程组、初始条件和边界条件设流体是黏性流体。一维可压缩黏性流动问题，在数学上可以用一维可压缩黏性流动方程组来描述。量纲为一的一维方程组为： （B.1）其中 （B.1a） （B.1b）其中和分别是量纲为一的密度、速度、压力和单位体积总能，为流体的黏性项。为普朗特数（此处公式中是有量纲量），为雷诺数，为比定容热容，为比定压热容，是量纲为一的量，称为气体绝热指数，为当地声速。求解区域为。取。初始条件：在时刻，其他物理量采用线性插值得到。边界条件：左边界处： （B.2）右边界处： （B.3）3.二阶精度迎风型差分格式一维方程组中的对流项采用的二阶精度迎风型差分格式： （B.4） （B.4a）其中。向量在第个特征方向上分量为： （B.4b） （B.4c） （B.4d） （B.4e） （B.4f） （B.4g）流通量矢量的非线性系数矩阵为： （B.5）非线性系数矩阵的特征值为： （B.6）非线性系数矩阵的右特征矢量为： （B.7） （B.8）一维方程组中的黏性项采用二阶精度中心差分格式。4.计算结果分析采用用语言和语言对一维可压缩黏性流动问题编制了计算程序，并对雷诺数的流动进行了计算，计算结果如图B.2和图 B.3所示。 图 B.3 语言程序得到的结果图 B.2 语言程序得到的结果图B.2和图 B.3是计算达到稳定后激波间断位置和密度、速度、压力和单位质量内能的分布。由上述计算结果中可以看出，采用二阶精度迎风型差分格式计算一维可压缩黏性流动问题得到的数值解和经典文献中的结果是完全一致的。计算结果表明，迎风型差分格式能够精确地捕捉激波间断，计算效果较好。由于本问题中黏性较大，所以计算得到的激波比较光滑，有一定的宽度。一维可压缩黏性流动问题的解是连续、光滑的。B-2 一维可压缩黏性流动问题的数值计算源程序1. 语言源程序/ UpwindTVD_1D.c/-/ 二阶迎风型差分格式求解一维可压缩黏性流动问题（语言版本） /-#include stdio.h#include math.h#define gama 1.4#define Tt 5.0#define im 201 /网格数/全局变量：double Q3im,Qold3im ; /Q: rou, rou*u, Edouble rouim,uim,pim,Tim,Eim,aim ;double Pr,Re,cv,cp,Ma,dx,dt ;/-void initial()double xl,xr,x;double ul,Tl,ur,Tr;/进出口的u，T值int i;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) ;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 ;return result ;/-double minmod(double x,double y)double result ;if(x*yfabs(y)result=y ;else if(fabs(x)fabs(y)result=x ;return result ;void Q2U()int i ;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 ;/-void UpwindTVD_1D_Solver()int i,l,k ;double dq3im,lamda3im,sigma3im,f3im ;double fl3,fwave3im,fr3,g3im,Gv3im ;double alfa3im,beta3im,Rn33im,R33im ;double ut,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 ;/计算左右特征向量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/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 ;fwaveli=0.0 ;gli=0.0 ;/计算alfafor(i=0; i=im-2; i+)for(l=0; l=2; l+)for(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)betali=(gli+1-gli)/alfali ;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 ) ;for(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 ;Qli=Qli-dt/dx*(frl-fll) ;/黏性项的处理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*(Ti+1-2.*Ti+Ti-1) +(Ti+1-Ti-1)*(Ti+1-Ti-1)/4.)*cp/Pr)/Re/(dx*dx) ;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) ;/-double error(double Q13im,double Q23im)int i,l ;double err,ddu ;err=1.0e-10 ;for(i=0; iim; i+)for(l=0; l3 ; l+)if(Q2li!=0.0)ddu=fabs(Q2li-Q1li)/Q2li) ;elseddu=fabs(Q2li-Q1li) ;if(errddu) err=ddu ;return err ;/-void output()double x,rou1,u1,E1,p1,T1 ; /临时变量int i ;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*Ma*Ma) ; fprintf(fp,%15.8e %15.8e %15.8e %15.8e %15.8en, x,rou1,u1,p1,cv*T1 ) ;fclose(fp) ;/-main()double t,err ;int n,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(%10d time: %16.9e error: %16.9en,n,t,err) ;output() ;output() ;printf(Program end !n) ;-2. 语言源程序! UpwindTVD_1D.f!-! 二阶迎风型差分格式求解一维可压缩黏性流动问题!（语言版本） !-program UPWIND_TVD_1Dimplicit real*8 (a-h,o-z) parameter(mx=201,Tt=5.0 )dimension Q(3,mx),Qold(3,mx) ! Q: rou, rou*u, E dimension rou(mx),u(mx),p(mx),T(mx),E(mx)real*8 Macommon/para_def/Sf,Re,Ma,pr,gama,cp,cv,dx,dt,im im=mxcall Initialize(Q)time=0.0 n=01 n=n+1 time=time+dtdo i=1,imdo l=1,3 Qold(l,i)=Q(l,i)end doend docall UpwindTVD_1D_Solver(Q)if(n/10000*10000.eq.n)then write(*,20)n,time,error(Q,Qold) call OutPut(Q)end if20 format(1x,i10, time=,e16.9, error=,e16.9,1x)if(time.lt.Tt) .and. (error(Q,Qold).gt.1.0e-8)then goto 1end ifcall OutPut(Q)write(*,*) Program end !end!- subroutine Initialize(Q)implicit real*8 (a-h,o-z)real*8 Macommon/para_def/Sf,Re,Ma,pr,gama,cp,cv,dx,dt,im dimension rou(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=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) do i=1,im rou(i)=rou0 xi=(i-1)*dx u(i)=(xi-xl)/(xr-xl)*(u1-u0)+u0 T(i)=(xi-xl)/(xr-xl)*(T1-T0)+T0 P(i)=rou(i)*T(i)/(gama*Ma*2) E(i)=rou(i)*(cv*T(i)+0.5*u(i)*2)end dodo i=1,im Q(1,i)=rou(i) Q(2,i)=rou(i)*u(i) Q(3,i)=E(i)end doend!-subroutine Q2U(Q,rou,u,p,T,E,a)implicit real*8 (a-h,o-z)real*8 Macommon/para_def/Sf,Re,Ma,pr,gama,cp,cv,dx,dt,imdimension Q(3,im),rou(im),u(im)dimension p(im),T(im),E(im),a(im)do i=1,im rou(i)=Q(1,i) u(i)=Q(2,i)/Q(1,i) E(i)=Q(3,i) T(i)=(E(i)/rou(i)-0.5*u(i)*2)/cv p(i)=rou(i)*T(i)/(gama*Ma*2) a(i)=dsqrt(T(i)/Maend doend!-subroutine UpwindTVD_1D_Solver(Q)implicit real*8 (a-h,o-z)real*8 Ma,lamda,minmodcommon/para_def/Sf,Re,Ma,pr,gama,cp,cv,dx,dt,imdimension Q(3,im),Qold(3,im),dq(3,im),rou(im),u(im),p(im),T(im)dimension E(im),a(im),f(3,im),fwave(3,im),fr(3),fl(3)dimension Gv(3,im),g(3,im),lamda(3,im),sigma(3,im)dimension alfa(3,im),Rn(3,3,im),R(3,3,im),beta(3,im)call Q2U(Q,rou,u,p,T,E,a)do i=1,im f(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)end dodo i=1,im-1 lamda(1,i)=(u(i)+u(i+1)/2.0 lamda(2,i)=(u(i)-a(i)+u(i+1)-a(i+1)/2.0 lamda(3,i)=(u(i)+a(i)+u(i+1)+a(i+1)/2.0end dodo i=1,im-1do l=1,3 dq(l,i)=(Q(l,i+1)-Q(l,i) sigma(l,i)=qz(lamda(l,i)/2.0 ! 定常情况end doend do do i=1,im-1 ut=(u(i)+u(i+1)/2.0 !ut 临时变量 at=(a(i)+a(i+1)/2.0 !at 临时变量 R(1,1,i)=1.0 R(1,2,i)=1.0 R(1,3,i)=1.0 R(2,1,i)=ut R(2,2,i)=ut-at R(2,3,i)=ut+at R(3,1,i)=ut*ut/2.0 R(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/at Rn(1,3,i)=-(gama-1.0)/at/at Rn(2,1,i)=ut/2.0/at+(gama-1.0)*ut*ut/4.0/at/at Rn(2,2,i)=-(gama-1.0)*ut/2./at/at-1.0/2.0/at Rn(2,3,i)=(gama-1.0)/2.0/at/at Rn(3,1,i)=-ut/2.0/at+(gama-1.0)*ut*ut/4.0/at/at Rn(3,2,i)=(1.0-gama)*ut/2./at/at+1.0/2.0/at Rn(3,3,i)=(gama-1.0)/2./at/atend do do i=1,imdo l=1,3 alfa(l,i)=0.0 fwave(l,i)=0.0 g(l,i)=0.0end doend dodo i=1,im-1do l=1,3do k=1,3 alfa(l,i)=alfa(l,i)+Rn(l,k,i)*dq(k,i)end doend doend do do i=2,im-1do l=1,3 g(l,i)=minmod(sigma(l,i)*alfa(l,i),sigma(l,i-1)*alfa(l,i-1)end doend dodo i=1,im-1do l=1,3 if(alfa(l,i).ne.0.0)then beta(l,i)=(g(l,i+1)-g(l,i)/alfa(l,i) else beta(l,i)=0.0 end if end doend dodo i=1,im-1 do l=1,3do k=1,3 fwave(l,i)=fwave(l,i)+0.5*R(l,k,i)*(g(k,i)+g(k,i+1) * -qz(lamda(k,i)+beta(k,i)*alfa(k,i) end doend doend do do i=2,im-1do l=1,3 fl(l)= 0.5*(f(l,i-1)+f(l,i) +fwave(l,i-1) fr(l)= 0.5*(f(l,i+1)+f(l,i) +fwave(l,i) Q(l,i)=Q(l,i)-dt/dx*(fr(l)-fl(l)end doend do! 处理黏性项 do i=2,im-1 Gv(1,i)=0.0 Gv(2,i)=(T(i)*(u(i+1)-2.*u(i)+u(i-1)+(T(i+1)-T(i-1)* &(u(i+1)-u(i-1)/4.)*4./3./re/dx*2 Gv(3,i)=(T(i)*(u(i+1)-2.*u(i)+u(i-1)*u(i)+(T(i+1)-T(i-1)* &(u(i+1)-u(i-1)*u(i)/4.+T(i)*(u(i+1)-u(i-1)*2/4.)*4./3. &+(T(i)*(T(i+1)-2.*T(i)+T(i-1)+(T(i+1)-T(i-1)*2/4.) &*cp/pr)/re/dx*2 end dodo i=2,im-1do l=1,3 Q(l,i)=Q(l,i)+Gv(l,i)*dtend doend do! 边界条件 Q(1,im)=Q(1,im-1