北航计算流体力学大作业代码.doc-资料库

qq_31936569-9557840-16359647605505794367.doc.pdf-第1页.png

第1页 / 共12页

qq_31936569-9557840-16359647605505794367.doc.pdf-第2页.png

第2页 / 共12页

qq_31936569-9557840-16359647605505794367.doc.pdf-第3页.png

第3页 / 共12页

qq_31936569-9557840-16359647605505794367.doc.pdf-第4页.png

第4页 / 共12页

qq_31936569-9557840-16359647605505794367.doc.pdf-第5页.png

第5页 / 共12页

qq_31936569-9557840-16359647605505794367.doc.pdf-第6页.png

第6页 / 共12页

qq_31936569-9557840-16359647605505794367.doc.pdf-第7页.png

第7页 / 共12页

qq_31936569-9557840-16359647605505794367.doc.pdf-第8页.png

第8页 / 共12页

! ! ! ! This code is for the CFD course. The mesh generation module, the time-marching processing module, and the post-processing module are available in this code. Students need to work out the module which evaluates the numerical advective flux. The flow case to be calculated is NACA0012 flow, for which with different specified boundary conditions the flow could be subsonic, transonic, or supersonic. !######################################################################! ! ! ! ! ! ! ! !######################################################################! ! ! !======================================================================! ! ! ! !======================================================================! The module definition in which contains the global variables and arrays. These variables and arrays can be "seen" in any subroutine ! when the module is included at the begining of the subroutine. ! ! ! ! ! MODULE MODGLOB ! !---The dimensions of the computational mesh; IN is the number of mesh !---points along streamwise; JN is that on the normal direction. INTEGER IN, JN ! !---NPASS is the index of the current time-marching step; NPASSM is the !---maximum time-marching loops. INTEGER NPASS, NPASSM ! !---The following parameters are in turn the static pressure at inlet, !---the density at inlet, the Mach number at inlet, the static pressure !---at outlet, and CFL number. REAL*8 SPINL, DNINL, VMINL, AGINL, SPOUT, CFL ! !---The following are thermodynamic parameters and their derivations, see !---their definitions in the code. REAL*8 GAMA, RGAS, GAMAM, GAMSM, GAMMH ! !---X and Y are mesh coordinates. !---XLNI,YLNI,SLNI are length vectors and the magnitudes of the length of the I-direction cell surfaces. !---XLNJ,YLNJ,SLNJ are length vectors and the magnitudes of the length of the J-direction cell surfaces. !---SARC are cell areas of the control volumes divided by mesh lines. REAL*8, ALLOCATABLE :: & X(:,:), Y(:,:),& XLNI(:,:), YLNI(:,:), SLNI(:,:), XLNJ(:,:), YLNJ(:,:), SLNJ(:,:), SARC(:,:)

! !---DENS,VELX,VELY,PRES are density, Vx, Vy, and static pressure. For !---cell-centered finite volume method as applied in this code, these !---flow variables are saved at the cell centers. Additionally, the pseudo !---cells are used to apply the boundary conditions, so that these arrays !---include the values at centers of pseudo cells. !---FLUX is the residual (sum of the advective fluxes). REAL*8, ALLOCATABLE :: & DENS(:,:), VELX(:,:), VELY(:,:), PRES(:,:), FLUX(:,:,:) ! END MODULE !====================== End of Module definition ======================= ! ! !========================== Main program begin ========================= PROGRAM MAIN USE MODGLOB IMPLICIT REAL*8(A-H,O-Z) ! PI = DACOS(-1.D0) !---Define the dimension of the computational mesh: IN = 281 JN = 61 !---Along the flow direction. !---On the normal direction. ! WRITE(*,*) 'Input the type of inflow:' WRITE(*,*) '1=supersonic, 2=subsonic without shock, 3=subsonic with attack angle' READ(*,*) NTASK ! !---Define the B.C. at inlet: SPINL = 50000.D0 DNINL = 1.22D0 IF(NTASK==1) THEN VMINL = 2.D0 AGINL = 0.D0 ELSE IF(NTASK==2) THEN VMINL = 0.7D0 AGINL = 0.D0 ELSE VMINL = 0.7D0 AGINL = 5.D0/180.D0*PI END IF ! !---Define the B.C. at outlet: SPOUT = 50000.D0

