一維黎曼問(wèn)題數(shù)值解與計(jì)算程序

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

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

3、總能。理想氣體狀態(tài)方程：p= -1 -1 E - J u2 v2（ A.3）初始條件：J = 1,5 = 0,山=1 ；= 0.125, u2 = 0, p2 = 0.1。邊界條件：X - -1和X=1處為自由輸出條件，UohUUnhUn/。3. 二階精度MacCormac差分格式n -2MacCormac兩步差分格式:2 nn - nUj -Uj -r fj - fjju；12、n nUj Ujn4 起)jf-f j 2(A.4)其中r二衛(wèi)。計(jì)算實(shí)踐表明，MacCormac兩步差分格式不能抑制激波附近非物Ax理振蕩。因此在計(jì)算激波時(shí)，必須采用人工黏性濾波方法:(A.5)需要引入一個(gè)與-nn

4、1nQ nnUi,j =Ui,j- U1,j -2嘔-UiJ,j為了在激波附近人工黏性起作用，而在光滑區(qū)域人工黏性為零,密度（或者壓力）相關(guān)的開(kāi)關(guān)函數(shù):Ipi-Pi-R-pjpi卅-Pi+Pi-Pe =(A.6)由式（A.6）可以看出，在光滑區(qū)域，密度變化很緩，因此 二值也很小；而在激波 附近密度變化很陡，V值就很大。帶有開(kāi)關(guān)函數(shù)的前置人工黏性濾波方法為:(A.7)_ 1Uinj =unj +羅日（ut - 2uinj + ulj其中參數(shù)往往需要通過(guò)實(shí)際試算來(lái)確定，也可采用線(xiàn)性近似方法得到:(A.8)由于聲速不會(huì)超過(guò)3,所以取|a|=3，在本計(jì)算中取=0.25。4. 計(jì)算結(jié)果分析計(jì)算分別采用標(biāo)

5、準(zhǔn)的C語(yǔ)言和Fortran77語(yǔ)言編寫(xiě)程序。計(jì)算中網(wǎng)格數(shù)取 1000，計(jì)算總時(shí)間為T(mén) =0.4。計(jì)算得到在T =0.4時(shí)刻的密度、速度和壓力分布如圖A.2 （ C語(yǔ)言計(jì)算結(jié)果）和圖A.3 （ Fortran77語(yǔ)言計(jì)算結(jié)果）所示。采用兩種不同語(yǔ)言編寫(xiě)程序所得到的計(jì)算結(jié)果完全吻合。從圖A.2和圖A.3中可以發(fā)現(xiàn)，MacCormac兩步差分格式能很好地捕捉激波, 計(jì)算得到的激波面很陡、很窄，計(jì)算激波精度是很高的。采用帶開(kāi)關(guān)函數(shù)的前置人工濾波法能消除激波附近的非物理振蕩，計(jì)算效果很好從圖A.2和圖A.3中可以看出通過(guò)激波后氣體的密度、 壓力和速度都是增加 的；在壓力分布中存在第二個(gè)臺(tái)階，表明在這里存

6、在一個(gè)接觸間斷，在接觸間斷 兩側(cè)壓力是有間斷的，而密度和速度是相等的。這個(gè)計(jì)算結(jié)果正確地反映了一維 Riema nn問(wèn)題的物理特性，并被激波管實(shí)驗(yàn)所驗(yàn)證。4 2a cin.I” Il.I 00.51.01.52.00.0 -0.0 - -d-Ii，I”I.0.D0.5101.52.DD6n-5 4 2 o.D.g.圖A.2采用C語(yǔ)言程序得到的一維Riema nn問(wèn)題密度、速度和壓力分布圖A.3采用Fortran77語(yǔ)言程序得到的一維Riema nn可題密度、速度和壓力分布A-2 一維Riemanr問(wèn)題數(shù)值計(jì)算源程序1. C語(yǔ)言源程序/ MacCormackID.cpp :定義控制臺(tái)應(yīng)用程序的入

7、口點(diǎn)。/*利用MacCormac差分格式求解一維激波管問(wèn)題（C語(yǔ)言版本）*/#in clude stdafx.h#in clude #in clude #in clude #define GAMA 1.4 氣體常數(shù)#define PI 3.141592654#defi ne L 2.0計(jì)算區(qū)域#define TT 0.4/ 總時(shí)間#define Sf 0.8/ 時(shí)間步長(zhǎng)因子#define J 1000/ 網(wǎng)格數(shù)/全局變量double UJ+23,UfJ+23,EfJ+23;/*計(jì)算時(shí)間步長(zhǎng)入口 ： U,當(dāng)前物理量，dx，網(wǎng)格寬度；返回 : 時(shí)間步長(zhǎng)。*/double CFL(double UJ

