[kar1209]

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:

1. Specify the boundary conditions, specify the guess values for pressure field and velocity field(for solving the non linearities)
2. fix the pressure at a cell to make the system of equations determinant
3. calculate the diffusion coefficients at the faces of the staggered control volumes
4. 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)
5. assign values to coefficients aP, aE, etc for the x momentum equation
6. repeat steps 4 and 5 for y momentum equation
7. 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
8. after leaving the gauss seidel loop, calculate residuals of a and y momentum equations and continuity equation using the obtained u and v velocities
9. calculate coefficients for pressure correction equation
10. solve the system of pressure correction equations using Gauss-Seidel with a similar converge criteria
11. update u, v and p and check for mass imbalance of all cells
12. 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.

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

[ztdep]

hope my code will give you some hints