CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Wiki > Solve Pressure Correction.f90 - Solution of pressure-correctio...

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

From CFD-Wiki

(Difference between revisions)
Jump to: navigation, search
 
Line 4: Line 4:
! solution of pressure-correction quation and correction U,V, and P
! solution of pressure-correction quation and correction U,V, and P
!Copyright (C) 2010  Michail Kiričkov
!Copyright (C) 2010  Michail Kiričkov
 +
!Copyright (C) 2016  Michail Kiričkov, Kaunas University for Technology
!This program is free software; you can redistribute it and/or
!This program is free software; you can redistribute it and/or
Line 22: Line 23:
Subroutine Solve_Pressure_Correction
Subroutine Solve_Pressure_Correction
include 'icomm_1.f90'
include 'icomm_1.f90'
-
 
+
integer test
-
      DpU(2:NXmax,2:NYmax) = 1./Ap(2:NXmax,2:NYmax,1)
+
    DpU(2:NXmax,2:NYmax) = 1./Ap(2:NXmax,2:NYmax,1)
-
      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  
!This is 1./AP coefficient involved in the derivation of the pressure equation  
!and correction velocities. Check for simple algorithm in dedicated books.
!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.
+
!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)
! in fact for this code APU(i,j) = AP(i,j,1)
! APV(i,j) = AP(i,j,2)
! APV(i,j) = AP(i,j,2)
! because he works with vector componnents style
! because he works with vector componnents style
-
 
+
!*****************************************************************************   
-
!****************************************************************************************   
+
! vertical  faces East - e
! vertical  faces East - e
   Do 101 I=1,NXmax
   Do 101 I=1,NXmax
   Do 101 J=1,NYmax-1
   Do 101 J=1,NYmax-1
-
 
+
    DXc_e = 0.
-
    If((i.ne.1).and.(i.ne.NXmax))then
+
S_e  = 0.
-
 
+
VOL_e = 0.
-
+
APU_e = 0.
-
      DXc_e = Xc(i+1,j+1) - Xc(i,j+1) !this is the distance between node P and node E
+
Ul_e  = 0.
-
 
+
DPl_e = 0. 
-
    S_e  =  Y(i  ,j+1) -  Y(i,j  ) !This is the east surface S_e
+
DPx_e = 0.
-
 
+
U_e  = 0.
 +
If(i.ne.NXmax)then
 +
      !this is the distance between node P and node E
 +
        DXc_e = Xc(i+1,j+1) - Xc(i,j+1)  
 +
        !This is the east surface S_e
 +
        S_e  =  Y(i  ,j+1) -  Y(i,j  )  
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)  ) !this is the iterpolation of 1/APU at east location
+
        !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
+
         APU_e =      0.5 * (  DpU(i,j+1)  +  DpU(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
+
        !this is the interpolation of velocity U on east location
-
 
+
         Ul_e  =      0.5 * (    F(i,j+1,1) +    F(i+1,j+1,1) )  
-
         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
+
        ! the same for interpolation of DPx
-
 
+
         DPl_e =      0.5 * ( DPx_c(i,j+1)  + DPx_c(i+1,j+1)  )  
-
     
+
        ! 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
+
         DPx_e =            ( F(i+1,j+1,4) - F(i,j+1,4) ) / DXc_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.
+
         U_e  = Ul_e  + APU_e * VOL_e * ( DPl_e - DPx_e)  
-
 
+
    !--------------------------------------------------------------     
-
+
Con_we(i,j) =  U_e  * S_e ! this is the convective flux.
End If
End If
-
+
      Con_we(1,j) = 0.
-
+
       Con_we(NXmax,j) = 0.
-
      Con_e(1,j) = 0.
+
101 continue
-
       Con_e(NXmax,j) = 0.
+
!********************************************************************
-
+
-
101 continue
+
-
 
+
-
!****************************************************************************************
+
! horisontal  faces
! horisontal  faces
   Do 102 I=1,NXmax-1
   Do 102 I=1,NXmax-1
   Do 102 J=1,NYmax
   Do 102 J=1,NYmax
-
 
+
    DYc_n = 0.
-
    If((j.ne.1).and.(j.ne.NYmax))then
+
    S_n  = 0.
-
 
+
    VOL_n = 0.
 +
        APV_n = 0.
 +
        Vl_n  = 0.
 +
        DPl_n = 0.
 +
        DPy_n = 0.
 +
  V_n  = 0.
 +
If((j.ne.NYmax))then
      DYc_n = Yc(i+1,j+1) - Yc(i+1,j)
      DYc_n = Yc(i+1,j+1) - Yc(i+1,j)
    S_n  =  X(i+1,j  ) -  X(i  ,j)
    S_n  =  X(i+1,j  ) -  X(i  ,j)
-
 
    VOL_n = DYc_n * S_n  
    VOL_n = DYc_n * S_n  
-
 
!-----------------------------------------------------------
!-----------------------------------------------------------
         APV_n =      0.5 * (  DpV(i+1,j)  +  DpV(i+1,j+1)  )
         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) )
         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)  )
         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
         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 )
   V_n  = Vl_n  + APV_n * VOL_n * ( DPl_n - DPy_n )
