一維黎曼問題數(shù)值解與計算程序

1、實用標(biāo)準(zhǔn)文案一維Riemann問題數(shù)值解與計算程序一維Riemann可題，即激波管問題，是一個典型的一維可壓縮無黏氣體動力 學(xué)問題，并有解析解。對它采用二階精度 MacCormac的步差分格式進(jìn)行數(shù)值求 解。同時，為了初學(xué)者入門和練習(xí)方便，這里給出了用C語言和Fortran7現(xiàn)寫的 計算一維Riemanrn可題的計算程序，供大家學(xué)習(xí)參考。A-1利用MacCormacK步差分格式求解一維Riemann、可題1. 一維 Riemand可題一維Riemann訶題實際上就是激波管問題。激波管是一根兩端封閉、內(nèi)部充滿氣體的直管，如圖A.1所示。在直管中由一薄膜將激波管隔開， 在薄膜兩側(cè)充A/Z/7Z/7

2、/Z/7/77/檢，%文檔邊界條件：x 1和x 1處為自由輸出條件，UoU1 , UN UN 1。3.二階精度MacCormacki分格式MacCormac I斷步差分格式:n1Uj2n11nuj-ujj2jnUj r1 n _uj2fjnfjnnfj11211 n _ fj2(A.4)其中r 0計算實踐表明，MacCormac新步差分格式不能抑制激波附近非物 x理振蕩。因此在計算激波時，必須采用人工黏性濾波方法：一n nUUui ,jui ,j為了在激波附近人工黏性起作用,1 n 。n n-ui 1,j 2ui,jui 1,j而在光滑區(qū)域人工黏性為零,(A.5)需要引入一個與密度（或者壓力）

3、相關(guān)的開關(guān)函數(shù):(A.6)由式（A.6）可以看出，在光滑區(qū)域，密度變化很緩，因此值也很??；而在激波附近密度變化很陡，值就很大。帶有開關(guān)函數(shù)的前置人工黏性濾波方法為：(A.7)(A.8)-nn 1n Q n nui,jui,j -ui 1,j 2ui,jui 1,j其中參數(shù)往往需要通過實際試算來確定，也可采用線性近似方法得到:t|a | 1 t|a|由于聲速不會超過3,所以取|a| 3,在本計算中取0.254.計算結(jié)果分析計算分別采用標(biāo)準(zhǔn)的 C語言和Fortran77g言編寫程序。計算中網(wǎng)格數(shù)取 1000,計算總時間為T 0.4。計算得到在T 0.4時刻的密度、速度和壓力分布 如圖A.2 （C語

4、言計算結(jié)果）和圖A.3 （ Fortran7M言計算結(jié)果）所示。采用兩 種不同語言編寫程序所得到的計算結(jié)果完全吻合。從圖A.2和圖A.3中可以發(fā)現(xiàn)，MacCormac新步差分格式能很好地捕捉激波， 計算得到的激波面很陡、很窄，計算激波精度是很高的。采用帶開關(guān)函數(shù)的前置人工濾波法能消除激波附近的非物理振蕩，計算效果很好。從圖A.2和圖A.3中可以看出通過激波后氣體的密度、壓力和速度都是增加的；在壓力分布中存在第二個臺階，表明在這里存在一個接觸間斷，在接觸間斷 兩側(cè)壓力是有間斷的，而密度和速度是相等的。這個計算結(jié)果正確地反映了一維 Riemann訶題的物理特性，并被激波管實驗所驗證。圖A.2采用C

5、語言程序得到的一維 Riemann題密度、速度和壓力分布圖A.3采用Fortran7旃言程序得到的一維 Riemann可題密度、速度和壓力分布A-2 一維Riemann題數(shù)值計算源程序1. C語言源程序/ MacCormackID.cpp :定義控制臺應(yīng)用程序的入口點。/*# 利用MacCormackE分格式求解一維激波管問題（C語言版本）*/#include "stdafX.h"#include <stdio.h>#include <stdlib.h>#include <math.h># define GAMA 1.4 氣體常數(shù)#def

