# Solve Pressure Correction.f90 - Solution of pressure-correction equation and correction of U,V and P

(Difference between revisions)
 Revision as of 12:59, 3 May 2010 (view source)Michail (Talk | contribs)← Older edit Revision as of 19:30, 15 August 2012 (view source)Michail (Talk | contribs) Newer edit → Line 26: Line 26: DpV(2:NXmax,2:NYmax) = 1./Ap(2:NXmax,2:NYmax,2) DpV(2:NXmax,2:NYmax) = 1./Ap(2:NXmax,2:NYmax,2) + !This is 1./AP coefficient involved in the derivation of the pressure equation + !and correction velocities. Check for simple algorithm in dedicated books. + !AP is the central coefficent there is one for U momentum equation --> APU and one for V-momentum APV. + ! in fact for this code APU(i,j) = AP(i,j,1) + ! APV(i,j) = AP(i,j,2) + ! because he works with vector componnents style !**************************************************************************************** !**************************************************************************************** Line 35: Line 41: - DXc_e = Xc(i+1,j+1) - Xc(i,j+1) + DXc_e = Xc(i+1,j+1) - Xc(i,j+1) !this is the distance between node P and node E - S_e  =  Y(i  ,j+1) -  Y(i,j  ) + + S_e  =  Y(i  ,j+1) -  Y(i,j  ) !This is the east surface S_e VOL_e = DXc_e * S_e VOL_e = DXc_e * S_e !------------------------------------------------------------- !------------------------------------------------------------- - APU_e =      0.5 * (  DpU(i,j+1)  +  DpU(i+1,j+1)  ) + APU_e =      0.5 * (  DpU(i,j+1)  +  DpU(i+1,j+1)  ) !this is the iterpolation of 1/APU at east location - Ul_e  =      0.5 * (    F(i,j+1,1) +    F(i+1,j+1,1) ) + Ul_e  =      0.5 * (    F(i,j+1,1) +    F(i+1,j+1,1) ) !this is the interpolation of velocity U on east location - DPl_e =      0.5 * ( DPx_c(i,j+1)  + DPx_c(i+1,j+1)  ) + DPl_e =      0.5 * ( DPx_c(i,j+1)  + DPx_c(i+1,j+1)  ) ! the same for interpolation of DPx + + DPx_e =            ( F(i+1,j+1,4) - F(i,j+1,4) ) / DXc_e ! this is the gradient of pressure at east location usint central difference - DPx_e =            ( F(i+1,j+1,4) - F(i,j+1,4) ) / DXc_e - U_e  = Ul_e  + APU_e * VOL_e * ( DPl_e - DPx_e) + U_e  = Ul_e  + APU_e * VOL_e * ( DPl_e - DPx_e) ! This is the interpolated velocity using at east location using the interpolation of Rhie and Chow !-------------------------------------------------------------- !-------------------------------------------------------------- - Con_e(i,j) =  U_e  * S_e + Con_e(i,j) =  U_e  * S_e ! this is the convective flux.

## Revision as of 19:30, 15 August 2012

