Non-staggered Cavity Simulation

!********NAVID SALARVAND*********! !EASTERN MEDITERRANEAN UNIVERSITY! PROGRAM CAVITY NON-STAGERRED IMPLICIT NONE INTEGER :: i, j, Count=0, k INTEGER, PARAMETER :: Nx = 100, Ny = 100 REAL(8) :: ResidualU, ResidualV, ResidualPc REAL(8), PARAMETER :: Ly = 1.0, Lx = 1.0 REAL(8), PARAMETER :: Dx = Lx/(Nx-2), Dy = Ly/(Ny-2) REAL(8), PARAMETER :: RE = 1000.0 REAL(8), PARAMETER :: epsilon = 1e-8 REAL(8), PARAMETER :: Ul = 1.0, AlphaU = 0.7, AlphaV= 0.7, AlphaP = 0.3 REAL(8), DIMENSION(Nx, Ny) :: X, Y REAL(8), DIMENSION(Nx, Ny) :: U=0.0, V=0.0, Pc=0.0 REAL(8), DIMENSION(Nx, Ny) :: ApU, ApV, Ap, An, Aw, As, Ae, SuU, SuV, Sp, & & Ape, Apn, Aps, Apw, P=0.0, ScU, ScV REAL(8), DIMENSION(Nx, Ny) :: Uold=0.0, Vold=0.0 REAL(8), DIMENSION(Nx, Ny) :: ApPc, AnPc, AwPc, AsPc, AePc, SuPc, SpPc REAL(8), DIMENSION(Nx, Ny) :: FluxU=0, FluxV=0, FluxUold=0, FluxVold=0 REAL(8), DIMENSION(Nx, Ny) :: FracE, FracW, FracN, FracS, Dn, Dw, Ds, De, & & DpU, DpV, PcfaceH, PcfaceV REAL(8), DIMENSION(Nx, Ny) :: PfaceH=0, PfaceV=0 REAL(8), DIMENSION(Nx, Ny) :: AtermU=0, StermU=0 , PtermU=0, OtermU=0 REAL(8), DIMENSION(Nx, Ny) :: AtermV=0, StermV=0 , PtermV=0, OtermV=0 CALL FILES CALL MESH(Nx, Ny, Dx, Dy, Lx, Ly, X, Y) CALL FRAC(Dx, Dy, Nx, Ny, X, Y, FracE, FracW, FracN, FracS) DO k=1,4000 Count = Count + 1 CALL PFACE(Nx, Ny, P, FracE, FracN, PfaceH, PfaceV) CALL BOUNDARIES(Ul, Nx, Ny, ApU, ApV, An, Aw, As, Ae, SuU, SuV, Sp, SuPc, ScU, ScV, & & AnPc, AwPc, AsPc, AePc) CALL COEFFICIENTS(X, Y, Nx, Ny, Dn, Dw, Ds, De, AlphaU, & & AlphaV, AlphaP, RE, Dx, Dy, & & ApU, ApV,Ap, An, Aw, As, Ae, SuU, SuV, Sp, & & ScU, ScV, Uold, Vold, Ape, Apn, & & Apw, Aps, ApPc, AnPc, AwPc, AsPc, AePc, SuPc, & & PfaceH, PfaceV, FluxU, FluxV, & & FracE, FracW, FracN, FracS) CALL SOLVER (Nx, Ny, ApU, An, Aw, As, Ae, SuU, U) CALL SOLVER (Nx, Ny, ApV, An, Aw, As, Ae, SuV, V) CALL FLUX(Nx, Ny, AlphaU, AlphaV, FracE, FracN, Dx, Dy, X, Y, Ape, Apn, & & FluxUold, FluxVold, ScU, ScV & & , PfaceH, PfaceV, FluxU, FluxV) Pc=0 PcfaceH=0 PcfaceV=0 DO i=1, 8 CALL SOLVER (Nx, Ny, ApPc, AnPc, AwPc, AsPc, AePc, SuPc, Pc) END DO CALL CORRECTION(Nx, Ny, PcfaceH, PcfaceV, FracE, FracN, & & Pc, Dn, Dw, Ds, De, FluxUold, FluxVold, Uold, Vold, & & FluxU, FluxV, U, V, P) CALL RESIDUALS(U, Aw, As, Ae, An, ResidualU, ApU, Nx, Ny, SuU) CALL RESIDUALS(V, Aw, As, Ae, An, ResidualV, ApV, Nx, Ny, SuV) CALL RESIDUALS(Pc, AwPc, AsPc, AePc, AnPc, ResidualPc, ApPc, Nx, Ny, SuPc)