-
        !-----------------------------------------------------------
+
    !-----------------------------------------------------------
-
       
+
Con_ns(i,j) =  V_n  * S_n  
-
Con_n(i,j) =  V_n  * S_n  
+
    End If
-
 
+
        Con_ns(i,1) = 0.
-
 
+
      Con_ns(i,NYmax) = 0.
-
End If
+
-
 
+
-
   
+
-
    Con_n(i,1) = 0.
+
-
      Con_n(i,NYmax) = 0.
+
-
 
+
  102 continue
  102 continue
-
 
-
!****************************************************************************************
 
   Summ = 0.
   Summ = 0.
-
 
Do 200 I=2,NXmax
Do 200 I=2,NXmax
Do 200 J=2,NYmax
Do 200 J=2,NYmax
-
 
S_e  =  Y(i  ,j) -  Y(i,j-1)
S_e  =  Y(i  ,j) -  Y(i,j-1)
S_n  =  X(i  ,j) -    X(i-1,j)
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
-
An(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
-
As(i,j) = 0.5 * (  DpV(i,j)  +  DpV(i,j-1)  )  * S_n * S_n
+
    Aw(i,j) = 0.5 * (  DpU(i,j)  +  DpU(i-1,j)  )  * S_e * S_e
-
Ae(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))
-
Aw(i,j) = 0.5 * (  DpU(i,j)  +  DpU(i-1,j)  )  * S_e * S_e
+
     Sp(i,j,3) = -1. * ( Con_ns(i-1,j) - Con_ns(i-1,j-1) &
-
+
                      + Con_we(i,j-1) - Con_we(i-1,j-1) )
-
Ap(i,j,3) = (An(i,j) + As(i,j) + Ae(i,j) + Aw(i,j))
+
    Summ = Summ + Sp(i,j,3) !  
-
   
+
-
     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
  200 continue
-
 
+
    write(*,*)'Sum',summ
-
write(*,*)'Sum',summ
+
    write(*,*)'solve PP'
-
!***********************************************************************************
+
-
!***********************************************************************************
+
-
write(*,*) Summ
+
-
 
+
-
write(*,*)'solve PP'
+
  niter = 0
  niter = 0
 +
20 continue
 +
        call Convergence_Criteria(3,Res_sum_before,niter)
 +
        niter= niter + 1
 +
        Call TDMA_1(3)
 +
call Convergence_Criteria(3,Res_sum_After,niter)
-
+
    if((abs(Res_sum_before-Res_sum_After).Ge.1.0E-10) &
-
20 call Convergence_Criteria(3,Res_sum_before)
+
    .and.(niter.le.1500).and.(abs(Res_sum_After).ge.1.0E-06))go to 20  
-
        niter= niter + 1
+
write(*,*)'Pressure correction Res_sum_before-Res_sum_After' &
-
+
    ,Res_sum_before-Res_sum_After,niter
-
Call TDMA_1(3)
+
! pressure correction
-
call Convergence_Criteria(3,Res_sum_After)
+
Do 301 I=2,NXmax
-
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  
+
Do 301 J=2,NYmax
-
write(*,*)'Res_sum_before-Res_sum_After',Res_sum_before-Res_sum_After,niter
+
        F(i,j,4) = F(i,j,4) + 0.2 *   F(i,j,3)
-
 
+
301 continue
-
   
+
-
 
+
-
+
-
!***********************************************************************************
+
-
!***********************************************************************************
+
!  velocities and pressure correction
!  velocities and pressure correction
-
 
+
Do 302 I=2,NXmax
-
Do 300 I=2,NXmax
+
Do 302 J=2,NYmax
-
Do 300 J=2,NYmax
+
-
+
DY = Y(i,j)-Y(i,j-1)
DY = Y(i,j)-Y(i,j-1)
DX = X(i,j)-X(i-1,j)
DX = X(i,j)-X(i-1,j)
-
 
PPe = 0.5 * ( F(i,j,3) + F(i+1,j,3) )
PPe = 0.5 * ( F(i,j,3) + F(i+1,j,3) )
PPw = 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) )
PPn = 0.5 * ( F(i,j,3) + F(i,j+1,3) )
PPs = 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,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,2) = F(i,j,2) + (PPs - PPn)  * DpV(i,j)  * Dx  
-
 
+
  302 continue
