CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > General Forums > Main CFD Forum

Patankar CFD FORTRAN 90 Code FVM

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree3Likes
  • 1 Post By siddiquesil
  • 1 Post By vesp
  • 1 Post By siddiquesil

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   February 6, 2019, 12:52
Default Patankar CFD FORTRAN 90 Code FVM
  #1
New Member
 
Numan Siddique
Join Date: Jan 2019
Location: India
Posts: 24
Rep Power: 4
siddiquesil is on a distinguished road
!===============================================
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
optimux likes this.

Last edited by siddiquesil; February 6, 2019 at 14:16.
siddiquesil is offline   Reply With Quote

Old   February 6, 2019, 12:55
Default
  #2
New Member
 
Numan Siddique
Join Date: Jan 2019
Location: India
Posts: 24
Rep Power: 4
siddiquesil is on a distinguished road
!################################################# #################
!================================================= =================
!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
!================================================= =================
siddiquesil is offline   Reply With Quote

Old   February 6, 2019, 12:57
Default
  #3
New Member
 
Numan Siddique
Join Date: Jan 2019
Location: India
Posts: 24
Rep Power: 4
siddiquesil is on a distinguished road
!================================================= =====================
!***************** 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

!================================================= =====================
siddiquesil is offline   Reply With Quote

Old   February 6, 2019, 12:59
Default
  #4
New Member
 
Numan Siddique
Join Date: Jan 2019
Location: India
Posts: 24
Rep Power: 4
siddiquesil is on a distinguished road
!===============================================
!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

!===============================================
siddiquesil is offline   Reply With Quote

Old   February 6, 2019, 13:03
Default
  #5
New Member
 
Numan Siddique
Join Date: Jan 2019
Location: India
Posts: 24
Rep Power: 4
siddiquesil is on a distinguished road
!================================================= =====================
!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

!================================================= =====================
siddiquesil is offline   Reply With Quote

Old   February 6, 2019, 13:04
Default
  #6
New Member
 
Numan Siddique
Join Date: Jan 2019
Location: India
Posts: 24
Rep Power: 4
siddiquesil is on a distinguished road
!================================================= =====================
!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

!================================================= =====================
siddiquesil is offline   Reply With Quote

Old   February 6, 2019, 13:06
Default
  #7
New Member
 
Numan Siddique
Join Date: Jan 2019
Location: India
Posts: 24
Rep Power: 4
siddiquesil is on a distinguished road
!================================================= =====================
! 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

!================================================= =====================
siddiquesil is offline   Reply With Quote

Old   February 6, 2019, 13:10
Default
  #8
New Member
 
Numan Siddique
Join Date: Jan 2019
Location: India
Posts: 24
Rep Power: 4
siddiquesil is on a distinguished road
!================================================= =====================
!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

!================================================= =====================
siddiquesil is offline   Reply With Quote

Old   February 6, 2019, 13:12
Default
  #9
New Member
 
Numan Siddique
Join Date: Jan 2019
Location: India
Posts: 24
Rep Power: 4
siddiquesil is on a distinguished road
!================================================= =====================
!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

!================================================= =====================
siddiquesil is offline   Reply With Quote

Old   February 6, 2019, 13:13
Default
  #10
New Member
 
Numan Siddique
Join Date: Jan 2019
Location: India
Posts: 24
Rep Power: 4
siddiquesil is on a distinguished road
!================================================= =====================
!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
siddiquesil is offline   Reply With Quote

Old   February 6, 2019, 13:17
Default
  #11
New Member
 
Numan Siddique
Join Date: Jan 2019
Location: India
Posts: 24
Rep Power: 4
siddiquesil is on a distinguished road
!================================================= =====================
!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 ************************
!================================================= =====================
siddiquesil is offline   Reply With Quote

Old   February 6, 2019, 13:41
Default
  #12
Member
 
Join Date: Aug 2018
Posts: 74
Rep Power: 4
vesp is on a distinguished road
wtf? what is the purpose of this post?
sbaffini likes this.
vesp is offline   Reply With Quote

Old   January 2, 2020, 17:44
Default
  #13
New Member
 
Hayder Mohammed
Join Date: Apr 2013
Posts: 11
Rep Power: 10
Hayder Mohammed is on a distinguished road
Send a message via MSN to Hayder Mohammed Send a message via Yahoo to Hayder Mohammed Send a message via Skype™ to Hayder Mohammed
Dear sir I need a fortran subroutine to solve melting and solidification by using enthalpy method ...finite volume or finite difference method ..thanks
Hayder Mohammed is offline   Reply With Quote

Old   January 2, 2020, 22:22
Default
  #14
New Member
 
Numan Siddique
Join Date: Jan 2019
Location: India
Posts: 24
Rep Power: 4
siddiquesil is on a distinguished road
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
Hayder Mohammed likes this.
siddiquesil is offline   Reply With Quote

Old   January 2, 2020, 23:01
Default
  #15
Super Moderator
 
Praveen. C
Join Date: Mar 2009
Location: Bangalore
Posts: 336
Blog Entries: 6
Rep Power: 15
praveen is on a distinguished road
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.
praveen is offline   Reply With Quote

Old   January 3, 2020, 02:20
Default
  #16
New Member
 
Numan Siddique
Join Date: Jan 2019
Location: India
Posts: 24
Rep Power: 4
siddiquesil is on a distinguished road
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?
siddiquesil is offline   Reply With Quote

Old   January 3, 2020, 04:18
Default
  #17
New Member
 
Hayder Mohammed
Join Date: Apr 2013
Posts: 11
Rep Power: 10
Hayder Mohammed is on a distinguished road
Send a message via MSN to Hayder Mohammed Send a message via Yahoo to Hayder Mohammed Send a message via Skype™ to Hayder Mohammed
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.
Hayder Mohammed is offline   Reply With Quote

Old   January 3, 2020, 04:25
Default
  #18
New Member
 
Numan Siddique
Join Date: Jan 2019
Location: India
Posts: 24
Rep Power: 4
siddiquesil is on a distinguished road
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.
siddiquesil is offline   Reply With Quote

Old   January 7, 2020, 20:25
Default
  #19
New Member
 
bokholef fouqd
Join Date: Jan 2020
Posts: 1
Rep Power: 0
amcic1996 is on a distinguished road
Quote:
Originally Posted by siddiquesil View Post
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 th pornjk porn800 redtubee reason?
Please go through this brief instruction before attempting anything with this code

Last edited by amcic1996; January 11, 2020 at 09:08.
amcic1996 is offline   Reply With Quote

Reply

Tags
code, fortran 90, fvm, patankar

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
CFD Salary CFD Main CFD Forum 17 January 3, 2017 17:09
A simple CFD code for teaching basic CFD? Christoph Lund Main CFD Forum 13 September 14, 2005 04:36
CFD code structure (F90) ma Main CFD Forum 4 January 10, 2005 20:47
OOP for CFD code Jongtae Kim Main CFD Forum 26 October 20, 2000 06:11
cfd job Dr. Don I anyanwu Main CFD Forum 20 May 17, 1999 15:13


All times are GMT -4. The time now is 15:32.