附录B一维可压NS方程B

附录B一维可压缩黏性流动问题的数值解法与计算程序一维可压缩黏性流动是气体动力学中最经典的黏性流动问题，对它采用迎风型TVD差分算法进行数值求解。同时，为了初学者入门和练习方便,这里给出了由C语言和Fortran77语言编写的、计算一维可压缩黏性流动问题的计算程序，供大家学习参考。B-1利用二阶迎风型TVD差分格式求解一维可压缩黏性流动问题1.一维可压缩黏性流动问题在两端开口管道中充满了可压缩黏性流体。当黏性流体以超声速从左向右运动时，一定会在管道中形成一道正激波，如图B.1所示。111up、、和222up、、分别为激波波前和波后的参数。该问题可简化为一维可压缩黏性流动问题。当数值解达到稳定时，在管道中可求解得到一道稳定的激波。2.基本方程组、初始条件和边界条件设流体是黏性流体。一维可压缩黏性流动问题，在数学上可以用一维可压缩黏性流动N-S方程组来描述。量纲为一的一维N-S方程组为：txxufs（B.1）其中122304,,343psuuuupsRexEEpusCKuuTRexPrRrxufs（B.1a）图B.1一维可压缩黏性流动问题示意图2221,,21,,,1pvpvvcTpECTuPrMkcULTRecaMcM（B.1b）其中up、、和E分别是量纲为一的密度、速度、压力和单位体积总能，s为流体的黏性项。Pr为普朗特数（此处公式中pck、、是有量纲量），Re为雷诺数，vc为比定容热容，pc为比定压热容，,pvcc是量纲为一的量，称为气体绝热指数，a为当地声速。求解区域为10x。取2,50,0.72,1.4,MRePrkT。初始条件：在0t时刻，,01x，其他物理量采用线性插值得到。边界条件：左边界0x处：0,0,0,1tutTt（B.2）右边界1x处：12222,0211,1121121,1111xxtxMutMTtMM（B.3）3.二阶精度迎风型TVD差分格式一维N-S方程组中的对流项采用Harten的二阶精度迎风型TVD差分格式：11122nnnnjjjjrtxsuuhh（B.4）111122212jjjjjhffRΦ（B.4a）其中trx。向量Φ在第l个特征方向上分量l为：111112222llllljjjjjjggQa（B.4b）11112222,llllljjjjjgminmod（B.4c）111212212,00,0lljjlljljjljgg（B.4d）1111222llljjjRu（B.4e）2,,0.050.51,2zzQzzz（B.4f）2112212121212lliililiQrQ非定常定常,,（B.4g）流通量矢量f的非线性Jiacobian系数矩阵-1=ARΛR为：222320103312232112uuuaauuufAuu（B.5）非线性Jiacobian系数矩阵A的特征值为：123123000000,,uuauaΛ（B.6）非线性Jiacobian系数矩阵A的右特征矢量-1、RR为：211112uuauauhuahuaR（B.7）222221222222211112111124222111124222uuaaauuuaaaaauuuaaaaaR（B.8）一维N-S方程组中的黏性项txs采用二阶精度中心差分格式。4.计算结果分析采用用C语言和Fortran77语言对一维可压缩黏性流动问题编制了计算程序，并对雷诺数50Re的流动进行了计算，计算结果如图B.2和图B.3所示。图B.2和图B.3是计算达到稳定后激波间断位置和密度()、速度()u、压力()p和单位质量内能()e的分布。由上述计算结果中可以看出，采用二阶精度迎风型TVD差分格式计算一维可压缩黏性流动问题得到的数值解和经典文献中的结果是完全一致的。计算结果表明，迎风型TVD差分格式能够精确地捕捉激波间断，计算效果较好。由于本问题中黏性较大(50)Re，所以计算得到的激波比较光滑，有一定的宽度。一维可压缩黏性流动问题的解是连续、光滑的。B-2一维可压缩黏性流动问题的数值计算源程序1.C语言源程序//UpwindTVD_1D.c//------------------------------------------------------------------------//二阶迎风型TVD差分格式求解一维可压缩黏性流动问题（C语言版本）图B.2Fortran77语言程序得到的结果图B.3C语言程序得到的结果//-------------------------------------------------------------------------#includestdio.h#includemath.h#definegama1.4#defineTt5.0#defineim201//网格数//全局变量：doubleQ[3][im],Qold[3][im];//Q:[rou,rou*u,E]doublerou[im],u[im],p[im],T[im],E[im],a[im];doublePr,Re,cv,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);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;rou[i]=1.0;u[i]=(x-xl)/(xr-xl)*(ur-ul)+ul;T[i]=(x-xl)/(xr-xl)*(Tr-Tl)+Tl;p[i]=rou[i]*T[i]/(gama*Ma*Ma);E[i]=rou[i]*(cv*T[i]+0.5*u[i]*u[i]);}for(i=0;iim;i++){Q[0][i]=rou[i];Q[1][i]=rou[i]*u[i];Q[2][i]=E[i];}}//-------------------------------------------------------------------doubleqz(doublex){doubleresult;if(fabs(x)0.1){result=fabs(x);}else{result=(x*x/0.1+0.1)/2.0;}returnresult;}//-------------------------------------------------------------------doubleminmod(doublex,doubley){doubleresult;if(x*y=0.0){result=0.0;}elseif(fabs(x)fabs(y)){result=y;}elseif(fabs(x)fabs(y)){result=x;}returnresult;}voidQ2U(){inti;for(i=0;iim;i++){rou[i]=Q[0][i];u[i]=Q[1][i]/rou[i];E[i]=Q[2][i];T[i]=(E[i]/rou[i]-0.5*u[i]*u[i])/cv;p[i]=rou[i]*T[i]/(gama*Ma*Ma);a[i]=sqrt(T[i])/Ma;}}//-------------------------------------------------------------------voidUpwindTVD_1D_Solver(){inti,l,k;doubledq[3][im],lamda[3][im],sigma[3][im],f[3][im];doublefl[3],fwave[3][im],fr[3],g[3][im],Gv[3][im];doublealfa[3][im],beta[3][im],Rn[3][3][im],R[3][3][im];doubleut,at;//临时变量Q2U();for(i=0;iim;i++){f[0][i]=rou[i]*u[i];f[1][i]=rou[i]*u[i]*u[i]+p[i];f[2][i]=u[i]*(E[i]+p[i]);}for(i=0;i=im-2;i++){lamda[0][i]=(u[i]+u[i+1])/2.0;lamda[1][i]=(u[i]-a[i]+u[i+1]-a[i+1])/2.0;lamda[2][i]=(u[i]+a[i]+u[i+1]+a[i+1])/2.0;}for(i=0;i=im-2;i++)for(l=0;l=2;l++){dq[l][i]=Q[l][i+1]-Q[l][i];sigma[l][i]=qz(lamda[l][i])/2.0;}//计算左右特征向量for(i=0;i=im-2;i++){ut=(u[i]+u[i+1])/2.0;at=(a[i]+a[i+1])/2.0;R[0][0][i]=1.0;R[0][1][i]=1.0;R[0][2][i]=1.0;R[1][0][i]=ut;R[1][1][i]=ut-at;R[1][2][i]=ut+at;R[2][0][i]=ut*ut/2.0;R[2][1][i]=at*at/(gama-1.0)+ut*ut/2.0-ut*at;R[2][2][i]=at*at/(gama-1.0)+ut*ut/2.0+ut*at;Rn[0][0][i]=1.0-(gama-1.0)*ut*ut/2.0/at/at;Rn[0][1][i]=(gama-1.0)*ut/at/at;Rn[0][2][i]=-(gama-1.0)/at/at;Rn[1][0][i]=ut/2.0/at+(gama-1.0)*ut*ut/4.0/at/at;Rn[1][1][i]=-(gama-1.0)*ut/2./at/at-1.0/2.0/at;Rn[1][2][i]=(gama-1.0)/2.0/at/at;Rn[2][0][i]=-ut/2.0/at+(gama-1.0)*ut*ut/4.0/at/at;Rn[2][1]