1、精選優(yōu)質(zhì)文檔-傾情為你奉上一維問題數(shù)值解與計(jì)算程序一維問題，即激波管問題，是一個(gè)典型的一維可壓縮無黏氣體動(dòng)力學(xué)問題，并有 解析解。對(duì)它采用二階精度兩步差分格式進(jìn)行數(shù)值求解。同時(shí)，為了初學(xué)者入門和練習(xí)方便,這里給出了用語言和編寫的計(jì)算一維問題的計(jì)算程序，供大家學(xué)習(xí)參考。A-1利用兩步差分格式求解一維問題1.一維問題圖A.1 激波管問題示意圖一維問題實(shí)際上就是激波管問題。激波管是一根兩端封閉、內(nèi)部充滿氣體的直管，如圖A.1所示。在直管中由一薄膜將激波管隔開，在薄膜兩側(cè)充有均勻理想氣體（可以是同一種氣體，也可以是不同種氣體），薄膜兩側(cè)氣體的壓力、密度不同。當(dāng)時(shí)，氣體處于靜止?fàn)顟B(tài)。當(dāng)時(shí)，薄膜瞬時(shí)突然破

2、裂，氣體從高壓端沖向低壓端，同時(shí)在管內(nèi)形成激波、稀疏波和接觸間斷等復(fù)雜波系。2.基本方程組、初始條件和邊界條件設(shè)氣體是理想氣體。一維問題在數(shù)學(xué)上可以用一維可壓縮無黏氣體方程組來描述。在直角坐標(biāo)系下量綱為一的一維方程組為： (A.1)其中 (A.2) 這里、分別是流體的密度、速度、壓力和單位體積總能。理想氣體狀態(tài)方程： （A.3）初始條件：；。邊界條件：和處為自由輸出條件，,。3.二階精度差分格式兩步差分格式： (A.4)其中。計(jì)算實(shí)踐表明，兩步差分格式不能抑制激波附近非物理振蕩。因此在計(jì)算激波時(shí)，必須采用人工黏性濾波方法： (A.5)為了在激波附近人工黏性起作用，而在光滑區(qū)域人工黏性為零，需要

3、引入一個(gè)與密度（或者壓力）相關(guān)的開關(guān)函數(shù)： (A.6)由式(A.6) 可以看出，在光滑區(qū)域，密度變化很緩，因此值也很小；而在激波附近密度變化很陡，值就很大。帶有開關(guān)函數(shù)的前置人工黏性濾波方法為： (A.7)其中參數(shù)往往需要通過實(shí)際試算來確定，也可采用線性近似方法得到：(A.8)由于聲速不會(huì)超過3，所以取，在本計(jì)算中取。4.計(jì)算結(jié)果分析計(jì)算分別采用標(biāo)準(zhǔn)的語言和語言編寫程序。計(jì)算中網(wǎng)格數(shù)取，計(jì)算總時(shí)間為。計(jì)算得到在時(shí)刻的密度、速度和壓力分布如圖A.2（語言計(jì)算結(jié)果）和圖A.3（語言計(jì)算結(jié)果）所示。采用兩種不同語言編寫程序所得到的計(jì)算結(jié)果完全吻合。從圖A.2和圖A.3中可以發(fā)現(xiàn)，兩步差分格式能很好地

4、捕捉激波，計(jì)算得到的激波面很陡、很窄，計(jì)算激波精度是很高的。采用帶開關(guān)函數(shù)的前置人工濾波法能消除激波附近的非物理振蕩，計(jì)算效果很好。從圖A.2和圖A.3中可以看出通過激波后氣體的密度、壓力和速度都是增加的；在壓力分布中存在第二個(gè)臺(tái)階，表明在這里存在一個(gè)接觸間斷，在接觸間斷兩側(cè)壓力是有間斷的，而密度和速度是相等的。這個(gè)計(jì)算結(jié)果正確地反映了一維問題的物理特性，并被激波管實(shí)驗(yàn)所驗(yàn)證。 圖A.2采用語言程序得到的一維問題密度、速度和壓力分布圖A.3采用語言程序得到的一維問題密度、速度和壓力分布A-2 一維問題數(shù)值計(jì)算源程序1. 語言源程序/ MacCormack1D.cpp : 定義控制臺(tái)應(yīng)用程序的入