1

WRITE(2,*), count, ResidualU WRITE(4,*), count, ResidualV WRITE(6,*), count, ResidualPc print*, count, ResidualU, ResidualV, ResidualPc IF ((ResidualU + ResidualV + ResidualPc)<epsilon) EXIT END DO CALL WRITES !**********SUBROUTINES & FUNCTIONS**********! CONTAINS !***SUBROUTINES***! SUBROUTINE FILES OPEN OPEN OPEN OPEN OPEN OPEN OPEN OPEN

(1, (2, (3, (4, (5, (6, (7, (8,

FILE="U_RESULTS.PLT") FILE="U_RESIDUAL.PLT") FILE="V_RESULTS.PLT") FILE="V_RESIDULAS.PLT") FILE="P_RESULTS.PLT") FILE="Pc_RESIDUALS.PLT") FILE="U.PLT") FILE="V.PLT")

WRITE(1,*) 'VARIABLES = "U(Y) at X=0.5", "Y"' WRITE(1,*) 'ZONE J=', Ny, ' F=POINT' WRITE(3,*) 'VARIABLES = "X", "V(X) at Y=0.5"' WRITE(3,*) 'ZONE I=', Nx, ' F=POINT' WRITE(2,*) 'VARIABLES = "Iteration", "U Residual"' WRITE(4,*) 'VARIABLES = "Iteration", "V Residual"' WRITE(6,*) 'VARIABLES = "Iteration", "P Correction Residual"' WRITE(5,*) 'VARIABLES = "X", "Y", "PRESSURE"' WRITE(5,*) 'ZONE I=', Nx, ' J=',Ny, ' F=POINT' WRITE(7,*) 'VARIABLES = "X", "Y", "VELOCITY"' WRITE(7,*) 'ZONE I=', Nx, ' J=',Ny, ' F=POINT' WRITE(8,*) 'VARIABLES = "X", "Y", "VELOCITY"' WRITE(8,*) 'ZONE I=', Nx, ' J=',Ny, ' F=POINT' END SUBROUTINE FILES SUBROUTINE WRITES DO j=1, Ny WRITE(1,*) U(Nx/2,j), Y(Nx/2,j) END DO DO i=1, Nx WRITE(3,*) X(i,Ny/2), V(i,Ny/2) END DO DO j=1, Ny DO i=1, Nx WRITE(5,*) X(i,j), Y(i,j), P(i,j) WRITE(7,*) X(i,j), Y(i,j), U(i,j) WRITE(8,*) X(i,j), Y(i,j), V(i,j) END DO END DO pause END SUBROUTINE WRITES

2

