附录B一维可压缩黏性流动问题的数值解法与计算程序.docx

一维可压缩黏性流动是气体动力学中最经典的黏性流动问题，对它采用迎风 型TVDg分算法进行数值求解同时，为了初学者入门和练习方便，这里给出了 由C语言和Fortran77语言编写的、计算一维可压缩黏性流动问题的计算程序，供 大家学习参考B-1利用二阶迎风型TVD差分格式求解一维可压缩黏性流动问题1. 一维可压缩黏性流动问题在两端开口管道中充满了可压缩黏性流体当黏性流体以超声速从左向右运 动时，一定会在管道中形成一道正激波，如图B.1所示p、u、p和p、u、p分 1 1 1 2 2 2别为激波波前和波后的参数该问题可简化为一维可压缩黏性流动问题当数值 解达到稳定时，在管道中可求解得到一道稳定的激波111 Pl A1「2 1 W1 u2 ——P\ - A外一X A图B.1 一维可压缩黏性流动问题示意图2. 基本方程组、初始条件和边界条件设流体是黏性流体一维可压缩黏性流动问题，在数学上可以用一维可压缩黏性流动N-S方程组来描述量纲为一的一维N-S方程组为:8u 8f 8s一 + 二=——8tdxdx（B.1）其中_ p -/ pu、r s)u =pu，f =pu 2 + p,s =sE«E + p )匕l S3 ^04日8 u业空+CJL四v 3 Re 8 x PrRr 8 x /（B.1a）PT p = ,y M 28砰—pU L 「Re 8—8—, c旦81=y （y -1） M 28c,y =〜，cVcxk8a = M8(B.lb)其中p、u、p和E分别是量纲为一的密度、速度、压力和单位体积总能，s为流体的黏性项。

Pr为普朗特数（此处公式中c「七、七是有量纲量），Re为雷诺数，c为比定容热容，c为比定压热容，c , c是量纲为一的量，称为气体绝热指数，V p p Va为当地声速求解区域为 0 < x <1取M = 2, Re = 50, Pr = 0.72, y= 1.4, k =日=T 初始条件：在t = 0时刻，p（x,0）=1，其他物理量采用线性插值得到B.2)边界条件：左边界 G = 0）处：p（0,t）=u（0,t）= T（0,t）=1右边界（x = 1）处:8p(x, t)8 xT (1, t )=x=12 + M 2u (1, t )=y-1 8土 M 2y - 1 8y-1)「y -1TTT J〔E +(y+1) MI(2 y—M 2 - ky +1 8(B.3)3. 二阶精度迎风型TVE差分格式其中一维N-S方程组中的对流项采用Harten的二阶精度迎风型TVE差分格式:hk j+2u n+1 = u n - rj j1 （h 广 2 f + fj + 22d s .+ Atd x-hnj -2 J\_ _ + R 妙j j+1 .1 .1k j+ 2 j+ 2)(B.4)(B.4a)r =冬。

向量垂在第i个特征方向上分量中i为: Ax* r中11 = gj + gj+1 -Q a +y j+2\l1j+2)j+(B.4b)gi = minmod bi ai ,bi ai j .1.1 .1.1\ J + 2 j+2 j - 2 j - 2 /g"+1 - g" a i '0,=(R-1)J+2|Z| >EQ (z )』b i.1i +2x A ui.1J + 2|z| 京'.,1k i+2 J4q[r ],2 k，+」J2(0.05<£ <0.5)一 rk i 2.1i+2非定常定常流通量矢量f的非线性Jiacobia系数矩阵A = R A R-1为:0 1A (u )=^f = - 3-Y u 2 (3-Y)ua u 2Y-2 ua2 a 2 3-Yk U 一”一1 「一1+ ^\0Y-1Y u非线性Jiacobian系数矩阵A的特征值为：(X 0 0)A = (j X 0k 0 02 X JX = u, X = u — a, X = u + a(B.4c)(B.4d)(B.4e)(B.4f)(B.4g)(B.5)(B.6)非线性血瓦4系数矩阵A的右特征矢量R、R-1为:( \勘1 1 1R = u u — a u + a1—u 2 h — ua h + uak 2 Jr G-Du 2(y - 1)uy-1'1 —2a 2a2a2u G-1)u 2(y - 1)u1y-1V_—1 —+— —里・12a 4a 22a 22a2a 2u(Y- 1)u 2(y - 1)u1y-1+ +——[2a 4a22a22a2a2 J(B.8)一维N-S方程组中的黏性项若也采用二阶精度中心差分格式。

4. 计算结果分析采用用C语言和Fortran?语言对一维可压缩黏性流动问题编制了计算程序，并对雷诺数Re = 50的流动进行了计算计算结果如图B.2和图B.3所示图B.2 Fortran?!语言程序得到的结果图B.3 C语言程序得到的结果图B.2和图B.3是计算达到稳定后激波间断位置和密度(p)、速度(u)、压 力(p)和单位质量内能(e)的分布由上述计算结果中可以看出，采用二阶精度迎 风型TVD差分格式计算一维可压缩黏性流动问题得到的数值解和经典文献中的 结果是完全一致的计算结果表明，迎风型TVE差分格式能够精确地捕捉激波间 断，计算效果较好由于本问题中黏性较大(Re = 50)，所以计算得到的激波比较 光滑，有一定的宽度一维可压缩黏性流动问题的解是连续、光滑的B-2 一维可压缩黏性流动问题的数值计算源程序1. C语言源程序// UpwindTVD_1D.c// //二阶迎风型TVD差分格式求解一维可压缩黏性流动问题(C语言版本)// #include "stdio.h”#include "math.h”#define gama 1.4#define Tt 5.0#define im 201 //网格数〃全局变量：double Q[3][im],Qold[3][im] ; //Q: [rou, rou*u, E]double rou[im],u[im],p[im],T[im],E[im],a[im];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 ;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; i0.1){result=fabs(x);}else{result=(x*x/0.1+0.1)/2.0 ;}return result ;}// double minmod(double x,double y){double result ;if(x*y<=0.0){result=0.0 ;}else if(fabs(x)>fabs(y)){result=y ;}else if(fabs(x)