|
Patankar CFD FORTRAN 90 Code FVM
!===============================================
This is a general purpose transient 2-d CFD code (FORTRAN 90)based on FVM. !************************************************* ************ Note to the users: Please go through this brief instruction before attempting anything with this code !************************************************* ************ 1. Most of the variables used in the code are listed in the file "list_of_var". However, some new variables are also used, typically for the USER portion. Those variables are explained just before opening the input file in subroutine GRID 2. The coordinate system used is right-handed. The systems are as: mode=1 : cartesian; mode=2: cylindrical(axisym); mode=3: polar 3. The user is reqd. to use only six subroutines almost. They are: a) Subroutine GRID: For specification of input data and initialization of arrays. It calls subroutine UGRID if one selects uniform gridding; and calls subroutine MESH if one selects non-uniform gridding. The selection is made by giving a value to a variable called NGRID introduced in subroutine GRID. Do not forget to specify following in subroutine grid: # mode: to specify coordinate system. For a system involving radial direction, also specify r(1). # lsolve: to select which variables you want to solve; eg. lsolve(4)=.true. will make eqn corresponding to nf=4(ie., energy eqn) to be solved. Note that nf=1,2,3,4,np,6 correspond to solving for u,v,pc,T,p,conc. in the code. !# ngrid: to specify choice of grid. For non-uniform grid, subroutine MESH is called. It takes input from a file "grid.dat", where proportional lengths of non-uniform grids are specified. However, for class-problems; generally uniform grid will suffice, for which specify NGRID =1 # no. of grids in coordinate dirns, domain size, time step(dt) etc. # eru,erv,ert: For convergence criteria # property values as reqd. # initialise all arrays, logical variables b) Subroutine START: For giving initial condns for variables (say u,v,T) for the problem c) Subroutine DENSE : For specifying density of the working medium. d) Subroutine BOUND: For giving boundary condns. for the problem e) Subroutine GAMSOR: For giving diffusion coeff and source term for all nfs (i.e., all eqns). As per symbols used in Patankar's book; the symbol 'ap' used in the code corresponds to 'Sp' and the symbol 'con' used in the code corresponds to 'Sc'. 'gam(:,: )' corresponds to the diffusion coefficient 'gamma' as per general 'phi' formulation f) Subroutine Printout: This subroutine will give the plot files which contain isotherms and velocity vectors plot compatible with matlab or tecplot (as you choose). There is a logical switch called ltecplot declared in the programme. If you declare ltecplot as true, then you will get matlab files only. You can run the ".m" files generated in MATLAB to see the outputs.However, you are encouraged to develop output format according to your own choice and needs. 4. For steady state soln; use a typically large value of time step(dt) as per instruction of the instructor 5. The code deals with a sample problem of natural convection in a rectangular cavity.The cavity is heated from two vertical sides by employing a uniform heat flux (symbolised as FLUX in the code). The top and bottom are insulated. Due to heating, there is a natural convection driven flow inside the liquid that develops with time. We are interested in studying the evolution of temperature field and velocity field inside the cavity. 6.Once you get the code, first compile it as it is. For compilation of the code in UNIX Systems, either compile in IBM m/c using the compilation command: xlf90 code.f OR complile in a COMPAQ machine using the compilation command f90 -fast code.90 Note: for compiling and running in IBM, you have to rename the file from code.f90 to code.f Once it is compiled, execute the file a.out by typing a.out in the command prompt of your terminal. You will soon see that in your screen certain figures are appearing. We are basically printing the values of u,v,T and their relative errors in every iteration to monitor convergence. Once it converges within a timestep, it goes to the next timestep. In this programme, I have kept timestep, i.e., dt=10s. The time domain I have kept as 510s, and in subroutine printout I have stated that the MATLAB output should be written every 50 timestep. Accordingly, after 500s, this programme will open a file out500.m You can run the file in MATLAB to see various plots. Note that the streamfunction is symbolised as pc in the output. Do not confuse it with the pressure correction that we are solving to get the fluid flow.The figure 1 will give you temperature contours, figure 2 will give you velocity vectors and figure 3 will give you streamfunctions.Once you start solving problems on your own, initially you will be solving steady state problems only. For that, use dt=1.0e30 (a typically large no.) and tlast=1.0e30. Change the name !of the .m file you want to have by altering the filename in subroutine printout. Note that your tlast (i.e. end of time domain) should be same as dt for a steady state problem |
!################################################# #################
!================================================= ================= !GLOBAL DECLARATIONS !================================================= ================= MODULE GLOBAL IMPLICIT NONE INTEGER,PARAMETER :: ni=40 INTEGER,PARAMETER :: nj=40 INTEGER,PARAMETER :: nij=40 INTEGER,PARAMETER :: nfmax=10 INTEGER,PARAMETER :: np=11 INTEGER,PARAMETER :: nfx3=nfmax+1 REAL,SAVE,DIMENSION(NI,NJ) :: u,v,pc,t,p,conc,rho,gam,con,ur,vr REAL,SAVE,DIMENSION(NI,NJ) :: aip,aim,ajp,ajm,ap,ap0 REAL,SAVE,DIMENSION(NI,NJ) :: du,dv REAL,SAVE,DIMENSION(NI,NJ) :: tO,uO,vO,told,uold,vold,clold,uhat,vhat REAL,SAVE,DIMENSION(NI,NJ) :: fle,flw,fln,fls REAL,SAVE,DIMENSION(NI,NJ,6) :: cofu,cofv,cofp,cof REAL,SAVE,DIMENSION(NI) :: x,xu,xdif,xcv,xcvs,fv,fvp,fx,fxm,xcvi,xcvip,ara REAL,SAVE,DIMENSION(NJ) :: y,yv,ydif,ycv,ycvs,ycvr,ycvrs,arx,arxj,arxjp REAL,SAVE,DIMENSION(NJ) :: r,rmn,sx,sxmn,fy,fym REAL,SAVE,DIMENSION(NIJ) :: pt,qt REAL,SAVE,DIMENSION(NFX3) :: relax REAL,SAVE,DIMENSION(NFX3) :: ntimes REAL,SAVE :: rhocon,ert,eru,erv,time,flux REAL,SAVE :: xl,yl,dt,tlast,dt2,time2,dt3,time3,dt4,time4,dt5,t ime5 REAL,SAVE :: smax,ssum,flow,diff,acof REAL,SAVE :: ts,tl,g,beta,tref,rhoref,amiu,constant1,constant2, constant3,tk,cp REAL,SAVE :: ax1,ay1,rhoinitial, tinitial REAL,SAVE :: deltmx,delumx,delvmx INTEGER,SAVE :: itll,NF,NRHO,NGAM,L1,L2,L3,M1,M2,M3,ngrid,ny1 INTEGER,SAVE :: IST,JST,ITER,ISTP,IPREF,JPREF,MODE,ITERL,NTIME CHARACTER(LEN=4),DIMENSION(NFX3) :: title LOGICAL,SAVE,DIMENSION(NFX3) :: LSOLVE,LPRINT,LBLK LOGICAL,SAVE :: LSTOP,lconv,ltecplot END MODULE GLOBAL !================================================= ================= |
!================================================= =====================
!***************** PROGRAM MAIN *********************************** !================================================= ===================== PROGRAM MAIN USE GLOBAL IMPLICIT NONE CALL GRID CALL geomet CALL START mainloop: DO WHILE(.NOT.lstop) CALL DENSE CALL BOUND CALL OLDVAL CALL COEFF END DO mainloop WRITE(*,*)" PROGRAM FINISHED. " END PROGRAM MAIN !================================================= ===================== |
!===============================================
!this subroutine calculates discretization eqn coeffs as per power law !************************************************* ************ subroutine diflow use global implicit none real :: Temp acof=diff if(flow.eq.0.) return temp=diff-abs(flow)*0.1 acof=0. If(temp.le.0.) return temp=temp/diff acof=diff*temp**5 end subroutine diflow !=============================================== |
!================================================= =====================
!THIS SUBROUTINE SOLVES DISCRETISATION EQUATIONS BY 'LINE BY LINE TDMA' !************************************************* ********************* SUBROUTINE SOLVE USE GLOBAL IMPLICIT NONE INTEGER :: I,J,II,JJ,N,NT,ISTF,JSTF,IT1,IT2,JT1,JT2 REAL :: BL,BLP,BLM,BLC,TEMP,DENOM REAL,DIMENSION(NI,NJ) ::f SELECT CASE(nf) CASE(1) f = u CASE(2) f = v CASE(3) f = pc CASE(4) f = t CASE(np) f = p CASE(6) f =conc END SELECT !---------------------------------------------------------------------- ISTF=IST-1 JSTF=JST-1 IT1=L2+IST IT2=L3+IST JT1=M2+JST JT2=M3+JST !---------------------------------------------------------------------- DO 999 NT=1,NTIMES(NF) N=NF !---------------------------------------------------------------------- IF(.NOT.LBLK(NF)) GO TO 10 PT(ISTF)=0. QT(ISTF)=0. DO 11 I=IST,L2 BL=0. BLP=0. BLM=0. BLC=0. DO 12 J=JST,M2 BL=BL+AP(I,J) IF(J.NE.M2) BL=BL-AJP(I,J) IF(J.NE.JST) BL=BL-AJM(I,J) BLP=BLP+AIP(I,J) BLM=BLM+AIM(I,J) BLC=BLC+CON(I,J)+AIP(I,J)*F(I+1,J)+AIM(I,J)*F(I-1,J) & & +AJP(I,J)*F(I,J+1)+AJM(I,J)*F(I,J-1)-AP(I,J)*F(I,J) 12 END DO DENOM=BL-PT(I-1)*BLM IF(ABS(DENOM/BL).LT.1.E-10) DENOM=1.E30 PT(I)=BLP/DENOM QT(I)=(BLC+BLM*QT(I-1))/DENOM 11 END DO BL=0. DO 13 II=IST,L2 I=IT1-II BL=BL*PT(I)+QT(I) DO 13 J=JST,M2 13 F(I,J)=F(I,J)+BL PT(JSTF)=0. QT(JSTF)=0. DO 21 J=JST,M2 BL=0. BLP=0. BLM=0. BLC=0. DO 22 I=IST,L2 BL=BL+AP(I,J) IF(I.NE.L2) BL=BL-AIP(I,J) IF(I.NE.IST) BL=BL-AIM(I,J) BLP=BLP+AJP(I,J) BLM=BLM+AJM(I,J) BLC=BLC+CON(I,J)+AIP(I,J)*F(I+1,J)+AIM(I,J)*F(I-1,J) & & +AJP(I,J)*F(I,J+1)+AJM(I,J)*F(I,J-1)-AP(I,J)*F(I,J) 22 END DO DENOM=BL-PT(J-1)*BLM IF(ABS(DENOM/BL).LT.1.E-10) DENOM=1.E30 PT(J)=BLP/DENOM QT(J)=(BLC+BLM*QT(J-1))/DENOM 21 END DO BL=0. DO 23 JJ=JST,M2 J=JT1-JJ BL=BL*PT(J)+QT(J) DO 23 I=IST,L2 23 F(I,J)=F(I,J)+BL 10 CONTINUE !---------------------------------------------------------------------- DO 90 J=JST,M2 PT(ISTF)=0. QT(ISTF)=F(ISTF,J) DO 70 I=IST,L2 DENOM=AP(I,J)-PT(I-1)*AIM(I,J) PT(I)=AIP(I,J)/DENOM TEMP=CON(I,J)+AJP(I,J)*F(I,J+1)+AJM(I,J)*F(I,J-1) QT(I)=(TEMP+AIM(I,J)*QT(I-1))/DENOM 70 END DO DO 80 II=IST,L2 I=IT1-II 80 F(I,J)=F(I+1,J)*PT(I)+QT(I) 90 END DO !---------------------------------------------------------------------- DO 190 JJ=JST,M3 J=JT2-JJ PT(ISTF)=0. QT(ISTF)=F(ISTF,J) DO 170 I=IST,L2 DENOM=AP(I,J)-PT(I-1)*AIM(I,J) PT(I)=AIP(I,J)/DENOM TEMP=CON(I,J)+AJP(I,J)*F(I,J+1)+AJM(I,J)*F(I,J-1) QT(I)=(TEMP+AIM(I,J)*QT(I-1))/DENOM 170 END DO DO 180 II=IST,L2 I=IT1-II 180 F(I,J)=F(I+1,J)*PT(I)+QT(I) 190 END DO !---------------------------------------------------------------------- DO 290 I=IST,L2 PT(JSTF)=0. QT(JSTF)=F(I,JSTF) DO 270 J=JST,M2 DENOM=AP(I,J)-PT(J-1)*AJM(I,J) PT(J)=AJP(I,J)/DENOM TEMP=CON(I,J)+AIP(I,J)*F(I+1,J)+AIM(I,J)*F(I-1,J) QT(J)=(TEMP+AJM(I,J)*QT(J-1))/DENOM 270 END DO DO 280 JJ=JST,M2 J=JT1-JJ 280 F(I,J)=F(I,J+1)*PT(J)+QT(J) 290 END DO !---------------------------------------------------------------------- DO 390 II=IST,L3 I=IT2-II PT(JSTF)=0. QT(JSTF)=F(I,JSTF) DO 370 J=JST,M2 DENOM=AP(I,J)-PT(J-1)*AJM(I,J) PT(J)=AJP(I,J)/DENOM TEMP=CON(I,J)+AIP(I,J)*F(I+1,J)+AIM(I,J)*F(I-1,J) QT(J)=(TEMP+AJM(I,J)*QT(J-1))/DENOM 370 END DO DO 380 JJ=JST,M2 J=JT1-JJ 380 F(I,J)=F(I,J+1)*PT(J)+QT(J) 390 END DO !---------------------------------------------------------------------- 999 END DO SELECT CASE(NF) CASE(1) u = f CASE(2) v = f CASE(3) pc = f CASE(4) t = f CASE(np) p = f CASE(6) conc = f END SELECT CALL RESET END SUBROUTINE SOLVE !================================================= ===================== |
!================================================= =====================
!THIS SUBROUTINE RESETS 'ap' AND 'con' TO zero !************************************************* ********************* SUBROUTINE RESET USE GLOBAL IMPLICIT NONE INTEGER :: I,J DO J=2,M2 DO I=2,L2 CON(I,J)=0. AP(I,J)=0. END DO END DO END SUBROUTINE RESET !================================================= ===================== |
!================================================= =====================
! THIS SUBROUTINE GENERATES THE CONTROL VOLUMES&VARIOUS RELATED GEOMETRICAL PARAMETERS !************************************************* ********************* SUBROUTINE GEOMET USE GLOBAL IMPLICIT NONE INTEGER :: I,J,IM4,IM5 1 FORMAT(//15X,'COMPUTATION IN CARTESIAN COORDINATES') 2 FORMAT(//15X,'COMPUTATION FOR AXISYMMETRIC SITUATION') 3 FORMAT(//15X,'COMPUTATION IN POLAR COORDINATES') 4 FORMAT(14X,40(1H*),//) !----------------------------------------------------------------------- NRHO=NP+1 NGAM=NRHO+1 L2=L1-1 L3=L2-1 M2=M1-1 M3=M2-1 X(1)=XU(2) DO 5 I=2,L2 5 X(I)=0.5*(XU(I+1)+XU(I)) X(L1)=XU(L1) Y(1)=YV(2) DO 10 J=2,M2 10 Y(J)=0.5*(YV(J+1)+YV(J)) Y(M1)=YV(M1) DO 15 I=2,L1 15 XDIF(I)=X(I)-X(I-1) DO 18 I=2,L2 18 XCV(I)=XU(I+1)-XU(I) DO 20 I=3,L2 20 XCVS(I)=XDIF(I) XCVS(3)=XCVS(3)+XDIF(2) XCVS(L2)=XCVS(L2)+XDIF(L1) DO 22 I=3,L3 XCVI(I)=0.5*XCV(I) 22 XCVIP(I)=XCVI(I) XCVIP(2)=XCV(2) XCVI(L2)=XCV(L2) DO 35 J=2,M1 35 YDIF(J)=Y(J)-Y(J-1) DO 40 J=2,M2 40 YCV(J)=YV(J+1)-YV(J) DO 45 J=3,M2 45 YCVS(J)=YDIF(J) YCVS(3)=YCVS(3)+YDIF(2) YCVS(M2)=YCVS(M2)+YDIF(M1) IF(MODE.NE.1) GO TO 55 DO 52 J=1,M1 RMN(J)=1.0 52 R(J)=1.0 GO TO 56 55 DO 50 J=2,M1 50 R(J)=R(J-1)+YDIF(J) RMN(2)=R(1) DO 60 J=3,M2 60 RMN(J)=RMN(J-1)+YCV(J-1) RMN(M1)=R(M1) 56 CONTINUE DO 57 J=1,M1 SX(J)=1. SXMN(J)=1. IF(MODE.NE.3) GO TO 57 SX(J)=R(J) IF(J.NE.1) SXMN(J)=RMN(J) 57 END DO DO 62 J=2,M2 YCVR(J)=R(J)*YCV(J) ARX(J)=YCVR(J) IF(MODE.NE.3) GO TO 62 ARX(J)=YCV(J) 62 END DO DO 64 J=4,M3 64 YCVRS(J)=0.5*(R(J)+R(J-1))*YDIF(J) YCVRS(3)=0.5*(R(3)+R(1))*YCVS(3) YCVRS(M2)=0.5*(R(M1)+R(M3))*YCVS(M2) IF(MODE.NE.2) GO TO 67 DO 65 J=3,M3 ARXJ(J)=0.25*(1.+RMN(J)/R(J))*ARX(J) 65 ARXJP(J)=ARX(J)-ARXJ(J) GO TO 68 67 DO 66 J=3,M3 ARXJ(J)=0.5*ARX(J) 66 ARXJP(J)=ARXJ(J) 68 ARXJP(2)=ARX(2) ARXJ(M2)=ARX(M2) DO 70 J=3,M3 FV(J)=ARXJP(J)/ARX(J) 70 FVP(J)=1.-FV(J) DO 85 I=3,L2 FX(I)=0.5*XCV(I-1)/XDIF(I) 85 FXM(I)=1.-FX(I) FX(2)=0. FXM(2)=1. FX(L1)=1. FXM(L1)=0. DO 90 J=3,M2 FY(J)=0.5*YCV(J-1)/YDIF(J) 90 FYM(J)=1.-FY(J) FY(2)=0. FYM(2)=1. FY(M1)=1. FYM(M1)=0. !---------------------------------------------------------------------- IM4=M1-3 IM5=M1-4 117 FORMAT(5E10.4) IF(MODE.EQ.1) PRINT 1 IF(MODE.EQ.2) PRINT 2 IF(MODE.EQ.3) PRINT 3 PRINT 4 open(unit=20,file="grids.dat") write(20,*)l1,m1 write(20,*)x write(20,*)y close(20) END SUBROUTINE GEOMET !================================================= ===================== |
!================================================= =====================
!THIS SUBROUTINE FORMS COEFFS. FOR DISCRETISATION EQNS. !************************************************* ********************* SUBROUTINE COEFF USE GLOBAL IMPLICIT NONE INTEGER :: I,J,K,N REAL :: ZERO,REL,FL,FLM,FLP,GM,GMM,VOL,APT,AREA REAL :: SXT,SXB,ARHO,TMX,UMX,VMX,DELT,DELU,DELV,BTIME,BBK CHARACTER(LEN=80) :: fname !**************** COEFFICIENTS FOR THE U EQUATION ************** ITERL=1 ZERO=0.0 lconv = .false. 99 continue CALL RESET !---------------------------------------------------------------------- NF=1 IF(.NOT.LSOLVE(NF)) GO TO 500 IST=3 JST=2 CALL GAMSOR REL=1.-RELAX(NF) DO 102 I=3,L2 FL=XCVI(I)*V(I,2)*RHO(I,1) FLM=XCVIP(I-1)*V(I-1,2)*RHO(I-1,1) FLOW=R(1)*(FL+FLM) DIFF=R(1)*(XCVI(I)*GAM(I,1)+XCVIP(I-1)*GAM(I-1,1))/YDIF(2) CALL DIFLOW 102 AJM(I,2)=ACOF+AMAX1(ZERO,FLOW) DO 103 J=2,M2 FLOW=ARX(J)*U(2,J)*RHO(1,J) DIFF=ARX(J)*GAM(1,J)/(XCV(2)*SX(J)) CALL DIFLOW AIM(3,J)=ACOF+AMAX1(ZERO,FLOW) DO 103 I=3,L2 IF(I.EQ.L2) GO TO 104 FL=U(I,J)*(FX(I)*RHO(I,J)+FXM(I)*RHO(I-1,J)) FLP=U(I+1,J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J) ) FLOW=ARX(J)*0.5*(FL+FLP) DIFF=ARX(J)*GAM(I,J)/(XCV(I)*SX(J)) GO TO 105 104 FLOW=ARX(J)*U(L1,J)*RHO(L1,J) DIFF=ARX(J)*GAM(L1,J)/(XCV(L2)*SX(J)) 105 CALL DIFLOW AIM(I+1,J)=ACOF+AMAX1(ZERO,FLOW) AIP(I,J)=AIM(I+1,J)-FLOW IF(J.EQ.M2) GO TO 106 FL=XCVI(I)*V(I,J+1)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*R HO(I,J)) FLM=XCVIP(I-1)*V(I-1,J+1)*(FY(J+1)*RHO(I-1,J+1)+FYM(J+1)* & & RHO(I-1,J)) GM=GAM(I,J)*GAM(I,J+1)/(YCV(J)*GAM(I,J+1)+YCV(J+1)*GAM(I,J)+ & & 1.0E-30)*XCVI(I) GMM=GAM(I-1,J)*GAM(I-1,J+1)/(YCV(J)*GAM(I-1,J+1)+YCV(J+1)* & & GAM(I-1,J)+1.E-30)*XCVIP(I-1) DIFF=RMN(J+1)*2.*(GM+GMM) GO TO 107 106 FL=XCVI(I)*V(I,M1)*RHO(I,M1) FLM=XCVIP(I-1)*V(I-1,M1)*RHO(I-1,M1) DIFF=R(M1)*(XCVI(I)*GAM(I,M1)+XCVIP(I-1)*GAM(I-1,M1))/YDIF(M1) 107 FLOW=RMN(J+1)*(FL+FLM) CALL DIFLOW AJM(I,J+1)=ACOF+AMAX1(ZERO,FLOW) AJP(I,J)=AJM(I,J+1)-FLOW VOL=YCVR(J)*XCVS(I) APT=(RHO(I,J)*XCVI(I)+RHO(I-1,J)*XCVIP(I-1)) & &/(XCVS(I)*DT) AP(I,J)=AP(I,J)-APT CON(I,J)=CON(I,J)+APT*UO(I,J) AP(I,J)=(-AP(I,J)*VOL+AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J)) & &/RELAX(NF) CON(I,J)=CON(I,J)*VOL+REL*AP(I,J)*U(I,J) DU(I,J)=VOL/(XDIF(I)*SX(J)) DU(I,J)=DU(I,J)/AP(I,J) 103 CONTINUE !---------------------------------------------------------------------- cofu(:,:,1) = con(:,: ) cofu(:,:,2) = aip(:,: ) cofu(:,:,3) = aim(:,: ) cofu(:,:,4) = ajp(:,: ) cofu(:,:,5) = ajm(:,: ) cofu(:,:,6) = ap(:,: ) !***************** COEFFICIENTS FOR THE V EQUATION *********** NF=2 CALL RESET IST=2 JST=3 CALL GAMSOR REL=1.-RELAX(NF) DO 202 I=2,L2 AREA=R(1)*XCV(I) FLOW=AREA*V(I,2)*RHO(I,1) DIFF=AREA*GAM(I,1)/YCV(2) CALL DIFLOW 202 AJM(I,3)=ACOF+AMAX1(ZERO,FLOW) DO 203 J=3,M2 FL=ARXJ(J)*U(2,J)*RHO(1,J) FLM=ARXJP(J-1)*U(2,J-1)*RHO(1,J-1) FLOW=FL+FLM DIFF=(ARXJ(J)*GAM(1,J)+ARXJP(J-1)*GAM(1,J-1))/(XDIF(2)*SXMN(J)) CALL DIFLOW AIM(2,J)=ACOF+AMAX1(ZERO,FLOW) DO 203 I=2,L2 IF(I.EQ.L2) GO TO 204 FL=ARXJ(J)*U(I+1,J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*R HO(I,J)) FLM=ARXJP(J-1)*U(I+1,J-1)*(FX(I+1)*RHO(I+1,J-1)+FXM(I+1)* & & RHO(I,J-1)) GM=GAM(I,J)*GAM(I+1,J)/(XCV(I)*GAM(I+1,J)+XCV(I+1)*GAM(I,J)+ & & 1.E-30)*ARXJ(J) GMM=GAM(I,J-1)*GAM(I+1,J-1)/(XCV(I)*GAM(I+1,J-1)+XCV(I+1)* & & GAM(I,J-1)+1.0E-30)*ARXJP(J-1) DIFF=2.*(GM+GMM)/SXMN(J) GO TO 205 204 FL=ARXJ(J)*U(L1,J)*RHO(L1,J) FLM=ARXJP(J-1)*U(L1,J-1)*RHO(L1,J-1) DIFF=(ARXJ(J)*GAM(L1,J)+ARXJP(J-1)*GAM(L1,J-1))/(XDIF(L1)*SXMN(J)) 205 FLOW=FL+FLM CALL DIFLOW AIM(I+1,J)=ACOF+AMAX1(ZERO,FLOW) AIP(I,J)=AIM(I+1,J)-FLOW IF(J.EQ.M2) GO TO 206 AREA=R(J)*XCV(I) FL=V(I,J)*(FY(J)*RHO(I,J)+FYM(J)*RHO(I,J-1))*RMN(J) FLP=V(I,J+1)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J) )*RMN(J+1) FLOW=(FV(J)*FL+FVP(J)*FLP)*XCV(I) DIFF=AREA*GAM(I,J)/YCV(J) GO TO 207 206 AREA=R(M1)*XCV(I) FLOW=AREA*V(I,M1)*RHO(I,M1) DIFF=AREA*GAM(I,M1)/YCV(M2) 207 CALL DIFLOW AJM(I,J+1)=ACOF+AMAX1(ZERO,FLOW) AJP(I,J)=AJM(I,J+1)-FLOW VOL=YCVRS(J)*XCV(I) SXT=SX(J) IF(J.EQ.M2) SXT=SX(M1) SXB=SX(J-1) IF(J.EQ.3) SXB=SX(1) APT=(ARXJ(J)*RHO(I,J)*0.5*(SXT+SXMN(J))+ARXJP(J-1)*RHO(I,J-1)* & &0.5*(SXB+SXMN(J)))/(YCVRS(J)*DT) AP(I,J)=AP(I,J)-APT CON(I,J)=CON(I,J)+APT*VO(I,J) AP(I,J)=(-AP(I,J)*VOL+AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J)) & &/RELAX(NF) CON(I,J)=CON(I,J)*VOL+REL*AP(I,J)*V(I,J) DV(I,J)=VOL/YDIF(J) DV(I,J)=DV(I,J)/AP(I,J) 203 CONTINUE !---------------------------------------------------------------------- cofv(:,:,1) = con(:,: ) cofv(:,:,2) = aip(:,: ) cofv(:,:,3) = aim(:,: ) cofv(:,:,4) = ajp(:,: ) cofv(:,:,5) = ajm(:,: ) cofv(:,:,6) = ap(:,: ) !******************* CALCULATE UHAT AND VHAT ******************** DO J=2,M2 DO I=3,L2 UHAT(I,J)=(COFU(I,J,2)*U(I+1,J)+COFU(I,J,3)*U(I-1,J)+COFU(I,J,4) & & *U(I,J+1)+COFU(I,J,5)*U(I,J-1)+COFU(I,J,1))/COFU(I,J,6) END DO END DO DO J=3,M2 DO I=2,L2 VHAT(I,J)=(COFV(I,J,2)*V(I+1,J)+COFV(I,J,3)*V(I-1,J)+COFV(I,J,4) & & *V(I,J+1)+COFV(I,J,5)*V(I,J-1)+COFV(I,J,1))/COFV(I,J,6) END DO END DO !******************** COEFFICIENTS FOR THE PRESSURE EQUATION ***** NF=3 CALL RESET IST=2 JST=2 CALL GAMSOR DO 410 J=2,M2 DO 410 I=2,L2 VOL=YCVR(J)*XCV(I) 410 CON(I,J)=CON(I,J)*VOL DO 402 I=2,L2 ARHO=R(1)*XCV(I)*RHO(I,1) CON(I,2)=CON(I,2)+ARHO*V(I,2) 402 AJM(I,2)=0. DO 403 J=2,M2 ARHO=ARX(J)*RHO(1,J) CON(2,J)=CON(2,J)+ARHO*U(2,J) AIM(2,J)=0. DO 403 I=2,L2 IF(I.EQ.L2) GO TO 404 ARHO=ARX(J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J)) FLOW=ARHO*UHAT(I+1,J) CON(I,J)=CON(I,J)-FLOW CON(I+1,J)=CON(I+1,J)+FLOW AIP(I,J)=ARHO*DU(I+1,J) AIM(I+1,J)=AIP(I,J) GO TO 405 404 ARHO=ARX(J)*RHO(L1,J) CON(I,J)=CON(I,J)-ARHO*U(L1,J) AIP(I,J)=0. 405 IF(J.EQ.M2) GO TO 406 ARHO=RMN(J+1)*XCV(I)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)* RHO(I,J)) FLOW=ARHO*VHAT(I,J+1) CON(I,J)=CON(I,J)-FLOW CON(I,J+1)=CON(I,J+1)+FLOW AJP(I,J)=ARHO*DV(I,J+1) AJM(I,J+1)=AJP(I,J) GO TO 407 406 ARHO=RMN(M1)*XCV(I)*RHO(I,M1) CON(I,J)=CON(I,J)-ARHO*V(I,M1) AJP(I,J)=0. 407 AP(I,J)=AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J) 403 continue DO 421 J=2,M2 DO 421 I=2,L2 AP(I,J)=AP(I,J)/RELAX(NP) 421 CON(I,J)=CON(I,J)+(1.0-RELAX(NP))*AP(I,J)*P(I,J) !---------------------------------------------------------------------- cofp(:,:,2) = aip(:,: ) cofp(:,:,3) = aim(:,: ) cofp(:,:,4) = ajp(:,: ) cofp(:,:,5) = ajm(:,: ) NF=NP CALL SOLVE !******************** COMPUTE U* AND V* ************************** NF=1 IST=3 JST=2 DO 550 K=1,6 DO 550 J=JST,M2 DO 550 I=IST,L2 550 COF(I,J,K)=COFU(I,J,K) !---------------------------------------------------------------------- con(:,: ) = cof(:,:,1) aip(:,: ) = cof(:,:,2) aim(:,: ) = cof(:,:,3) ajp(:,: ) = cof(:,:,4) ajm(:,: ) = cof(:,:,5) ap(:,: ) = cof(:,:,6) !---------------------------------------------------------------------- DO 551 J=JST,M2 DO 551 I=IST,L2 551 CON(I,J)=CON(I,J)+DU(I,J)*AP(I,J)*(P(I-1,J)-P(I,J)) CALL SOLVE NF=2 IST=2 JST=3 DO 552 K=1,6 DO 552 J=JST,M2 DO 552 I=IST,L2 552 COF(I,J,K)=COFV(I,J,K) !---------------------------------------------------------------------- con(:,: ) = cof(:,:,1) aip(:,: ) = cof(:,:,2) aim(:,: ) = cof(:,:,3) ajp(:,: ) = cof(:,:,4) ajm(:,: ) = cof(:,:,5) ap(:,: ) = cof(:,:,6) !---------------------------------------------------------------------- DO 553 J=JST,M2 DO 553 I=IST,L2 553 CON(I,J)=CON(I,J)+DV(I,J)*AP(I,J)*(P(I,J-1)-P(I,J)) CALL SOLVE !*************** COEFFICIENTS FOR THE PRESSURE CORRECTION EQUATION *** NF=3 CALL RESET IST=2 JST=2 DO 554 K=2,5 DO 554 J=JST,M2 DO 554 I=IST,L2 554 COF(I,J,K)=COFP(I,J,K) !---------------------------------------------------------------------- aip(:,: ) = cof(:,:,2) aim(:,: ) = cof(:,:,3) ajp(:,: ) = cof(:,:,4) ajm(:,: ) = cof(:,:,5) !---------------------------------------------------------------------- CALL GAMSOR SMAX=0. SSUM=0. DO 510 J=2,M2 DO 510 I=2,L2 VOL=YCVR(J)*XCV(I) 510 CON(I,J)=CON(I,J)*VOL DO 502 I=2,L2 ARHO=R(1)*XCV(I)*RHO(I,1) CON(I,2)=CON(I,2)+ARHO*V(I,2) 502 END DO DO 503 J=2,M2 ARHO=ARX(J)*RHO(1,J) CON(2,J)=CON(2,J)+ARHO*U(2,J) DO 503 I=2,L2 IF(I.EQ.L2) GO TO 504 ARHO=ARX(J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J)) FLOW=ARHO*U(I+1,J) CON(I,J)=CON(I,J)-FLOW CON(I+1,J)=CON(I+1,J)+FLOW GO TO 505 504 ARHO=ARX(J)*RHO(L1,J) CON(I,J)=CON(I,J)-ARHO*U(L1,J) 505 IF(J.EQ.M2) GO TO 506 ARHO=RMN(J+1)*XCV(I)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)* RHO(I,J)) FLOW=ARHO*V(I,J+1) CON(I,J)=CON(I,J)-FLOW CON(I,J+1)=CON(I,J+1)+FLOW GO TO 507 506 ARHO=RMN(M1)*XCV(I)*RHO(I,M1) CON(I,J)=CON(I,J)-ARHO*V(I,M1) 507 AP(I,J)=AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J) PC(I,J)=0. SMAX=AMAX1(SMAX,ABS(CON(I,J))) SSUM=SSUM+CON(I,J) 503 CONTINUE CALL SOLVE !*************** COME HERE TO CORRECT THE VELOCITIES ************ DO 521 J=2,M2 DO 521 I=2,L2 IF(I.NE.2) U(I,J)=U(I,J)+DU(I,J)*(PC(I-1,J)-PC(I,J)) IF(J.NE.2) V(I,J)=V(I,J)+DV(I,J)*(PC(I,J-1)-PC(I,J)) 521 CONTINUE 500 CONTINUE !************** COEFFICIENTS FOR OTHER EQUATIONS **************** CALL RESET IST=2 JST=2 DO 600 N=4,NFMAX NF=N IF(.NOT.LSOLVE(NF)) GO TO 600 CALL GAMSOR REL=1.-RELAX(NF) DO 602 I=2,L2 AREA=R(1)*XCV(I) FLOW=AREA*V(I,2)*RHO(I,1) DIFF=AREA*GAM(I,1)/YDIF(2) CALL DIFLOW 602 AJM(I,2)=ACOF+AMAX1(ZERO,FLOW) DO 603 J=2,M2 FLOW=ARX(J)*U(2,J)*RHO(1,J) DIFF=ARX(J)*GAM(1,J)/(XDIF(2)*SX(J)) CALL DIFLOW AIM(2,J)=ACOF+AMAX1(ZERO,FLOW) DO 603 I=2,L2 IF(I.EQ.L2) GO TO 604 FLOW=ARX(J)*U(I+1,J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)* RHO(I,J)) DIFF=ARX(J)*2.*GAM(I,J)*GAM(I+1,J)/((XCV(I)*GAM(I+1,J)+ & & XCV(I+1)*GAM(I,J)+1.0E-30)*SX(J)) GO TO 605 604 FLOW=ARX(J)*U(L1,J)*RHO(L1,J) DIFF=ARX(J)*GAM(L1,J)/(XDIF(L1)*SX(J)) 605 CALL DIFLOW AIM(I+1,J)=ACOF+AMAX1(ZERO,FLOW) AIP(I,J)=AIM(I+1,J)-FLOW AREA=RMN(J+1)*XCV(I) IF(J.EQ.M2) GO TO 606 FLOW=AREA*V(I,J+1)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RH O(I,J)) DIFF=AREA*2.*GAM(I,J)*GAM(I,J+1)/(YCV(J)*GAM(I,J+1)+ & & YCV(J+1)*GAM(I,J)+1.0E-30) GO TO 607 606 FLOW=AREA*V(I,M1)*RHO(I,M1) DIFF=AREA*GAM(I,M1)/YDIF(M1) 607 CALL DIFLOW AJM(I,J+1)=ACOF+AMAX1(ZERO,FLOW) AJP(I,J)=AJM(I,J+1)-FLOW 603 CONTINUE !---------------------------------------------------------------------- DO 349 I=1,L1 DO 349 J=1,M1 VOL=YCVR(J)*XCV(I) APT=RHO(I,J)/DT AP(I,J)=AP(I,J)-APT if(nf.eq.4) CON(I,J)=CON(I,J)+APT*TO(I,J) AP(I,J)=(-AP(I,J)*VOL+AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J)) & &/RELAX(NF) if(nf.eq.4) CON(I,J)=CON(I,J)*VOL+REL*AP(I,J)*t(I,J) if(nf.eq.np) CON(I,J)=CON(I,J)*VOL+REL*AP(I,J)*p(I,J) !---------------------------------------------------------------------- if(nf.eq.4) TOLD(i,j)=T(i,j) !---------------------------------------------------------------------- 349 CONTINUE CALL SOLVE !---------------------------------------------------------------------- 600 END DO !*****************---------------CONVERGENCE CRITERIA---------------------------- tmx = maxval(t) umx = maxval(u) vmx = maxval(v) delt = maxval(abs(t-told)) delu = maxval(abs(u-uold)) delv = maxval(abs(v-vold)) if(tmx.gt.0.0) then deltmx = delt/tmx else deltmx=0.0 endif if(umx.gt.0.0) then delumx = delu/umx else delumx=0.0 endif if(vmx.gt.0.0) then delvmx = delv/vmx else delvmx=0.0 endif if(time.lt.120.0) then WRITE(*,*)'\n',iterl,'tmx=',tmx,'delt=',deltmx,'um x=',umx,'delu=',delumx,'vmx=',vmx,'delv=',delvmx endif if(iterl.ge.istp) then if((delt.lt.0.15).and.(delu.lt.0.005)) then lconv = .true. endif endif if(lconv) goto 59 IF(DELTMX.GT.ERT) GOTO 51 IF(DELUMX.GT.ERU) GOTO 51 IF(DELVMX.GT.ERV) GOTO 51 write(*,*)'converged' goto 59 51 IF(ITERL.EQ.ISTP) THEN WRITE(*,*)'exceeded max. no. of iteration' WRITE(*,*)'\n',iterl,'tmx=',tmx,'delt=',deltmx,'um x=',umx,'delu=',delumx,'vmx=',vmx,'delv=',delvmx GOTO 59 STOP ENDIF 40 FORMAT(I5,6F10.6,2X,2F12.4) ITERL=ITERL+1 DO 765 I=2,L1 DO 765 J=2,M1 UOLD(I,J)=U(I,J) VOLD(I,J)=V(I,J) 765 CONTINUE CALL BOUND GOTO 99 59 ITLL=ITLL+1 BTIME=TIME+DT write(*,*)'time=',time,'end of inner iterations' fname="velocity.dat" open(UNIT=10,file=fname) write(10,*)u write(10,*)v close(10) TIME=TIME+DT !plotting every 50 timestep IF(MOD(INT(TIME/DT),50).eq.0) THEN CALL INTERPOLATION ENDIF CALL ENRBAL BBK=AMOD(TIME,real(NTIME)) IF(BBK.LE.0.00001)THEN ITER=ITER+1 CALL BOUND END IF IF(TIME.GE.TLAST) LSTOP=.TRUE. DO 775 I=2,L1 DO 775 J=2,M1 UOLD(I,J)=U(I,J) VOLD(I,J)=V(I,J) 775 CONTINUE END SUBROUTINE COEFF !================================================= ===================== |
!================================================= =====================
!this subroutine stores initial values for the coming timestep !************************************************* ********************* subroutine oldval use global implicit none integer :: I,j do i=1,L1 do j=1,m1 to(i,j)=t(i,j) uo(i,j)=u(i,j) vo(i,j)=v(i,j) end do end do end subroutine oldval !================================================= ===================== |
!================================================= =====================
!this subroutine generates uniform grid !************************************************* ********************* subroutine ugrid use global implicit none integer :: I,j real :: Dx,dy xu(2)=0. Dx=xl/float(l1-2) do 1 i=3,l1 1 xu(i)=xu(i-1)+dx yv(2)=0. Dy=yl/float(m1-2) do 2 j=3,m1 2 yv(j)=yv(j-1)+dy end subroutine ugrid |
!================================================= =====================
!THIS SUBROTINE CHECKS OVERALL ENERGY BALANCE !************************************************* ********************* SUBROUTINE ENRBAL USE GLOBAL IMPLICIT NONE END SUBROUTINE ENRBAL !================================================= ===================== !THIS SUBROUTINE GIVES INTERPOLATED VELOCITIES OF PROBLEM VARIABLES !************************************************* ********************* SUBROUTINE INTERPOLATION USE GLOBAL do j=1,m1 u(1,j)=u(2,j) end do do i=1,l1 v(i,1)=v(i,2) end do do i=2,l2 do j=1,m1 ur(i,j)=0.5*(u(i+1,j)+u(i,j)) end do end do do i=1,l1 do j=2,m2 vr(i,j)=0.5*(v(i,j+1)+v(i,j)) end do end do CALL PRINTOUT END SUBROUTINE INTERPOLATION !================================================= ===================== !THIS SUBROUTINE DIRECTLY GIVES MATLAB PLOTS OF PROBLEM VARIABLES !************************************************* ********************* SUBROUTINE PRINTOUT USE GLOBAL character::filename*20,fname*20,fname1*20,fname2*2 0 integer:: ileunit,ile,ile1 !print*,'************test',MOD(INT(TIME/DT),20) !plotting every 50 timestep IF(MOD(INT(TIME/DT),50).ne.0) return !CALCULATE THE STREAM FUNCTION--------------------------------------- pc(2,2)=0. DO I=2,L1 IF(I.NE.2) pc(I,2)=pc(I-1,2)-RHO(I-1,1)*V(I-1,2)*R(1)*XCV(I-1) DO J=3,M1 RHOM=FX(I)*RHO(I,J-1)+FXM(I)*RHO(I-1,J-1) pc(I,J)=pc(I,J-1)+RHOM*U(I,J-1)*ARX(J-1) ENDDO ENDDO !if(time.le.2500) goto700 !goto 711 if(int(time).eq.500) filename='out500.m' open(unit=55,file=FileName) write(55,FMT='(A3,<l1>(E15.8,1X),A2)') "x=[" , x , "];" write(55,FMT='(A3,<m1>(E15.8,1X),A2)') "y=[" , y , "];" write(55,*) "t=[" do j=1,m1 write(55,FMT='(<l1>(E15.8,1X),A1)') t(:,j), ";" end do write(55,*) "];" write(55,*) "pc=[" do j=1,m1 write(55,FMT='(<l1>(E15.8,1X),A1)') pc(:,j), ";" end do write(55,*) "];" write(55,*) "u=[" do j=1,m1 write(55,FMT='(<l1>(E15.8,1X),A1)') ur(:,j), ";" end do write(55,*) "];" write(55,*) "v=[" do j=1,m1 write(55,FMT='(<l1>(E15.8,1X),A1)') vr(:,j), ";" end do write(55,*) "];" write(55,*) "figure(1)" write(55,*) "contour(x,y,t);" write(55,*) "figure(2)" write(55,*) "quiver(x,y,u,v)" write(55,*) "figure(3)" write(55,*) "contour(x,y,pc)" write(55,*) "title('Pictures at time=",time,"seconds');" close(55) ! 711 continue if(ltecplot) then !techplot plotting ileunit=int(time) ile=int(time) write(fname,1000)ileunit 1000 format('temp.',i5.5) open(UNIT=ileunit,FILE=fname) !directs to an output file write(ileunit,*) "VARIABLES = 'X','Y','T'" write(ileunit,*) "ZONE I=",l1,"J=",m1 do j=1,m1 do i=1,l1 write(ileunit,702)(x(i)),(y(j)),(t(i,j)) 702 format(1x,f8.4,2x,f8.4,2x,f8.4) enddo enddo close(ileunit) !techplot plotting of velocity write(fname1,1001)ile 1001 format('vel.',i5.5) OPEN(UNIT=ile,FILE=fname1) write(ile,*) "VARIABLES = 'X','Y','U','V'" write(ile,*) "ZONE I=",l1,"J=",m1 do j=1,m1 do i=1,l1 write(ile,701)(x(i)),(y(j)),(ur(i,j)),(vr(i,j)) 701 format(1x,f8.4,2x,f8.4,2x,f8.4,2x,f8.4) ! 701 format((1x,4f10.2),A1) enddo enddo close(ile) ile1=int(time) write(fname2,1002)ile1 1002 format('stream.',i5.5) open(UNIT=ile1,FILE=fname2) !directs to an output file write(ile1,*) "VARIABLES = 'X','Y','PC'" write(ile1,*) "ZONE I=",l1,"J=",m1 do j=1,m1 do i=1,l1 write(ile1,703)(x(i)),(y(j)),(pc(i,j)) 703 format(1x,f8.4,2x,f8.4,2x,f8.4) enddo enddo close(ile1) endif END SUBROUTINE PRINTOUT !================================================= ===================== !================= USER PART BEGINS HERE ========================== !================================================= ===================== !================================================= ===================== !THIS SUBROUTINE INITIALISES DIFFERENT VARIABLES AND READS INPUT DATA !************************************************* ********************* SUBROUTINE GRID USE GLOBAL IMPLICIT NONE INTEGER :: I,J IPREF=1 JPREF=1 MODE=1 RELAX =0.6 LSOLVE = .FALSE. LSOLVE(1)=.TRUE. LSOLVE(2)=.TRUE. LSOLVE(3)=.TRUE. LSOLVE(4)=.TRUE. LSOLVE(11)=.TRUE. LSTOP=.FALSE. DO I=1,NFX3 NTIMES(I)=1 LBLK(I)=.TRUE. END DO TIME=0. TITLE(1)='UVEL' TITLE(2)='VVEL' TITLE(4)='TEMP' TITLE(6)='CONC' deltmx = 1.0 !================================================= ======================= !=========== meaning of symbols: !This is a sample data file for free convection in a rect. cavity !Symbol meanings are: !xl : x-direction domain length !yl : y-direction domain length !l1 :no. of grids in x-dirn !m1 :no. of grids in y-dirn !ngrid : switch for grid type; if ngrid=1 then uniform grid, else non-uniform !tk : thermal conductivity !g : acceleration due to gravity !beta : volumetric expansion coefficient !cp : specific heat at const. press. !tref : ref. temp. for buoyancy source term !amiu : dynamic viscosity !rhocon : density !flux : heat flux at boundary !relax(nf) : relaxation factor for a particular nf !dt : time step size to start with !tlast : total time domain !time2 : time till the time step value of dt is valid !dt2 : time step value after time t2; similarly dt3..... !eru,erv,ert: convergence tolerances for u,v,t (relative error) !istp : maximum no. of iterations in each timestep !ntimes(nf) : no. of times corr. to which eqn corr. to a particular nf is solved !tinitial :initial value of temp.for a transient problem !rhoinitial :initial value of density for a transient problem !***************** INPUT DATA HERE *************************** XL=0.22 ; YL=0.11 ; L1=40;M1=40 NTIME=1 ; NGRID=1 ; ISTP=500 tk=0.62 ; cp=4179.0 rhocon=997.0 ;RHOINITIAL=997.0 DT=10.0 ; TLAST=510.0 ERU=0.00001; ERV=0.00001 ; ERT=0.000001 tinitial=23.0 ;TREF=23.0; BETA=0.0002761 ; AMIU=0.000855 ;FLUX=500.0 g=9.81; relax(1)=0.5; relax(2)=0.5; relax(3)=0.5; relax(4)=0.5 ltecplot=.false. !---------------------------------------------------------------------- ITER=0 R(1)=0.0 CALL ugrid !************ CON,AP,U,V,RHO,PC AND P ARRAYS ARE INITIALIZED HERE ** DO J=1,M1 DO I=1,L1 U(I,J)=0. V(I,J)=0. PC(I,J)=0. CONC(i,j)=0. P(I,J)=0. RHO(I,J)= rhocon DU(I,J)=0. DV(I,J)=0. UOLD(I,J)=0. VOLD(I,J)=0. AIP(I,J)=0. AIM(I,J)=0. AJM(I,J)=0. AJP(I,J)=0. CON(I,J)=0. AP(I,J)=0. AP0(I,J)=0. told(i,j)=0. END DO END DO cofu = 0.0 cofv = 0.0 cof = 0.0 cofp = 0.0 END SUBROUTINE grid !====================================================================== !THIS SUBROUTINE GIVES INITIAL CONDITIONS FOR THE PROBLEM !************************************************* ********************* SUBROUTINE START USE GLOBAL IMPLICIT NONE INTEGER :: I,J !real:;temp DO I=1,L1 DO J=1,M1 U(I,J)=0.0 V(I,J)=0.0 T(I,J)=TINITIAL TOLD(I,J)=tinitial rho(i,j)=rhocon END DO END DO END SUBROUTINE start !================================================= ===================== !THIS SUBROUTINE CALCULATES THE DENSITY !************************************************* ********************* SUBROUTINE dense USE GLOBAL IMPLICIT NONE INTEGER :: I,J DO I=1,L1 DO J=1,M1 rho(i,j)=rhoinitial END DO END DO END SUBROUTINE DENSE !================================================= ===================== !THIS SUBROUTINE GIVES BOUNDARY CONDITIONS FOR THE PROBLEM !************************************************* ********************* SUBROUTINE bound USE GLOBAL IMPLICIT NONE INTEGER::I,J !U-VELOCITY BOUNDARY CONDITION DO I=2,L1 DO J=1,M1 !TOP & BOTTOM u(i,m1)=0.0 U(I,1)=0.0 !SIDES U(2,J)=0.0 U(L1,J)=0.0 END DO END DO DO I=1,L1 DO J=2,M1 V(1,J)=0.0 V(L1,J)=0.0 V(I,2)=0.0 V(I,M1)=0.0 END DO END DO !TEMPERATURE BOUNDARY CONDITION DO J=1,M1 !LEFT T(1,J)=T(2,J)+FLUX*XDIF(2)/TK !SIDES T(L1,J)=T(L2,J)+(FLUX*XDIF(L1)/TK) ENDDO DO I=1,L1 !BOTTOM AND TOP T(I,1)=T(I,2) t(i,m1)=t(i,m2) END DO END SUBROUTINE bound !================================================= ===================== !THIS SUBROUTINE GIVES DIFF. COEFF. & SOURCE TERMS FOR DISCRETISED EQNS. !************************************************* ********************* SUBROUTINE gamsor USE GLOBAL IMPLICIT NONE INTEGER :: I,J REAL:: TAV DO I=1,L1 DO J=1,M1 IF(nf.EQ.4) gam(i,j)=tk/cp IF(NF/=4) GAM(I,J)=AMIU END DO END DO IF(NF==2) THEN DO I=2,L2 DO J=3,M2 Tav=FY(J)*T(I,J)+FYM(J)*T(I,J-1) CON(I,J)=CON(I,J)+rhoinitial*g*BETA*(Tav-tref) END DO END DO ENDIF if(nf.eq.4) then do i=1,l1 gam(i,m1)=0.0 enddo endif END SUBROUTINE gamsor !================================================= ===================== !******************** END OF THE CODE ************************ !================================================= ===================== |
wtf? what is the purpose of this post?
|
Dear sir I need a fortran subroutine to solve melting and solidification by using enthalpy method ...finite volume or finite difference method ..thanks
|
You can do melting and solidification using this code.
U just need to add source term and change the boundary condition here in this code. The rest all are same. Read Patankar Book to understand this code. best of luck |
Please do not post code like this. It is not a good way to share code.
You can put the code here (assuming you have the rights to do it) https://www.cfd-online.com/Wiki/Source_code_archive and make a new post and link to it. |
Ok.
Thank you sir. I will do it Sir, I have converted this code for spherical coordinate. Work title"Flow over a sphere under moderate Reynolds number" My U-velocity is converging but V-velocity is not converging. What would be the reason? |
which code do you mean sir ?....actually I want to simulate a heat transfer through wall subjected to a solar radiation and heat convection.
the wall consists of three layers : the outer layer is a mortar followed by a layer of brick +phase change material and the inner layer is a mortar. I want to do this simulation using FORTRAN 90 code using enthalpy -finite volume method so if you help me I'll be very grateful to you. |
See, first u try to simulate a cavity filled with phased change material.
Upper and lower side of the cavity are insulated. Left side temperature is higher than the melting point temp and right side is kept at zero degree. |
Quote:
|
Hi Sir,
Do you know how to convert Fluent mesh file format to Fortran FVM use? |
| All times are GMT -4. The time now is 22:45. |