SUBROUTINE MESH(Nx, Ny, Dx, Dy, Lx, Ly, X, Y) IMPLICIT NONE INTEGER :: i, j INTEGER, INTENT(IN) :: Nx, Ny REAL(8), INTENT(IN) :: Ly, Lx, Dx, Dy REAL(8), DIMENSION(Nx,Ny), INTENT(OUT) :: X, Y X(1,:) = 0 X(Nx,:) = Lx Y(:,1) = 0 Y(:,Ny) = Ly DO i=2, Nx-1 X(i,:) = Dx/2 + (i-2)*Dx END DO DO j=2, Ny-1 Y(:,j) = Dy/2 + (j-2)*Dy END DO END SUBROUTINE MESH SUBROUTINE FRAC(Dx, Dy, Nx, Ny, X, Y, FracE, FracW, FracN, FracS) IMPLICIT NONE INTEGER :: i, j, Nx, Ny REAL(8), INTENT(IN) :: Dx, Dy REAL(8), INTENT(IN), DIMENSION(Nx, Ny) :: X, Y REAL(8), INTENT(OUT), DIMENSION(Nx, Ny) :: FracE, FracW, FracN, FracS DO i=2, Nx-1 FracE(i,2:Ny-1) = Dx/(2*(X(i+1,2:Ny-1)-X(i,2:Ny-1))) FracW(i,2:Ny-1) = Dx/(2*(X(i,2:Ny-1)-X(i-1,2:Ny-1))) END DO DO j=2, Ny-1 FracN(2:Nx-1,j) = Dy/(2*(Y(2:Nx-1,j+1)-Y(2:Nx-1,j))) FracS(2:Nx-1,j) = Dy/(2*(Y(2:Nx-1,j)-Y(2:Nx-1,j-1))) END DO END SUBROUTINE FRAC SUBROUTINE BOUNDARIES(Ul, Nx, Ny, ApU, ApV, An, Aw, As, Ae, SuU, SuV, Sp, SuPc, ScU,& & ScV, AnPc, AwPc, AsPc, AePc) IMPLICIT NONE INTEGER :: i, j INTEGER, INTENT(IN) :: Nx, Ny REAL(8), INTENT(IN) :: Ul REAL(8), DIMENSION(Nx, Ny), INTENT(OUT) :: Apu, Apv, An, Aw, As, Ae, SuU, SuV, & & SuPc, Sp, ScU, ScV REAL(8), DIMENSION(Nx, Ny), INTENT(OUT) :: AnPc, AwPc, AsPc, AePc SuU = 0.0; SuV =0.0; Sp = 0.0; SuPc = 0.0 ScU = 0.0; ScV = 0.0 An = 0.0; Aw =0.0; As = 0.0; Ae = 0.0 AnPc = 0.0; AwPc =0.0; AsPc = 0.0; AePc = 0.0 !west boundary Ap(1,2:Ny-1) = 1.0 ApU(1,2:Ny-1) = 1.0 !U ApV(1,2:Ny-1) = 1.0 !V ApPc(1,2:Ny-1) = 1.0 !Pc !east boundary Ap(Nx,2:Ny-1) = 1.0 ApU(Nx,2:Ny-1) = 1.0 !U ApV(Nx,2:Ny-1) = 1.0 !V ApPc(Nx,2:Ny-1) = 1.0 !Pc !north boundary Ap(:,Ny) = 1.0 ApU(:,Ny) = 1.0 !U

3

SuU(:,Ny) = Ul ApV(:,Ny) = 1.0 !V ApPc(:,Ny) = 1.0 !Pc !south boundary Ap(:,1) = 1.0 ApU(:,1) = 1.0 !U ApV(:,1) = 1.0 !V ApPc(:,1) = 1.0 !Pc END SUBROUTINE BOUNDARIES SUBROUTINE COEFFICIENTS(X, Y, Nx, Ny, Dn, Dw, Ds, De, AlphaU, & & AlphaV, AlphaP, RE, Dx, Dy, & & ApU, ApV,Ap, An, Aw, As, Ae, SuU, SuV, Sp, & & ScU, ScV, Uold, Vold, Ape, Apn, & & Apw, Aps, ApPc, AnPc, AwPc, AsPc, AePc, SuPc, & & PfaceH, PfaceV, FluxU, FluxV, & & FracE, FracW, FracN, FracS) IMPLICIT NONE INTEGER :: i, j INTEGER, INTENT(IN) :: Nx, Ny REAL(8), DIMENSION(Nx,Ny) :: Ab REAL(8), INTENT(IN) :: AlphaU, AlphaV, AlphaP, RE, Dx, Dy REAL(8), DIMENSION(Nx,Ny), INTENT(IN) :: X, Y, PfaceH, PfaceV, Uold, Vold REAL(8), DIMENSION(Nx,Ny), INTENT(IN) :: Sp, ScU, ScV REAL(8), DIMENSION(Nx,Ny), INTENT(IN) :: FluxU, FluxV, FracE, FracW, FracN, FracS REAL(8), DIMENSION(Nx,Ny), INTENT(OUT) :: ApU, ApV, Ap, An, Aw, As, Ae, SuU, SuV REAL(8), DIMENSION(Nx,Ny), INTENT(OUT) :: Ape, Apn, Apw, Aps REAL(8), DIMENSION(Nx,Ny), INTENT(OUT) :: ApPc, AnPc, AwPc, AsPc, AePc, SuPc REAL(8), DIMENSION(Nx,Ny), INTENT(OUT) :: Dn, Dw, Ds, De Dn = 0; Dw = 0; Ds = 0; De = 0; Ab = 0 DO i=2, Nx-1 DO j=2, Ny-1 !General coefficients for V and U An(i,j) = (1./RE)*Dx/(Y(i,j+1)-Y(i,j)) Aw(i,j) = (1./RE)*Dy/(X(i,j)-X(i-1,j)) As(i,j) = (1./RE)*Dx/(Y(i,j)-Y(i,j-1)) Ae(i,j) = (1./RE)*Dy/(X(i+1,j)-X(i,j))