-
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
Do 501 I=1,NXmaxC
-
 
       F(i,1,4) = F(i,2,4)            !   
       F(i,1,4) = F(i,2,4)            !   
       F(i,NYmaxC,4) = F(i,NYmaxC-1,4)!   
       F(i,NYmaxC,4) = F(i,NYmaxC-1,4)!   
-
 
-
 
  501 continue
  501 continue
  Do 502 J=1,NYmaxC
  Do 502 J=1,NYmaxC
-
 
       F(1,j,4) = F(2,j,4)            !   
       F(1,j,4) = F(2,j,4)            !   
       F(NXmaxC,j,4) = F(NXmaxC-1,j,4) !   
       F(NXmaxC,j,4) = F(NXmaxC-1,j,4) !   
-
 
-
 
  502 continue
  502 continue
-
 
  F(:,:,4) = F(:,:,4) - F(3,4,4)  
  F(:,:,4) = F(:,:,4) - F(3,4,4)  
-
!***********************************************************************************
+
999 continue
-
 
+
-
 
+
Return
Return

Latest revision as of 15:03, 19 May 2016


!Sample program for solving Lid-driven cavity flow test using SIMPLE-algorithm
! solution of pressure-correction quation and correction U,V, and P
!Copyright (C) 2010  Michail Kiričkov
!Copyright (C) 2016  Michail Kiričkov, Kaunas University for Technology

!This program is free software; you can redistribute it and/or
!modify it under the terms of the GNU General Public License
!as published by the Free Software Foundation; either version 2
!of the License, or (at your option) any later version.

!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'
integer test
    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
 	     DXc_e = 0.
		 S_e   = 0.
		 VOL_e = 0.
		 APU_e = 0.
		 Ul_e  = 0.
		 DPl_e = 0.   
		 DPx_e = 0.
		 U_e   = 0.
	If(i.ne.NXmax)then
 	     !this is the distance between node P and node E
         DXc_e = Xc(i+1,j+1) - Xc(i,j+1) 
         !This is the east surface S_e
         S_e   =  Y(i  ,j+1) -  Y(i,j  ) 
		 VOL_e = DXc_e * S_e 
	!-------------------------------------------------------------
        !this is the iterpolation of 1/APU at east location         
         APU_e =       0.5 * (   DpU(i,j+1)   +   DpU(i+1,j+1)   ) 
        !this is the interpolation of velocity U on east location
         Ul_e  =       0.5 * (     F(i,j+1,1) +     F(i+1,j+1,1) ) 
        ! the same for interpolation of DPx
         DPl_e =       0.5 * ( DPx_c(i,j+1)   + DPx_c(i+1,j+1)   ) 
        ! 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 
        ! This is the interpolated velocity using at east location 
        !         using the interpolation of Rhie and Chow 	
         U_e   = Ul_e  + APU_e * VOL_e * ( DPl_e - DPx_e) 
    !--------------------------------------------------------------    
		 Con_we(i,j) =   U_e  * S_e ! this is the convective flux.
	End If
 	     Con_we(1,j) = 0.
  	     Con_we(NXmax,j) = 0.
101 continue
!********************************************************************
! horisontal  faces
   Do 102 I=1,NXmax-1
   Do 102 J=1,NYmax
 	     DYc_n = 0.
	     S_n   = 0.
	     VOL_n = 0.
         APV_n = 0.
         Vl_n  = 0.
         DPl_n = 0.
         DPy_n = 0.
   		 V_n   = 0.
	If((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_ns(i,j) =   V_n  * S_n 
    End If
         Con_ns(i,1) = 0.
 	     Con_ns(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_ns(i-1,j) - Con_ns(i-1,j-1) &
                      + Con_we(i,j-1) - Con_we(i-1,j-1) )
    Summ = Summ + Sp(i,j,3) ! 
 200 continue
    write(*,*)'Sum',summ
    write(*,*)'solve PP'
 niter = 0
20 continue
        call Convergence_Criteria(3,Res_sum_before,niter)
        niter= niter + 1
        Call TDMA_1(3)
		call Convergence_Criteria(3,Res_sum_After,niter)

    if((abs(Res_sum_before-Res_sum_After).Ge.1.0E-10) &
    .and.(niter.le.1500).and.(abs(Res_sum_After).ge.1.0E-06))go to 20		 
	write(*,*)'Pressure correction Res_sum_before-Res_sum_After' &
    ,Res_sum_before-Res_sum_After,niter
!  pressure correction	
	Do 301 I=2,NXmax
	Do 301 J=2,NYmax
        F(i,j,4) = F(i,j,4) + 0.2 *   F(i,j,3) 
 301 continue
!  velocities and pressure correction
	Do 302 I=2,NXmax
	Do 302 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 
 302 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) 
999 continue

Return
End

My wiki