5、口點(diǎn)。/*-*利用差分格式求解一維激波管問題（語言版本）* -*/#include stdafx.h #include #include #include #define GAMA 1.4/氣體常數(shù)#define PI 3.#define 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+23,double dx) int i

6、; double maxvel,p,u,vel; maxvel=1e-100; for(i=1;imaxvel)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/(

7、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的定義見Euler方程組。-*/void U2E(d

8、ouble 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;/*-一維差分格式求解器入口: U， 上一時(shí)刻的U矢量，Uf、Ef，臨時(shí)變量， dx，網(wǎng)格寬度，dt, 時(shí)間步長(zhǎng)；出口: U， 計(jì)算得到的當(dāng)前時(shí)刻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/d

9、x; 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;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=Uik-r*(Efi+1k-Efik

10、); /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é)果, 用數(shù)據(jù)格式畫圖入口: U， 當(dāng)前時(shí)刻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=

11、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(TTT) dt=CFL(U,dx); T+=dt; printf(T=%10g dt=%10gn,T,dt); MacCormack_1DSolver(U,Uf,Ef,dx,dt); bou

12、nd(U,dx); Output(U,dx);-2. 語言源程序! MacCormack1D.for -!利用差分格式求解一維激波管問題（語言版本） -*/ program MacCormack1D implicit double precision (a-h,o-z) parameter (M=1000) common /G_def/ GAMA,PI,J,JJ,dL,TT,Sf dimension U(0:M+1,0:2),Uf(0:M+1,0:2) dimension Ef(0:M+1,0:2)!氣體常數(shù) GAMA=1.4 PI=3.!網(wǎng)格數(shù) J=M!計(jì)算區(qū)域 dL=2.0!總時(shí)間 TT=

13、0.4!時(shí)間步長(zhǎng)因子 Sf=0.8 call 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 1 call Output(U,dx) end!-!計(jì)算時(shí)間步長(zhǎng)!入口: U， 當(dāng)前物理量，dx, 網(wǎng)格寬度；!返回: 時(shí)間步長(zhǎng)。!- double precision function CFL(U,dx) implicit double precision (a-h,o-z) common

14、 /G_def/ GAMA,PI,J,JJ,dL,TT,Sf dimension U(0:J+1,0:2) dmaxvel=1e-10 do 10 i=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=vel10 continue CFL=Sf*dx/dmaxvel end !-!初始化!入口: 無；!出口: U, 已經(jīng)給定的初始值，dx，網(wǎng)格寬度。!- subroutine Init(U,dx) implici

15、t 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.0 u1=0 v1=0 p1=1.0 rou2=0.125 u2=0 v2=0 p2=0.1 dx=dL/J do 20 i=0,J/2 U(i,0)=rou1 U(i,1)=rou1*u1 U(i,2)=p1/(GAMA-1)+0.5*rou1*u1*u120 continue do 21 i=J/2+1,J+1 U(i,0)=rou2 U(i,1)=rou2*u2 U(i,2)=p2/(GA

16、MA-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,2 U(0,k)=U(1,k)30 continue!右邊界 do 31 k=0,2 U(J+1,k)=U(J,k)31 continue end!-!根據(jù)U計(jì)算E!入口: U，當(dāng)前U矢量

17、；!出口: E，計(jì)算得到的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,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)

18、+p)*uu40 continue end!-!一維差分格式求解器!入口: 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/dx dnu=0.25 do 60 i=1,J do 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)!人工黏性項(xiàng) 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,J U(i,k)=Ef(i,k)61 continue call U2E(U,Ef,