https://www.cfd-online.com/W/index.php?title=Solve_UV.f90_-_Solution_of_the_momentum_equations_for_U_and_V&feed=atom&action=history
Solve UV.f90 - Solution of the momentum equations for U and V - Revision history
2024-03-28T08:14:24Z
Revision history for this page on the wiki
MediaWiki 1.16.5
https://www.cfd-online.com/W/index.php?title=Solve_UV.f90_-_Solution_of_the_momentum_equations_for_U_and_V&diff=23616&oldid=prev
Michail at 14:58, 19 May 2016
2016-05-19T14:58:47Z
<p></p>
<a href="https://www.cfd-online.com/W/index.php?title=Solve_UV.f90_-_Solution_of_the_momentum_equations_for_U_and_V&diff=23616&oldid=10342">Show changes</a>
Michail
https://www.cfd-online.com/W/index.php?title=Solve_UV.f90_-_Solution_of_the_momentum_equations_for_U_and_V&diff=10342&oldid=prev
Michail at 13:34, 4 May 2010
2010-05-04T13:34:01Z
<p></p>
<table style="background-color: white; color:black;">
<col class='diff-marker' />
<col class='diff-content' />
<col class='diff-marker' />
<col class='diff-content' />
<tr valign='top'>
<td colspan='2' style="background-color: white; color:black;">← Older revision</td>
<td colspan='2' style="background-color: white; color:black;">Revision as of 13:34, 4 May 2010</td>
</tr><tr><td colspan="2" class="diff-lineno">Line 53:</td>
<td colspan="2" class="diff-lineno">Line 53:</td></tr>
<tr><td class='diff-marker'> </td><td style="background: #eee; color:black; font-size: smaller;"><div>! upwind differencing (all other will be included into the source term) </div></td><td class='diff-marker'> </td><td style="background: #eee; color:black; font-size: smaller;"><div>! upwind differencing (all other will be included into the source term) </div></td></tr>
<tr><td class='diff-marker'> </td><td style="background: #eee; color:black; font-size: smaller;"></td><td class='diff-marker'> </td><td style="background: #eee; color:black; font-size: smaller;"></td></tr>
<tr><td class='diff-marker'>-</td><td style="background: #ffa; color:black; font-size: smaller;"><div><del class="diffchange diffchange-inline"> </del>Conv_w = Area_w * ( F(i,j,1) + F(i-1,j ,1) ) * 0.5 </div></td><td class='diff-marker'>+</td><td style="background: #cfc; color:black; font-size: smaller;"><div><ins class="diffchange diffchange-inline"> </ins>Conv_w = Area_w * ( F(i,j,1) + F(i-1,j ,1) ) * 0.5 </div></td></tr>
<tr><td class='diff-marker'>-</td><td style="background: #ffa; color:black; font-size: smaller;"><div><del class="diffchange diffchange-inline"> </del>Conv_e = Area_e * ( F(i,j,1) + F(i+1,j ,1) ) * 0.5</div></td><td class='diff-marker'>+</td><td style="background: #cfc; color:black; font-size: smaller;"><div><ins class="diffchange diffchange-inline"> </ins>Conv_e = Area_e * ( F(i,j,1) + F(i+1,j ,1) ) * 0.5</div></td></tr>
<tr><td class='diff-marker'> </td><td style="background: #eee; color:black; font-size: smaller;"><div> Conv_s = Area_s * ( F(i,j,2) + F(i ,j-1,2) ) * 0.5</div></td><td class='diff-marker'> </td><td style="background: #eee; color:black; font-size: smaller;"><div> Conv_s = Area_s * ( F(i,j,2) + F(i ,j-1,2) ) * 0.5</div></td></tr>
<tr><td class='diff-marker'> </td><td style="background: #eee; color:black; font-size: smaller;"><div> Conv_n = Area_n * ( F(i,j,2) + F(i ,j+1,2) ) * 0.5</div></td><td class='diff-marker'> </td><td style="background: #eee; color:black; font-size: smaller;"><div> Conv_n = Area_n * ( F(i,j,2) + F(i ,j+1,2) ) * 0.5</div></td></tr>
<tr><td colspan="2" class="diff-lineno">Line 60:</td>
<td colspan="2" class="diff-lineno">Line 60:</td></tr>
<tr><td class='diff-marker'> </td><td style="background: #eee; color:black; font-size: smaller;"><div> if(i.eq.2 )Conv_w = 0.</div></td><td class='diff-marker'> </td><td style="background: #eee; color:black; font-size: smaller;"><div> if(i.eq.2 )Conv_w = 0.</div></td></tr>
<tr><td class='diff-marker'> </td><td style="background: #eee; color:black; font-size: smaller;"><div> if(i.eq.NXmax)Conv_e = 0.</div></td><td class='diff-marker'> </td><td style="background: #eee; color:black; font-size: smaller;"><div> if(i.eq.NXmax)Conv_e = 0.</div></td></tr>
<tr><td class='diff-marker'>-</td><td style="background: #ffa; color:black; font-size: smaller;"><div><del class="diffchange diffchange-inline"> </del>if(j.eq.2 )Conv_s = 0.</div></td><td class='diff-marker'>+</td><td style="background: #cfc; color:black; font-size: smaller;"><div><ins class="diffchange diffchange-inline"> </ins>if(j.eq.2 )Conv_s = 0.</div></td></tr>
<tr><td class='diff-marker'> </td><td style="background: #eee; color:black; font-size: smaller;"><div> if(j.eq.NYmax)Conv_n = 0.</div></td><td class='diff-marker'> </td><td style="background: #eee; color:black; font-size: smaller;"><div> if(j.eq.NYmax)Conv_n = 0.</div></td></tr>
<tr><td class='diff-marker'> </td><td style="background: #eee; color:black; font-size: smaller;"></td><td class='diff-marker'> </td><td style="background: #eee; color:black; font-size: smaller;"></td></tr>
</table>
Michail
https://www.cfd-online.com/W/index.php?title=Solve_UV.f90_-_Solution_of_the_momentum_equations_for_U_and_V&diff=10312&oldid=prev
Michail: New page: <pre> !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 !This progra...
2010-05-03T08:50:54Z
<p>New page: <pre> !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 !This progra...</p>
<p><b>New page</b></p><div><pre><br />
<br />
!Sample program for solving Lid-driven cavity flow test using SIMPLE-algorithm<br />
! solution of momentum equation for U and V modul<br />
!Copyright (C) 2010 Michail Kiričkov<br />
<br />
!This program is free software; you can redistribute it and/or<br />
!modify it under the terms of the GNU General Public License<br />
!as published by the Free Software Foundation; either version 2<br />
!of the License, or (at your option) any later version.<br />
<br />
!This program is distributed in the hope that it will be useful,<br />
!but WITHOUT ANY WARRANTY; without even the implied warranty of<br />
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the<br />
!GNU General Public License for more details.<br />
<br />
!You should have received a copy of the GNU General Public License<br />
!along with this program; if not, write to the Free Software<br />
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.<br />
<br />
!**********************************************************************<br />
Subroutine Solve_UV<br />
include 'icomm_1.f90'<br />
<br />
<br />
! calculation of fluxes<br />
! all geometry has rectangular 2D notation<br />
<br />
Do 100 I= 2,NXmax<br />
Do 100 J= 2,NYmax<br />
<br />
Gam_e = ( Gam(i+1,j ) + Gam(i ,j ) ) * 0.5<br />
Gam_w = ( Gam(i-1,j ) + Gam(i ,j ) ) * 0.5<br />
Gam_s = ( Gam(i ,j-1) + Gam(i ,j ) ) * 0.5<br />
Gam_n = ( Gam(i ,j+1) + Gam(i ,j ) ) * 0.5<br />
<br />
<br />
!----------------------------------------------------<br />
Area_w = Y(i-1,j)-Y(i-1,j-1) <br />
Area_e = Y(i ,j)-Y(i ,j-1) <br />
<br />
Area_s = X(i,j-1)-X(i-1,j-1) <br />
Area_n = X(i,j )-X(i-1,j ) <br />
!----------------------------------------------------<br />
<br />
Del_w = Xc(i ,j)-Xc(i-1,j) <br />
Del_e = Xc(i+1,j)-Xc(i ,j) <br />
<br />
Del_s = Yc(i,j )-Yc(i,j-1) <br />
Del_n = Yc(i,j+1)-Yc(i,j ) <br />
!----------------------------------------------------<br />
<br />
! upwind differencing (all other will be included into the source term) <br />
<br />
Conv_w = Area_w * ( F(i,j,1) + F(i-1,j ,1) ) * 0.5 <br />
Conv_e = Area_e * ( F(i,j,1) + F(i+1,j ,1) ) * 0.5<br />
Conv_s = Area_s * ( F(i,j,2) + F(i ,j-1,2) ) * 0.5<br />
Conv_n = Area_n * ( F(i,j,2) + F(i ,j+1,2) ) * 0.5<br />
<br />
if(i.eq.2 )Conv_w = 0.<br />
if(i.eq.NXmax)Conv_e = 0.<br />
if(j.eq.2 )Conv_s = 0.<br />
if(j.eq.NYmax)Conv_n = 0.<br />
<br />
Diff_e = Area_e * Gam_e / Del_e<br />
Diff_w = Area_w * Gam_w / Del_w <br />
Diff_s = Area_s * Gam_s / Del_s<br />
Diff_n = Area_n * Gam_n / Del_n <br />
<br />
<br />
Aw(i,j) = Diff_w + max( Conv_w,0.)<br />
Ae(i,j) = Diff_e + max(-1.* Conv_e,0.)<br />
As(i,j) = Diff_s + max( Conv_s,0.)<br />
An(i,j) = Diff_n + max(-1.* Conv_n,0.)<br />
<br />
Ap(i,j,1:2)= Aw(i,j) + Ae(i,j) + An(i,j) + As(i,j) <br />
<br />
Sp(i,j,1:2)= 0.<br />
<br />
!-------------------------------- HLPA SCHEME----------------------------<br />
! go to 600 ! (now HLPA is "off")<br />
<br />
DO 500 nf=1,2<br />
<br />
! Subroutine HLPA(Uw,Fww,Fw,Fp,Fe,Delta_f)<br />
<br />
if( (i.GT.2).AND.(i.LT.NXmax-0).and.(j.GT.2).AND.(j.LT.NYmax-0) ) then <br />
<br />
!------------------ w face -------------------<br />
Fww = F(i-2,j,nf)<br />
Fw = F(i-1,j,nf)<br />
Fp = F(i ,j,nf)<br />
Fe = F(i+1,j,nf)<br />
<br />
call HLPA(Conv_w,Fww,Fw,Fp,Fe,Delta_f)<br />
<br />
Sp(i,j,nf) = Sp(i,j,nf) + Conv_w * Delta_f<br />
<br />
!------------------ e face--------------------<br />
<br />
Fww = F(i-1,j,nf)<br />
Fw = F(i ,j,nf)<br />
Fp = F(i+1,j,nf)<br />
Fe = F(i+2,j,nf)<br />
<br />
call HLPA(Conv_e,Fww,Fw,Fp,Fe,Delta_f)<br />
<br />
Sp(i,j,nf) = Sp(i,j,nf) + Conv_e * Delta_f * (-1.)<br />
<br />
!------------------ s face--------------------<br />
Fww = F(i ,j-2,nf)<br />
Fw = F(i ,j-1,nf)<br />
Fp = F(i ,j ,nf)<br />
Fe = F(i ,j+1,nf)<br />
<br />
call HLPA(Conv_s,Fww,Fw,Fp,Fe,Delta_f)<br />
<br />
Sp(i,j,nf) = Sp(i,j,nf) + Conv_s * Delta_f<br />
<br />
!------------------ n face--------------------<br />
<br />
Fww = F(i ,j-1,nf)<br />
Fw = F(i ,j ,nf)<br />
Fp = F(i ,j+1,nf)<br />
Fe = F(i ,j+2,nf)<br />
<br />
call HLPA(Conv_n,Fww,Fw,Fp,Fe,Delta_f)<br />
<br />
Sp(i,j,nf) = Sp(i,j,nf) + Conv_n * Delta_f *(-1.)<br />
<br />
<br />
end if <br />
<br />
500 continue <br />
<br />
<br />
600 continue <br />
<br />
<br />
!------------------------------------------------------------------------<br />
<br />
100 continue ! coefficient cycle<br />
<br />
!----------------------------- pressure gradient ------------------------<br />
<br />
Do 200 I= 2,NXmax<br />
Do 200 J= 2,NYmax<br />
<br />
DX = X(i,j) - X(i-1,j)<br />
DY = Y(i,j) - Y(i,j-1)<br />
<br />
VOL = DX * DY<br />
<br />
PE = ( F(i,j,4) + F(i+1,j,4) ) * 0.5<br />
PW = ( F(i,j,4) + F(i-1,j,4) ) * 0.5<br />
PN = ( F(i,j,4) + F(i,j+1,4) ) * 0.5<br />
PS = ( F(i,j,4) + F(i,j-1,4) ) * 0.5<br />
<br />
DPx_c(i,j) = (PE-PW)/DX<br />
DPy_c(i,j) = (PN-PS)/DY<br />
<br />
Sp(i,j,1) = Sp(i,j,1) - DPx_c(i,j) * VOL<br />
Sp(i,j,2) = Sp(i,j,2) - DPy_c(i,j) * VOL<br />
<br />
200 continue<br />
<br />
!---------------------------- under-relaxation ---------------------------------<br />
<br />
Alfa = 0.85<br />
Urf = 1. / Alfa<br />
<br />
<br />
Ap(1:NXmaxC,1:NYmaxC,1) = Ap(1:NXmaxC,1:NYmaxC,1) * Urf<br />
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 <br />
<br />
Ap(1:NXmaxC,1:NYmaxC,2) = Ap(1:NXmaxC,1:NYmaxC,2) * Urf<br />
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 <br />
<br />
!---------------------------------------------------------------------------------<br />
<br />
!******************************************************************<br />
<br />
niter = 0<br />
write(*,*)'solve U'<br />
call Convergence_Criteria(1,Res_sum_before)<br />
<br />
10 continue<br />
niter= niter + 1<br />
Call TDMA_1(1)<br />
call Convergence_Criteria(1,Res_sum_After)<br />
If((abs(Res_sum_before-Res_sum_After).Ge.0.00000000001).and.niter.le.20)then<br />
<br />
Res_sum_before = Res_sum_After<br />
go to 10 <br />
End if<br />
<br />
<br />
<br />
niter = 0<br />
write(*,*)'solve V'<br />
call Convergence_Criteria(2,Res_sum_before)<br />
<br />
20 continue<br />
<br />
niter= niter + 1<br />
Call TDMA_1(2)<br />
call Convergence_Criteria(2,Res_sum_After)<br />
If((abs(Res_sum_before-Res_sum_After).Ge.0.00000000001).and.niter.le.20)then<br />
<br />
Res_sum_before = Res_sum_After<br />
go to 20 <br />
End if<br />
<br />
!***********************************************************************<br />
<br />
!---------------------------------------------------------------------------------<br />
Return<br />
End<br />
<br />
</pre></div>
Michail