```
!Sample program for solving Lid-driven cavity flow test using SIMPLE-algorithm
! solution of pressure-correction quation and correction U,V, and P

!This program is free software; you can redistribute it and/or
!modify it under the terms of the GNU General Public License

!This program is distributed in the hope that it will be useful,
!but WITHOUT ANY WARRANTY; without even the implied warranty of
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!GNU General Public License for more details.

!You should have received a copy of the GNU General Public License
!along with this program; if not, write to the Free Software
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.

!**********************************************************************
Subroutine Solve_Pressure_Correction
include 'icomm_1.f90'

DpU(2:NXmax,2:NYmax) = 1./Ap(2:NXmax,2:NYmax,1)
DpV(2:NXmax,2:NYmax) = 1./Ap(2:NXmax,2:NYmax,2)

!This is 1./AP coefficient involved in the derivation of the pressure equation
!and correction velocities. Check for simple algorithm in dedicated books.
!AP is the central coefficent there is one for U momentum equation --> APU and one for V-momentum APV.
! in fact for this code APU(i,j) = AP(i,j,1)
! APV(i,j) = AP(i,j,2)
! because he works with vector componnents style

!****************************************************************************************
! vertical  faces East - e
Do 101 I=1,NXmax
Do 101 J=1,NYmax-1

If((i.ne.1).and.(i.ne.NXmax))then

DXc_e = Xc(i+1,j+1) - Xc(i,j+1) !this is the distance between node P and node E

S_e   =  Y(i  ,j+1) -  Y(i,j  ) !This is the east surface S_e

VOL_e = DXc_e * S_e

!-------------------------------------------------------------
APU_e =       0.5 * (   DpU(i,j+1)   +   DpU(i+1,j+1)   ) !this is the iterpolation of 1/APU at east location
Ul_e  =       0.5 * (     F(i,j+1,1) +     F(i+1,j+1,1) ) !this is the interpolation of velocity U on east location
DPl_e =       0.5 * ( DPx_c(i,j+1)   + DPx_c(i+1,j+1)   ) ! the same for interpolation of DPx

DPx_e =             ( F(i+1,j+1,4) - F(i,j+1,4) ) / DXc_e ! this is the gradient of pressure at east location usint central difference

U_e   = Ul_e  + APU_e * VOL_e * ( DPl_e - DPx_e) ! This is the interpolated velocity using at east location using the interpolation of Rhie and Chow

!--------------------------------------------------------------
Con_e(i,j) =   U_e  * S_e ! this is the convective flux.

End If

Con_e(1,j) = 0.
Con_e(NXmax,j) = 0.

101 continue

!****************************************************************************************
! horisontal  faces
Do 102 I=1,NXmax-1
Do 102 J=1,NYmax

If((j.ne.1).and.(j.ne.NYmax))then

DYc_n = Yc(i+1,j+1) - Yc(i+1,j)
S_n   =  X(i+1,j  ) -  X(i  ,j)

VOL_n = DYc_n * S_n

!-----------------------------------------------------------
APV_n =       0.5 * (   DpV(i+1,j)   +   DpV(i+1,j+1)   )
Vl_n  =       0.5 * (     F(i+1,j,2) +     F(i+1,j+1,2) )
DPl_n =       0.5 * ( DPy_c(i+1,j)   + DPy_c(i+1,j+1)   )

DPy_n = ( F(i+1,j+1,4) - F(i+1,j,4) ) / DYc_n

V_n   = Vl_n  + APV_n * VOL_n * ( DPl_n - DPy_n )
!-----------------------------------------------------------

Con_n(i,j) =   V_n  * S_n

End If

Con_n(i,1) = 0.
Con_n(i,NYmax) = 0.

102 continue

!****************************************************************************************
Summ = 0.

Do 200 I=2,NXmax
Do 200 J=2,NYmax

S_e   =  Y(i  ,j) -  Y(i,j-1)
S_n   =  X(i  ,j) -    X(i-1,j)

An(i,j) = 0.5 * (   DpV(i,j)   +   DpV(i,j+1)   )  * S_n * S_n
As(i,j) = 0.5 * (   DpV(i,j)   +   DpV(i,j-1)   )  * S_n * S_n
Ae(i,j) = 0.5 * (   DpU(i,j)   +   DpU(i+1,j)   )  * S_e * S_e
Aw(i,j) = 0.5 * (   DpU(i,j)   +   DpU(i-1,j)   )  * S_e * S_e

Ap(i,j,3) = (An(i,j) + As(i,j) + Ae(i,j) +	Aw(i,j))

Sp(i,j,3) = -1. * ( Con_n(i-1,j) - Con_n(i-1,j-1) + Con_e(i,j-1) - Con_e(i-1,j-1) )

Summ = Summ + Sp(i,j,3) !
200 continue

write(*,*)'Sum',summ
!***********************************************************************************
!***********************************************************************************
write(*,*) Summ

write(*,*)'solve PP'
niter = 0

20	call Convergence_Criteria(3,Res_sum_before)
niter= niter + 1

Call TDMA_1(3)
call Convergence_Criteria(3,Res_sum_After)
If((abs(Res_sum_before-Res_sum_After).Ge.0.0000001).and.(niter.le.500).and.(abs(Res_sum_After).ge.0.0001))go to 20
write(*,*)'Res_sum_before-Res_sum_After',Res_sum_before-Res_sum_After,niter

!***********************************************************************************
!***********************************************************************************
!  velocities and pressure correction

Do 300 I=2,NXmax
Do 300 J=2,NYmax

DY = Y(i,j)-Y(i,j-1)
DX = X(i,j)-X(i-1,j)

PPe = 0.5 * ( F(i,j,3) + F(i+1,j,3) )
PPw = 0.5 * ( F(i,j,3) + F(i-1,j,3) )
PPn = 0.5 * ( F(i,j,3) + F(i,j+1,3) )
PPs = 0.5 * ( F(i,j,3) + F(i,j-1,3) )

F(i,j,1) = F(i,j,1) + (PPw - PPe)  * DpU(i,j)  * Dy
F(i,j,2) = F(i,j,2) + (PPs - PPn)  * DpV(i,j)  * Dx

F(i,j,4) = F(i,j,4) + 0.1 *   F(i,j,3)

300 continue

!***********************************************************************************
! correct convective fluxes  through vertical East faces
Do 701 I=1,NXmax
Do 701 J=1,NYmax-1

If((i.ne.1).and.(i.ne.NXmax))then

Con_e(i,j) =  Con_e(i,j) + Ae(i,j+1) * (F(i+1,j+1,3)-F(i,j+1,3) )

end if

if((i.eq.1).or.(i.eq.NXmax))Con_e(i,j)=0.

701 continue

! correct  convective fluxes through horisontal  North faces
Do 702 I=1,NXmax-1
Do 702 J=1,NYmax

If((j.ne.1).and.(j.ne.NYmax))then

Con_n(i,j) = Con_n(i,j) + An(i+1,j) * (F(i+1,j+1,3)-F(i+1,j,3) )

end if

if((j.eq.1).or.(j.eq.NYmax))Con_n(i,j)=0.

702 continue

!***********************************************************************************
Do 501 I=1,NXmaxC

F(i,1,4) = F(i,2,4)            !
F(i,NYmaxC,4) = F(i,NYmaxC-1,4)!

501 continue

Do 502 J=1,NYmaxC

F(1,j,4) = F(2,j,4)             !
F(NXmaxC,j,4) = F(NXmaxC-1,j,4) !

502 continue

F(:,:,4) = F(:,:,4) - F(3,4,4)
!***********************************************************************************

Return
End

```