附录B一维可压NS方程B.docxVIP

附录B一维可压NS方程B 附录B一维可压NS方程B PAGE 附录B一维可压NS方程B 附录B一维可压缩黏性流动问题的数值解法与计算程序 一维可压缩黏性流动是气体动力学中最经典的黏性流动问题，对它采用迎风型差分算法进行数值求解。同时，为了初学者入门和练习方便,这里给出了由语言和语言编写的、计算一维可压缩黏性流动问题的计算程序，供大家学习参考。 B-1利用二阶迎风型差分格式求解一维可压缩黏性流动问题 1.一维可压缩黏性流动问题 在两端开口管道中充满了可压缩黏性流体。当黏性流体以超声速从左向右运动时，一定会在管道中形成一道正激波，如图所示。和分别为激波波前和波后的参数。该问题可简化为一维可压缩黏性流动问题。当数值解达到稳定时，在管道中可求解得到一道稳定的激波。 图一维可压缩黏性 图一维可压缩黏性流动问题示意图 2.基本方程组、初始条件和边界条件 设流体是黏性流体。一维可压缩黏性流动问题，在数学上可以用一维可压缩黏性流动方程组来描述。量纲为一的一维方程组为： （） 其中 （） （） 其中和分别是量纲为一的密度、速度、压力和单位体积总能，为流体的黏性项。为普朗特数（此处公式中是有量纲量），为雷诺数，为比定容热容，为比定压热容，是量纲为一的量，称为气体绝热指数，为当地声速。 求解区域为。取。 初始条件：在时刻，，其他物理量采用线性插值得到。 边界条件：左边界处：（） 右边界处： （） 3.二阶精度迎风型差分格式 一维方程组中的对流项采用的二阶精度迎风型差分格式： （） （） 其中。向量在第个特征方向上分量为： （） （） （） （） （） （） 流通量矢量的非线性系数矩阵为： （） 非线性系数矩阵的特征值为： （） 非线性系数矩阵的右特征矢量为： （） （） 一维方程组中的黏性项采用二阶精度中心差分格式。 4.计算结果分析 采用用语言和语言对一维可压缩黏性流动问题编制了计算程序，并对雷诺数的流动进行了计算，计算结果如图和图所示。 图语言程序得到的结果 图语言程序得到的结果 图语言程序得到的结果 图和图是计算达到稳定后激波间断位置和密度、速度、压力和单位质量内能的分布。由上述计算结果中可以看出，采用二阶精度迎风型差分格式计算一维可压缩黏性流动问题得到的数值解和经典文献中的结果是完全一致的。计算结果表明，迎风型差分格式能够精确地捕捉激波间断，计算效果较好。由于本问题中黏性较大，所以计算得到的激波比较光滑，有一定的宽度。一维可压缩黏性流动问题的解是连续、光滑的。 B-2一维可压缩黏性流动问题的数值计算源程序 1.语言源程序 at/at; Rn[1][2][i]=/at/at; Rn[2][0][i]=-ut/at+*ut*ut/at/at; Rn[2][1][i]=*ut/2./at/at+at; Rn[2][2][i]=/2./at/at; } for(i=0;i=im-1;i++) for(l=0;l=2;l++) { alfa[l][i]=; fwave[l][i]=; g[l][i]=; } u[i]+u[i-1])+(T[i+1]-T[i-1]) *(u[i+1]-u[i-1])/4.)*4./3./Re/(dx*dx); Gv[2][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]) *(u[i+1]-u[i-1])/4.)*4./3.+(T[i]*(T[i+1]-2.*T[i]+T[i-1]) +(T[i+1]-T[i-1])*(T[i+1]-T[i-1])/4.)*cp/Pr)/Re/(dx*dx); } for(i=1;i=im-2;i++) for(l=0;l=2;l++) { Q[l][i]=Q[l][i]+Gv[l][i]*dt; } 言源程序 ! ! !二阶迎风型差分格式求解一维可压缩黏性流动问题 !（语言版本） ! programUPWIND_TVD_1D implicitreal*8(a-h,o-z) parameter(mx=201,Tt= dimensionQ(3,mx),Qold(3,mx)!Q:[rou,rou*u,E] dimensionrou(mx),u(mx),p(mx),T(mx),E(mx) real*8Ma common/para_def/Sf,Re,Ma,pr,gama,cp,cv,dx,dt,im im=mx callInitialize(Q) time= n=0 1n=n+1 time=time+dt doi=1,im dol