Solve Pressure Correction.f90 - Solution of pressure-correction equation and correction of U,V and P
From CFD-Wiki
(Difference between revisions)
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) | |
- | + | 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. | |
- | + | S_e = 0. | |
- | + | VOL_e = 0. | |
- | + | APU_e = 0. | |
- | DXc_e = Xc(i+1,j+1) - Xc(i,j+1) ! | + | 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 | VOL_e = DXc_e * S_e | ||
- | |||
!------------------------------------------------------------- | !------------------------------------------------------------- | ||
- | APU_e = 0.5 * ( DpU(i,j+1) + DpU(i+1,j+1) ) !this is the | + | !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) ) ! | + | 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 | + | !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 ! | + | ! 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) | + | 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 | End If | ||
- | + | Con_we(1,j) = 0. | |
- | + | Con_we(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. | |
- | + | 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 | |
- | + | End If | |
- | + | Con_ns(i,1) = 0. | |
- | + | Con_ns(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 | |
- | + | 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) ! | |
- | + | ||
- | Sp(i,j,3) = -1. * ( | + | |
- | + | ||
- | + | ||
200 continue | 200 continue | ||
- | + | write(*,*)'Sum',summ | |
- | + | write(*,*)'solve PP' | |
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | 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) & | |
- | + | .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 | |
- | 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 | + | Do 302 J=2,NYmax |
- | Do | + | |
- | + | ||
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 | |
- | + | !********************************************************************** | |
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | ! | + | |
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