! !---Define the CFL number: CFL = 0.8D0 ! !---Define the maximum time-marching loops: NPASSM = 10000 ! ! ! ! ! ! ALLOCATE( X(IN,JN), Y(IN,JN) ) ALLOCATE( XLNI(IN,JN-1), YLNI(IN,JN-1), SLNI(IN,JN-1),& XLNJ(IN-1,JN), YLNJ(IN-1,JN), SLNJ(IN-1,JN), SARC(IN-1,JN-1) ) ALLOCATE( DENS(-1:IN+1,-1:JN+1), VELX(-1:IN+1,-1:JN+1), VELY(-1:IN+1,-1:JN+1),& PRES(-1:IN+1,-1:JN+1), FLUX(4,IN-1,JN-1) ) = 287.06D0 RGAS GAMA = 1.4D0 GAMAM = GAMA - 1.D0 GAMSM = GAMA/GAMAM GAMMH = GAMAM/2.D0 CALL READMESH CALL INIFLOW DO NPASS=1,NPASSM CALL DEFBCON CALL RESIDUAL CALL UPDATEQ END DO CALL DEFBCON CALL POSTOUT STOP END !============================ End of MAIN ============================== ! !======================================================================! ! !======================================================================! Generate the computational mesh of the nozzle. ! SUBROUTINE READMESH USE MODGLOB IMPLICIT REAL*8(A-H,O-Z) DIMENSION ITMP(50), JTMP(50) ! !------------------------- Read in the mesh --------------------------

OPEN(61,FILE='ingrid.dat',FORM='UNFORMATTED') READ(61) NBKM READ(61)((ITMP(N),JTMP(N)),N=1,NBKM) DO N=1,NBKM READ(61) ((X(I,J),I=1,IN),J=1,JN),& ((Y(I,J),I=1,IN),J=1,JN) END DO CLOSE(61) ! OPEN(71,FILE='pltgrid.plt') WRITE(71,*) 'TITLE="NONAME"' WRITE(71,*) 'VARIABLES="X","Y"' WRITE(71,*) 'ZONE T="NONAME", I=',IN,',J=',JN,', F=POINT' DO J=1,JN DO I=1,IN WRITE(71,*) X(I,J), Y(I,J) END DO END DO CLOSE(71) ! !------------ Calculate the surface length vectors and ----------------- !------------ the areas of 2D control volumes divided ----------------- !------------ by mesh lines. ----------------- 21 ! 22 ! DO 21 J=1,JN-1 DO 21 I=1,IN XLNI(I,J) = Y(I,J+1) - Y(I,J YLNI(I,J) = X(I,J SLNI(I,J) = DSQRT( XLNI(I,J)*XLNI(I,J) + YLNI(I,J)*YLNI(I,J) ) ) ) - X(I,J+1) CONTINUE DO 22 J=1,JN DO 22 I=1,IN-1 XLNJ(I,J) = Y(I YLNJ(I,J) = X(I+1,J) - X(I SLNJ(I,J) = DSQRT( XLNJ(I,J)*XLNJ(I,J) + YLNJ(I,J)*YLNJ(I,J) ) ,J) - Y(I+1,J) ,J) CONTINUE DO 23 J=1,JN-1 DO 23 I=1,IN-1 I1 = I + 1 J1 = J + 1 SARC(I,J) = 0.5D0*( DABS( ( X(I ,J ) - X(I1,J ) )*Y(I1,J1)& +( X(I1,J ) - X(I1,J1) )*Y(I ,J )& +( X(I1,J1) - X(I ,J ) )*Y(I1,J ) )&

+ DABS( ( X(I ,J ) - X(I1,J1) )*Y(I ,J1)& +( X(I1,J1) - X(I ,J1) )*Y(I ,J )& +( X(I ,J1) - X(I ,J ) )*Y(I1,J1) ) ) 23 ! CONTINUE RETURN END !========================== End of READMESH ============================ ! !======================================================================! ! ! !======================================================================! The following method generates a mean flow field. Generate the initial flow field. ! ! SUBROUTINE INIFLOW USE MODGLOB IMPLICIT REAL*8(A-H,O-Z) DNINI = DNINL SPINI = SPINL VTINI = VMINL*DSQRT(GAMA*SPINL/DNINL) DO 11 J=1,JN-1 DO 11 I=1,IN-1 DENS(I,J) = DNINI VELX(I,J) = VTINI*DCOS(AGINL) VELY(I,J) = VTINI*DSIN(AGINL) PRES(I,J) = SPINI CONTINUE RETURN END ! ! 11 ! !=========================== End of INIFLOW ============================ ! !======================================================================! ! !======================================================================! Define the boundary conditions. ! SUBROUTINE DEFBCON USE MODGLOB IMPLICIT REAL*8(A-H,O-Z) ! !---Define the B.C. at far field boundary(J=JN): TEMP = 1.D0 + GAMMH*VMINL**2 TPINL = SPINL*TEMP**GAMSM STINL = SPINL/DNINL/RGAS