+ + + +

max(0.0,-FluxV(i,j)*Dx) max(FluxU(i-1,j)*Dy,0.0) max(FluxV(i,j-1)*Dx,0.0) max(0.0,-FluxU(i,j)*Dy)

Ab(i,j) = Dx*FluxV(i,j) - Dx*FluxV(i,j-1) + Dy*FluxU(i,j) - Dy*FluxU(i-1,j) Ap(i,j) = An(i,j) + Aw(i,j) + As(i,j) + Ae(i,j) + Ab(i,j) - & & Sp(i,j)*Dx*Dy !U coefficient and source term ApU(i,j) = Ap(i,j)/AlphaU SuU(i,j) = -Dy*(PfaceH(i,j)-PfaceH(i-1,j)) & & + ((1-AlphaU)/AlphaU)*Ap(i,j)*Uold(i,j) + ScU(i,j)*Dx*Dy !V coefficient and source term ApV(i,j) = Ap(i,j)/AlphaV SuV(i,j) = -Dx*(PfaceV(i,j)-PfaceV(i,j-1)) & & + ((1-AlphaV)/AlphaV)*Ap(i,j)*Vold(i,j) + ScV(i,j)*Dx*Dy END DO END DO DO i=2, Nx-1 DO j=2, Ny-1 !Cell faces coeffiicients Ape(i,j) = INTERPOLATION(FracE(i,j), An(i+1,j)+ & & Aw(i+1,j)+As(i+1,j)+Ae(i+1,j), & & An(i,j)+Aw(i,j)+As(i,j)+Ae(i,j)) & & - (X(i+1,j)-X(i,j))*Dy* & & INTERPOLATION(FracE(i,j), Sp(i+1,j), Sp(i,j)) Apn(i,j) = INTERPOLATION(FracN(i,j), An(i,j+1)+ & & Aw(i,j+1)+As(i,j+1)+Ae(i,j+1), & & An(i,j)+Aw(i,j)+As(i,j)+Ae(i,j)) & & - (Y(i,j+1)-Y(i,j))*Dx* & & INTERPOLATION(FracN(i,j), Sp(i,j+1), Sp(i,j))

4

Apw(i,j) = INTERPOLATION(FracW(i,j), An(i-1,j)+ & & Aw(i-1,j)+As(i-1,j)+Ae(i-1,j), & & An(i,j)+Aw(i,j)+As(i,j)+Ae(i,j)) & & - (X(i,j)-X(i-1,j))*Dy* & & INTERPOLATION(FracW(i,j), Sp(i-1,j), Sp(i,j)) Aps(i,j) = INTERPOLATION(FracS(i,j), An(i,j-1)+ & & Aw(i,j-1)+As(i,j-1)+Ae(i,j-1), & & An(i,j)+Aw(i,j)+As(i,j)+Ae(i,j)) & & - (Y(i,j)-Y(i,j-1))*Dx* & & INTERPOLATION(FracS(i,j), Sp(i,j-1), Sp(i,j)) END DO END DO Apw(2,2:Ny-1) = Ap(1,2:Ny-1) Aps(2:Nx-1,2) = Ap(2:Nx-1,1) Ape(Nx-1,2:Ny-1) = Ap(Nx,2:Ny-1) Apn(2:Nx-1,Ny-1) = Ap(2:Nx-1,Ny) DO i=2, Nx-1 DO j=2, Ny-1 !P correction coeffiicients Dn(i,j) = AlphaV*Dx/Apn(i,j) Dw(i,j) = AlphaU*Dy/Apw(i,j) Ds(i,j) = AlphaV*Dx/Aps(i,j) De(i,j) = AlphaU*Dy/Ape(i,j) AnPc(i,j) AwPc(i,j) AsPc(i,j) AePc(i,j)

