附錄F有限體積法計(jì)算方腔流(F)

實(shí)用文案附錄F 二維不可壓縮黏性流體方腔流動(dòng)問(wèn)題的有限體積算法與計(jì)算程序二維方腔流動(dòng)問(wèn)題是一個(gè)不可壓縮黏性流動(dòng)中典型流動(dòng)。 雖然目前尚不能求得它的解析解，但是它常被用來(lái)作為檢驗(yàn)各種數(shù)值算法計(jì)算精度和可靠性的算例。文獻(xiàn)中幾乎大多數(shù)算法都對(duì)它進(jìn)行過(guò)計(jì)算。 在本算例中采用有限體積算法三階迎風(fēng)型QUICK離散格式對(duì)它進(jìn)行數(shù)值求解。同時(shí)，為了初學(xué)者入門(mén)和練習(xí)方便,這里給出了用C語(yǔ)言和Fortran77語(yǔ)言編寫(xiě)的計(jì)算二維不可壓縮黏性方腔流動(dòng)問(wèn)題計(jì)算程序，供大家學(xué)習(xí)參考。F-1利用有限體積算法三階迎風(fēng)型 QUICK離散格式求解二維不可壓縮黏性流體方腔流動(dòng)問(wèn)題1.二維不可壓縮黏性流體方腔流動(dòng)問(wèn)題二維不可壓縮黏性流體方腔流動(dòng)（ cavity flow）：有一正方形腔室，其量綱為一的寬度為1.0，里面充滿靜止的不可壓縮黏性流體，方腔內(nèi)初始時(shí)刻壓力和密度為p=1.0、=1.0它周?chē)诿妫ㄗ笥冶诿婧?/p>

圖F.1二維不可壓縮黏性方腔流問(wèn)題示意圖底面）固定不動(dòng)，上壁面以量綱為一的速度1.0沿著上壁面方向自左向右運(yùn)動(dòng)（圖F.1）?；痉匠探M、初始條件和邊界條件設(shè)流體是黏性流體。二維方腔流動(dòng)問(wèn)題在數(shù)學(xué)上可以由二維不可壓縮黏性流動(dòng)N-S方程組來(lái)表示，把它寫(xiě)成通用變量的微分方程組形式，有：圖F.1二維不可壓縮黏性方uv腔流動(dòng)問(wèn)題示意圖（F.1）Sxyxxyy標(biāo)準(zhǔn)實(shí)用文案其中u為變量 在水平x方向的流速，v為 在垂直y方向的流速， 為黏度，S為源項(xiàng)。源項(xiàng)中不僅包含壓力梯度項(xiàng)，也包含時(shí)間導(dǎo)數(shù)項(xiàng)。初始條件：方腔上壁面以量綱為一的速度 1.0沿著上壁面方向自左向右運(yùn)動(dòng)。邊界條件：流動(dòng)速度u、v均可采用無(wú)滑移邊界條件， 壓力p采用自由輸出邊界條件。3.計(jì)算網(wǎng)格劃分和控制體單元與節(jié)點(diǎn)定義采用交錯(cuò)網(wǎng)格，圖 F.2和圖F.3是計(jì)算網(wǎng)格、控制體單元和節(jié)點(diǎn)示意圖。圖F.2方腔流動(dòng)計(jì)算網(wǎng)格、圖F.3計(jì)算采用的交錯(cuò)網(wǎng)格示意圖控制體單元和節(jié)點(diǎn)示意圖節(jié)點(diǎn)P所在主控制單元如圖 F.2中有陰影部分所示。在 x方向與節(jié)點(diǎn)P相鄰的節(jié)點(diǎn)為W和E，在y方向與節(jié)點(diǎn)P相鄰的節(jié)點(diǎn)為S和N，主控制單元界面分別為w、s、e、n。壓力p和速度u、v分別在三套不同網(wǎng)格中如圖 F.3中有陰影部分所示。4.有限體積算法三階迎風(fēng)型 QUICK離散格式標(biāo)準(zhǔn)實(shí)用文案對(duì)方程（F.1）在圖F.2所示節(jié)點(diǎn)P所在控制體單元內(nèi)積分，有：udVVyvdVVxdVVydVSdV（F.2）VxxyV由于二維不可壓縮黏性流體方腔流動(dòng)是二維問(wèn)題，因此控制體單元體積V僅是面積，而它的邊界是長(zhǎng)度。設(shè)AwAey,AsAnx，利用Gauss定理，可將方程（F.2）改寫(xiě)成如下有限體積算法離散格式：uAeuAwvAnvAsAAAAS（F.3）xyxexwynys對(duì)上式中、、、采用一階向前差分近似，則有：xexwynysxy

E P,e x xN P,n y y

