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

1、word附录B 一维可压缩黏性流动问题的数值解法与计算程序一维可压缩黏性流动是气体动力学中最经典的黏性流动问题，对它采用迎风型差分算法进展数值求解。同时，为了初学者入门和练习方便,这里给出了由语言和语言编写的、计算一维可压缩黏性流动问题的计算程序，供大家学习参考。B-1 利用二阶迎风型差分格式求解一维可压缩黏性流动问题在两端开口管道中充满了可压缩黏性流体。当黏性流体以超声速从左向右运动时，一定会在管道中形成一道正激波，如图B.1所示。和分别为激波波前和波后的参数。该问题可简化为一维可压缩黏性流动问题。当数值解达到稳定时，在管道中可求解得到一道稳定的激波。2.根本方程组、初始条件和边界条件设流体是黏性流体。一维可压缩黏性流动问题，在数学上可以用一维可压缩黏性流动方程组来描述。量纲为一的一维方程组为： B.1其中 B.1aB.1b其中和分别是量纲为一的密度、速度、压力和单位体积总能，为流体的黏性项。为普朗特数此处公式中是有量纲量，为雷诺数，为比定容热容，为比定压热容，是量纲为一的量，称为气体绝热指数，为当地声速。求解区域为。取。初始条件：在时刻，其他物理量采用线性插值得到。边界条件：左边界

2、处： B.2右边界处：B.3差分格式一维方程组中的对流项采用的二阶精度迎风型差分格式： B.4B.4a其中。向量在第个特征方向上分量为：B.4b B.4c B.4d B.4e B.4f B.4g流通量矢量的非线性系数矩阵为： B.5非线性系数矩阵的特征值为： B.6非线性系数矩阵的右特征矢量为： B.7 B.8一维方程组中的黏性项采用二阶精度中心差分格式。采用用语言和语言对一维可压缩黏性流动问题编制了计算程序，并对雷诺数的流动进展了计算，计算结果如图B.2和图 B.3所示。图 B.3 语言程序得到的结果图 B.2 语言程序得到的结果、速度、压力和单位质量内能的分布。由上述计算结果中可以看出，采用二阶精度迎风型差分格式计算一维可压缩黏性流动问题得到的数值解和经典文献中的结果是完全一致的。计算结果明确，迎风型差分格式能够准确地捕捉激波连续，计算效果较好。由于本问题中黏性较大，所以计算得到的激波比拟光滑，有一定的宽度。一维可压缩黏性流动问题的解是连续、光滑的。B-2 一维可压缩黏性流动问题的数值计算源程序1.语言源程序/-/ 二阶迎风型差分格式求解一维可压缩黏性流动问题语言版本 /-#inc

3、lude stdio.h#include math.h#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

4、)/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

5、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

6、 ;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 ;

《附录B一维可压N-S方程(B)》由会员大米分享，可在线阅读，更多相关《附录B一维可压N-S方程(B)》请在金锄头文库上搜索。