《附录C 二维斜激波刚壁反射》由会员分享，可在线阅读，更多相关《附录C 二维斜激波刚壁反射（23页珍藏版）》请在金锄头文库上搜索。
1、附录C 二维斜激波在平面刚壁上反射问题的数值解与计算程序斜激波在平面刚壁上反射问题是具有解析解的二维可压缩无黏流动问题。对 它采用具有二阶精度Lax-Wendrof两步差分格式进行数值求解。同时，为了初学 者入门和练习方便,这里给出了由C语言和Fortran77吾言编写的计算斜激波在平 面刚壁上反射问题的计算程序，供大家学习参考。C-1利用二阶 Lax-Wendroff两步差分格式求解斜激波在平面刚性壁面上反射问题1.二维斜激波在平面刚性壁面上反射问题设一斜激波以29。入射角投向一平面刚性壁面，并与此平面刚性壁面相撞, 并在刚性壁面上反射，形成入射波和反射波，如图C.1所示。图C.1二维斜激波
2、在平面刚性壁面上反射问题示意图2基本方程组、初始条件和边界条件设气体为理想气体。斜激波在平面刚性壁面上反射问题在数学上可以用二维可压缩无黏流动Eule方程组来描述。量纲为一的二维Eule方程组为:其中:血+f+电=0dtdxdy(C.1)f pu 、f pv、u =pu，/ =pu 2 + p,g =puv(C.1a)pvpuvpv 2 + p1 E J1 (E + p)u J1 (E + p)v J这里p、u、v、p和E分别是密度、X和y方向的速度分量、压力和单位体积的总能。理想气体状态方程:C.1b)(Y -l)pe = (丫一 1) E-2 p(u2 + v2)初始条件：在t二0时，斜激
3、波刚好和平面刚性壁面相遇，根据激波马赫数和 入射角 = 29，可计算出初始流场分布如下:C.2a)左边界(x = 0, y = 0 -1.0)为均匀来流条件:u 2.90, v = 0.0, p = 1.0, p = 0.714291 1 1 1上边界(x 0 - 4.0, y 1.0)为均匀来流条件:u 2.61934, v - 0.506322 2 (C.2b) p 1.69997, p 1.5281922边界条件：右边界(x 4.0, y 1.0)为自由输出条件：C.2c)du Qv 0dx dyC.2d)下边界为刚性壁面，在水平刚性壁面上法向速度v 0w3.二阶精度Lax-Wendro
4、两步差分格式斜激波在平面刚性壁面上反射问题是一个二维无黏流动问题，可以采用算子 分裂法把它分裂成两个一维流动问题。分裂后对每个一维流动问题采用二阶精度 Lax-Wendro两步差分格式进行数值计算。为了尽可能地减少人为因素，可采用 如下分裂计算步骤：f 1 At L (At)L | 1 At12 Jyx (2 丿Un+1 LxC.3)其中L (At)是方程式(C. 1)在x方向的差分算子:xC.3a)更 + df - 0dt dx其中Ly(M是方程式QD在y方向的差分算子:du+ dtdg (u )dxC.3b)对每个一为无黏流动采用二阶精度Lax-Wendrof两步差分格式:+丄1j+丄22
5、Un+1 = Unj j(C.4)计算实践表明，二阶精度Lax-Wendro两步差分格式不能抑制激波附近的非物理振荡。因此在计算激波时，必须采用带开关函数的前置人工黏性滤波方法:U n = Un + 耳 9 Un 2U+ U ）i, j i, j 2i+1, j i, ji1, j(C.5)参数耳往往需要通过实际试算来确定，也可采用线性近似方法得到:Ip.- P.p pii 1|p.+1 - P.+p pii19 =(C.5a)r =Ax1 Ax丿(C.5b)由于声速不超过3,所以取I a I= 3，r = 0.25 。4.计算结果分析计算分别采用标准的C语言和Fortran77吾言编写程序。
6、计算网格数取100x400，计算时间为T = 3.0。计算得到的在T = 3.0时刻的密度等值线如图C.2（C语言计算结果）和图C.3（ Fortran77吾言计算结果）所示。图C.4（ C语言计算结果）和图C.5（ Fortran77吾言计算结果）分别是在T = 3.0时刻，在y = 0.5位置上，沿x方向的压力分布。对图C.2、图C.3、图C.4和图C.5进行比较， 可以发现两个程序的计算结果是完全吻合的。从图C.4和图C.5中可以发现，压力分布存在两个压力台阶，第一个台阶是 通过入射激波后压力增加，而第二个台阶是通过反射激波后压力再一次增加。图C.2采用C语言程序计 算得到的密度等值线图
7、-C.3 -图C.3采用Fortran77吾言程序计算得到的密度等值线图图C.4采用C语言程序计算得到的在t = 3.0时刻、y = 0.5位置沿x方向的压力分布图C.5采用Fortran77吾言程序计算得到的t二3.0时刻、y = 0.5位置沿x方向的压力分布C-2二维斜激波在平面刚性壁面上反射问题的数值计算源程序1. C语言源程序/ Lax-Wendroff2D.cpp : 定义控制台应用程序的入口点。/*利用 Lax-Wendroff 差分格式求解二维平面激波反射问题( C 语言版本) * */#include stdafx.h#include #include #include #de
8、fine GAMA 1.4/气体常数#define PI 3.141592654#define MIN(x,y) (x)(y)?(x):(y)#define Lx 4.0/计算区域#define Ly 1.0#define TT 3.0/总时间#define Sf 0.8/时间步长因子#define Jx 400/网格数#define Jy 100/全局变量double UJx+2Jy+24,UfJx+2Jy+24,EFfJx+2Jy+24;/*计算时间步长入口： U,当前物理量，dx、dy,网格宽度；返回: 时间步长。*/double CFL(double UJx+2Jy+24,double
9、 dx,double dy)int i,j;double maxvel,p,u,v,vel;maxvel=1e-100;for(i=1;i=Jx;i+)for(j=1;jmaxvel)maxvel=vel;return Sf*MIN(dx,dy)/maxvel;/*初始化入口: 无；出口： U,已经给定的初始值，dx、dy,网格宽度。*/void Init(double UJx+2Jy+24,double & dx,double & dy)/初始化int i,j;double rou1=1.0,u1=2.9,v1=0,p1=0.71429; /初始条件double rou2=1.69997,u
10、2=2.61934,v2=-0.50632,p2=1.52819;double theta,lx;dx=Lx/Jx;dy=Ly/Jy;theta=29*PI/180;入射激波角度 29 for(j=0;j=Jy+1;j+)for(i=0;i=Jx+1;i+)lx=(1-j*dy)/tan(theta);if(i*dx=lx)Uij0=rou1;Uij1=rou1*u1;Uij2=rou1*v1;Uij3=p1/(GAMA-1)+rou1*(u1*u1+v1*v1)/2;elseUij0=rou2;Uij1=rou2*u2;Uij2=rou2*v2;Uij3=p2/(GAMA-1)+rou2*(
11、u2*u2+v2*v2)/2;/*边界条件入口: dx、dy， 网格宽度； 出口: U， 已经给定边界。 */ void bound(double UJx+2Jy+24,double dx,double dy) int i,j,k;double rou1=1.0, u1=2.9, v1=0,p1=0.71429; /初始条件double rou2=1.69997,u2=2.61934,v2=-0.50632,p2=1.52819;/左边界for(j=0;j=Jy+1;j+)U0j0=rou1;U0j1=rou1*u1;U0j2=rou1*v1;U0j3=p1/(GAMA-1)+rou1*(u1
12、*u1+v1*v1)/2;/右边界 for(j=0;j=Jy+1;j+)for(k=0;k4;k+)UJx+1jk=UJxjk;/上边界 for(i=0;i=Jx+1;i+)UiJy+10=rou2;UiJy+11=rou2*u2;UiJy+12=rou2*v2;UiJy+13=p2/(GAMA-1)+rou2*(u2*u2+v2*v2)/2; /下边界 for(i=0;i=Jx+1;i+)Ui00=Ui10;Ui01=Ui11;Ui02=0;Ui03=Ui13;/*根据U计算E入口： U,当前U矢量；出口： E,计算得到的E矢量，U、E,定义见Euler方程组。*/void U2E(double U4,double E4)double u,v,p; u=U1/U0; v=U2/U0;p=(GAMA-1)*(U3-0.5*(U1*U1+U2*U2)/U0); E0=U1;E1=U0*u*u+p;E2=U0*u*v; E3=(U3+p)*u;/*根据 U 计算 F入口： U，当前U矢量；出口： E，计算得到的F矢量，U、F定义见Euler方程组。*/void U2F(double U4,double F4)/根据 U 计算 Fdouble u,v,p; u=U1/U0; v=U2/U0;p