Wwx（F.4）Ssy同時(shí)記：FeueAe,FwuwAw（F.5）FnvnAn,FsvsAsDeeAe,DwwAw,xPExPW（F.6）nAn,DssAsDnyPNyPS則可由式（F.2）寫(xiě)成：FeeFwwFnnFssDeEPDwPW（F.7）DnNPDsPSSxy式中P、E、W、N、S、De、Dw、Dn、Ds都是控制體單元內(nèi)節(jié)點(diǎn)上的已知量，如果利用差分計(jì)算得到控制體單元邊界上的流通量Fee、Fww、Fnn、Fss，就可以求出節(jié)點(diǎn)上未知量P、E、W、N、S。標(biāo)準(zhǔn)實(shí)用文案圖F.4三階迎風(fēng)型 QUICK離散格式示意圖為了便于討論，現(xiàn)對(duì)一維對(duì)流擴(kuò)散方程的三階迎風(fēng)型 QUICK離散格式進(jìn)行分析：在三階迎風(fēng)型QUICK離散格式中，計(jì)算主控制單元界面上流動(dòng)量 需要取主控制單元界面兩側(cè) 3個(gè)節(jié)點(diǎn)處的流動(dòng)量值進(jìn)行插值計(jì)算得到， 其中兩個(gè)節(jié)點(diǎn)位于界面緊鄰的兩側(cè)，第三個(gè)節(jié)點(diǎn)位于迎風(fēng)一側(cè)較遠(yuǎn)鄰點(diǎn)，如圖 F.4所示。當(dāng)ue 0,uw 0時(shí)，通過(guò)WW、W和P三個(gè)節(jié)點(diǎn)值擬合曲線來(lái)計(jì)算主控制單元左側(cè)界面參數(shù) w。通過(guò)節(jié)點(diǎn)W、P和E三個(gè)節(jié)點(diǎn)值擬合曲線來(lái)計(jì)算主控制單元右側(cè)界面參數(shù) e。當(dāng)ue 0,uw 0，則分別通過(guò)節(jié)點(diǎn)W、P、E和P、E、EE三個(gè)節(jié)點(diǎn)值計(jì)算主控制單元左、右兩側(cè)界面參數(shù) w和e。根據(jù)上述計(jì)算原則，可以得到界面參數(shù) w計(jì)算公式如下：當(dāng)uw0時(shí)，界面參數(shù)w計(jì)算公式為：636（F.8a）ww8wW8wP8wWW當(dāng)ue0時(shí)，界面參數(shù)w計(jì)算公式為：636（F.8b）e8P8E8W對(duì)于一維無(wú)源項(xiàng)一維對(duì)流擴(kuò)散方程三階迎風(fēng)型QUICK離散格式：當(dāng)ue0,uw0時(shí)，三階迎風(fēng)型QUICK離散格式為：標(biāo)準(zhǔn)實(shí)用文案aP P aWWaEEaWWWW其中aWDw3Fw1Fe88aD3FEe8eaWW1Fw8aPaWaEaWWFeFw同理，若ue0,uw0，三階迎風(fēng)型QUICK離散格式為：aPPaWWaEEaEEEE其中aWDw3Fw,8aEDe6Fe1Fw,88aWW1Fe8aPaWaEaEEFeFw

F.9）F.9a）F.10）F.10a）將兩種流動(dòng)方向離散方程（ F.9）和（F.10）合并后，可得到統(tǒng)一的一維對(duì)流擴(kuò)散方程三階迎風(fēng)型 QUICK離散格式：aPPaWWaEEaEEEEaEEEE（F.11）其中aWDw6wFw1eFe31wFw888aWW1wFw8aEDe3eFe61eFe11wFw（F.11a）888aEE11eFe8aPaWaEaWWaEEFeFw標(biāo)準(zhǔn)實(shí)用文案式中當(dāng)Fw時(shí)，w；01當(dāng)Fw時(shí)，w；00（F.11b）當(dāng)Fe時(shí)，e；01當(dāng)Fe時(shí)，e。00同理，可以得到帶有源項(xiàng)的二維對(duì)流擴(kuò)散方程三階迎風(fēng)型QUICK離散格式為：aPPaWWaEEaEEEEaEEEE（F.12）aSSaNNaSSSSaNNNNSV其中S為有限體積算法中源項(xiàng)平均值。式中各個(gè)系數(shù)為：aWDw6wFw1eFe31wFw888aWW1wFw8361aEDe8eFe81eFe81wFw1aEE1eFe86131aSDssFsnFnsFs888aSS1sFs8aNDn3nFn61nFn11sFs888aNN11nFn8aPaWaEaWWaEEaSaNaSSaNNFeFwFnFs（F.12a）式中當(dāng)Fk0時(shí)，當(dāng)Fk0時(shí)，