= = = =

Dx*Dn(i,j) Dy*Dw(i,j) Dx*Ds(i,j) Dy*De(i,j)

AsPc(i,2) = 0 AePc(Nx-1,j) = 0 AnPc(i,Ny-1) = 0 AwPc(2,j) = 0 ApPc(i,j) = AnPc(i,j) + AwPc(i,j) + AsPc(i,j) + AePc(i,j) SuPc(i,j) = Dy*FluxU(i-1,j) - Dy*FluxU(i,j) + Dx*FluxV(i,j-1) - Dx*FluxV(i,j) END DO END DO END SUBROUTINE COEFFICIENTS SUBROUTINE FLUX(Nx, Ny, AlphaU, AlphaV, FracE, FracN, Dx, Dy, X, Y, & & Ape, Apn, FluxUold, FluxVold, ScU, ScV & & , PfaceH, PfaceV, FluxU, FluxV) IMPLICIT NONE INTEGER :: i, j INTEGER, INTENT(IN) :: Nx, Ny REAL(8), INTENT(IN) :: AlphaU, AlphaV, Dx, Dy REAL(8), DIMENSION(Nx,Ny), INTENT(IN) :: X, Y, Ape, Apn, FluxUold, FluxVold, ScU, ScV REAL(8), DIMENSION(Nx,Ny), INTENT(IN) :: PfaceH, PfaceV, FracE, FracN REAL(8), DIMENSION(Nx,Ny), INTENT(OUT) :: FluxU, FluxV REAL(8), DIMENSION(Nx,Ny) :: AtermU, StermU , PtermU, OtermU REAL(8), DIMENSION(Nx,Ny) :: AtermV, StermV , PtermV, OtermV AtermU=0; StermU=0; PtermU=0; OtermU=0 AtermV=0; StermV=0; PtermV=0; OtermV=0 DO i=2, Nx-2 DO j=2, Ny-1 AtermU (i,j) = FracE(i,j)*Ap(i+1,j)*U(i+1,j) + (1-FracE(i,j))*Ap(i,j)*U(i,j) StermU(i,j) = AlphaU*((FracE(i,j)*ScU(i+1,j) + & & (1-FracE(i,j))*ScU(i,j))*(X(i+1,j)-X(i,j))*Dy & & - FracE(i,j)*ScU(i+1,j)*Dx*Dy - & & (1-FracE(i,j))*ScU(i,j)*Dx*Dy) PtermU(i,j) = AlphaU*(-Dy*(P(i+1,j)-P(i,j)) + & & FracE(i,j)*Dy*(PfaceH(i+1,j)-PfaceH(i,j)) &

5

& + (1-FracE(i,j))*Dy*(PfaceH(i,j)-PfaceH(i-1,j))) OtermU(i,j) = (1-AlphaU)*(Ape(i,j)*FluxUold(i,j) - & & FracE(i,j)*Ap(i+1,j)*Uold(i+1,j) - & & (1-FracE(i,j))*Ap(i,j)*Uold(i,j)) FluxU(i,j) = (1/Ape(i,j)) * (AtermU(i,j) + StermU(i,j) + & & PtermU(i,j) + OtermU(i,j)) END DO END DO DO i=2, Nx-1 DO j=2, Ny-2 AtermV (i,j) = FracN(i,j)*Ap(i,j+1)*V(i,j+1) + (1-FracN(i,j))*Ap(i,j)*V(i,j) StermV(i,j) = AlphaV*((FracN(i,j)*ScV(i,j+1) + & & (1-FracN(i,j))*ScV(i,j))*(Y(i,j+1)-Y(i,j))*Dx & & - FracN(i,j)*ScV(i,j+1)*Dx*Dy - & & (1-FracN(i,j))*ScV(i,j)*Dx*Dy) PtermV(i,j) = AlphaV*(-Dx*(P(i,j+1)-P(i,j)) + & & FracN(i,j)*Dx*(PfaceV(i,j+1)-PfaceV(i,j)) & & + (1-FracN(i,j))*Dx*(PfaceV(i,j)-PfaceV(i,j-1))) OtermV(i,j) = (1-AlphaV)*(Apn(i,j)*FluxVold(i,j) - & & FracN(i,j)*Ap(i,j+1)*Vold(i,j+1) - & & (1-FracN(i,j))*Ap(i,j)*Vold(i,j)) FluxV(i,j) = (1/Apn(i,j)) * (AtermV(i,j) + StermV(i,j) + & & PtermV(i,j) + OtermV(i,j)) END DO END DO END SUBROUTINE FLUX SUBROUTINE PFACE(Nx, Ny, P, FracE, FracN, PfaceH, PfaceV) IMPLICIT NONE INTEGER :: i,j INTEGER, INTENT(IN) :: Nx, Ny REAL(8), DIMENSION(Nx,Ny), INTENT(IN) :: P, FracE, FracN REAL(8), DIMENSION(Nx,Ny), INTENT(OUT) :: PfaceH, PfaceV DO i=2, Nx-2 DO j=2, Ny-1 PfaceH(i,j) = INTERPOLATION(FracE(i,j), P(i+1,j), P(i,j)) END DO END DO DO j=2, Ny-2 DO i=2, Nx-1 PfaceV(i,j) = INTERPOLATION(FracN(i,j), P(i,j+1), P(i,j)) END DO END DO PfaceH(1,2:Ny-1) = P(1,2:Ny-1) PfaceH(Nx-1,2:Ny-1) = P(Nx,2:Ny-1) PfaceV(2:Nx-1,1) = P(2:Nx-1,1) PfaceV(2:Nx-1,Ny-1) = P(2:Nx-1,Ny) END SUBROUTINE PFACE SUBROUTINE CORRECTION(Nx, Ny, PcfaceH, PcfaceV, FracE, FracN, & & Pc, Dn, Dw, Ds, De, FluxUold, FluxVold, Uold, Vold, & & FluxU, FluxV, U, V, p) IMPLICIT NONE

