| kar1209 |
November 18, 2017 07:36 |
SIMPLE algorithm for lid-driven cavity not converging
Hi,
I am intending to solve incompressible 2D steady flow inside a lid-driven square cavity. I am using SIMPLE algorithm on a staggered grid. I have written the code in FORTRAN and the steps in my code are as follows:
- Specify the boundary conditions, specify the guess values for pressure field and velocity field(for solving the non linearities)
- fix the pressure at a cell to make the system of equations determinant
- calculate the diffusion coefficients at the faces of the staggered control volumes
- calculate the convective flow rates at the faces of staggered control volumes using the average velocities of (P_cell - E_cell, P_cell - W_ cell, etc)
- assign values to coefficients aP, aE, etc for the x momentum equation
- repeat steps 4 and 5 for y momentum equation
- solve x momentum equation and then y momentum equation(Gauss seidel convergence criteria is given in terms of max(error) where error is a 2D array containing the error (abs(u_new - u_old)) value of each cell
- after leaving the gauss seidel loop, calculate residuals of a and y momentum equations and continuity equation using the obtained u and v velocities
- calculate coefficients for pressure correction equation
- solve the system of pressure correction equations using Gauss-Seidel with a similar converge criteria
- update u, v and p and check for mass imbalance of all cells
- use the updated u, v, p for calculating momentum equation coefficients and repeat
Please find the code attached(the code attached is for Reynolds number 10). As you can see in the code my Gauss Seidel tolerance values for convergence are 1E-9 for u and v velocities and 1E-8 for pressure correction. Also my under-relaxation factors (URFs) are 0.5 for u and v and 0.01 for pressure.
When I run the code, the first outer iteration completes and in the second iteration, the Gauss Seidel for y momentum equation does not converge. I have been trying different URFs but still the problem exists.
It would be nice if I could get some suggestions.
Code:
program main
implicit none
integer :: imax, jmax, i, j, iter = 0, it
!u_temp and v_temp are the matrices with 2 less columns and 2 less rows than u and v respectively
!u_temp and v_temp temporarily store output from gauss-seidel
double precision, allocatable :: u(:,:), u_temp(:,:), u_aP(:,:), u_aE(:,:), u_aW(:,:), u_aN(:,:), u_aS(:,:), &
& u_b(:,:), u_error(:,:)
double precision, allocatable :: v(:,:), v_temp(:,:), v_aP(:,:), v_aE(:,:), v_aW(:,:), v_aN(:,:), v_aS(:,:), &
& v_b(:,:), v_error(:,:)
double precision, allocatable :: p(:,:), pprime(:,:), aP(:,:), aE(:,:), aW(:,:), aN(:,:), aS(:,:), b(:,:), psi(:,:), &
& d_u(:,:), d_v(:,:), p_temp(:,:), p_error(:,:)
double precision :: Fe, Fw, Fn, Fs, De, Dw, Dn, Ds
double precision :: H = 1.0, W = 1.0, del_x, del_y, rho = 1, mu = 0.1
double precision :: U_W = 0.0, U_E = 0.0, V_S = 0.0, V_N = 0.0, U_S, U_N, V_W, V_E
double precision :: alpha_u = 0.5, alpha_v = 0.5, alpha_p = 0.01
double precision :: u_residual, v_residual, c_residual, tolerance = 1e-4, temp, r
write(*,*),'Enter the number of cells (cells for scalars) required in the x direction'
read(*,*), imax
write(*,*),'Enter the number of cells (cells for scalars) required in the y direction'
read(*,*), jmax
del_x = W/imax
del_y = H/jmax
!diffusion coefficients for calculating coefficients in momentum equation
De = mu*del_y/del_x
Dw = mu*del_y/del_x
Dn = mu*del_x/del_y
Ds = mu*del_x/del_y
allocate(u(imax+1,jmax), u_temp(2:imax,jmax), u_aP(2:imax,jmax), u_aE(2:imax,jmax), u_aW(2:imax,jmax), &
& u_aN(2:imax,jmax), u_aS(2:imax,jmax), u_b(2:imax,jmax), d_u(imax+1,jmax), u_error(2:imax,jmax))
allocate(v(imax,jmax+1), v_temp(imax,2:jmax), v_aP(imax,2:jmax), v_aE(imax,2:jmax), v_aW(imax,2:jmax), &
& v_aN(imax,2:jmax), v_aS(imax,2:jmax), v_b(imax,2:jmax), d_v(imax,jmax+1), v_error(imax,2:jmax))
allocate(p(imax,jmax), pprime(imax,jmax), aP(imax,jmax), aE(imax,jmax), aW(imax,jmax), aN(imax,jmax), &
& aS(imax,jmax), b(imax,jmax), psi(imax,jmax), p_temp(imax,jmax), p_error(imax,jmax))
!initial guess values - p*, u_g, v_g, p'
p = 0.0
pprime = 0.0
u = 0.5
v = 0.5
!setting the p' value zero at cell 1
pprime(1,1) = 0.0
!boundary conditions falling at the cell centers of u and v control volumes
do j = 1, jmax
u(1,j) = U_W
u(imax+1,j) = U_E
end do
do i = 1, imax
v(i,1) = V_S
v(i,jmax+1) = V_N
end do
!u BCs which fall on 'v' cell centers and v BCs which fall on 'u' cell centers
U_S = 0.0
U_N = 1.0
V_W = 0.0
V_E = 0.0
10 FORMAT('Iteration',T12, 'u-velocity @ (',I3,',',I3,')')
20 FORMAT(I8,T12,F12.6)
30 FORMAT(/'Iteration',T15,'c-residual',T30,'u-residual',T45,'v-residual'/)
40 FORMAT(I5,T15,e12.6,T30,e12.6,T45,e12.6)
50 FORMAT(/'Iteration',T10,'aP',T20,'aE',T30,'aW',T40,'aN',T50,'aS'/)
60 FORMAT(I5,T10,f10.6,T20,f10.6,T30,f10.6,T40,f10.6,T50,f10.6)
70 FORMAT(/I6,T7,I10,T20,'u-error:'e15.6,T55,'u =',T58,f12.6/)
75 FORMAT(/I6,T7,I10,T20,'v-error:'e15.6,T55,'u =',T58,f12.6/)
80 FORMAT(/I6,T7,I10,T20,'p''-error:'e15.6,T55,'u =',T58,f12.6/)
90 FORMAT(1000(f15.6,','))
open(4, FILE = 'iter vs u_vel.txt')
write(4,10),imax/2,jmax/2
open(5, FILE = 'aP_etc_u-vel.xlsx')
do while(iter<1000)
iter = iter+1
!coefficients for u momentum equation over u control volumes
do j = 1, jmax
do i = 2, imax
!calculating convective fluxes at the faces of the 'u' control volumes
Fe = rho*del_y*(u(i+1,j) + u(i,j))/2
Fw = rho*del_y*(u(i,j) + u(i-1,j))/2
Fs = rho*del_x*(v(i-1,j) + v(i,j))/2
Fn = rho*del_x*(v(i-1,j+1) + v(i,j+1))/2
!adding the pressure gradient terms to b
if (j == 1) then
u_aE(i,j) = De + max(-Fe,0.0)
u_aW(i,j) = Dw + max(Fw,0.0)
u_aN(i,j) = Dn + max(-Fn,0.0)
u_aS(i,j) = 0
u_aP(i,j) = u_aE(i,j) + u_aW(i,j) + u_aN(i,j) + Fe - Fw + Fn + 2*Ds
u_b(i,j) = (2*Ds + Fs)*U_S + (p(i-1,j) - p(i,j))*del_y
else if (j == jmax) then
u_aE(i,j) = De + max(-Fe,0.0)
u_aW(i,j) = Dw + max(Fw,0.0)
u_aN(i,j) = 0
u_aS(i,j) = Ds + max(Fs,0.0)
u_aP(i,j) = u_aE(i,j) + u_aW(i,j) + u_aS(i,j) + Fe - Fw - Fs + 2*Dn
u_b(i,j) = (2*Dn - Fn)*U_N + (p(i-1,j) - p(i,j))*del_y
else
u_aE(i,j) = De + max(-Fe,0.0)
u_aW(i,j) = Dw + max(Fw,0.0)
u_aN(i,j) = Dn + max(-Fn,0.0)
u_aS(i,j) = Ds + max(Fs,0.0)
u_aP(i,j) = u_aE(i,j) + u_aW(i,j) + u_aN(i,j) + u_aS(i,j) + Fe - Fw + Fn - Fs
u_b(i,j) = (p(i-1,j) - p(i,j))*del_y
end if
end do
end do
!coefficients for v* momentum equation over staggered control volumes for v
do j = 2, jmax
do i = 1, imax
!calculating convective fluxes at the faces of the 'v' control volumes
Fe = rho*del_y*(u(i+1,j) + u(i+1,j-1))/2
Fw = rho*del_y*(u(i,j) + u(i,j-1))/2
Fs = rho*del_x*(v(i,j) + v(i,j-1))/2
Fn = rho*del_x*(v(i,j) + v(i,j+1))/2
!adding the pressure gradient terms to b
if (i == 1) then
v_aE(i,j) = De + max(-Fe,0.0)
v_aW(i,j) = 0
v_aN(i,j) = Dn + max(-Fn,0.0)
v_aS(i,j) = Ds + max(Fs,0.0)
v_aP(i,j) = v_aE(i,j) + v_aN(i,j) + v_aS(i,j) + Fe + Fn - Fs + 2*Dw
v_b(i,j) = (2*Dw + Fw)*V_W + (p(i,j-1) - p(i,j))*del_x
else if (i == imax) then
v_aE(i,j) = 0
v_aW(i,j) = Dw + max(Fw,0.0)
v_aN(i,j) = Dn + max(-Fn,0.0)
v_aS(i,j) = Ds + max(Fs,0.0)
v_aP(i,j) = v_aW(i,j) + v_aN(i,j) + v_aS(i,j) - Fw + Fn - Fs + 2*De
v_b(i,j) = (2*De - Fe)*V_E + (p(i,j-1) - p(i,j))*del_x
else
v_aE(i,j) = De + max(-Fe,0.0)
v_aW(i,j) = Dw + max(Fw,0.0)
v_aN(i,j) = Dn + max(-Fn,0.0)
v_aS(i,j) = Ds + max(Fs,0.0)
v_aP(i,j) = v_aE(i,j) + v_aW(i,j) + v_aN(i,j) + v_aS(i,j) + Fe - Fw + Fn - Fs
v_b(i,j) = (p(i,j-1) - p(i,j))*del_x
end if
end do
end do
!solving the u* momentum equation
it = 0
do
it = it + 1
do j = 1, jmax
do i = 2, imax
u_temp(i,j) = u(i,j)
end do
end do
do j = 1, jmax
do i = 2, imax
u(i,j) = (u_aE(i,j)*u(i+1,j) + u_aW(i,j)*u(i-1,j) + u_aN(i,j)*u(i,j+1) + u_aS(i,j)*u(i,j-1) + &
& u_b(i,j) + (1 - alpha_u)*u_aP(i,j)*u(i,j)/alpha_u)*alpha_u/u_aP(i,j)
end do
end do
do j = 1, jmax
do i = 2, imax
u_error(i,j) = abs(u(i,j) - u_temp(i,j))
end do
end do
r = maxval(u_error)
write(*,70),iter,it,r,u((imax/2),(jmax/2))
if(r < 1e-9) exit
end do
!solving the v* momentum equation
it = 0
do
it = it + 1
do j = 2, jmax
do i = i, imax
v_temp(i,j) = v(i,j)
end do
end do
do j = 2, jmax
do i = 1, imax
v(i,j) = (v_aE(i,j)*v(i+1,j) + v_aW(i,j)*v(i-1,j) + v_aN(i,j)*v(i,j+1) + v_aS(i,j)*v(i,j-1) + &
& v_b(i,j) + (1 - alpha_v)*v_aP(i,j)*v(i,j)/alpha_v)*alpha_v/v_aP(i,j)
end do
end do
do j = 2, jmax
do i = 1, imax
v_error(i,j) = abs(v(i,j) - v_temp(i,j))
end do
end do
r = maxval(v_error)
write(*,75),iter,it,r,u((imax/2),(jmax/2))
if(r < 1e-9) exit
end do
!calculating u residual
temp = 0
do j = 1, jmax
do i = 2, imax
temp = temp + abs(u_aP(i,j)*u(i,j) - (u_aE(i,j)*u(i+1,j) + u_aW(i,j)*u(i-1,j) + &
& u_aN(i,j)*u(i,j+1) + u_aS(i,j)*u(i,j-1) + u_b(i,j)))
end do
end do
u_residual = temp/((imax-1)*jmax)
!calculating v residual
temp = 0
do j = 2, jmax
do i = 1, imax
temp = temp + abs(v_aP(i,j)*v(i,j) - (v_aE(i,j)*v(i+1,j) + v_aW(i,j)*v(i-1,j) + &
& v_aN(i,j)*v(i,j+1) + v_aS(i,j)*v(i,j-1) + v_b(i,j)))
end do
end do
v_residual = temp/(imax*(jmax-1))
!coefficients for pressure correction equation
do j = 1, jmax
do i = 1, imax
!de = d_u(i+1,j)
!dw = d_u(i,j)
!dn = d_v(i,j+1)
!ds = d_v(i,j)
if(i == imax) then
d_u(i+1,j) = 0
else
d_u(i+1,j) = del_y/u_aP(i+1,j)
end if
if (i == 1) then
d_u(i,j) = 0
else
d_u(i,j) = del_y/u_aP(i,j)
end if
if(j == jmax) then
d_v(i,j+1) = 0
else
d_v(i,j+1) = del_x/v_aP(i,j+1)
end if
if(j == 1) then
d_v(i,j) = 0
else
d_v(i,j) = del_x/v_aP(i,j)
end if
aE(i,j) = rho*d_u(i+1,j)*del_y
aW(i,j) = rho*d_u(i,j)*del_y
aN(i,j) = rho*d_v(i,j+1)*del_x
aS(i,j) = rho*d_v(i,j)*del_x
aP(i,j) = aE(i,j) + aW(i,j) + aN(i,j) + aS(i,j)
b(i,j) = rho*(-u(i+1,j)*del_y + u(i,j)*del_y - v(i,j+1)*del_x + v(i,j)*del_x)
end do
end do
!storing u_aP, u_aE, u_aW, u_aN, u_aS for a u cell in the center
if(mod(iter,10) == 1)write(5,50)
write(5,60),iter,u_aP(imax/2,jmax/2), u_aE(imax/2,jmax/2), u_aW(imax/2,jmax/2), &
& u_aN(imax/2,jmax/2), u_aS(imax/2,jmax/2)
!calculating continuity residual
temp = 0
do j = 1, jmax
do i = 1, imax
temp = temp + abs(b(i,j))
end do
end do
c_residual = temp/(imax*jmax)
if(mod(iter,10) == 1) write(*,30)
write(*,40),iter, c_residual, u_residual, v_residual
if((u_residual + v_residual + c_residual) < tolerance) then
write(*,*),'the solution is converged'
exit
end if
!solving pressure correction equation
it = 0
do
it = it + 1
do j = 1, jmax
do i = 1, imax
p_temp(i,j) = pprime(i,j)
end do
end do
do j = 1, jmax
do i = 1, imax
if(i == 1 .AND. j == 1) cycle
pprime(i,j) = (aE(i,j)*pprime(i+1,j) + aW(i,j)*pprime(i-1,j) + aN(i,j)*pprime(i,j+1) + &
& aS(i,j)*pprime(i,j-1) + b(i,j))/aP(i,j)
end do
end do
temp = 0
do j = 1, jmax
do i = 1, imax
p_error(i,j) = abs(pprime(i,j) - p_temp(i,j))
end do
end do
r = maxval(p_error)
write(*,80),iter,it,r,u((imax/2),(jmax/2))
if(r < 1e-8) exit
end do
!correcting pressure
do j = 1, jmax
do i = 1, imax
p(i,j) = p(i,j) + alpha_p*pprime(i,j)
end do
end do
!correcting u velocities
do j = 1, jmax
do i = 2, imax
u(i,j) = u(i,j) + d_u(i,j)*(pprime(i-1,j) - pprime(i,j))
end do
end do
!correcting v velocities
do j = 2, jmax
do i = 1, imax
v(i,j) = v(i,j) + d_v(i,j)*(pprime(i,j-1) - pprime(i,j))
end do
end do
if(mod(iter,100) == 1)write(4,20),iter,u(imax/2,jmax/2)
b = 0
temp = 0
do j = 1, jmax
do i = 1, imax
b(i,j) = -u(i+1,j) + u(i,j) - v(i,j+1) + v(i,j)
temp = temp + b(i,j)
end do
end do
temp = temp/(imax*jmax)
write(*,*)'Mass imbalance after velocity correction is',temp
end do
close(4)
close(5)
write(*,*),'the u velocities at cell centers after',iter,'iterations are:'
open(1, FILE = 'u.vel.txt')
do j = 1, jmax
write(*,'(1000f10.6)'),(u(i,j), i = 1, imax+1)
write(1,'(1000f10.6)'),(u(i,j), i = 1, imax+1)
end do
close(1)
write(*,*),'the v velocities at cell centers after',iter,'iterations are:'
open(2, FILE = 'v.vel.xlsx')
do j = 1, jmax+1
write(*,90),(v(i,j), i = 1, imax)
write(2,90),(v(i,j), i = 1, imax)
end do
close(2)
write(*,*),'the pressure at cell centers after',iter,'iterations are:'
open(3, FILE = 'presure.txt')
do j = 1, jmax
write(*,'(1000f10.6)'),(p(i,j), i = 1, imax)
write(3,'(1000f10.6)'),(p(i,j), i = 1, imax)
end do
close(3)
deallocate(u, u_temp, u_aP, u_aE, u_aW, u_aN, u_aS, u_b, d_u)
deallocate(v, v_temp, v_aP, v_aE, v_aW, v_aN, v_aS, v_b, d_v)
deallocate(p, pprime, aP, aE, aW, aN, aS, b, psi, p_temp)
end program
|
