高等流体力学课程论文(附计算程序)

研究生课程考核试卷（适用于课程论文、提交报告）设52/1.110/ms的气体以smV/1.4的速度以零攻角定常绕流长度为1Lm的大平板，试用数值解讨论边界层内的流动规律。解：（1）根据问题的特性，忽略质量力，边界层内水平方向速度为u,竖直方向速度为v,则该流动的控制方程为：C.E．0yvxuM.E.xpyuxuyuvxuu1)(2222ypyvxvyvvxvu1)(2222选特征参数ue，L，将方程组中各变量无量纲化，使各无量纲化变量的取值范围为（0，1）或取1的数量级记作（～1）做变换：x*=x/L∈（0，1）；y*=x/δ=y/AL∈（0，1）u*=u/ue=u/V∞∈（0，1）；v*=v/vmax=v/Bue=v/BV∞∈（0，1）p*=p/ρv∞2（～1）代入方程组，并比较各自的数量级关系，有0][***yvABxuLV1AB*1在连续性方程中去掉任何一项都缺乏物理依据和工程事实，故应考虑两项具有相同数量级，即A=B。]1//[)(**2**222**22******2xpyuALVxuLVLVyuvABxuuLV1AB*1Re1*1Re1*21A*11两个粘性扩散项中，第一项具有很低的数量级，可以去掉，由于边界层中粘性作用不容忽视，故应保留，则有A2=1/Re]11//[)(**2**222**22******2ypAByvALVxvLVLBVyvvABxvuLBV1AB*1Re1*1Re1*21A*1AB1*1对流项与粘性扩散项具有相同的数量级，而压力梯度具有高得多的数量级。去掉各方程中具有明显较低数量级的项，得到简化的近似方程组，即边界层方程为：C.E．0yvxuM.E.xpyuyuvxuu1220yp从上式可以看出，在任一过流断面上，边界层内各点的压力与其外边界上的势流压力pe相等，即)()(),(xpxpyxpeCvpVpee222121这里2222eeeeuvuV则0)2121(11122xuuuvpxxpxpeeee边界层方程进一步简化为：C.E．0yvxuM.E.22yuyuvxuu其定解的边界条件为：Uuyvuy:0,0:0（2）采用无量纲相似性解法引入流函数ψ（x，y），将未知数由两个化为一个令xvyu,做相似性变换引入无因次变量,()()2/yyfgxgxUxU()()fufgUgxUUfyyY1()2uUfUfxxx()()uUUffyygx222()()uUUffyygxg()()2UgvfgUffxxx代入方程，有：21()()22()eUgUUufUffFffxxgxg化简得：0fffB.C.00,(0)0(0)0eyuuff0)0(0,00fvy1)()(,''fUUfuy（3）将边值问题变为初值问题设(0)fA（A为非零常数）做变换，令31A，则有1133()()()()FfAfFA112333[()]ffFAAFAF23()fAFAF43()fAFAF即有0FFF边界条件相应变换为(0)0;(0)0;(0)1FFF由前一个方程可得23321()()1[]()fAFAF根据GEAR算法，编写出求解程序，用于求A、和边界层厚度的程序如下：PROGRAMMAINEXTERNALF,JACOBIDIMENSIONY(3,8),S(3,10),YM(3),D(3),P(3,3),Z(3,549900)DIMENSIONS02(3),ER(3),T(59900),IIS(3),JJS(3)DOUBLEPRECISIONY,S,YM,D,P,Z,S02,ER,T,A,B,H,HMIN,HMAXDOUBLEPRECISIONX(1001),DELT(1001),Y_I(549900),U_Y(549900),V_Y(549900)DATAY/0D0,0D0,1D0,21*0.0/A=0.0B=6.0WRITE(*,*)H=0.0001HMIN=0.000001HMAX=0.0001EPS=1.D-7WRITE(*,*)N=59900M=3CALLGGEAR(A,B,H,HMIN,HMAX,EPS,M,N,Y,T,Z,KF,&F,JACOBI,D,P,S,S02,YM,ER,IIS,JJS)Y=0D0Y(3,1)=1d0/Z(2,N)**1.5CALLGGEAR(A,B,H,HMIN,HMAX,EPS,M,N,Y,T,Z,KF,&F,JACOBI,D,P,S,S02,YM,ER,IIS,JJS)WRITE(*,20)KF20FORMAT(1X,'KF=',I4)DOI=1,NY_D=Z(2,I)ETA_D=T(I)IF(Y_D.GE.0.9999)EXITENDDOOPEN(UNIT=11,FILE='DELT.TXT')OPEN(UNIT=12,FILE='UV_Y.TXT')DOI=1,1001X(I)=1D-3*(I-1)DELT(I)=ETA_D*DSQRT(2*1.1D-5*X(I)/4.1D0)WRITE(11,110)X(I),DELT(I)ENDDO110FORMAT(1X,2F16.8)DOI=1,NY_I(I)=T(I)*DSQRT(2D0*1.1D-5*0.5D0/4.1D0)U_Y(I)=4.1D0*Z(2,I)V_Y(I)=DSQRT(4.1D0*1.1D-5/2D0/0.5D0)*(T(I)*Z(2,I)-Z(1,I))WRITE(12,111)Y_I(I),U_Y(I),V_Y(I)ENDDO111FORMAT(1X,3F16.8)WRITE(*,*)WRITE(*,30)30FORMAT(7X,'T',14X,'Y(1)',11X,'Y(2)',11X,'Y(3)')OPEN(UNIT=10,FILE='RESULT.TXT')DO40I=1,N!WRITE(*,50)T(I),Z(1,I),Z(2,I),Z(3,I)WRITE(10,50)T(I),Z(1,I),Z(2,I),Z(3,I)40CONTINUE50FORMAT(1X,F10.6,5X,3D15.6)WRITE(*,*)ENDSUBROUTINEJACOBI(T,Y,P,M)IMPLICITNONEINTEGER::M,I,JDOUBLEPRECISION::TDOUBLEPRECISION,DIMENSION(M)::Y,Y0DOUBLEPRECISION,DIMENSION(M,M)::PINTEGER,DIMENSION(M,M)::EDOUBLEPRECISION::H0DOUBLEPRECISION,DIMENSION(M)::Y1,Y2,&Y3,Y4,Y5DOUBLEPRECISION,DIMENSION(M)::FF1,FF2,&FF3,FF4,FF5H0=0.000001E=0D0Y0=Y!****************单Ì£¤位?矩?阵¨®***************DOI=1,ME(I,I)=1ENDDO!**********利¤?用®?五?点Ì?公?式º?数ºy值¦Ì微¡é分¤?*********DOI=1,MY1=Y0-3*H0*E(I,:)Y2=Y0-2*H0*E(I,:)Y3=Y0-H0*E(I,:)Y4=Y0Y5=Y0+H0*E(I,:)CALLF(T,Y1,M,FF1)CALLF(T,Y2,M,FF2)CALLF(T,Y3,M,FF3)CALLF(T,Y4,M,FF4)CALLF(T,Y5,M,FF5)P(:,I)=(-FF1+6*FF2-18*FF3+&10*FF4+3*FF5)/12/H0ENDDOENDSUBROUTINEJACOBISUBROUTINEBRINV(A,N,L,IS,JS)DIMENSIONA(N,N),IS(N),JS(N)DOUBLEPRECISIONA,T,DL=1DO100K=1,ND=0.0DO10I=K,NDO10J=K,NIF(ABS(A(I,J)).GT.D)THEND=ABS(A(I,J))IS(K)=IJS(K)=JENDIF10CONTINUEIF(D+1.0.EQ.1.0)THENL=0WRITE(*,20)RETURNENDIF20FORMAT(1X,'ERR**NOTINV')DO30J=1,NT=A(K,J)A(K,J)=A(IS(K),J)A(IS(K),J)=T30CONTINUEDO40I=1,NT=A(I,K)A(I,K)=A(I,JS(K))A(I,JS(K))=T40CONTINUEA(K,K)=1/A(K,K)DO50J=1,NIF(J.NE.K)THENA(K,J)=A(K,J)*A(K,K)ENDIF50CONTINUEDO70I=1,NIF(I.NE.K)THENDO60J=1,NIF(J.NE.K)THENA(I,J)=A(I,J)-A(I,K)*A(K,J)ENDIF60CONTINUEENDIF70CONTINUEDO80I=1,NIF(I.NE.K)THENA(I,K)=-A(I,K)*A(K,K)ENDIF80CONTINUE100CONTINUEDO130K=N,1,-1DO110J=1,NT=A(K,J)A(K,J)=A(JS(K),J)A(JS(K),J)=T110CONTINUEDO120I=1,NT=A(I,K)A(I,K)=A(I,IS(K))A(I,IS(K))=T120CONTINUE130CONTINUERETURNENDSUBROUTINEGGEAR(A,B,H,HMIN,HMAX,EPS,M,N,Y,T,Z,KF,*F,JACOBI,D,P,S,S02,YM,ER,IIS,JJS)DIMENSIONT(N),Z(M,N),P(M,M),Y(M,8),D(M),S(M,10)DIMENSIONYM(M),S02(M),ER(M),AA(7),PP(7,3),IIS(M),JJS(M)DOUBLEPRECISIONT,Z,D,P,Y,S,YM,S02,ER,AA,PP,A,B,H,*HMIN,HMAX,HW,HD,RM,T0,TD,R,DD,PR1,PR2,PR3,RRDATAPP/2.0,4.5,7.333,10.42,13.7,17.15,1.0,*3.0,6.0,9.167,12.5,15.98,1.0,1.0,*1.0,1.0,0.5,0.1667,0.04133,0.008267,1.0/AA(2)=-1.0JT=0NN=0KF=1T0=ANQ=1CALLF(T0,Y,M,D)DO10I=1,M10Y(I,2)=D(I)*HHW=HK=2DO15I=1,M15YM(I)=1.020IRT=1KF=1NN=NN+1T(NN)=T0DO30I=1,M30Z(I,NN)=Y(I,1)IF(T0.GE.B)RETURNIF(NN.EQ.N)RETURNDO40I=1,MDO40J=1,K40S(I,J)=Y(I,J)HD=HWIF(H.NE.HD)THENRM=H/HDIRT1=1CALLG75(HMIN,HD,RM,HMAX,Y,S,H,IDB,K,M)ENDIFNQD=NQTD=