6

INTEGER :: i,j INTEGER, INTENT(IN) :: Nx, Ny REAL(8), DIMENSION(Nx,Ny) :: PcfaceH, PcfaceV, DpU, DpV REAL(8), DIMENSION(Nx,Ny), INTENT(IN) :: FracE, FracN REAL(8), DIMENSION(Nx,Ny), INTENT(INOUT) :: Pc REAL(8), DIMENSION(Nx,Ny), INTENT(IN) :: Dn, Dw, Ds, De REAL(8), DIMENSION(Nx,Ny), INTENT(OUT) :: FluxUold, FluxVold, Uold, Vold REAL(8), DIMENSION(Nx,Ny), INTENT(INOUT) :: FluxU, FluxV REAL(8), DIMENSION(Nx,Ny), INTENT(INOUT) :: P REAL(8), DIMENSION(Nx,Ny), INTENT(INOUT) :: U, V !Border pressure-correction extrapolation Pc(1,2:Ny-1) = Pc(2,2:Ny-1) - (Pc(3,2:Ny-1)-Pc(2,2:Ny-1))* & & (X(2,2:Ny-1)-X(1,2:Ny-1))/(X(3,2:Ny-1)-X(2,2:Ny-1)) Pc(Nx,2:Ny-1) = Pc(Nx-1,2:Ny-1) + (Pc(Nx-1,2:Ny-1)-Pc(Nx-2,2:Ny-1))* & & (X(Nx,2:Ny-1)-X(Nx-1,2:Ny-1))/(X(Nx-1,2:Ny-1)-X(Nx-2,2:Ny-1)) Pc(1:Nx,1) = Pc(1:Nx,2) - (Pc(1:Nx,3)-Pc(1:Nx,2))* & & (Y(1:Nx,2)-Y(1:Nx,1))/(Y(1:Nx,3)-Y(1:Nx,2)) Pc(1:Nx,Ny) = Pc(1:Nx,Ny-1) + (Pc(1:Nx,Ny-1)-Pc(1:Nx,Ny-2))* & & (Y(1:Nx,Ny)-Y(1:Nx,Ny-1))/(Y(1:Nx,Ny-1)-Y(1:Nx,Ny-2)) DO i=2, Nx-2 DO j=2, Ny-1 PcfaceH(i,j) = INTERPOLATION(FracE(i,j), Pc(i+1,j), Pc(i,j)) END DO END DO DO j=2, Ny-2 DO i=2, Nx-1 PcfaceV(i,j) = INTERPOLATION(FracN(i,j), Pc(i,j+1), Pc(i,j)) END DO END DO PcfaceH(1,2:Ny-1) = Pc(1,2:Ny-1) PcfaceH(Nx-1,2:Ny-1) = Pc(Nx,2:Ny-1) PcfaceV(2:Nx-1,1) = Pc(2:Nx-1,1) PcfaceV(2:Nx-1,Ny-1) = Pc(2:Nx-1,Ny) DO i=2, Nx-1 DO j=2, Ny-1 DpU(i,j) = AlphaU*Dy/Ap(i,j) DpV(i,j) = AlphaV*Dx/Ap(i,j) U(i,j) = U(i,j) + DpU(i,j)*(PcfaceH(i-1,j)-PcfaceH(i,j)) V(i,j) = V(i,j) + Dpv(i,j)*(PcfaceV(i,j-1)-PcfaceV(i,j)) END DO END DO DO i=2, Nx-1 DO j=2, Ny-2 FluxV(i,j) = FluxV(i,j) + Dn(i,j)*(Pc(i,j) - Pc(i,j+1)) END DO END DO DO i=2, Nx-2 DO j=2, Ny-1 FluxU(i,j) = FluxU(i,j) + De(i,j)*(Pc(i,j) - Pc(i+1,j)) END DO END DO