6、ine PI 3.141592654# define L 2.0計算區(qū)域# define TT 0.4/總時間# define Sf 0.8時間步長因子# define J 1000網(wǎng)格數(shù)/全局變量double UJ+23,UfJ+23,EfJ+23;/*計算時間步長入口： U ,當(dāng)前物理量，dx,網(wǎng)格寬度； 返回：時間步長。*/double CFL(double UJ+23,double dx)int i;double maxvel,p,u,vel;maxvel=1e-100;for(i=1;i<=J;i+)u=Ui1/Ui0;p=(GAMA-1)*(Ui2-0.5*Ui0*u*u);

7、vel=sqrt(GAMA*p/Ui0)+fabs(u);if(vel>maxvel)maxvel=vel;return Sf*dx/maxvel;/*初始化入口 ： 無；出口： U,已經(jīng)給定的初始值，dx, 網(wǎng)格寬度。*/void Init(double UJ+23,double & dx)初始條件int i;double rou1=1.0 ,u1=0.0,p1=1.0; /double rou2=0.125,u2=0.0,p2=0.1;dx=L/J;for(i=0;i<=J/2;i+)Ui0=rou1;Ui1=rou1*u1;Ui2=p1/(GAMA-1)+rou1*u

8、1*u1/2;for(i=J/2+1;i<=J+1;i+)Ui0=rou2;Ui1=rou2*u2;Ui2=p2/(GAMA-1)+rou2*u2*u2/2;/*邊界條件入口： dx ,網(wǎng)格寬度；出口： U ,已經(jīng)給定的邊界。*/void bound(double UJ+23,double dx)int k;/左邊界for(k=0;k<3;k+)U0k=U1k;/右邊界for(k=0;k<3;k+)UJ+1k=UJk;/*根據(jù)U計算E入口： U ,當(dāng)前U矢量；出口： E,計算得到的E矢量，U、E的定義見Euler方程組。*/void U2E(double U3,double

9、E3)double u,p;u=U1/U0;p=(GAMA-1)*(U2-0.5*U1*U1/U0);E0=U1;E1=U0*u*u+p;E2=(U2+p)*u;/*一維MacCormack分格式求解器入口： U , 上一時刻的 U矢量，Uf、Ef,臨時變量，dx ,網(wǎng)格寬度，dt,時間步長； 出口： U ,計算得到的當(dāng)前時刻 U矢量。*/void MacCormack_1DSolver(double UJ+23,doubleEfJ+23,double dx,double dt) int i,k;double r,nu,q;r=dt/dx;nu=0.25;for(i=1;i<=J;i+)

10、 q=fabs(fabs(Ui+10-Ui0)-fabs(Ui0-Ui-1。) /(fabs(Ui+10-Ui0)+fabs(Ui0-Ui-10)+1e-100); for(k=0;k<3;k+)Efik=Uik+0.5*nu*q*(Ui+1k-2*Uik+Ui-1k); for(k=0;k<3;k+) for(i=1;i<=J;i+)Uik=Efik;for(i=0;i<=J+1;i+)U2E(Ui,Efi);for(i=0;i<=J;i+) for(k=0;k<3;k+)Ufik=Uik-r*(Efi+1k-Efik); U(n+1/2)(i+1/2)f

11、or(i=0;i<=J;i+)U2E(Ufi,Efi); E(n+1/2)(i+1/2)for(i=1;i<=J;i+)for(k=0;k<3;k+)Uik=0.5*(Uik+Ufik)-0.5*r*(Efik-Efi-1k); U(n+1)(i) /*UfJ+23,double開關(guān)函數(shù)人工黏性項輸出結(jié)果，用Origin數(shù)據(jù)格式畫圖入口： U,當(dāng)前時刻U矢量，dx,網(wǎng)格寬度；出口 ：無。*/void Output(double UJ+23,double dx) int i;FILE *fp;double rou,u,p;fp=fopen("result.txt&qu

12、ot;,"w");for(i=0;i<=J+1;i+)rou=Ui0;u=Ui1/rou;p=(GAMA-1)*(Ui2-0.5*Ui0*u*u);fprintf(fp,"%20f%20.10e%20.10e%20.10e%20.10en",i*dx,rou,u,p,Ui2); fclose(fp);/*主函數(shù)入口 ： 無； 出口 ：無。*/void main() double T,dx,dt;Init(U,dx);T=0;while(T<TT)dt=CFL(U,dx);T+=dt;printf("T=%10g dt=%10gn&q

13、uot;,T,dt);MacCormack_1DSolver(U,Uf,Ef,dx,dt);bound(U,dx);Output(U,dx);2. Fortran77s言源程序! MacCormackID.forFortran77§言版本)!利用MacCormacki分格式求解一維激波管問題（*/program MacCormackIDimplicit double precision (a-h,o-z)parameter (M=1000)common /G_def/ GAMA,PI,J,JJ,dL,TT,Sfdimension U(0:M+1,0:2),Uf(0:M+1,0:2)d

