一维黎曼问题数值解与计算程序.docx
- 文档编号:10570374
- 上传时间:2023-02-21
- 格式:DOCX
- 页数:15
- 大小:100.63KB
一维黎曼问题数值解与计算程序.docx
《一维黎曼问题数值解与计算程序.docx》由会员分享,可在线阅读,更多相关《一维黎曼问题数值解与计算程序.docx(15页珍藏版)》请在冰豆网上搜索。
一维黎曼问题数值解与计算程序
一维Riemann问题数值解与计算程序
一维Rieman"可题,即激波管问题,是一个典型的一维可压缩无黏气体动力
学问题,并有解析解。
对它采用二阶精度MacCormac两步差分格式进行数值求解。
同时,为了初学者入门和练习方便,这里给出了用C语言和Fortran77编写的计算一维Riemanri可题的计算程序,供大家学习参考。
A-1利用MacCormac两步差分格式求解一维Riemann问题
1.一维Riemanr问题
一维Riemand可题实际上就是激波管问题。
激波管是一根两端封闭、内部充
满气体的直管,如图A.1所示。
在直管中由一薄膜将激波管隔开,在薄膜两侧充
有均匀理想气体(可以是同一种气体,也可以是不同种气体),薄膜两侧气体的压力、密度不同。
当t乞0时,气体处于静止状态。
当t=0时,薄膜瞬时突然破裂,气体从高压端冲向低压端,同时在管内形成激波、稀疏波和接触间断等复杂波系。
2.基本方程组、初始条件和边界条件
//////////////////////////////////////
图A.1激波管问题示意图
设气体是理想气体。
一维Riemanri可题在数学上可以用一维可压缩无黏气体
Eulei方程组来描述。
在直角坐标系下量纲为一的一维Euler方程组为:
■^+f=0,—1兰x兰1(A.1)
.:
t:
x
rP、
「Pu、
其中
u=
Pu
f=
Pu2+p
(A.2)
l(E+p)u』 这里「、u、p、E分别是流体的密度、速度、压力和单位体积总能。 理想气体 状态方程: p=-1◎-1E-Ju2v2(A.3) 初始条件: J=1,5=0,山=1;=0.125,u2=0,p2=0.1。 边界条件: X--1和X=1处为自由输出条件,UohU—UnhUn/。 3.二阶精度MacCormac差分格式 n-2 MacCormac两步差分格式: 2nn-n Uj-Uj-rfj-fjj 其中r二卫。 计算实践表明,MacCormac两步差分格式不能抑制激波附近非物 Ax 理振荡。 因此在计算激波时,必须采用人工黏性滤波方法: (A.5) 需要引入一个与 -nn1nQnn Ui,j=Ui,j-U「1,j-2嘔-UiJ,j 为了在激波附近人工黏性起作用,而在光滑区域人工黏性为零, 密度(或者压力)相关的开关函数: 由式(A.6)可以看出,在光滑区域,密度变化很缓,因此二值也很小;而在激波附近密度变化很陡,V值就很大。 带有开关函数的前置人工黏性滤波方法为: (A.7) _1 Uinj=unj+罗日(ut-2uinj+ulj 其中参数往往需要通过实际试算来确定,也可采用线性近似方法得到: (A.8) 由于声速不会超过3,所以取|a|=3,在本计算中取=0.25。 4.计算结果分析 计算分别采用标准的C语言和Fortran77语言编写程序。 计算中网格数取1000,计算总时间为T=0.4。 计算得到在T=0.4时刻的密度、速度和压力分布 如图A.2(C语言计算结果)和图A.3(Fortran77语言计算结果)所示。 采用两 种不同语言编写程序所得到的计算结果完全吻合。 从图A.2和图A.3中可以发现,MacCormac两步差分格式能很好地捕捉激波,计算得到的激波面很陡、很窄,计算激波精度是很高的。 采用带开关函数的前置 人工滤波法能消除激波附近的非物理振荡,计算效果很好 从图A.2和图A.3中可以看出通过激波后气体的密度、压力和速度都是增加的;在压力分布中存在第二个台阶,表明在这里存在一个接触间断,在接触间断两侧压力是有间断的,而密度和速度是相等的。 这个计算结果正确地反映了一维Riemann问题的物理特性,并被激波管实验所验证。 A-2一维Riemanr问题数值计算源程序 1.C语言源程序 //MacCormackID.cpp: 定义控制台应用程序的入口点。 // /* *利用MacCormac差分格式求解一维激波管问题(C语言版本) * */ 〃#include"stdafx.h" #include #include #include #definePI3.141592654 #defineL2.0〃计算区域 #defineTT0.4//总时间 #defineSf0.8//时间步长因子 #defineJ1000//网格数 //全局变量 doubleU[J+2][3],Uf[J+2][3],Ef[J+2][3]; /* 计算时间步长 入口: U,当前物理量,dx,网格宽度; 返回: 时间步长。 */ doubleCFL(doubleU[J+2][3],doubledx) { inti; doublemaxvel,p,u,vel; maxvel=1e-100; for(i=1;i<=J;i++) { u=U[i][1]/U[i][0]; p=(GAMA-1)*(U[i][2]-0.5*U[i][0]*u*u); vel=sqrt(GAMA*p/U[i][0])+fabs(u);if(vel>maxvel)maxvel=vel; } returnSf*dx/maxvel; } /* 初始化 入口: 无; 出口: U,已经给定的初始值, dx,网格宽度。 */ voidInit(doubleU[J+2][3],double&dx) { inti; doublerou1=1.0,u1=0.0,p1=1.0;//初始条件 doublerou2=0.125,u2=0.0,p2=0.1; dx=L/J; for(i=0;i<=J/2;i++) { U[i][0]=rou1; U[i][1]=rou1*u1;U[i][2]=p1/(GAMA-1)+rou1*u1*u1/2; } for(i=J/2+1;i<=J+1;i++) { U[i][0]=rou2; U[i][1]=rou2*u2; U[i][2]=p2/(GAMA-1)+rou2*u2*u2/2; } } /* 边界条件 入口: dx,网格宽度;出口: U,已经给定的边界。 */voidbound(doubleU[J+2][3],doubledx){ intk; //左边界for(k=0;k<3;k++)U[0][k]=U[1][k]; //右边界 for(k=0;k<3;k++)U[J+1][k]=U[J][k]; } /* 根据U计算E入口: U,当前U矢量;出口: E,计算得到的E矢量, U、E的定义见Euler方程组。 */voidU2E(doubleU[3],doubleE[3]){ doubleu,p; u=U[1]/U[0];p=(GAMA-1)*(U[2]-0.5*U[1]*U[1]/U[0]); E[0]=U[1]; E[1]=U[0]*u*u+p;E[2]=(U[2]+p)*u; } /* 一维MacCormac差分格式求解器 入口: u,上一时刻的u矢量,Uf、Ef,临时变量, dx,网格宽度,dt,时间步长; 出口: u,计算得到的当前时刻u矢量。 */ Ef[J+2][3],double voidMacCormack_1DSolver(doubleu[J+2][3],doubleuf[J+2][3],doubledx,doubledt) { inti,k; doubler,nu,q; r=dt/dx; nu=0.25; for(i=1;i<=J;i++) { q=fabs(fabs(u[i+1][0]-u[i][0])-fabs(u[i][0]-u[i-1][0]))/(fabs(u[i+1][0]-u[i][0])+fabs(u[i][0]-u[i-1][0])+1e-100);//开关函数 for(k=0;k<3;k++) Ef[i][k]=u[i][k]+0.5*nu*q*(u[i+1][k]-2*u[i][k]+u[i-1][k]);//人工黏性项} for(k=0;k<3;k++)for(i=1;i<=J;i++)u[i][k]=Ef[i][k]; for(i=0;i<=J+1;i++)u2E(u[i],Ef[i]); for(i=0;i<=J;i++) for(k=0;k<3;k++)uf[i][k]=u[i][k]-r*(Ef[i+1][k]-Ef[i][k]);//u(n+1/2)(i+1/2) for(i=0;i<=J;i++)u2E(uf[i],Ef[i]);//E(n+1/2)(i+1/2) for(i=1;i<=J;i++) for(k=0;k<3;k++)u[i][k]=0.5*(u[i][k]+uf[i][k])-0.5*r*(Ef[i][k]-Ef[i-1][k]);//u(n+1)(i) } /* 输出结果,用Origin数据格式画图 入口: u,当前时刻u矢量,dx,网格宽度; 出口: 无。 */ voidOutput(doubleU[J+2][3],doubledx) { inti; FILE*fp; doublerou,u,p; fp=fopen("result.txt","w"); for(i=0;i<=J+1;i++) { rou=U[i][0]; u=U[i][1]/rou; p=(GAMA-1)*(U[i][2]-0.5*U[i][0]*u*u); fprintf(fp,"%20f%20.10e%20.10e%20.10e%20.10e\n",i*dx,rou,u,p,U[i][2]);}fclose(fp); } /* 主函数入口: 无;出口: 无。 */ voidmain() { doubleT,dx,dt; Init(U,dx); T=0;while(T {dt=CFL(U,dx); T+=dt; printf("T=%10gdt=%10g\n",T,dt); MacCormack_1DSolver(U,Uf,Ef,dx,dt); bound(U,dx); } Output(U,dx); } 2.Fortran7语言源程序 ! MacCormack1D.for ! 利用MacCormac差分格式求解一维激波管问题(Fortran77语言版本)*/programMacCormack1Dimplicitdoubleprecision(a-h,o-z)parameter(M=1000)common/G_def/GAMA,PI,J,JJ,dL,TT,Sf dimensionU(0: M+1,0: 2),Uf(0: M+1,0: 2)dimensionEf(0: M+1,0: 2) ! 气体常数 GAMA=1.4PI=3.1415926 ! 网格数 J=M ! 计算区域 dL=2.0 ! 总时间 TT=0.4 ! 时间步长因子 Sf=0.8 callInit(U,dx) T=0 1dt=CFL(U,dx) T=T+dtwrite(*,*)'T=',T,'dt=',dtcallMacCormack_1D_Solver(U,Uf,Ef,dx,dt)callbound(U,dx) if(T.lt.TT)goto1 callOutput(U,dx)end ! ! 计算时间步长 ! 入口: U,当前物理量,dx,网格宽度; ! 返回: 时间步长。 ! doubleprecisionfunctionCFL(U,dx)implicitdoubleprecision(a-h,o-z)common/G_def/GAMA,PI,J,JJ,dL,TT,SfdimensionU(0: J+1,0: 2) dmaxvel=1e-10 do10i=1,J uu=U(i,1)/U(i,0)p=(GAMA-1)*(U(i,2)-0.5*U(i,0)*uu*uu)vel=dsqrt(GAMA*p/U(i,0))+dabs(uu)if(vel.gt.dmaxvel)dmaxvel=vel 10continue CFL=Sf*dx/dmaxvel end ! ! 初始化 ! 入口: 无; ! 出口: U,已经给定的初始值,dx,网格宽度。 ! subroutineInit(U,dx)implicitdoubleprecision(a-h,o-z)common/G_def/GAMA,PI,J,JJ,dL,TT,SfdimensionU(0: J+1,0: 2) ! 初始条件 rou1=1.0 u1=0 v1=0 p1=1.0 rou2=0.125 u2=0 v2=0 p2=0.1 dx=dL/J do20i=0,J/2 U(i,0)=rou1U(i,1)=rou1*u1 U(i,2)=p1/(GAMA-1)+0.5*rou1*u1*u1 20continue do21i=J/2+1,J+1 U(i,0)=rou2 U(i,1)=rou2*u2 U(i,2)=p2/(GAMA-1)+0.5*rou2*u2*u221continue end ! ! 边界条件 ! 入口: dx,网格宽度; ! 出口: U,已经给定边界。 ! subroutinebound(U,dx)implicitdoubleprecision(a-h,o-z)common/G_def/GAMA,PI,J,JJ,dL,TT,SfdimensionU(0: J+1,0: 2) ! 左边界 do30k=0,2 U(0,k)=U(1,k) 30continue ! 右边界 do31k=0,2 U(J+1,k)=U(J,k) 31continue end ! ! 根据U计算E ! 入口: U,当前U矢量; ! 出口: E,计算得到的E矢量, ! U、E定义见Euler方程组。 ! subroutineU2E(U,E,is,in)implicitdoubleprecision(a-h,o-z)common/G_def/GAMA,PI,J,JJ,dL,TT,SfdimensionU(0: J+1,0: 2),E(0: J+1,0: 2) do40i=is,in uu=U(i,1)/U(i,0)p=(GAMA-1)*(U(i,2) $-0.5*U(i,1)*U(i,1)/U(i,0)) E(i,0)=U(i,1) E(i,1)=U(i,0)*uu*uu+p E(i,2)=(U(i,2)+p)*uu 40continue end ! ! 一维MacCormac差分格式求解器 ! 入口: U,上一时刻U矢量, ! Uf、Ef,临时变量, ! dx,网格宽度,dt,,时间步长; ! 出口: U,计算得到得当前时刻U矢量。 ! subroutineMacCormack_1D_Solver(U,Uf,Ef,dx,dt)implicitdoubleprecision(a-h,o-z)common/G_def/GAMA,PI,J,JJ,dL,TT,SfdimensionU(0: J+1,0: 2),Uf(0: J+1,0: 2)dimensionEf(0: J+1,0: 2) r=dt/dx dnu=0.25 do60i=1,J do60k=0,2 ! 开关函数 q=dabs(dabs(U(i+1,0)-U(i,0))-dabs(U(i,0)-U(i-1,0))) $/(dabs(U(i+1,0)-U(i,0))+dabs(U(i,0)-U(i-1,0))+1e-10)! 人工黏性项 Ef(i,k)=U(i,k)+0.5*dnu*q*(U(i+1,k)-2*U(i,k)+U(i-1,k)) 60continue do61k=0,2 do61i=1,J U(i,k)=Ef(i,k) 61continue callU2E(U,Ef,0,J+1) do63i=0,J do63k=0,2 ! U(n+1/2)(i+1/2) Uf(i,k)=U(i,k)-r*(Ef(i+1,k)-Ef(i,k)) 63continue ! E(n+1/2)(i+1/2) callU2E(Uf,Ef,0,J) do64i=1,J do64k=0,2 ! U(n+1)(i) U(i,k)=0.5*(U(i,k)+Uf(i,k))-0.5*r*(Ef(i,k)-Ef(i-1,k)) 64continue end ! ! 输出结果,用Origin数据格式画图 ! 入口: U,当前时刻U矢量, ! dx,网格宽度; ! 出口: 无。 ! subroutineOutput(U,dx) implicitdoubleprecision(a-h,o-z)common/G_def/GAMA,PI,J,JJ,dL,TT,SfdimensionU(0: J+1,0: 2) open(1,file='result.txt',status='unknown') do80i=0,J+1 rou=U(i,0) uu=U(i,1)/roup=(GAMA-1)*(U(i,2)-0.5*U(i,0)*uu*uu)write(1,81)i*dx,rou,uu,p,U(i,2) 80continue close (1) 81format(D20.10,D20.10,D20.10,D20.10,D20.10)end
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 黎曼 问题 数值 计算 程序