Solve UV.f90 - Solution of the momentum equations for U and V
From CFD-Wiki
(Difference between revisions)
Line 3: | Line 3: | ||
!Sample program for solving Lid-driven cavity flow test using SIMPLE-algorithm | !Sample program for solving Lid-driven cavity flow test using SIMPLE-algorithm | ||
! solution of momentum equation for U and V modul | ! solution of momentum equation for U and V modul | ||
- | !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_UV | Subroutine Solve_UV | ||
include 'icomm_1.f90' | include 'icomm_1.f90' | ||
- | + | Istoch_nar = 0.; | |
- | + | ||
! calculation of fluxes | ! calculation of fluxes | ||
! all geometry has rectangular 2D notation | ! all geometry has rectangular 2D notation | ||
- | |||
Do 100 I= 2,NXmax | Do 100 I= 2,NXmax | ||
Do 100 J= 2,NYmax | Do 100 J= 2,NYmax | ||
- | |||
Gam_e = ( Gam(i+1,j ) + Gam(i ,j ) ) * 0.5 | Gam_e = ( Gam(i+1,j ) + Gam(i ,j ) ) * 0.5 | ||
Gam_w = ( Gam(i-1,j ) + Gam(i ,j ) ) * 0.5 | Gam_w = ( Gam(i-1,j ) + Gam(i ,j ) ) * 0.5 | ||
Gam_s = ( Gam(i ,j-1) + Gam(i ,j ) ) * 0.5 | Gam_s = ( Gam(i ,j-1) + Gam(i ,j ) ) * 0.5 | ||
Gam_n = ( Gam(i ,j+1) + Gam(i ,j ) ) * 0.5 | Gam_n = ( Gam(i ,j+1) + Gam(i ,j ) ) * 0.5 | ||
- | |||
- | |||
!---------------------------------------------------- | !---------------------------------------------------- | ||
Area_w = Y(i-1,j)-Y(i-1,j-1) | Area_w = Y(i-1,j)-Y(i-1,j-1) | ||
Area_e = Y(i ,j)-Y(i ,j-1) | Area_e = Y(i ,j)-Y(i ,j-1) | ||
- | |||
Area_s = X(i,j-1)-X(i-1,j-1) | Area_s = X(i,j-1)-X(i-1,j-1) | ||
Area_n = X(i,j )-X(i-1,j ) | Area_n = X(i,j )-X(i-1,j ) | ||
!---------------------------------------------------- | !---------------------------------------------------- | ||
- | |||
Del_w = Xc(i ,j)-Xc(i-1,j) | Del_w = Xc(i ,j)-Xc(i-1,j) | ||
Del_e = Xc(i+1,j)-Xc(i ,j) | Del_e = Xc(i+1,j)-Xc(i ,j) | ||
- | |||
Del_s = Yc(i,j )-Yc(i,j-1) | Del_s = Yc(i,j )-Yc(i,j-1) | ||
Del_n = Yc(i,j+1)-Yc(i,j ) | Del_n = Yc(i,j+1)-Yc(i,j ) | ||
!---------------------------------------------------- | !---------------------------------------------------- | ||
- | |||
! upwind differencing (all other will be included into the source term) | ! upwind differencing (all other will be included into the source term) | ||
- | + | Conv_w = Area_w * ( F(i,j,1) + F(i-1,j ,1) ) * 0.5 | |
- | + | Conv_e = Area_e * ( F(i,j,1) + F(i+1,j ,1) ) * 0.5 | |
- | + | ||
Conv_s = Area_s * ( F(i,j,2) + F(i ,j-1,2) ) * 0.5 | Conv_s = Area_s * ( F(i,j,2) + F(i ,j-1,2) ) * 0.5 | ||
Conv_n = Area_n * ( F(i,j,2) + F(i ,j+1,2) ) * 0.5 | Conv_n = Area_n * ( F(i,j,2) + F(i ,j+1,2) ) * 0.5 | ||
- | |||
if(i.eq.2 )Conv_w = 0. | if(i.eq.2 )Conv_w = 0. | ||
if(i.eq.NXmax)Conv_e = 0. | if(i.eq.NXmax)Conv_e = 0. | ||
- | + | if(j.eq.2 )Conv_s = 0. | |
if(j.eq.NYmax)Conv_n = 0. | if(j.eq.NYmax)Conv_n = 0. | ||
- | |||
Diff_e = Area_e * Gam_e / Del_e | Diff_e = Area_e * Gam_e / Del_e | ||
Diff_w = Area_w * Gam_w / Del_w | Diff_w = Area_w * Gam_w / Del_w | ||
Diff_s = Area_s * Gam_s / Del_s | Diff_s = Area_s * Gam_s / Del_s | ||
Diff_n = Area_n * Gam_n / Del_n | Diff_n = Area_n * Gam_n / Del_n | ||
- | |||
- | |||
Aw(i,j) = Diff_w + max( Conv_w,0.) | Aw(i,j) = Diff_w + max( Conv_w,0.) | ||
Ae(i,j) = Diff_e + max(-1.* Conv_e,0.) | Ae(i,j) = Diff_e + max(-1.* Conv_e,0.) | ||
As(i,j) = Diff_s + max( Conv_s,0.) | As(i,j) = Diff_s + max( Conv_s,0.) | ||
An(i,j) = Diff_n + max(-1.* Conv_n,0.) | An(i,j) = Diff_n + max(-1.* Conv_n,0.) | ||
- | + | Check_Flux(i,j) = Conv_e - Conv_w + Conv_s - Conv_n | |
- | Ap(i,j,1:2)= Aw(i,j) + Ae(i,j) + An(i,j) + As(i,j) | + | Ap(i,j,1:2)= Aw(i,j) + Ae(i,j) + An(i,j) + As(i,j)! + Check_flux(i,j) |
- | + | ||
Sp(i,j,1:2)= 0. | Sp(i,j,1:2)= 0. | ||
- | + | Istoch_nar_abs = Istoch_nar_abs + ABS(Check_Flux(i,j)) | |
+ | Istoch_nar_whole = Istoch_nar_whole + (Check_Flux(i,j)) | ||
!-------------------------------- HLPA SCHEME---------------------------- | !-------------------------------- HLPA SCHEME---------------------------- | ||
- | + | go to 600 ! (now HLPA is "off") | |
- | + | ||
DO 500 nf=1,2 | DO 500 nf=1,2 | ||
- | |||
! Subroutine HLPA(Uw,Fww,Fw,Fp,Fe,Delta_f) | ! Subroutine HLPA(Uw,Fww,Fw,Fp,Fe,Delta_f) | ||
- | |||
if( (i.GT.2).AND.(i.LT.NXmax-0).and.(j.GT.2).AND.(j.LT.NYmax-0) ) then | if( (i.GT.2).AND.(i.LT.NXmax-0).and.(j.GT.2).AND.(j.LT.NYmax-0) ) then | ||
- | |||
!------------------ w face ------------------- | !------------------ w face ------------------- | ||
Fww = F(i-2,j,nf) | Fww = F(i-2,j,nf) | ||
Line 92: | Line 75: | ||
Fp = F(i ,j,nf) | Fp = F(i ,j,nf) | ||
Fe = F(i+1,j,nf) | Fe = F(i+1,j,nf) | ||
- | |||
call HLPA(Conv_w,Fww,Fw,Fp,Fe,Delta_f) | call HLPA(Conv_w,Fww,Fw,Fp,Fe,Delta_f) | ||
- | |||
Sp(i,j,nf) = Sp(i,j,nf) + Conv_w * Delta_f | Sp(i,j,nf) = Sp(i,j,nf) + Conv_w * Delta_f | ||
- | |||
!------------------ e face-------------------- | !------------------ e face-------------------- | ||
- | |||
Fww = F(i-1,j,nf) | Fww = F(i-1,j,nf) | ||
Fw = F(i ,j,nf) | Fw = F(i ,j,nf) | ||
Fp = F(i+1,j,nf) | Fp = F(i+1,j,nf) | ||
Fe = F(i+2,j,nf) | Fe = F(i+2,j,nf) | ||
- | |||
call HLPA(Conv_e,Fww,Fw,Fp,Fe,Delta_f) | call HLPA(Conv_e,Fww,Fw,Fp,Fe,Delta_f) | ||
- | |||
Sp(i,j,nf) = Sp(i,j,nf) + Conv_e * Delta_f * (-1.) | Sp(i,j,nf) = Sp(i,j,nf) + Conv_e * Delta_f * (-1.) | ||
- | |||
!------------------ s face-------------------- | !------------------ s face-------------------- | ||
Fww = F(i ,j-2,nf) | Fww = F(i ,j-2,nf) | ||
Line 113: | Line 89: | ||
Fp = F(i ,j ,nf) | Fp = F(i ,j ,nf) | ||
Fe = F(i ,j+1,nf) | Fe = F(i ,j+1,nf) | ||
- | |||
call HLPA(Conv_s,Fww,Fw,Fp,Fe,Delta_f) | call HLPA(Conv_s,Fww,Fw,Fp,Fe,Delta_f) | ||
- | |||
Sp(i,j,nf) = Sp(i,j,nf) + Conv_s * Delta_f | Sp(i,j,nf) = Sp(i,j,nf) + Conv_s * Delta_f | ||
- | |||
!------------------ n face-------------------- | !------------------ n face-------------------- | ||
- | |||
Fww = F(i ,j-1,nf) | Fww = F(i ,j-1,nf) | ||
Fw = F(i ,j ,nf) | Fw = F(i ,j ,nf) | ||
Fp = F(i ,j+1,nf) | Fp = F(i ,j+1,nf) | ||
Fe = F(i ,j+2,nf) | Fe = F(i ,j+2,nf) | ||
- | |||
call HLPA(Conv_n,Fww,Fw,Fp,Fe,Delta_f) | call HLPA(Conv_n,Fww,Fw,Fp,Fe,Delta_f) | ||
- | |||
Sp(i,j,nf) = Sp(i,j,nf) + Conv_n * Delta_f *(-1.) | Sp(i,j,nf) = Sp(i,j,nf) + Conv_n * Delta_f *(-1.) | ||
- | |||
- | |||
end if | end if | ||
- | |||
500 continue | 500 continue | ||
- | |||
- | |||
600 continue | 600 continue | ||
- | |||
- | |||
!------------------------------------------------------------------------ | !------------------------------------------------------------------------ | ||
- | |||
100 continue ! coefficient cycle | 100 continue ! coefficient cycle | ||
- | |||
!----------------------------- pressure gradient ------------------------ | !----------------------------- pressure gradient ------------------------ | ||
- | |||
Do 200 I= 2,NXmax | Do 200 I= 2,NXmax | ||
Do 200 J= 2,NYmax | Do 200 J= 2,NYmax | ||
- | |||
DX = X(i,j) - X(i-1,j) | DX = X(i,j) - X(i-1,j) | ||
DY = Y(i,j) - Y(i,j-1) | DY = Y(i,j) - Y(i,j-1) | ||
- | |||
VOL = DX * DY | VOL = DX * DY | ||
- | |||
PE = ( F(i,j,4) + F(i+1,j,4) ) * 0.5 | PE = ( F(i,j,4) + F(i+1,j,4) ) * 0.5 | ||
PW = ( F(i,j,4) + F(i-1,j,4) ) * 0.5 | PW = ( F(i,j,4) + F(i-1,j,4) ) * 0.5 | ||
PN = ( F(i,j,4) + F(i,j+1,4) ) * 0.5 | PN = ( F(i,j,4) + F(i,j+1,4) ) * 0.5 | ||
PS = ( F(i,j,4) + F(i,j-1,4) ) * 0.5 | PS = ( F(i,j,4) + F(i,j-1,4) ) * 0.5 | ||
- | |||
DPx_c(i,j) = (PE-PW)/DX | DPx_c(i,j) = (PE-PW)/DX | ||
DPy_c(i,j) = (PN-PS)/DY | DPy_c(i,j) = (PN-PS)/DY | ||
- | |||
Sp(i,j,1) = Sp(i,j,1) - DPx_c(i,j) * VOL | Sp(i,j,1) = Sp(i,j,1) - DPx_c(i,j) * VOL | ||
Sp(i,j,2) = Sp(i,j,2) - DPy_c(i,j) * VOL | Sp(i,j,2) = Sp(i,j,2) - DPy_c(i,j) * VOL | ||
- | |||
200 continue | 200 continue | ||
- | |||
!---------------------------- under-relaxation --------------------------------- | !---------------------------- under-relaxation --------------------------------- | ||
- | + | Alfa = 0.8 | |
- | Alfa = 0. | + | |
Urf = 1. / Alfa | Urf = 1. / Alfa | ||
- | |||
- | |||
Ap(1:NXmaxC,1:NYmaxC,1) = Ap(1:NXmaxC,1:NYmaxC,1) * Urf | Ap(1:NXmaxC,1:NYmaxC,1) = Ap(1:NXmaxC,1:NYmaxC,1) * Urf | ||
- | Sp(1:NXmaxC,1:NYmaxC,1) = Sp(1:NXmaxC,1:NYmaxC,1) + (1. - Alfa )* Ap(1:NXmaxC,1:NYmaxC,1)*F(1:NXmaxC,1:NYmaxC,1) ! / Alfa | + | Sp(1:NXmaxC,1:NYmaxC,1) = Sp(1:NXmaxC,1:NYmaxC,1) + (1. - Alfa )* & |
- | + | Ap(1:NXmaxC,1:NYmaxC,1)*F(1:NXmaxC,1:NYmaxC,1) ! / Alfa | |
Ap(1:NXmaxC,1:NYmaxC,2) = Ap(1:NXmaxC,1:NYmaxC,2) * Urf | Ap(1:NXmaxC,1:NYmaxC,2) = Ap(1:NXmaxC,1:NYmaxC,2) * Urf | ||
- | Sp(1:NXmaxC,1:NYmaxC,2) = Sp(1:NXmaxC,1:NYmaxC,2) + (1. - Alfa )* Ap(1:NXmaxC,1:NYmaxC,2)*F(1:NXmaxC,1:NYmaxC,2) ! / Alfa | + | Sp(1:NXmaxC,1:NYmaxC,2) = Sp(1:NXmaxC,1:NYmaxC,2) + (1. - Alfa )* & |
- | + | Ap(1:NXmaxC,1:NYmaxC,2)*F(1:NXmaxC,1:NYmaxC,2) ! / Alfa | |
!--------------------------------------------------------------------------------- | !--------------------------------------------------------------------------------- | ||
- | |||
- | |||
- | |||
niter = 0 | niter = 0 | ||
write(*,*)'solve U' | write(*,*)'solve U' | ||
- | call Convergence_Criteria(1,Res_sum_before) | + | call Convergence_Criteria(1,Res_sum_before,niter) |
- | + | ||
10 continue | 10 continue | ||
niter= niter + 1 | niter= niter + 1 | ||
Call TDMA_1(1) | Call TDMA_1(1) | ||
- | call Convergence_Criteria(1,Res_sum_After) | + | call Convergence_Criteria(1,Res_sum_After,niter) |
- | If((abs(Res_sum_before-Res_sum_After).Ge. | + | If((abs(Res_sum_before-Res_sum_After).Ge.1.00E-10).and.niter.le.550)then |
- | + | ||
Res_sum_before = Res_sum_After | Res_sum_before = Res_sum_After | ||
go to 10 | go to 10 | ||
End if | End if | ||
- | |||
- | |||
niter = 0 | niter = 0 | ||
write(*,*)'solve V' | write(*,*)'solve V' | ||
- | call Convergence_Criteria(2,Res_sum_before) | + | call Convergence_Criteria(2,Res_sum_before,niter) |
- | + | ||
20 continue | 20 continue | ||
- | |||
niter= niter + 1 | niter= niter + 1 | ||
Call TDMA_1(2) | Call TDMA_1(2) | ||
- | call Convergence_Criteria(2,Res_sum_After) | + | call Convergence_Criteria(2,Res_sum_After,niter) |
- | If((abs(Res_sum_before-Res_sum_After).Ge. | + | If((abs(Res_sum_before-Res_sum_After).Ge.1.00E-10).and.niter.le.550)then |
- | + | ||
Res_sum_before = Res_sum_After | Res_sum_before = Res_sum_After | ||
go to 20 | go to 20 | ||
End if | End if | ||
- | |||
- | |||
- | |||
!--------------------------------------------------------------------------------- | !--------------------------------------------------------------------------------- | ||
Return | Return |
Latest revision as of 14:58, 19 May 2016
!Sample program for solving Lid-driven cavity flow test using SIMPLE-algorithm ! solution of momentum equation for U and V modul ! 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_UV include 'icomm_1.f90' Istoch_nar = 0.; ! calculation of fluxes ! all geometry has rectangular 2D notation Do 100 I= 2,NXmax Do 100 J= 2,NYmax Gam_e = ( Gam(i+1,j ) + Gam(i ,j ) ) * 0.5 Gam_w = ( Gam(i-1,j ) + Gam(i ,j ) ) * 0.5 Gam_s = ( Gam(i ,j-1) + Gam(i ,j ) ) * 0.5 Gam_n = ( Gam(i ,j+1) + Gam(i ,j ) ) * 0.5 !---------------------------------------------------- Area_w = Y(i-1,j)-Y(i-1,j-1) Area_e = Y(i ,j)-Y(i ,j-1) Area_s = X(i,j-1)-X(i-1,j-1) Area_n = X(i,j )-X(i-1,j ) !---------------------------------------------------- Del_w = Xc(i ,j)-Xc(i-1,j) Del_e = Xc(i+1,j)-Xc(i ,j) Del_s = Yc(i,j )-Yc(i,j-1) Del_n = Yc(i,j+1)-Yc(i,j ) !---------------------------------------------------- ! upwind differencing (all other will be included into the source term) Conv_w = Area_w * ( F(i,j,1) + F(i-1,j ,1) ) * 0.5 Conv_e = Area_e * ( F(i,j,1) + F(i+1,j ,1) ) * 0.5 Conv_s = Area_s * ( F(i,j,2) + F(i ,j-1,2) ) * 0.5 Conv_n = Area_n * ( F(i,j,2) + F(i ,j+1,2) ) * 0.5 if(i.eq.2 )Conv_w = 0. if(i.eq.NXmax)Conv_e = 0. if(j.eq.2 )Conv_s = 0. if(j.eq.NYmax)Conv_n = 0. Diff_e = Area_e * Gam_e / Del_e Diff_w = Area_w * Gam_w / Del_w Diff_s = Area_s * Gam_s / Del_s Diff_n = Area_n * Gam_n / Del_n Aw(i,j) = Diff_w + max( Conv_w,0.) Ae(i,j) = Diff_e + max(-1.* Conv_e,0.) As(i,j) = Diff_s + max( Conv_s,0.) An(i,j) = Diff_n + max(-1.* Conv_n,0.) Check_Flux(i,j) = Conv_e - Conv_w + Conv_s - Conv_n Ap(i,j,1:2)= Aw(i,j) + Ae(i,j) + An(i,j) + As(i,j)! + Check_flux(i,j) Sp(i,j,1:2)= 0. Istoch_nar_abs = Istoch_nar_abs + ABS(Check_Flux(i,j)) Istoch_nar_whole = Istoch_nar_whole + (Check_Flux(i,j)) !-------------------------------- HLPA SCHEME---------------------------- go to 600 ! (now HLPA is "off") DO 500 nf=1,2 ! Subroutine HLPA(Uw,Fww,Fw,Fp,Fe,Delta_f) if( (i.GT.2).AND.(i.LT.NXmax-0).and.(j.GT.2).AND.(j.LT.NYmax-0) ) then !------------------ w face ------------------- Fww = F(i-2,j,nf) Fw = F(i-1,j,nf) Fp = F(i ,j,nf) Fe = F(i+1,j,nf) call HLPA(Conv_w,Fww,Fw,Fp,Fe,Delta_f) Sp(i,j,nf) = Sp(i,j,nf) + Conv_w * Delta_f !------------------ e face-------------------- Fww = F(i-1,j,nf) Fw = F(i ,j,nf) Fp = F(i+1,j,nf) Fe = F(i+2,j,nf) call HLPA(Conv_e,Fww,Fw,Fp,Fe,Delta_f) Sp(i,j,nf) = Sp(i,j,nf) + Conv_e * Delta_f * (-1.) !------------------ s face-------------------- Fww = F(i ,j-2,nf) Fw = F(i ,j-1,nf) Fp = F(i ,j ,nf) Fe = F(i ,j+1,nf) call HLPA(Conv_s,Fww,Fw,Fp,Fe,Delta_f) Sp(i,j,nf) = Sp(i,j,nf) + Conv_s * Delta_f !------------------ n face-------------------- Fww = F(i ,j-1,nf) Fw = F(i ,j ,nf) Fp = F(i ,j+1,nf) Fe = F(i ,j+2,nf) call HLPA(Conv_n,Fww,Fw,Fp,Fe,Delta_f) Sp(i,j,nf) = Sp(i,j,nf) + Conv_n * Delta_f *(-1.) end if 500 continue 600 continue !------------------------------------------------------------------------ 100 continue ! coefficient cycle !----------------------------- pressure gradient ------------------------ Do 200 I= 2,NXmax Do 200 J= 2,NYmax DX = X(i,j) - X(i-1,j) DY = Y(i,j) - Y(i,j-1) VOL = DX * DY PE = ( F(i,j,4) + F(i+1,j,4) ) * 0.5 PW = ( F(i,j,4) + F(i-1,j,4) ) * 0.5 PN = ( F(i,j,4) + F(i,j+1,4) ) * 0.5 PS = ( F(i,j,4) + F(i,j-1,4) ) * 0.5 DPx_c(i,j) = (PE-PW)/DX DPy_c(i,j) = (PN-PS)/DY Sp(i,j,1) = Sp(i,j,1) - DPx_c(i,j) * VOL Sp(i,j,2) = Sp(i,j,2) - DPy_c(i,j) * VOL 200 continue !---------------------------- under-relaxation --------------------------------- Alfa = 0.8 Urf = 1. / Alfa Ap(1:NXmaxC,1:NYmaxC,1) = Ap(1:NXmaxC,1:NYmaxC,1) * Urf Sp(1:NXmaxC,1:NYmaxC,1) = Sp(1:NXmaxC,1:NYmaxC,1) + (1. - Alfa )* & Ap(1:NXmaxC,1:NYmaxC,1)*F(1:NXmaxC,1:NYmaxC,1) ! / Alfa Ap(1:NXmaxC,1:NYmaxC,2) = Ap(1:NXmaxC,1:NYmaxC,2) * Urf Sp(1:NXmaxC,1:NYmaxC,2) = Sp(1:NXmaxC,1:NYmaxC,2) + (1. - Alfa )* & Ap(1:NXmaxC,1:NYmaxC,2)*F(1:NXmaxC,1:NYmaxC,2) ! / Alfa !--------------------------------------------------------------------------------- niter = 0 write(*,*)'solve U' call Convergence_Criteria(1,Res_sum_before,niter) 10 continue niter= niter + 1 Call TDMA_1(1) call Convergence_Criteria(1,Res_sum_After,niter) If((abs(Res_sum_before-Res_sum_After).Ge.1.00E-10).and.niter.le.550)then Res_sum_before = Res_sum_After go to 10 End if niter = 0 write(*,*)'solve V' call Convergence_Criteria(2,Res_sum_before,niter) 20 continue niter= niter + 1 Call TDMA_1(2) call Convergence_Criteria(2,Res_sum_After,niter) If((abs(Res_sum_before-Res_sum_After).Ge.1.00E-10).and.niter.le.550)then Res_sum_before = Res_sum_After go to 20 End if !--------------------------------------------------------------------------------- Return End