14、imension Ef(0:M+1,0:2)!氣體常數(shù)GAMA=1.4PI=3.1415926!網(wǎng)格數(shù)J=M!計算區(qū)域dL=2.0!總時間TT=0.4!時間步長因子Sf=0.8call Init(U,dx)T=01 dt=CFL(U,dx)T=T+dtwrite(*,*)'T=',T,'dt=',dtcall MacCormack_1D_Solver(U,Uf,Ef,dx,dt) call bound(U,dx)if(T.lt.TT)goto 1call Output(U,dx)end計算時間步長入口： U,當(dāng)前物理量，dx,網(wǎng)格寬度;返回：時間步長。doubl

15、e precision function CFL(U,dx) implicit double precision (a-h,o-z) common /G_def/ GAMA,PI,J,JJ,dL,TT,Sf dimension U(0:J+1,0:2)dmaxvel=1e-10do 10 i=1,Juu=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=vel10 continueCFL=Sf*dx/dmaxvelend!初始化

16、!入口 ：無；!出口： U,已經(jīng)給定的初始值，dx ,網(wǎng)格寬度。!subroutine Init(U,dx)implicit double precision (a-h,o-z) common /G_def/ GAMA,PI,J,JJ,dL,TT,Sfdimension U(0:J+1,0:2)!初始條件rou1=1.0u1=0v1=0p1=1.0rou2=0.125u2=0v2=0p2=0.1dx=dL/Jdo 20 i=0,J/2U(i,0)=rou1U(i,1)=rou1*u1U(i,2)=p1/(GAMA-1)+0.5*rou1*u1*u120 continuedo 21 i=J/2+

17、1,J+1U(i,0)=rou2U(i,1)=rou2*u2U(i,2)=p2/(GAMA-1)+0.5*rou2*u2*u221 continue end!邊界條件!入口： dx ,網(wǎng)格寬度；!出口： U ,已經(jīng)給定邊界。!subroutine bound(U,dx)implicit double precision (a-h,o-z) common /G_def/ GAMA,PI,J,JJ,dL,TT,Sf dimension U(0:J+1,0:2)!左邊界do 30 k=0,2U(0,k)=U(1,k)30 continue!右邊界do 31 k=0,2U(J+1,k尸U(J,k)31

18、 continueend根據(jù)U計算E入口： U ,當(dāng)前U矢量；出口：E ,計算得到的E矢量，U 、E定義見 Euler方程組。subroutine U2E(U,E,is,in)implicit double precision (a-h,o-z) common /G_def/ GAMA,PI,J,JJ,dL,TT,Sf dimension U(0:J+1,0:2),E(0:J+1,0:2)do 40 i=is,inuu=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*

19、uu+pE(i,2)=(U(i,2)+p)*uu40 continueend一維MacCormackE分格式求解器入口： U , 上一日刻U矢量，Uf 、Ef,臨時變量，dx ,網(wǎng)格寬度，dt,時間步長；出口： U ,計算得到得當(dāng)前時刻 U矢量。subroutine MacCormack_1D_Solver(U,Uf,Ef,dx,dt) implicit double precision (a-h,o-z)common /G_def/ GAMA,PI,J,JJ,dL,TT,Sfdimension U(0:J+1,0:2),Uf(0:J+1,0:2)dimension Ef(0:J+1,0:2)r=dt/dxdnu=0.25do 60 i=1,Jdo 60 k=0,2!開關(guān)函數(shù)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) 60 continue do 61 k=0,2 do 61 i=1,JU(i,k尸Ef(i,k