Fortran语言-涡量流函数法中心差分格式的二维方腔顶盖驱动计算

涡量流函数法中心差分格式的二维方腔顶盖驱动计算本题是关于粘性流体方腔顶盖驱动的问题。采用涡量流函数法，在均分网格下，用中心差分格式进行计算，结果与文献中所采用的其他方法和格式进行比较，认为中心差分格式符合计算的精度，但是其明显缺点是计算过程不稳定。1、参数的无量纲化令X=其中，由涡量的定义式：QUOTEX=xH;Y=yH;U=u将速度的导数项无量纲化，并将上述定义式代入，即ω=令utHΩ=再由流函数的定义式：u=∂ψ/∂y同理，令ψ可以得到流函数的无量纲表达式为：Ψ=流函数边界条件X,Y∈ABΨ=0X,YX,Y∈CDΨ=0

涡量的边界条件（Thom公式）X,YX,YX,YX,Y∈ADΩ2、方程的无量纲处理过程2.1、非守恒型方程的处理（1）将以上假设的各式代入到非守恒型方程组：得到以下无量纲形式的方程：①涡量控制方程∂②流函数控制方程∂（2）采用有限差分法离散非守恒型涡量的无量纲控制方程：R整理得到：2+简化为：4+写成一般形式是:a其中ap=4；aWaS=2.2、守恒型方程的处理（1）将假设的各个无量纲量代入守恒型涡量方程得到其无量纲形式为：Re[（2）采用有限容积法离散守恒型涡量的无量纲控制方程：对于扩散项，有s==[对于对流项，有s对于界面上的取值，采用如下形式：Ue≥0时，Uw≥0时，Vn≥0时，ΩVs≥0时，Ωs其中，界面上涡量的插值采用中心差分格式，令fe+=ffn+=所以，Ωe=ΩP+关于界面上速度的插值，可以用上一层次的流函数来表示，形式如下：Ue=ΨN+ΨVn=Ψ由以上各式，整理得到守恒型涡量控制方程无量纲形式的离散方程为：R+R=综上，可整理得到（同位的均分网格）：41-+写成一般形式为：b其中，bbbbb2.3、对于流函数方程，离散可得到：2整理并写成一般形式是：c其中，cp=4；cW=1；cE并且，采用有限容积法离散的结果与采用这里的方法得到的结果是相同的。2.4、将压力的泊松方程进行无量纲化并离散令P代入压力泊松方程，得到其无量纲形式为：∂边界上的压力满足N-S方程，同样将其无量纲化，并与上式组成压力的封闭方程组，X方向上：∂PY方向上：∂P将以上三式所组成的方程组离散，得到：内点压力计算式为：2-左边界压力为：1+(右边界压力为：1-(上边界压力为：1下边界压力为：1+(由此可得到压力的封闭方程组，通过求得的速度场来计算压力场。3、程序编写思路：11、定义计算区域，并进行网格划分，得到网格步长，并且设定初场；2、更新涡量方程的系数（第一次更新时使用流函数Ψ0解涡量方程得到Ω用ΩP用新得到的流函数更新涡量的边界值判断收敛与否否是输出结果，进行比较 4、计算结果4.1、两种离散方式的计算时间的比较：（EPS=1E-6,i代表松弛因子）RE有限容积法与有限差分法计算时间的比较/sRE网格数网格数100(i=0.3)1000(i=0.1)有限容积法FVM有限差分法FDM有限容积法FVM有限差分法FDM20×2037.0287.481973.202132.5640×40509.17517.584415.804713.2080×806889.286906.0914359.8015536.744.2、雷诺数为100时各物理量场的分布图形（1）有限差分法的情况下：A流函数等值线（图1a/b/c）N=20×20N=40×40N=80×80B涡量等值线（图2a/b/c）N=20×20N=40×40N=80×80C压力等值线（图3a/b/c）N=20×20N=40×40N=80×80D中心线的速度分布（图4a/b）（2）有限容积法的情况下：A流函数等值线（图5a/b/c）N=20×20N=40×40N=80×80B涡量等值线（图6a/b/c）N=20×20N=40×40N=80×80C压力等值线（图7a/b/c）N=20×20N=40×40N=80×80D中心线的速度分布（图8a/b）4.3、雷诺数为1000时各物理量场的分布图形（1）有限差分法的情况下：A流函数等值线（图9a/b/c）N=20×20N=40×40N=80×80B涡量等值线（图10a/b/c）N=20×20N=40×40N=80×80C压力等值线（图11a/b/c）N=20×20N=40×40N=80×80D中心线的速度分布（图12a/b/c）（2）有限容积法的情况下：A流函数等值线（图13a/b/c）N=20×20N=40×40N=80×80B涡量等值线（图14a/b/c）N=20×20N=40×40N=80×80C压力等值线（图15a/b/c）N=20×20N=40×40N=80×80D中心线的速度分布（图16a/b/c）4.4、文献[1]上的结果（1）雷诺数100时候的流函数和涡量图形（图17a/b）（2）雷诺数1000时的流函数和涡量图形（图18a/b）5、讨论与小结在图4a/b、图8a/b、图12a/b/c、图16a/b/c中，将计算方腔的中心线速度分布与U.GHIA等人用CSI-MG方法在129×129网格数情况下得到的结果[1]相比较，可以看出本题的计算结果与基准解符合得比较好，说明计算方法具有合适的精度。但是，本题在整个计算过程中，所遇到的关于稳定性及收敛性方面的问题较多，正如文献[9]所讲，在涡量流函数的计算中，由于涡量控制方程是对流占优扩散的，使用中心差分的离散方法不能得到很好的计算结果。在低雷诺数RE=100时，计算很容易就得到符合基准解的结果。从中心线上的速度分布可以看到，结果与基准解吻合。但是，在较高雷诺数1000时，由于中心差分格式是条件稳定[10]，这造成了计算的不稳定，并且非常容易发散。关于涡量的上边界条件，若按照文献[9]上的形式:Ω=2/Y来做，则得到的速度会出现越界，但是这时速度的图形却更接近高雷诺数时候的基准解，且涡量和流函数也更加接近高雷诺数时候的基准解；若采用题目的边界条件，则得出的速度均有物理意义，且涡量及流函数的图形比较符合低雷诺数的解。这说明涡量边界条件的确定对于计算的准确性有很大的影响。（如图19a/b）a.采用文献上的涡量边界在b.采用文献上的涡量边界在RE=100时候的涡量等值线图RE=100时候的流函数等值线图本题在计算雷诺数RE=1000的时候，过程是明显不稳定的。将松弛因子取得很小（0.1），造成收敛性变差，计算很费时。原因之一是采用高斯赛德尔迭代，且内迭代次数取200次，而外迭代次数在结果收敛时（80×80网格）竟达到115832次之多。本题程序只适合用于计算雷诺数较低且网格数较密的情况。参考文献[1]UGhia,KNGhiaandCTShin.High-ResolutionsforincompressibleflowusingtheNavier-Stokesequationsandamultigridmethod.JournalofComputationalPhysics.1982,48:387-411[2]E.Erturk1.Numericalsolutionsof2-DsteadyincompressibledrivencavityflowathighReynoldsnumbers.JournalforNumericalMethodsinFluids.2005,48:747–774[3]Yih-FerngPeng,Yuo-HsienShiau,RobertR.Hwang.Transitionina2-Dlid-drivencavityflow.Computers&Fluids.2003;32:337–352[4]PatankarSV.Numericalheattransferandfluidflow.NewYork:McGraw-Hill;1980.[5]Charles-HenriBruneau,MazenSaad.The2Dlid-drivencavityproblemrevisited.Computers&Fluids.2006;35:326–348[6]李明秀,陶文铨,王秋旺.LATTICE-BOLTZMANN方法及其在顶盖驱动流数值模拟中的应用.工程热物理[J].2001,3.22(2):223-225[7]周高领,陈斌,宇波.不同形状空腔顶盖驱动流的数值模拟.中国工程热物理学会学术会议论文.[8]陈庆光.张永建.钱宝光.QUICK和乘方格式在顶盖驱动方腔流动数值计算中的比较[J]-山东科技大学学报(自然科学版)2002,3.21(1):8-11[9]吴晓冬,陈立亮.流函数-涡量法的二维方腔流数值模拟.中国铸造装备与技术[J].2007,1:36-38[10]王贤钢,俞茂铮,陶文铨.一种稳定性与精度协调一致的二阶差分格式.西安交通大学学报.2005,11.39(11):1194-1198源程序：（1）有限差分法源程序：!THISPROGRAMISWRITTENTOSOLVEA2-DDRIVENFLOWPROBLEMINASQUARE!CAVITYUSINGTHEVORTICITY-STREAMFUNCTIONMETHOD.!L1X方向的节点数!M1y方向的节点数!DXX方向的步长!DYY方向的步长!PHAI流函数!OMEGA涡量!RE雷诺数!ITS迭代次数!ERR1，ERR2前后迭代步各点涡量之差PROGRAMMAINIMPLICITNONEREALDX,DYREALERR1,ERR2,ERR3INTEGER,PARAMETER::L1=81,M1=81REAL::RE=100.0INTEGERi,j,KINTEGERITS!计算涡量时的跌代数INTEGERITSS!计算压力时的跌代数INTEGERN!内迭代次数INTEGER::ITSMAX=5000!最大跌代数REAL::EPS=1E-6!迭代精度REAL::ALPHA1=0.3!松弛因子REAL::ALPHA2=0.3REALU(L1,M1),V(L1,M1)!速度分量REALP(L1,M1)!压力分布REALXL(L1,M1),YM(L1,M1)!存放每个点的横坐标和纵坐标REALPHAI(L1,M1)REALOMEGA(L1,M1)REALAW(L1,M1),AE(L1,M1),AS(L1,M1),AN(L1,M1)!涡量方程中的各项系数REALOME_K(L1,M1)!存放每一迭代步得到的涡量值REALPHAI_K(L1,M1)!存放每一迭代步得到的流函数值REALP_K(L1,M1)!存放压力泊松方程迭代时上一步的值DOUBLEPRECISIONTIME_START,TIME_END!*CALLCPU_TIME(TIME_START)CALLMESHING(L1,M1,DX,DY,PHAI,OMEGA,U,V,p,XL,YM)ERR1=1.0ERR2=1.0ITS=0DOWHILE(ERR1>=EPS.OR.ERR2>=EPS.OR.ITS<ITSMAX)DOI=1,L1DOJ=1,M1OME_K(I,J)=OMEGA(I,J)!用OME_K储存前一迭代步的计算结果 PHAI_K(I,J)=PHAI(I,J)!用PHAI_K储存前一迭代步的计算结果ENDDOENDDODOI=2,L1-1!涡量方程系数的更新DOJ=2,M1-1AW(I,J)=1.0+RE*U(I,J)*DX/2.0AE(I,J)=1.0-RE*U(I,J)*DX/2.0AS(I,J)=1.0+RE*V(I,J)*DY/2.0AN(I,J)=1.0-RE*V(I,J)*DY/2.0ENDDOENDDODON=1,200DOJ=M1-1,2,-1!解涡量方程DOI=2,L1-1OMEGA(i,j)=(AW(i,j)*OMEGA(i-1,j)+AE(i,j)*OMEGA(i+1,j)& +AS(i,j)*OMEGA(i,j-1)+AN(i,j)*OMEGA(i,j+1))/4.0OMEGA(I,J)=OME_K(I,J)+ALPHA1*(OMEGA(I,J)-OME_K(I,J))ENDDOENDDOENDDOWRITE(*,*)'OMEGA(3,11)=',OMEGA(3,11)DON=1,200DOJ=M1-1,2,-1!解流函数方程 DOI=2,L1-1PHAI(i,j)=(PHAI(i-1,j)+PHAI(i+1,j)& +PHAI(i,j-1)+PHAI(i,j+1)-OMEGA(i,j)*DX**2)/4.0PHAI(I,J)=PHAI_K(I,J)+ALPHA2*(PHAI(I,J)-PHAI_K(I,J))ENDDOENDDOENDDOWRITE(*,*)'PHAI(3,11)=',PHAI(3,11) DOI=2,L1-1!由流函数计算速度分布 DOJ=2,M1-1 U(I,J)=(PHAI(I,J+1)-PHAI(I,J))/DY+(2.0*PHAI_K(I,J)-PHAI_K(I,J+1)-PHAI_K(I,J-1))/(2.0*DY) V(I,J)=(PHAI(I-1,J)-PHAI(I,J))/DY+(2.0*PHAI_K(I,J)-PHAI_K(I+1,J)-PHAI_K(I-1,J))/(2.0*DX) ENDDO ENDDODOj=1,M1!左、右边界的涡量值OMEGA(1,J)=2.0*(PHAI(2,J)-PHAI(1,J))/DX**2OMEGA(L1,J)=2.0*(PHAI(L1-1,J)-PHAI(L1,J))/DX**2ENDDODOi=2,L1-1!上、下边界的涡量值OMEGA(i,M1)=2.0/DY OMEGA(i,1)=2.0*(PHAI(i,2)-PHAI(i,1))/DY**2ENDDO OMEGA(1,1)=(OMEGA(1,2)+OMEGA(2,1))/2.0 OMEGA(1,M1)=(OMEGA(1,M1-1)+OMEGA(2,M1))/2.0 OMEGA(L1,1)=(OMEGA(L1-1,1)+OMEGA(2,L1))/2.0 OMEGA(L1,M1)=(OMEGA(L1,M1-1)+OMEGA(L1-1,M1))/2.0ITS=ITS+1 ERR1=MAXVAL(ABS(OME_K-OMEGA)) ERR2=MAXVAL(ABS(PHAI_K-PHAI)) WRITE(*,*)'ERR1=',ERR1 WRITE(*,*)'ERR2=',ERR2 WRITE(*,*)'ITS=',ITS ENDDO!* ERR3=1.0ITSS=0!用计算得到的速度求解内点压力场分布DOWHILE(ERR3>=EPS.AND.ITSS<=ITSMAX) DOI=1,L1 DOJ=1,M1 P_K(I,J)=P(I,J) ENDDO ENDDO DOI=2,L1-1 DOJ=2,M1-1 P(I,J)=(P(I-1,J)+P(I+1,J)+P(I,J+1)+P(I,J-1)-0.25*((U(I+1,J)-U(I-1,J))*(V(I,J+1)-V(I,J-1))-&(V(I+1,J)-V(I-1,J))*(U(I,J+1)-U(I,J-1))))/4.0 P(I,J)=P_K(I,J)+0.3*(P(I,J)-P_K(I,J)) ENDDO ENDDO!* DOJ=2,M1-1!压力场的左、右边界值 P(1,J)=4.0/3.0*P(2,J)-1.0/3.0*P(3,J)-(2.0/(3.0*RE))*(-2.0*U(2,J)+U(3,J))/D P(L1,J)=4.0/3.0*P(L1-1,J)-1.0/3.0*P(L1-2,J)+(2.0/(3.0*RE))*(-2.0*U(L1-1,J)+U(L1-2,J))/DX ENDDO DOI=2,L1-1!压力场的上、下边界值 P(I,M1)=(4.0/3.0)*P(I,M1-1)-(1.0/3.0)*P(I,M1-2)+(2.0/(3.0*RE*DY))*(-2.0*V(I,M1-1)+V(I,M1-2)) P(I,1)=4.0/3.0*P(I,M1-1)-1.0/3.0*P(I,M1-2)-(2.0/(3.0*RE))*(-2.0*V(I,2)+V(I,3))/DX ENDDO!角点上的压力进行加权平均 P(1,1)=0.5*P(1,2)+0.5*P(2,1) P(1,M1)=0.5*P(1,M1-1)+0.5*P(2,M1) P(L1,1)=0.5*P(L1,2)+0.5*P(L1-1,1) P(L1,M1)=0.5*P(L1,M1-1)+0.5*P(L1-1,M1)ITSS=ITSS+1 ERR3=MAXVAL(ABS(P_K-P)) ENDDOCALLCPU_TIME(TIME_END) WRITE(*,*)'总共运行时间(S)',TIME_END-TIME_START!* OPEN(1,FILE='PHAI.DAT') OPEN(2,FILE='OMEGA.DAT') OPEN(3,FILE='V.DAT') OPEN(4,FILE='U.DAT') OPEN(5,FILE='PRESSURE.DAT') WRITE(1,*)'TITLE="PHAI"' WRITE(1,*)'VARIABLE="X","Y","PHAI"' WRITE(1,*)'ZONET="S"','I=',L1,'J=',M1,'C=BLACK' WRITE(2,*)'TITLE="OMEGA"' WRITE(2,*)'VARIABLE="X","Y","OMEGA"' WRITE(2,*)'ZONET="S"','I=',L1,'J=',M1,'C=BLACK' WRITE(3,*)'TITLE="V"' WRITE(3,*)'VARIABLE="X","Y","V"' WRITE(3,*)'ZONET="S"','I=',L1,'J=',M1,'C=BLACK' WRITE(4,*)'TITLE="U"' WRITE(4,*)'VARIABLE="X","Y","U"' WRITE(4,*)'ZONET="S"','I=',L1,'J=',M1,'C=BLACK' WRITE(5,*)'TITLE="PRESSURE"' WRITE(5,*)'VARIABLE="X","Y","PRESSURE"' WRITE(5,*)'ZONET="S"','I=',L1,'J=',M1,'C=BLACK'DOI=1,L1DOJ=1,M1WRITE(1,*)XL(I,J),YM(I,J),PHAI(I,J)WRITE(2,*)XL(I,J),YM(I,J),OMEGA(I,J)WRITE(3,*)XL(I,J),YM(I,J),V(I,J)WRITE(4,*)XL(I,J),YM(I,J),U(I,J)WRITE(5,*)XL(I,J),YM(I,J),P(I,J)ENDDOENDDO CLOSE(1) CLOSE(2) CLOSE(3) CLOSE(4) CLOSE(5)OPEN(6,FILE='V_X.TXT') OPEN(7,FILE='U_Y.TXT') DOI=1,L1 WRITE(6,*)V(I,(M1+1)/2.0) ENDDO DOJ=1,M1 WRITE(7,*)U((L1+1)/2.0,J) ENDDO CLOSE(6) CLOSE(7)!* CONTAINS !网格划分并且设定初场SUBROUTINEMESHING(L,M,X,Y,PHI,OME,U_X,V_Y,P0,XL,YM) INTEGERL,M REALX,Y,PHI(L,M),OME(L,M),U_X(L,M),V_Y(L,M),p0(L,M),XL(L,M),YM(L,M) X=1.0/(L-1) Y=1.0/(M-1) xL(1,:)=0 xL(L,:)=1 YM(:,1)=0 YM(:,M)=1 DOI=1,L DOJ=1,M XL(I,J)=XL(1,J)+(I-1)*X YM(I,J)=YM(I,1)+(J-1)*Y ENDDO ENDDO DOi=1,L!流函数、压力初场 DOj=1,M PHI(i,j)=0 P0(I,J)=1.0 ENDDO ENDDO DOi=1,L!涡量初场 DOj=1,M-1 OME(i,j)=0 ENDDO ENDDO DOi=1,L OME(i,M)=2.0/Y ENDDO DOI=1,L!速度初场 U_X(I,1)=0 V_Y(I,1)=0 U_X(I,M)=1.0 V_Y(I,M)=0 ENDDO DOJ=2,M-1 DOI=1,L U_X(I,J)=0 V_Y(I,J)=0ENDDO ENDDO ENDSUBROUTINE ENDPROGRAM（2）有限容积法源程序：!THISPROGRAMISWRITTENTOSOLVEA2-DDRIVENFLOWPROBLEMINASQUARE!CAVITYUSINGTHEVORTICITY-STREAMFUNCTIONMETHOD.!L1X方向的节点数!M1y方向的节点数!DXX方向的步长!DYY方向的步长!PHAI流函数!OMEGA涡量!RE雷诺数!ITS迭代次数!ERR1，ERR2前后迭代步各点涡量之差PROGRAMMAINIMPLICITNONEREALDX,DYREALERR1,ERR2INTEGER,PARAMETER::L1=21,M1=21INTEGER::RE=100INTEGERi,jINTEGERITS!跌代数INTEGER::ITSMAX=1000!最大跌代数REAL::EPS=1E-6!迭代精度REALN!内迭代次数REAL::ALPHA=0.1!松弛因子REALPHAI(L1,M1)REALOMEGA(L1,M1)REALAP(L1,M1),AW(L1,M1),AE(L1,M1),AS(L1,M1),AN(L1,M1)!一般形式的涡量方程中的各项系数REALOME_K(L1,M1)!存放每一迭代步得到的涡量值REALPHAI_K(L1,M1)!存放每一迭代步得到的流函数值REALCP,CW,CE,CS,CN!一般形式的流函数方程中的各项系数DOUBLEPRECISIONTIME_START,TIME_ENDCALLMESHING(L1,M1,DX,DY,PHAI,OMEGA) CP=2.0/DX**2+2.0/DY**2CW=1.0/DX**2CE=1.0/DX**2CS=1.0/DY**2CN=1.0/DY**2ERR1=1ERR2=1ITS=0 CALLCPU_TIME(TIME_START) CALLMESHING(L1,M1,DX,DY,PHAI,OMEGA,U,V,p,XL,YM) ERR1=1.0 ERR2=1.0 ITS=0 DOWHILE(ERR1>=EPS.OR.ERR2>=EPS.AND.ITS<ITSMAX)DOI=1,L1 DOJ=1,M1 OME_K(I,J)=OMEGA(I,J)!用OME_K储存前一迭代步的计算结果 PHAI_K(I,J)=PHAI(I,J)!用PHAI_K储存前一迭代步的计算结果ENDDOENDDO DOI=2,L1-1!涡量方程系数的更新 DOJ=2,M1-1AW(I,J)=1.0+RE*(U(i-1,j)+U(i,j))*DX/4.0AE(I,J)=1.0-RE*(U(i+1,j)+U(i,j))*DX/4.0AS(I,J)=1.0+RE*(V(i,j-1)+V(i,j))*DY/4.0AN(I,J)=1.0-RE*(V(i,j+1)+V(i,j))*DY/4.0ENDDO ENDDODON=1,200!内迭代次数取200DOi=2,L1-1!解涡量方程DOj=2,M1-1OMEGA(i,j)=(AW(i,j)*OMEGA(i-1,j)+AE(i,j)*OMEGA(i+1,j)& +AS(i,j)*OMEGA(i,j-1)+AN(i,j)*OMEGA(i,j+1))/AP(i,j)ENDDOENDDOENDDODON=1,200DOi=2,L1-1!解流函数方程DOj=2,M1-1PHAI(i,j)=(CW*PHAI(i-1,j)+CE*PHAI(i+1,j)& +CS*PHAI(i,j-1)+CN*PHAI(i,j+1)-OMEGA(i,j))/CPENDDOENDDOENDDOITS=ITS+1ERR1=MAXVAL(ABS(OME_K-OMEGA)) ERR2=MAXVAL(ABS(PHAI_K-PHAI))DOj=1,M1!用亚松弛处理左、右边界的涡量值OMEGA(1,J)=2*(PHAI(2,J)-PHAI(1,J))/DX**2OMEGA(L1,J)=2*(PHAI(L1-1,J)-PHAI(L1,J))/DX**2 OMEGA(1,J)=OME_K(1,j)+ALPHA*(OMEGA(1,J)-OME_K(1,J)) OMEGA(L1,J)=OME_K(L1,J)+ALPHA*(OMEGA(L1,J)-OME_K(L1,J))ENDDODOi=2,L1-1!用亚松弛处理上、下边界的涡量值OMEGA(i,1)=2*(PHAI(i,2)-PHAI(i,1))/DY**2OMEGA(i,M1)=2*(PHAI(i,M1-1)-PHAI(i,M1))/DY**2OMEGA(i,1)=OME_K(i,1)+ALPHA*(OMEGA(i,1)-OME_K(i,1))OMEGA(i,M1)=OME_K(i,M1)+ALPHA*(OMEGA(i,M1)-OME_K(i,M1))ENDDO ENDDO!* ERR3=1.0ITSS=0!用计算得到的速度求解内点压力场分布DOWHILE(ERR3>=EPS.AND.ITSS<=ITSMAX) DOI=1,L1 DOJ=1,M1 P_K(I,J)=P(I,J) ENDDO ENDDO DOI=2,L1-1 DOJ=2,M1-1 P(I,J)=(P(I-1,J)+P(I+1,J)+P(I,J+1)+P(I,J-1)-0.25*((U(I+1,J)-U(I-1,J))*(V(I,J+1)-V(I,J-1))-&(V(I+1,J)-V(I-1,J))*(U(I,J+1)-U(I,J-1))))/4.0 P(I,J)=P_K(I,J)+0.3*(P(I,J)-P_K(I,J)) ENDDO ENDDO!* DOJ=2,M1-1!压力场的左、右边界值 P(1,J)=4.0/3.0*P(2,J)-1.0/3.0*P(3,J)-(2.0/(3.0*RE))*(-2.0*U(2,J)+U(3,J))/D P(L1,J)=4.0/3.0*P(L1-1,J)-1.0/3.0*P(L1-2,J)+(2.0/(3.0*RE))*(-2.0*U(L1-1,J)+U(L1-2,J))/DX ENDDO DOI=2,L1-1!压力场的上、下边界值 P(I,M1)=(4.0/3.0)*P(I,M1-1)-(1.0/3.0)*P(I,M1-2)+(2.0/(3.0*RE*DY))*(-2.0*V(I,M1-1)+V(I,M1-2)) P(I,1)=4.0/3.0*P(I,M1-1)-1.0/3.0*P(I,M1-2)-(2.0/(3.0*RE))*(-2.0*V(I,2)+V(I,3))/DX ENDDO!角点上的压力进行加权平均 P(1,1)=0.5*P(1,2)+0.5*P(2,1) P(1,M1)=0.5*P(1,M1-1)+0.5*P(2,M1) P(L1,1)=0.5*P(L1,2)+