TTINL = STINL*TEMP VTINL = VMINL*DSQRT(GAMA*RGAS*STINL) DO 11 I=1,IN-1 SNX = XLNJ(I,JN)/SLNJ(I,JN) SNY = YLNJ(I,JN)/SLNJ(I,JN) VNFFB = VELX(I,JN-1)*SNX + VELY(I,JN-1)*SNY VTFFB = DSQRT(VELX(I,JN-1)**2+VELY(I,JN-1)**2) SNFFB = DSQRT(GAMA*PRES(I,JN-1)/DENS(I,JN-1)) VMCHB = VTFFB/SNFFB IF(VNFFB<0.D0) THEN IF(VMCHB>=1.D0) THEN DENS(I,JN) = DNINL VELX(I,JN) = VTINL*DCOS(AGINL) VELY(I,JN) = VTINL*DSIN(AGINL) PRES(I,JN) = SPINL ELSE SPI = DMIN1( PRES(I,JN-1), TPINL ) TMP = (TPINL/SPI)**(1.D0/GAMSM) VMI = DSQRT((TMP-1.D0)/GAMMH) STI = TTINL/TMP VTI = VMI*DSQRT(GAMA*RGAS*STI) DENS(I,JN) = SPI/STI/RGAS VELX(I,JN) = VTI*DCOS(AGINL) VELY(I,JN) = VTI*DSIN(AGINL) PRES(I,JN) = SPI END IF ELSE IF(VMCHB>=1.D0) THEN PRES(I,JN) = PRES(I,JN-1) ELSE PRES(I,JN) = SPOUT END IF DENS(I,JN) = DENS(I,JN-1) VELX(I,JN) = VELX(I,JN-1) VELY(I,JN) = VELY(I,JN-1) END IF CONTINUE 11 ! !---Define the wall B.C. at the airfoil surface(J=1): DO 21 I=61,220 VNORM = VELX(I,1)*XLNJ(I,1) + VELY(I,1)*YLNJ(I,1) VTEMP = 2.D0*VNORM/SLNJ(I,1)/SLNJ(I,1) DENS(I,0) = DENS(I,1) VELX(I,0) = VELX(I,1) - VTEMP*XLNJ(I,1)

VELY(I,0) = VELY(I,1) - VTEMP*YLNJ(I,1) PRES(I,0) = PRES(I,1) CONTINUE 21 ! !---Define the internal B.C.(J=1) : !---The first layer: DO 31 I=1,60 DENS(I,0) = DENS(281-I,1) VELX(I,0) = VELX(281-I,1) VELY(I,0) = VELY(281-I,1) PRES(I,0) = PRES(281-I,1) CONTINUE 31 ! !----The second layer: DO 32 I=1,60 DENS(I,-1) = DENS(281-I,2) VELX(I,-1) = VELX(281-I,2) VELY(I,-1) = VELY(281-I,2) PRES(I,-1) = PRES(281-I,2) CONTINUE 32 ! !---Define the internal B.C.(J=1) : !---The first layer: DO 41 I=221,IN-1 DENS(I,0) = DENS(281-I,1) VELX(I,0) = VELX(281-I,1) VELY(I,0) = VELY(281-I,1) PRES(I,0) = PRES(281-I,1) CONTINUE 41 ! !----The second layer: DO 42 I=221,IN-1 DENS(I,-1) = DENS(281-I,2) VELX(I,-1) = VELX(281-I,2) VELY(I,-1) = VELY(281-I,2) PRES(I,-1) = PRES(281-I,2) CONTINUE 42 ! !---Define the outlet B.C.(I=1) : DO 51 J=1,JN-1 VTBCD = DSQRT(VELX(1,J)**2+VELY(1,J)**2) SNBCD = DSQRT(GAMA*PRES(1,J)/DENS(1,J)) VMCHB = VTBCD/SNBCD IF(VMCHB>=1.D0) THEN

PRES(0,J) = PRES(1,J) ELSE PRES(0,J) = SPOUT END IF DENS(0,J) = DENS(1,J) VELX(0,J) = VELX(1,J) VELY(0,J) = VELY(1,J) CONTINUE 51 ! !---Define the outlet B.C.(I=IN) : DO 52 J=1,JN-1 VTBCD = DSQRT(VELX(IN-1,J)**2+VELY(IN-1,J)**2) SNBCD = DSQRT(GAMA*PRES(IN-1,J)/DENS(IN-1,J)) VMCHB = VTBCD/SNBCD IF(VMCHB>=1.D0) THEN PRES(IN,J) = PRES(IN-1,J) ELSE PRES(IN,J) = SPOUT END IF DENS(IN,J) = DENS(IN-1,J) VELX(IN,J) = VELX(IN-1,J) VELY(IN,J) = VELY(IN-1,J) 52 ! CONTINUE RETURN END !=========================== End of DEFBCON ============================ ! !======================================================================! ! !======================================================================! Calculate the residual which is the sum of numerical flux. ! SUBROUTINE RESIDUAL USE MODGLOB IMPLICIT REAL*8(A-H,O-Z) ! ! RETURN END !=========================== End of ADVFLUX ============================ ! !======================================================================! ! ! !======================================================================! Update the conservative variables at time level n+1 using a simple one-step time marching scheme. ! !

资料库

北航计算流体力学大作业代码.doc

相关推荐

课程资源

热门标签

最新资料