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

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

2、想氣體（可以是同一種氣體，也可以是不同種氣體），薄膜兩側(cè)氣體的壓力、密度不同。當(dāng) t 0 時，氣體處于靜止?fàn)顟B(tài)。當(dāng) t 0 時，薄膜瞬時突然破裂，氣體從高壓端沖向低壓端，同時在管內(nèi)形成激波、稀疏波和接觸間斷等復(fù)雜波系。2. 基本方程組、初始條件和邊界條件圖 A.1 激波管問題示意圖設(shè)氣體是理想氣體。一維Riemann問題在數(shù)學(xué)上可以用一維可壓縮無黏氣體Euler方程組來描述。在直角坐標(biāo)系下量綱為一的一維Euler方程組為：uf,1x1(A.1)tx0u其中uu , fu 2p(A.2)E(E p)u這里 、 u 、 p 、 E 分別是流體的密度、速度、壓力和單位體積總能。理想氣體狀態(tài)方程：p1

3、e1E1u2v2（A.3）2初始條件： 1 1, u10, p1 1 ； 20.125, u20, p20.1。供參考邊界條件： x1 和 x1 處為自由輸出條件， u0u1 , uNuN 1 。3. 二階精度 MacCormack差分格式MacCormack兩步差分格式：n1unjr f jnf jn 1uj2(A.4)n 11nn 11n 1n122f j2uj2u juj2r f j 1其中 rt 。計算實踐表明， MacCormack兩步差分格式不能抑制激波附近非物x理振蕩。因此在計算激波時，必須采用人工黏性濾波方法：uin,juin, j1uin1, j2uin, juin 1, j

4、(A.5)2為了在激波附近人工黏性起作用， 而在光滑區(qū)域人工黏性為零， 需要引入一個與密度（或者壓力）相關(guān)的開關(guān)函數(shù)：i1iii1(A.6)i1iii1由式 (A.6) 可以看出，在光滑區(qū)域，密度變化很緩，因此 值也很??；而在激波附近密度變化很陡， 值就很大。帶有開關(guān)函數(shù)的前置人工黏性濾波方法為：uin,juin, j1uin1, j 2uin, j uin 1, j(A.7)2其中參數(shù)往往需要通過實際試算來確定，也可采用線性近似方法得到：t | a | 1t | a |(A.8)xx由于聲速不會超過3，所以取 | a |3 ，在本計算中取0.25 。4. 計算結(jié)果分析計算分別采用標(biāo)準(zhǔn)的 C

5、語言和 Fortran77語言編寫程序。計算中網(wǎng)格數(shù)取 1000 ，計算總時間為 T 0.4。計算得到在 T 0.4 時刻的密度、速度和壓力分布如圖 A.2 （ C 語言計算結(jié)果）和圖 A.3 （ Fortran77語言計算結(jié)果）所示。采用兩種不同語言編寫程序所得到的計算結(jié)果完全吻合。從圖 A.2 和圖 A.3 中可以發(fā)現(xiàn)，MacCormack兩步差分格式能很好地捕捉激波，計算得到的激波面很陡、 很窄，計算激波精度是很高的。 采用帶開關(guān)函數(shù)的前置供參考人工濾波法能消除激波附近的非物理振蕩，計算效果很好。從圖 A.2 和圖 A.3 中可以看出通過激波后氣體的密度、壓力和速度都是增加的；在壓力分布

6、中存在第二個臺階，表明在這里存在一個接觸間斷，在接觸間斷兩側(cè)壓力是有間斷的， 而密度和速度是相等的。 這個計算結(jié)果正確地反映了一維Riemann問題的物理特性，并被激波管實驗所驗證。圖A.2采 用 C 語 言 程 序 得 到 的 一 維圖 A.3采用 Fortran77語言程序得到的一維Riemann問題密度、速度和壓力分布Riemann問題密度、速度和壓力分布A-2 一維 Riemann問題數(shù)值計算源程序1.C 語言源程序/ MacCormack1D.cpp :定義控制臺應(yīng)用程序的入口點。/*-* 利用 MacCormack差分格式求解一維激波管問題（C 語言版本）*-*/#include

