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

一維問(wèn)題數(shù)值解與計(jì)算程序一維問(wèn)題，即激波管問(wèn)題，是一個(gè)典型的一維可壓縮無(wú)黏氣體動(dòng)力學(xué)問(wèn)題，并有解析解。對(duì)它采用二階精度兩步差分格式進(jìn)行數(shù)值求解。同時(shí)，為了初學(xué)者入門和練習(xí)方便,這里給出了用語(yǔ)言和編寫的計(jì)算一維問(wèn)題的計(jì)算程序，供大家學(xué)習(xí)參考。A-1利用兩步差分格式求解一維問(wèn)題1.一維問(wèn)題圖A.1激波管問(wèn)題示意圖一維問(wèn)題實(shí)際上就是激波管問(wèn)題。激波管是一根兩端封閉、內(nèi)部充滿氣體的直管，如圖A.1所示。在直管中由一薄膜將激波管隔開，在薄膜兩側(cè)充有均勻理想氣體〔可以是同一種氣體，也可以是不同種氣體〕，薄膜兩側(cè)氣體的壓力、密度不同。當(dāng)時(shí)，氣體處于靜止?fàn)顟B(tài)。當(dāng)時(shí)，薄膜瞬時(shí)突然破裂，氣體從高壓端沖向低壓端，同時(shí)在管內(nèi)形成激波、稀疏波和接觸間斷等復(fù)雜波系。圖A.1激波管問(wèn)題示意圖2.根本方程組、初始條件和邊界條件設(shè)氣體是理想氣體。一維問(wèn)題在數(shù)學(xué)上可以用一維可壓縮無(wú)黏氣體方程組來(lái)描述。在直角坐標(biāo)系下量綱為一的一維方程組為：(A.1)其中()這里、、、分別是流體的密度、速度、壓力和單位體積總能。理想氣體狀態(tài)方程：〔A.3〕初始條件：；。邊界條件：和處為自由輸出條件，,。3.二階精度差分格式兩步差分格式：()其中。計(jì)算實(shí)踐說(shuō)明，兩步差分格式不能抑制激波附近非物理振蕩。因此在計(jì)算激波時(shí)，必須采用人工黏性濾波方法：()為了在激波附近人工黏性起作用，而在光滑區(qū)域人工黏性為零，需要引入一個(gè)與密度〔或者壓力〕相關(guān)的開關(guān)函數(shù)：()由式()可以看出，在光滑區(qū)域，密度變化很緩，因此值也很小；而在激波附近密度變化很陡，值就很大。帶有開關(guān)函數(shù)的前置人工黏性濾波方法為：()其中參數(shù)往往需要通過(guò)實(shí)際試算來(lái)確定，也可采用線性近似方法得到：()由于聲速不會(huì)超過(guò)3，所以取，在本計(jì)算中取。4.計(jì)算結(jié)果分析計(jì)算分別采用標(biāo)準(zhǔn)的語(yǔ)言和語(yǔ)言編寫程序。計(jì)算中網(wǎng)格數(shù)取，計(jì)算總時(shí)間為。計(jì)算得到在時(shí)刻的密度、速度和壓力分布如圖A.2〔語(yǔ)言計(jì)算結(jié)果〕和圖A.3〔語(yǔ)言計(jì)算結(jié)果〕所示。采用兩種不同語(yǔ)言編寫程序所得到的計(jì)算結(jié)果完全吻合。從圖A.2和圖A.3中可以發(fā)現(xiàn)，兩步差分格式能很好地捕捉激波，計(jì)算得到的激波面很陡、很窄，計(jì)算激波精度是很高的。采用帶開關(guān)函數(shù)的前置人工濾波法能消除激波附近的非物理振蕩，計(jì)算效果很好。從圖A.2和圖A.3中可以看出通過(guò)激波后氣體的密度、壓力和速度都是增加的；在壓力分布中存在第二個(gè)臺(tái)階，說(shuō)明在這里存在一個(gè)接觸間斷，在接觸間斷兩側(cè)壓力是有間斷的，而密度和速度是相等的。這個(gè)計(jì)算結(jié)果正確地反映了一維問(wèn)題的物理特性，并被激波管實(shí)驗(yàn)所驗(yàn)證。圖A.圖A.2采用語(yǔ)言程序得到的一維問(wèn)題密度、速度和壓力分布圖采用語(yǔ)言程序得到的一維問(wèn)題密度、速度和壓力分布A-2一維問(wèn)題數(shù)值計(jì)算源程序1.語(yǔ)言源程序//MacCormack1D.cpp:定義控制臺(tái)應(yīng)用程序的入口點(diǎn)。///*-----------------------------------------------------------------------------------------------------*利用差分格式求解一維激波管問(wèn)題〔語(yǔ)言版本〕*-------------------------------------------------------------------------------------------------------*///#include"stdafx.h"#include<stdio.h>#include<stdlib.h>#include<math.h>#defineJ1000//網(wǎng)格數(shù)//全局變量doubleU[J+2][3],Uf[J+2][3],Ef[J+2][3];/*-------------------------------------------------------計(jì)算時(shí)間步長(zhǎng)入口:U，當(dāng)前物理量，dx，網(wǎng)格寬度；返回:時(shí)間步長(zhǎng)。---------------------------------------------------------*/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;}/*-------------------------------------------------------初始化入口:無(wú)；出口:U，已經(jīng)給定的初始值，dx,網(wǎng)格寬度。---------------------------------------------------------*/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，網(wǎng)格寬度；出口:U，已經(jīng)給定的邊界。---------------------------------------------------------*/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];}/*-------------------------------------------------------根據(jù)U計(jì)算E入口:U，當(dāng)前U矢量；出口:E，計(jì)算得到的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;}/*-------------------------------------------------------一維差分格式求解器入口:U，上一時(shí)刻的U矢量，Uf、Ef，臨時(shí)變量，dx，網(wǎng)格寬度，dt,時(shí)間步長(zhǎng)；出口:U，計(jì)算得到的當(dāng)前時(shí)刻U矢量。---------------------------------------------------------*/voidMacCormack_1DSolver(doubleU[J+2][3],doubleUf[J+2][3],doubleEf[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);//開關(guān)函數(shù)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]);//人工黏性項(xiàng)}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)}/*-------------------------------------------------------輸出結(jié)果,用數(shù)據(jù)格式畫圖入口:U，當(dāng)前時(shí)刻U矢量，dx，網(wǎng)格寬度；出口:無(wú)。---------------------------------------------------------*/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);}/*-------------------------------------------------------主函數(shù)入口:無(wú)；出口:無(wú)。---------------------------------------------------------*/voidmain(){doubleT,dx,dt;Init(U,dx);T=0;while(T<TT){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.語(yǔ)言源程序!MacCormack1D.for------------------------------------------------------------------------------------------------------------!利用差分格式求解一維激波管問(wèn)題〔語(yǔ)言版本〕--------------------------------------------------------------------------------------------------------------*/programMacCormack1Dimplicitdoubleprecision(a-h,o-z)parameter(M=1000)common/G_def/GAMA,PI,J,JJ,dL,TT,SfdimensionU(0:M+1,0:2),Uf(0:M+1,0:2)dimensionEf(0:M+1,0:2)!氣體常數(shù)!網(wǎng)格數(shù)J=M!計(jì)算區(qū)域!總時(shí)間!時(shí)間步長(zhǎng)因子callInit(U,dx)T=01dt=CFL(U,dx)T=T+dtwrite(*,*)'T=',T,'dt=',dtcallMacCormack_1D_Solver(U,Uf,Ef,dx,dt) callbound(U,dx)if(T.lt.TT)goto1callOutput(U,dx)end!-------------------------------------------------------!計(jì)算時(shí)間步長(zhǎng)!入口:U，當(dāng)前物理量，dx,網(wǎng)格寬度；!返回:時(shí)間步長(zhǎng)。!-------------------------------------------------------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-10do10i=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=vel10continueCFL=Sf*dx/dmaxvelend!-------------------------------------------------------!初始化!入口:無(wú)；!出口:U,已經(jīng)給定的初始值，dx，網(wǎng)格寬度。!-------------------------------------------------------subroutineInit(U,dx)implicitdoubleprecision(a-h,o-z)common/G_def/GAMA,PI,J,JJ,dL,TT,SfdimensionU(0:J+1,0:2)!初始條件u1=0v1=0u2=0v2=0dx=dL/Jdo20i=0,J/2U(i,0)=rou1U(i,1)=rou1*u1U(i,2)=p1/(GAMA-1)+0.5*rou1*u1*u120continuedo21i=J/2+1,J+1U(i,0)=rou2U(i,1)=rou2*u2U(i,2)=p2/(GAMA-1)+0.5*rou2*u2*u221continueend!-------------------------------------------------------!邊界條件!入口:dx，網(wǎng)格寬度；!出口:U，已經(jīng)給定邊界。!-------------------------------------------------------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,2U(J+1,k)=U(J,k)31continueend!-------------------------------------------------------!根據(jù)U計(jì)算E!入口:U，當(dāng)前U矢量；!出口:E，計(jì)算得到的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,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)*uu40continueend!-------------------------------------------------------!一維差分格式求解器!入口:U，上一時(shí)刻U矢量，!Uf、Ef，臨時(shí)變量，!dx，網(wǎng)格寬度，dt,，時(shí)間步長(zhǎng)；!出口:U，計(jì)算得到得當(dāng)前時(shí)刻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/dxdo60i=1,Jdo60k=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))60continuedo61k=0,2do61i=1,JU(i,k)=Ef(i,k)61continuecallU2E(U,Ef,0,J+1)do63i=0,Jdo63k=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,Jdo64k=0,2!U(n+1)(i)U(i,k)=0.5*(U(i,k)+Uf(i,k))