kk

1；0。 （F.12b）k w,e,s,n源項(xiàng)S為：up）S（F.13tx標(biāo)準(zhǔn)實(shí)用文案nn1S離散若把u表示tn時(shí)刻動(dòng)量，u表示tn1時(shí)刻動(dòng)量，則可以得到源項(xiàng)格式為：n1nSdVuPuPxypepwy（F.14）tV最后，得到有限體積算法二維對(duì)流擴(kuò)散方程三階迎風(fēng)型QUICK離散格式：aPuPnaWuWnaEuEnaSuSnaNuNnn1n（F.15）uPuPxynnypepwt式中系數(shù)ak為一階迎風(fēng)格式中各對(duì)應(yīng)系數(shù)。5.計(jì)算結(jié)果分析利用三階迎風(fēng)型QUICK離散格式和相應(yīng)的初始條件和邊界條件，求解二維不可壓縮黏性流體方腔流動(dòng)問(wèn)題。圖 F.5是不同雷諾數(shù)Re條件下采用三階迎風(fēng)型QUICK離散格式得到的二維不可壓縮黏性流體方腔流動(dòng)的計(jì)算結(jié)果。計(jì)算結(jié)果和文獻(xiàn)中其他高精度算法得到的計(jì)算結(jié)果進(jìn)行了比較，兩者計(jì)算結(jié)果十分吻合，能把方腔下壁面兩個(gè)底角附近二次小渦清晰地計(jì)算出來(lái)。 這表明有限體積算法三階迎風(fēng)型 QUICK離散格式具有相當(dāng)高的計(jì)算精度。Re=100 Re=1000標(biāo)準(zhǔn)實(shí)用文案Re=5000 Re=10000圖F.5不同雷諾數(shù) Re條件下采用三階迎風(fēng)型 QUICK離散格式計(jì)算二維不可壓縮黏性方腔流動(dòng)的計(jì)算結(jié)果從圖F.5中可以看出：二維不可壓縮黏性流體方腔流動(dòng)的中心大渦并不在中心位置，方腔內(nèi)流動(dòng)也并不對(duì)稱(chēng)。這是因?yàn)?，方腔上壁面以量綱為一的速度 1.0沿著上壁面方向自左向右運(yùn)動(dòng)時(shí)，在方腔上壁面兩側(cè)的兩個(gè)頂角處不再滿足邊界條件，這是一個(gè)帶奇性的方腔流動(dòng)。計(jì)算結(jié)果表明，方腔流動(dòng)和雷諾數(shù)有關(guān)，隨著雷諾數(shù) Re增加，計(jì)算精度在降低。當(dāng)雷諾數(shù) Re較低時(shí)，方腔下壁面的兩個(gè)底角附近的二次小渦十分清晰，隨著雷諾數(shù)Re的增加，二次小渦變得越來(lái)越模糊。由于三階迎風(fēng)型 QUICK離散格式計(jì)算精度較高，因此三階迎風(fēng)型 QUICK離散格式計(jì)算效果一般要比一階迎風(fēng)型離散格式相對(duì)來(lái)說(shuō)好些。標(biāo)準(zhǔn)實(shí)用文案F-2 二維不可壓縮黏性方腔流動(dòng)問(wèn)題數(shù)值計(jì)算源程序C語(yǔ)言源程序//fvm_quick_cavity.cpp/*----------------------------------------------------------------------------!利用有限體積算法三階迎風(fēng)型 QUICK離散格式和!人工壓縮算法求解方腔流動(dòng)問(wèn)題（ C語(yǔ)言版本）-------------------------------------------------------------------------------*/#include<stdio.h>#include<stdlib.h>#include<math.h>#defineMX100#defineMY100#defineRe100.00#definedt0.0005#definec22.25Doubleu[MX+1][MY+2],v[MX+2][MY+1],p[MX+2][MY+2],un[MX+1][MY+2],vn[MX+2][MY+1],pn[MX+2][MY+2];標(biāo)準(zhǔn)實(shí)用文案全局變量/*-------------------------------------------------------定義兩個(gè)要用到的函數(shù)--------------------------------------------------------*/doublemax(doublea,doubleb){if(a<b)returnb;elsereturna;}doublealfa(doublex){if(x>=0)return1.0;elsereturn0.0;}/*------------------------------------------------------------------------標(biāo)準(zhǔn)實(shí)用文案初始化入口:無(wú)；出口：u、v、p、dx、dy，初始速度、壓強(qiáng)和空間步長(zhǎng)。---------------------------------------------------------------------------*/voidinit(doubleu[MX+1][MY+2],doublev[MX+2][MY+1],doublep[MX+2][MY+2],double&dx,double&dy){inti,j;dx=1.0/MX;dy=1.0/MY;for(i=0;i<=MX;i++){for(j=0;j<=MY+1;j++){u[i][j]=0.0;if(j==MY+1)u[i][j]=4.0/3.0;if(j==MY)u[i][j]=2.0/3.0;}}for(i=0;i<=MX+1;i++)for(j=0;j<=MY;j++)v[i][j]=0.0;標(biāo)準(zhǔn)實(shí)用文案for(i=0;i<=MX+1;i++)for(j=0;j<=MY+1;j++)p[i][j]=1.0;}/*-----------------------------------------------------------------------------------------------一階迎風(fēng)型離散格式二維的三階迎風(fēng)型 QUICK離散格式為 9點(diǎn)格式，因此有兩層邊界網(wǎng)格需要處理，本程序采用一階迎風(fēng)型離散格式處理內(nèi)層，用物理邊界條件處理外層。入口:u、v、p、dx、dy、i、j，當(dāng)前速度、壓強(qiáng)，空間步長(zhǎng)和網(wǎng)格節(jié)點(diǎn)編號(hào)；出口：un，新的x方向速度。-----------------------------------------------------------------------------------------------*/void upwind_u(double u[MX+1][MY+2],double v[MX+2][MY+1],doublep[MX+2][MY+2],doubleun[MX+1][MY+2],doubledx,doubledy,inti,intj){doubleaw,ae,as,an,df,ap,miu;miu=1.0/Re;aw=miu+max(0.5*(u[i-1][ j]+u[i][j])*dy,0.0);標(biāo)準(zhǔn)實(shí)用文案ae=miu+max(0.0,-0.5*(u[i][ j]+u[i+1][j])*dy);as=miu+max(0.5*(v[i][ j-1]+v[i+1][ j-1])*dx,0.0);an=miu+max(0.0,-0.5*(v[i][j]+v[i+1][ j])*dx);df=0.5*(u[i+1][j]-u[i-1][j])*dy+0.5*(v[i][ j]+v[i+1][ j]-v[i][j-1]-v[i+1][j-1])*dx;ap=aw+ae+as+an+df;un[i][ j]=u[i][ j] +dt/dx/dy*(-ap*u[i][j]+aw*u[i-1][j]+ae*u[i+1][j]+as*u[i][ j-1]+an*u[i][j+1])-dt*(p[i+1][j]-p[i][j])/dx;}/*------------------------------------------------------------------------------------------------入口:u、v、p、dx、dy、i、j，當(dāng)前速度、壓強(qiáng)，空間步長(zhǎng)和網(wǎng)格節(jié)點(diǎn)編號(hào)；出口：vn，新的y方向速度。-----------------------------------------------------------------------------------------------------*/void upwind_v(double u[MX+1][MY+2],double v[MX+2][MY+1],doublep[MX+2][MY+2],doublevn[MX+2][MY+1],doubledx,doubledy,inti,intj){doubleaw,ae,as,an,df,ap,miu;miu=1.0/Re;aw=miu+max(0.5*(u[i-1][ j]+u[i-1][j+1])*dy,0.0);標(biāo)準(zhǔn)實(shí)用文案ae=miu+max(0.0,-0.5*(u[i][j]+u[i][ j+1])*dy);as=miu+max(0.5*(v[i][ j-1]+v[i][ j])*dx,0.0);an=miu+max(0.0,-0.5*(v[i][j]+v[i][ j+1])*dx);df=0.5*(u[i][ j]+u[i][ j+1]-u[i-1][ j]-u[i-1][j+1])*dy+0.5*(v[i][ j+1]-v[i][j-1])*dx;ap=aw+ae+as+an+df;vn[i][j]=v[i][ j] +dt/dx/dy*(-ap*v[i][j]+aw*v[i-1][ j]+ae*v[i+1][j]+as*v[i][ j-1]+an*v[i][ j+1])-dt*(p[i][j+1]-p[i][j])/dy;}/*----------------------------------------------------------------------三階迎風(fēng)型QUICK離散格式入口:u、v、p、dx、dy，當(dāng)前速度、壓強(qiáng)，空間步長(zhǎng)；出口：un、vn，新的速度。-----------------------------------------------------------------------*/void quick(double u[MX+1][MY+2],double v[MX+2][MY+1],doublep[MX+2][MY+2],doubleun[MX+1][MY+2],doublevn[MX+2][MY+1],doubledx,doubledy){doublemiu,fw,fe,fs,fn,df,aw,ae,as,an,aww,aee,ass,ann,ap;inti,j;miu=1.0/Re;標(biāo)準(zhǔn)實(shí)用文案for(i=2;i<=MX-2;i++){for(j=2;j<=MY-1;j++){fw=0.5*(u[i-1][j]+u[i][ j])*dy;fe=0.5*(u[i][j]+u[i+1][ j])*dy;fs=0.5*(v[i][j-1]+v[i+1][j-1])*dx;fn=0.5*(v[i][j]+v[i+1][j])*dx;df=fe-fw+fn-fs;aw=miu+0.750*alfa(fw)*fw+0.125*alfa(fe)*fe+0.375*(1.0-alfa(fw))*fw;aww=-0.125*alfa(fw)*fw;ae=miu-0.375*alfa(fe)*fe-0.750*(1.0-alfa(fe))*fe-0.125*(1.0-alfa(fw))*fw;aee=0.125*(1.0-alfa(fe))*fe;as=miu+0.750*alfa(fs)*fs+0.125*alfa(fn)*fn+0.375*(1.0-alfa(fs))*fs;ass=-0.125*alfa(fs)*fs;an=miu-0.375*alfa(fn)*fn-0.750*(1.0-alfa(fn))*fn-0.125*(1.0-alfa(fs))*fs;ann=0.125*(1.0-alfa(fn))*fn;ap=aw+ae+as+an+aww+aee+ass+ann+df;//aw、ae、as、an...均為有限體積算法中各項(xiàng)系數(shù)，詳見(jiàn)前文三階迎風(fēng)型 QUICK離散格式。un[i][j]=u[i][j]+標(biāo)準(zhǔn)實(shí)用文案dt/dx/dy*(-ap*u[i][ j]+aw*u[i-1][j]+ae*u[i+1][j]+as*u[i][ j-1]+an*u[i][ j+1]+aww*u[i-2][j]+aee*u[i+2][ j]+ass*u[i][ j-2]+ann*u[i][j+2])-dt*(p[i+1][j]-p[i][j])/dx;}}//------------------------------------------------------------------------------j=1;for(i=2;i<=MX-2;i++)upwind_u(u,v,p,un,dx,dy,i,j);j=MY;for(i=2;i<=MX-2;i++)upwind_u(u,v,p,un,dx,dy,i,j);i=1;for(j=1;j<=MY;j++)upwind_u(u,v,p,un,dx,dy,i,j);i=MX-1;for(j=1;j<=MY;j++)upwind_u(u,v,p,un,dx,dy,i,j);//內(nèi)層邊界由一階迎風(fēng)型離散格式得到 -----------------------------------------//--------------------------------------------------------------------------------for(i=1;i<=MX-1;i++)標(biāo)準(zhǔn)實(shí)用文案{un[i][0]=-un[i][1];un[i][MY+1]=2.0-un[i][MY];}for(j=0;j<=MY+1;j++){un[0][j]=0.0;un[MX][j]=0.0;}外層邊界條件按物理邊界條件給出//-------------------------------------------------------------------------------for(i=2;i<=MX-1;i++){for(j=2;j<=MY-2;j++){fw=0.5*(u[i-1][j]+u[i-1][j+1])*dy;fe=0.5*(u[i][j]+u[i][j+1])*dy;fs=0.5*(v[i][j-1]+v[i][ j])*dx;fn=0.5*(v[i][j]+v[i][ j+1])*dx;df=fe-fw+fn-fs;aw=miu+0.750*alfa(fw)*fw+0.125*alfa(fe)*fe+0.375*(1.0-alfa(fw))*fw;aww=-0.125*alfa(fw)*fw;標(biāo)準(zhǔn)實(shí)用文案ae=miu-0.375*alfa(fe)*fe-0.750*(1.0-alfa(fe))*fe-0.125*(1.0-alfa(fw))*fw;aee=0.125*(1.0-alfa(fe))*fe;as=miu+0.750*alfa(fs)*fs+0.125*alfa(fn)*fn+0.375*(1.0-alfa(fs))*fs;ass=-0.125*alfa(fs)*fs;an=miu-0.375*alfa(fn)*fn-0.750*(1.0-alfa(fn))*fn-0.125*(1.0-alfa(fs))*fs;ann=0.125*(1.0-alfa(fn))*fn;ap=aw+ae+as+an+aww+aee+ass+ann+df;vn[i][j]=v[i][j] +dt/dx/dy*(-ap*v[i][j]+aw*v[i-1][ j]+ae*v[i+1][j]+as*v[i][ j-1]+an*v[i][j+1]+aww*v[i-2][j]+aee*v[i+2][j]+ass*v[i][ j-2]+ann*v[i][ j+2])-dt*(p[i][j+1]-p[i][ j])/dy;}}//-----------------------------------------------------------------------------j=1;for(i=2;i<=MX-1;i++)upwind_v(u,v,p,vn,dx,dy,i,j);j=MY-1;for(i=2;i<=MX-1;i++)upwind_v(u,v,p,vn,dx,dy,i,j);標(biāo)準(zhǔn)實(shí)用文案i=1;for(j=1;j<=MY-1;j++)upwind_v(u,v,p,vn,dx,dy,i,j);i=MX;for(j=1;j<=MY-1;j++)upwind_v(u,v,p,vn,dx,dy,i,j);//----------------------------------------------------------------------------for(i=1;i<=MX;i++){vn[i][0]=0.0;vn[i][MY]=0.0;}for(j=0;j<=MY;j++){vn[0][j]=-vn[1][j];vn[MX+1][j]=-vn[MX][j];}}//----------------------------------------------------------------------------/*----------------------------------------------------------------------------標(biāo)準(zhǔn)實(shí)用文案更新壓強(qiáng)入口:un、vn、p、dx、dy，新的速度，當(dāng)前壓強(qiáng)，空間步長(zhǎng)；出口:pn，新的壓強(qiáng)。-----------------------------------------------------------------------------*/void getp(double un[MX+1][MY+2],double vn[MX+2][MY+1],doublep[MX+2][MY+2],doublepn[MX+2][MY+2],doubledx,doubledy){inti,j;for(i=1;i<=MX;i++)for(j=1;j<=MY;j++)pn[i][ j]=p[i][ j]-dt*c2/dx*(un[i][ j]-un[i-1][ j]+vn[i][ j]-vn[i][ j-1]);for(i=1;i<=MX;i++){pn[i][0]=pn[i][1];pn[i][MY+1]=pn[i][MY];}for(j=0;j<=MY+1;j++){pn[0][ j]=pn[1][j];pn[MX+1][ j]=pn[MX][j];}標(biāo)準(zhǔn)實(shí)用文案}/*------------------------------------------------------主程序-------------------------------------------------------*/voidmain(){doubledx,dy,err,value,uo,vo,po;inti,j,step;init(u,v,p,dx,dy);err=100.0;step=0;while(err>1e-4&&step<1e6)//err<1e-5 ，定常解判據(jù)； step，限制迭代次數(shù){printf("step=%d\n",step);step++;err=0.0;quick(u,v,p,un,vn,dx,dy);getp(un,vn,p,pn,dx,dy);//-------------------------------------------------------for(i=0;i<=MX;i++){標(biāo)準(zhǔn)實(shí)用文案for(j=0;j<=MY+1;j++){value=fabs(un[i][j]-u[i][ j])/dt;if(value>err) err=value;u[i][j]=un[i][j];}}for(i=0;i<=MX+1;i++){for(j=0;j<=MY;j++){value=fabs(vn[i][ j]-v[i][j])/dt;if(value>err) err=value;v[i][j]=vn[i][j];}}for(i=0;i<=MX+1;i++){for(j=0;j<=MY+1;j++){value=fabs(pn[i][j]-p[i][ j])/c2/dt;if(value>err) err=value;標(biāo)準(zhǔn)實(shí)用文案p[i][j]=pn[i][ j];}}printf("err=%le\n",err);//--------------------------------------------------------}輸出結(jié)果，用Tecplot數(shù)據(jù)格式畫(huà)圖FILE*fp;fp=fopen("Result.plt","w");fprintf(fp,"variables=x,y,u,v,p\nzonei=%d,j=%d,f=point\n",MX+1,MY+1);for(i=0;i<=MX;i++){for(j=0;j<=MY;j++){uo=0.5*(u[i][j]+u[i][j+1]);vo=0.5*(v[i][j]+v[i+1][j]);po=0.25*(p[i][j]+p[i+1][j]+p[i][j+1]+p[i+1][j+1]);fprintf(fp,"%20.10e%20.10e%20.10e%20.10e%20.10e\n",i*dx,j*dy,uo,vo,po);}}fclose(fp);標(biāo)準(zhǔn)實(shí)用文案----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------Fortran77語(yǔ)言源程序————————————————————————————————————!利用有限體積算法三階迎風(fēng)型 QUICK離散格式和人工壓縮算法求解方腔流動(dòng)問(wèn)題（Fortran77語(yǔ)言版本）————————————————————————————————————programQUICK_cavityparameter(mx=101,my=101)implicitdoubleprecision(a-h,o-z)dimensionu(mx,my+1),v(mx+1,my),p(mx+1,my+1)dimensionun(mx,my+1),vn(mx+1,my),pn(mx+1,my+1)common/ini/u,v,pc2=2.25標(biāo)準(zhǔn)實(shí)用文案re=100.0dt=0.0005dx=1.0/float(mx-1)dy=1.0/float(my-1)!----------------------------------------------------------------------------------------! u、v、p為t時(shí)刻值，un、vn、pn為t+1 時(shí)刻值，! mx、my為最大網(wǎng)格數(shù)， c2為虛擬壓縮系數(shù)的平方， re為雷諾數(shù)。!-----------------------------------------------------------------------------------------num=0err=100.00!nun，計(jì)數(shù)器；err，判斷人工壓縮法求解收斂的標(biāo)準(zhǔn)callinitial調(diào)入初始條件，以下為人工壓縮算法求解err=0.0callquick(u,v,p,un,vn,mx,my,dx,dy,dt,re)!QUICK離散格式求解動(dòng)量方程，得到 un、vncallcalp(p,un,vn,pn,mx,my,dx,dy,dt,c2)求壓強(qiáng)pncallcheck(u,v,p,un,vn,pn,mx,my,dx,dy,dt,c2,err)校驗(yàn)流場(chǎng)信息，判斷是否收斂，同時(shí)更新u、v、pwrite(*,*)'error=',err標(biāo)準(zhǔn)實(shí)用文案num=num+1write(*,*)num屏幕跟蹤輸出enddocalloutput(u,v,p,mx,my,dx,dy)輸出結(jié)果文件End!!subroutineinitial初始化流場(chǎng)parameter(mx=101,my=101)doubleprecisionu(mx,my+1),v(mx+1,my),p(mx+1,my+1)common/ini/u,v,pdoi=1,mx+1doj=1,my+1p(i,j)=1.0enddoenddodoi=1,mxdoj=1,my+1u(i,j)=0.0標(biāo)準(zhǔn)實(shí)用文案enddoenddodoi=1,mx+1doj=1,myv(i,j)=0.0enddoenddoendsubroutine!!subroutinequick(u,v,p,un,vn,mx,my,dx,dy,dt,re)以QUICK格式離散動(dòng)量方程implicitdoubleprecision(a-h,o-z)dimensionu(mx,my+1),v(mx+1,my),p(mx+1,my+1),un(mx,my+1),vn(mx+1,my)doubleprecisionmiumiu=1.0/re! 以 下 求 解 x 方 向 速 度un----------------------------------------------------------------------------doi=3,mx-2標(biāo)準(zhǔn)實(shí)用文案doj=3,my-1fw=0.5*(u(i-1,j)+u(i,j))*dyfe=0.5*(u(i,j)+u(i+1,j))*dyfs=0.5*(v(i,j-1)+v(i+1,j-1))*dxfn=0.5*(v(i,j)+v(i+1,j))*dxdf=fe-fw+fn-fsaw=miu+0.750*alfa(fw)*fw+0.125*alfa(fe)*fe+0.375*(1.0-alfa(fw))*fwaww=-0.125*alfa(fw)*fwae=miu-0.375*alfa(fe)*fe-0.750*(1.0-alfa(fe))*fe-0.125*(1.0-alfa(fw))*fwaee=0.125*(1.0-alfa(fe))*feas=miu+0.750*alfa(fs)*fs+0.125*alfa(fn)*fn+0.375*(1.0-alfa(fs))*fsass=-0.125*alfa(fs)*fsan=miu-0.375*alfa(fn)*fn-0.750*(1.0-alfa(fn))*fn-0.125*(1.0-alfa(fs))*fsann=0.125*(1.0-alfa(fn))*fnap=aw+ae+as+an+aww+aee+ass+ann+df!aw、ae、as、an...均為有限體積算法中各項(xiàng)系數(shù)，詳見(jiàn)前文三階迎風(fēng)型 QUICK離散格式un(i,j)=u(i,j)+dt/dx/dy*(-ap*u(i,j)+aw*u(i-1,j)+ae*u(i+1,j) &+as*u(i,j-1)+an*u(i,j+1)+aww*u(i-2,j)+aee*u(i+2,j) &+ass*u(i,j-2)+ann*u(i,j+2))-dt*(p(i+1,j)-p(i,j))/dxenddoenddo標(biāo)準(zhǔn)實(shí)用文案!-------------------------------------------------------------------------------------------j=2doi=3,mx-2callupbound_u(u,v,p,un,mx,my,dx,dy,dt,re,i,j)enddoj=mydoi=3,mx-2callupbound_u(u,v,p,un,mx,my,dx,dy,dt,re,i,j)enddoi=2doj=2,mycallupbound_u(u,v,p,un,mx,my,dx,dy,dt,re,i,j)enddoi=mx-1doj=2,mycallupbound_u(u,v,p,un,mx,my,dx,dy,dt,re,i,j)enddo!內(nèi)層邊界由一階迎風(fēng)型離散格式得到 ----------------------------------------------------!-------------------------------------------------------------------------------------------doi=2,mx-1標(biāo)準(zhǔn)實(shí)用文案un(i,1)=-un(i,2)un(i,my+1)=2.0-un(i,my)enddodoj=1,my+1un(1,j)=0.0un(mx,j)=0.0enddo外層邊界條件按物理邊界條件給出!-----------------------------------------------------------------------------以下求解y方向速度vn-----------------------------------------------doi=3,mx-1doj=3,my-2fw=0.5*(u(i-1,j)+u(i-1,j+1))*dyfe=0.5*(u(i,j)+u(i,j+1))*dyfs=0.5*(v(i,j-1)+v(i,j))*dxfn=0.5*(v(i,j)+v(i,j+1))*dxdf=fe-fw+fn-fsaw=miu+0.750*alfa(fw)*fw+0.125*alfa(fe)*fe+0.375*(1.0-alfa(fw))*fwaww=-0.125*alfa(fw)*fwae=miu-0.375*alfa(fe)*fe-0.750*(1.0-alfa(fe))*fe-0.125*(1.0-alfa(fw))*fwaee=0.125*(1.0-alfa(fe))*fe標(biāo)準(zhǔn)實(shí)用文案as=miu+0.750*alfa(fs)*fs+0.125*alfa(fn)*fn+0.375*(1.0-alfa(fs))*fsass=-0.125*alfa(fs)*fsan=miu-0.375*alfa(fn)*fn-0.750*(1.0-alfa(fn))*fn-0.125*(1.0-alfa(fs))*fsann=0.125*(1.0-alfa(fn))*fnap=aw+ae+as+an+aww+aee+ass+ann+dfvn(i,j)=v(i,j)+dt/dx/dy*(-ap*v(i,j)+aw*v(i-1,j)+ae*v(i+1,j) &+as*v(i,j-1)+an*v(i,j+1)+aww*v(i-2,j)+aee*v(i+2,j) &+ass*v(i,j-2)+ann*v(i,j+2))-dt*(p(i,j+1)-p(i,j))/dyenddoenddo!-----------------------------------------------------------------------------j=2doi=3,mx-1callupbound_v(u,v,p,vn,mx,my,dx,dy,dt,re,i,j)enddoj=my-1doi=3,mx-1callupbound_v(u,v,p,vn,mx,my,dx,dy,dt,re,i,j)enddoi=2doj=2,my-1callupbound_v(u,v,p,vn,mx,my,dx,dy,dt,re,i,j)標(biāo)準(zhǔn)實(shí)用文案enddoi=mxdoj=2,my-1callupbound_v(u,v,p,vn,mx,my,dx,dy,dt,re,i,j)enddo!----------------------------------------------------------------------------doi=2,mxvn(i,1)=0.0vn(i,my)=0.0enddodoj=1,myvn(1,j)=-vn(2,j)vn(mx+1,j)=-vn(mx,j)enddo!----------------------------------------------------------------------------Endsubroutine!!functionalfa(x)函數(shù),正1負(fù)0doubleprecisionalfa,x標(biāo)準(zhǔn)實(shí)用文案alfa=1.0elsealfa=0.0endifend!!subroutineupbound_u(u,v,p,un,mx,my,dx,dy,dt,re,i,j)以一階迎風(fēng)型離散格式得到內(nèi)層邊界un的值implicitdoubleprecision(a-h,o-z)dimensionu(mx,my+1),v(mx+1,my),p(mx+1,my+1),un(mx,my+1)doubleprecisionmiumiu=1.0/reaw=miu+max(0.5*(u(i-1,j)+u(i,j))*dy,0.0)ae=miu+max(0.0,-0.5*(u(i,j)+u(i+1,j))*dy)as=miu+max(0.5*(v(i,j-1)+v(i+1,j-1))*dx,0.0)an=miu+max(0.0,-0.5*(v(i,j)+v(i+1,j))*dx)df=0.5*(u(i+1,j)-u(i-1,j))*dy+0.5*(v(i,j)+v(i+1,j)-v(i,j-1)-v(i+1,j-1))*dxap=aw+ae+as+an+dfun(i,j)=u(i,j)+dt/dx/dy*(-ap*u(i,j)+aw*u(i-1,j)+ae*u(i+1,j) &+as*u(i,j-1)+an*u(i,j+1))-dt*(p(i+1,j)-p(i,j))/dxEndsubroutine標(biāo)準(zhǔn)實(shí)用文案!!subroutineupbound_v(u,v,p,vn,mx,my,dx,dy,dt,re,i,j)!以一階迎風(fēng)型離散格式得到內(nèi)層邊界 vn值implicitdoubleprecision(a-h,o-z)dimensionu(mx,my+1),v(mx+1,my),p(mx+1,my+1),vn(mx+1,my)doubleprecisionmiumiu=1.0/reaw=miu+max(0.5*(u(i-1,j)+u(i-1,j+1))*dy,0.0)ae=miu+max(0.0,-0.5*(u(i