8、+23,double dx)int i;double maxvel,p,u,vel;maxvel=1e-100;for(i=1;imaxvel)maxvel=vel;return Sf*dx/maxvel;/*初始化入口 ： 無(wú)；出口 ： 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; U

9、i2=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(k=0;k3;k+)U0k=U1k;/右邊界for(k=0;k3;k+)UJ+1k=UJk;/*根據(jù) U 計(jì)算 E 入口 ： U ， 當(dāng)前 U 矢量； 出口 ： E， 計(jì)算得到的 E 矢量，U、 E 的定義見(jiàn) Euler

10、 方程組。 */ 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;/*一維MacCormac差分格式求解器入口 ： u,上一時(shí)刻的 u矢量，Uf、Ef,臨時(shí)變量，dx，網(wǎng)格寬度，dt,時(shí)間步長(zhǎng)；出口 ： u, 計(jì)算得到的當(dāng)前時(shí)刻 u 矢量。 */EfJ+23,doublevoid MacCormack_1DSolver(double uJ+23,double ufJ+23,double dx,double dt)int i,

11、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); / 開(kāi)關(guān)函數(shù)for(k=0;k3;k+)Efik=uik+0.5*nu*q*(ui+1k-2*uik+ui-1k);/ 人工黏性項(xiàng) for(k=0;k3;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;k3;k+) ufik=ui

12、k-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;k3;k+) uik=0.5*(uik+ufik)-0.5*r*(Efik-Efi-1k); /u(n+1)(i)/*輸出結(jié)果 , 用 Origin 數(shù)據(jù)格式畫(huà)圖入口 ： u, 當(dāng)前時(shí)刻 u 矢量, dx, 網(wǎng)格寬度；出口 ： 無(wú)。*/void Output(double UJ+23,double dx)int i;FILE *fp;double rou,u,p;fp=fopen(resul

13、t.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.10e%20.10en,i*dx,rou,u,p,Ui2); fclose(fp);/*主函數(shù) 入口 : 無(wú)； 出口 : 無(wú)。*/void main()double T,dx,dt;Init(U,dx);T=0; while(TTT) dt=CFL(U,dx);T+=dt;printf(T=%10g dt=%10gn,T,dt);MacCormack_1DSolver(U,Uf,Ef,

14、dx,dt);bound(U,dx);Output(U,dx);2. Fortran7語(yǔ)言源程序! MacCormack1D.for!利用MacCormac差分格式求解一維激波管問(wèn)題(Fortran77語(yǔ)言版本) */ program MacCormack1D implicit 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.4 PI=

15、3.1415926! 網(wǎng)格數(shù)J=M! 計(jì)算區(qū)域dL=2.0! 總時(shí)間TT=0.4! 時(shí)間步長(zhǎng)因子Sf=0.8call Init(U,dx)T=01 dt=CFL(U,dx)T=T+dt write(*,*)T=,T,dt=,dt call MacCormack_1D_Solver(U,Uf,Ef,dx,dt) call bound(U,dx)if(T.lt.TT)goto 1call Output(U,dx) end! 計(jì)算時(shí)間步長(zhǎng)! 入口 : U ， 當(dāng)前物理量， dx, 網(wǎng)格寬度；! 返回 : 時(shí)間步長(zhǎng)。!double precision function CFL(U,dx) implic

16、it 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!初始化!入口 : 無(wú)；!出口 ： U,已經(jīng)給定的初始值，dx,網(wǎng)格寬度。!su

17、broutine Init(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)! 初始條件rou1=1.0u1=0v1=0p1=1.0rou2=0.125u2=0v2=0p2=0.1dx=dL/Jdo 20 i=0,J/2U(i,0)=rou1 U(i,1)=rou1*u1U(i,2)=p1/(GAMA-1)+0.5*rou1*u1*u120 continuedo 21 i=J/2+1,J+1U(i,0)=rou2U(i,1)=rou2*u2U

18、(i,2)=p2/(GAMA-1)+0.5*rou2*u2*u2 21 continueend!邊界條件!入口 ： 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 continueend! 根據(jù) U 計(jì)算 E!

19、 入口 ： U ，當(dāng)前 U 矢量；! 出口 ： E ，計(jì)算得到的 E 矢量，!U 、E 定義見(jiàn) 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*uu

20、+pE(i,2)=(U(i,2)+p)*uu40 continueend! 一維MacCormac差分格式求解器! 入口 : U ， 上一時(shí)刻 U 矢量，!Uf 、 Ef ，臨時(shí)變量，!dx，網(wǎng)格寬度，dt,時(shí)間步長(zhǎng)；! 出口 : U ， 計(jì)算得到得當(dāng)前時(shí)刻 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,Sf dimension 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! 開(kāi)關(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) ! 人工黏性項(xiàng)Ef(i,k)=U(i,k)+0.5*dnu*q*(U(i+1,k)-2*U(i,k)+U