7、"stdafx.h"#include <stdio.h>#include <stdlib.h>#include <math.h>#define GAMA 1.4/ 氣體常數(shù)#define PI3.141592654#define L2.0/ 計算區(qū)域供參考#define TT 0.4/ 總時間#define Sf 0.8/ 時間步長因子#define J1000/網(wǎng)格數(shù)/全局變量double UJ+23,UfJ+23,EfJ+23;/*-計算時間步長入口 : U ，當(dāng)前物理量，dx，網(wǎng)格寬度；返回 : 時間步長。-*/double CFL

8、(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);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=

9、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*u1*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

10、(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 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, 時間

11、步長；出口 : U ， 計算得到的當(dāng)前時刻U 矢量。-*/void MacCormack_1DSolver(double UJ+23,double UfJ+23,double EfJ+23,double dx,double dt)int i,k;double r,nu,q;r=dt/dx;nu=0.25;for(i=1;i<=J;i+)q=fabs(fabs(Ui+10-Ui0)-fabs(Ui0-Ui-10)/(fabs(Ui+10-Ui0)+fabs(Ui0-Ui-10)+1e-100); /開關(guān)函數(shù)for(k=0;k<3;k+)Efik=Uik+0.5*nu*q*(Ui+1k

12、-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)for(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-1

13、k); /U(n+1)(i)/*-輸出結(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","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.10

14、e%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=%10gdt=%10gn",T,dt);MacCormack_1DSolver(U,Uf,Ef,dx,dt);bound(U,dx);Output(U,dx);-供參考2.Fortran77語言源程序! MacCormack1D.for-! 利用 MacCormack差分格式求解

15、一維激波管問題（ Fortran77語言版本）-*/program MacCormack1Dimplicit 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)dimension 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=

16、T+dtwrite(*,*)'T=',T,'dt=',dtcall MacCormack_1D_Solver(U,Uf,Ef,dx,dt)call bound(U,dx)call Output(U,dx)end!-! 計算時間步長! 入口 : U ， 當(dāng)前物理量， dx, 網(wǎng)格寬度；! 返回 : 時間步長。!-供參考double precision function CFL(U,dx)implicit double precision (a-h,o-z)common /G_def/ GAMA,PI,J,JJ,dL,TT,Sfdimension U(0:J+1,0

17、: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)10 continue CFL=Sf*dx/dmaxvel end!-! 初始化! 入口: 無；! 出口 : 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

18、:2)! 初始條件 rou1=1.0 u1=0 v1=0 p1=1.0 rou2=0.125 u2=0 v2=0 p2=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+1,J+1 U(i,0)=rou2供參考U(i,1)=rou2*u2U(i,2)=p2/(GAMA-1)+0.5*rou2*u2*u221 continue end!-! 邊界條件! 入口 : dx ，網(wǎng)格寬度；! 出口 : U ， 已經(jīng)給定邊界。!-subroutine

19、 bound(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)! 左邊界do 30 k=0,2U(0,k)=U(1,k)30 continue! 右邊界do 31 k=0,2 U(J+1,k)=U(J,k)31 continue end!-! 根據(jù)U計算E! 入口 : U，當(dāng)前 U 矢量；! 出口 : E，計算得到的E 矢量，! U、E 定義見 Euler 方程組。!-subroutine U2E(U,E,is,in)implicit double

20、 precision (a-h,o-z)common /G_def/ GAMA,PI,J,JJ,dL,TT,Sfdimension 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*uu+pE(i,2)=(U(i,2)+p)*uu40continue供參考end!-! 一維 MacCormack差分格式求解器!入口 : U， 上一時刻U 矢量，! Uf 、 Ef ，臨時變量，! dx，網(wǎng)格寬

21、度， 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)