7

DO i=1, Nx DO j=1, Ny P(i,j) = P(i,j) + AlphaP*Pc(i,j) Uold(i,j) = U(i,j) Vold(i,j) = V(i,j) FluxVold(i,j) = FluxV(i,j) FluxUold(i,j) = FluxU(i,j) END DO END DO !Border pressure extrapolation P(1,2:Ny-1) = P(2,2:Ny-1) - (P(3,2:Ny-1)-P(2,2:Ny-1))* & & (X(2,2:Ny-1)-X(1,2:Ny-1))/(X(3,2:Ny-1)-X(2,2:Ny-1)) P(Nx,2:Ny-1) = P(Nx-1,2:Ny-1) + (P(Nx-1,2:Ny-1)-P(Nx-2,2:Ny-1))* & & (X(Nx,2:Ny-1)-X(Nx-1,2:Ny-1))/(X(Nx-1,2:Ny-1)-X(Nx-2,2:Ny-1)) P(1:Nx,1) = P(1:Nx,2) - (P(1:Nx,3)-P(1:Nx,2))* & & (Y(1:Nx,2)-Y(1:Nx,1))/(Y(1:Nx,3)-Y(1:Nx,2)) P(1:Nx,Ny) = P(1:Nx,Ny-1) + (P(1:Nx,Ny-1)-P(1:Nx,Ny-2))* & & (Y(1:Nx,Ny)-Y(1:Nx,Ny-1))/(Y(1:Nx,Ny-1)-Y(1:Nx,Ny-2)) END SUBROUTINE CORRECTION SUBROUTINE TDMA (a, b, c, d, N, ans) IMPLICIT NONE INTEGER :: k INTEGER, INTENT(IN) :: N REAL(8), DIMENSION(N) :: dp, cp REAL(8), DIMENSION(N), INTENT(IN) :: a, b, c, d REAL(8), DIMENSION(N), INTENT(OUT) :: ans dp(1)=-a(1)/d(1) cp(1)=c(1)/d(1) DO k=2, N dp(k) = -a(k)/(d(k)+b(k)*dp(k-1)) cp(k) = (-b(k)*cp(k-1)+c(k))/(d(k)+b(k)*dp(k-1)) END DO ans(N) = cp(N) DO k=N-1,1,-1 ans(k) = dp(k)*ans(k+1) + cp(k) END DO END SUBROUTINE TDMA SUBROUTINE H_SWEEPER (Nx, Ny, Ap, An, Aw, As, Ae, Su, T) IMPLICIT NONE INTEGER :: i, j INTEGER, INTENT(IN) :: Nx, Ny REAL(8), DIMENSION(Ny) :: a, b, c, d, ans REAL(8), DIMENSION(Nx,Ny), INTENT(IN) :: Ap, An, Aw, As, Ae, Su REAL(8), DIMENSION(Nx,Ny), INTENT(INOUT) :: T DO i=1, Nx DO j=1, Ny IF (i==1) THEN c(j) = Ae(i,j)*T(i+1,j) + Su(i,j) ELSE IF (i==Nx) THEN c(j) = Aw(i,j)*T(i-1,j) + Su(i,j) ELSE c(j) = Aw(i,j)*T(i-1,j) + Ae(i,j)*T(i+1,j) + Su(i,j) END IF a(j) = -An(i,j)

8

b(j) = -As(i,j) d(j) = Ap(i,j) END DO CALL TDMA (a, b, c, d, Ny, ans) DO j=1,Ny T(i,j) = ans(j) END DO T(i,1:Ny) = ans(1:Ny) END DO END SUBROUTINE H_SWEEPER SUBROUTINE V_SWEEPER (Nx, Ny, Ap, An, Aw, As, Ae, Su, T) IMPLICIT NONE INTEGER :: i, j INTEGER, INTENT(IN) :: Nx, Ny REAL(8), DIMENSION(Nx) :: a, b, c, d, ans REAL(8), DIMENSION(Nx,Ny), INTENT(IN) :: Ap, An, Aw, As, Ae, Su REAL(8), DIMENSION(Nx,Ny), INTENT(INOUT) :: T DO j=1, Ny DO i=1, Nx IF (j==1) THEN c(i) = An(i,j)*T(i,j+1) + Su(i,j) ELSE IF (j==Ny) THEN c(i) = As(i,j)*T(i,j-1) + + Su(i,j) ELSE c(i) = As(i,j)*T(i,j-1) + An(i,j)*T(i,j+1) + Su(i,j) END IF a(i) = -Ae(i,j) b(i) = -Aw(i,j) d(i) = Ap(i,j) END DO CALL TDMA (a, b, c, d, Nx, ans) DO I=1,Nx T(i,j) = ans(i) END DO T(1:Nx,j) = ans(1:Nx) END DO END SUBROUTINE V_SWEEPER SUBROUTINE SOLVER (Nx, Ny, Ap, An, Aw, As, Ae, Su, T) IMPLICIT NONE INTEGER :: i, j INTEGER, INTENT(IN) :: Nx, Ny REAL(8), DIMENSION(Nx,Ny), INTENT(IN) :: Ap, An, Aw, As, Ae, Su REAL(8), DIMENSION(Nx,Ny), INTENT(INOUT) :: T CALL V_SWEEPER (Nx, Ny, Ap, An, Aw, As, Ae, Su, T) CALL H_SWEEPER (Nx, Ny, Ap, An, Aw, As, Ae, Su, T) END SUBROUTINE SOLVER SUBROUTINE RESIDUALS(T, Aw, As, Ae, An, Residual, Ap, Nx, Ny, Su) IMPLICIT NONE INTEGER :: i, j INTEGER, INTENT(IN) :: Nx, Ny REAL(8), DIMENSION(Nx,Ny), INTENT(IN) :: Ap, An, Aw, As, Ae, Su REAL(8), DIMENSION(Nx,Ny), INTENT(IN) :: T REAL(8), INTENT(OUT) :: RESIDUAL Residual = 0.0 DO j=2, Ny-1 DO i=2, Nx-1 Residual = abs(Aw(i,j)*T(i-1,j) + Ae(i,j)*T(i+1,j) + & & An(i,j)*T(i,j+1) + As(i,j)*T(i,j-1) + Su(i,j) - & & Ap(i,j)*T(i,j)) + Residual END DO END DO END SUBROUTINE RESIDUALS

9

!***FUNCTIONS***! REAL(8) FUNCTION INTERPOLATION(F, M, P) IMPLICIT NONE REAL(8), INTENT(IN) :: F, M, P INTERPOLATION = F*M + (1-F)*P END FUNCTION INTERPOLATION END PROGRAM CAVITY NON-STAGERRED

10