CFD Online Discussion Forums - Pressure-projection scheme for a lid-driven cavity

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- Main CFD Forum (https://www.cfd-online.com/Forums/main/)

- - Pressure-projection scheme for a lid-driven cavity (https://www.cfd-online.com/Forums/main/164977-pressure-projection-scheme-lid-driven-cavity.html)

Pressure-projection scheme for a lid-driven cavity

Hi everyone

I was mixed by one problem for some time and wonder if there is someone who can cast some light on my head.

I'm trying to implement a pressure-projection code for the case of a lid-driven cavity. But I think my code is wrong and I don't know why. The results I get are false. (you can see the attached picture) I don't know if it is a boundary condition problem, if something is missing in my method, if the parameters are wrong or if it's simply a programming error (I have checked many times :D). I'm almost sure it is a simple issue and if someone who is a bit experienced with pressure projection method see this, it could be easy to point out the issue.

I'm a CFD begginer so I apologize by advance if it is a trivial issue, but I have been stucked for many hours and I really want this program to work :)

You can find my program attached in this topic, if someone find an error or some missing calculus, I will be really grateful :)

Thank you in advance !
See you.

My code works as follows :

First step : Intermediate velocity

Momentum :
$\frac{\vec{u}^*-\vec{u}^n}{\Delta t} + [(\vec{u}.\nabla)\vec{u}]^n = [\nu \nabla^2 \vec{u}]^n$
Intermediate velocity :
$\vec{u}^* = \vec{u}^n + \Delta t [-(\vec{u}.\nabla)\vec{u}+\nu\nabla^2\vec{u}]^n$

Second step : Fractional step, solution of the poisson equation

1. Compute the divergence of the velocity field (to have RHS)

$\frac{\vec{u}^{n+1}-\vec{u}^*}{\Delta t} = -\frac{1}{\rho}\nabla p$
and $\nabla \vec{u}^{n+1}\approx 0$ so

$\nabla.\vec{p}^{n+1}=\frac{\rho}{\Delta t}(\nabla . \vec{u}^*)=RHS$

Before solving this equation, we need to compute the Right Hand Side.

$RHS_{i,j}=\frac{\rho}{\Delta t}(\frac{u^*_{i,j}-u^*_{i-1,j}}{\Delta x}+\frac{v^*_{i,j}-v^*_{i,j-1}}{\Delta y})$

2. Solve the poisson equation

After computing the RHS, we need to solve this poisson equation iteratively. I use a S.O.R method :

$\frac{\omega p^{n+1}_{i-1,j}-2p^{n+1}_{i,j}+2(1-\omega)p^n_{i,j}+\omega p^{n+1}_{i+1,j}}{\Delta x^2} + \frac{\omega p^{n+1}_{i,j-1}-2p^{n+1}_{i,j}+2(1-\omega)p^n_{i,j}+\omega p^{n+1}_{i,j+1}}{\Delta y^2}$
$= \omega RHS_{i,j}$

So

$\underbrace{\frac{\omega}{\Delta x^2}}_{a_w}p^{n+1}_{i-1,j} + \underbrace{\frac{\omega}{\Delta x^2}}_{a_e}p^{n+1}_{i+1,j} +\underbrace{\frac{\omega}{\Delta y^2}}_{a_n}p^{n+1}_{i,j+1} +\underbrace{\frac{\omega}{\Delta y^2}}_{a_s}p^{n+1}_{i-1,j} \underbrace{-2[\frac{1}{\Delta x^2}+\frac{1}{\Delta y^2}]}_{a_p} p^{n+1}_{i,j}$
$+ 2(1-\omega)[\frac{1}{\Delta x^2}+\frac{1}{\Delta y^2}]p^{n}_{i,j} = \omega RHS$

So we compute the pressure this way :

$p^{n+1}_{i,j} = (1-\omega)p^n_{i,j}+\frac{1}{a_p}[a_s p^n_{i,j-1}+a_n p^n_{i,j+1}+a_w p^n_{i-1,j}+a_e p^n_{i+1,j}-\omega RHS_{i,j}]$

Third step : Velocity field update

$\vec{u}^{n+1}=\vec{u}^*-\frac{\Delta t}{\rho}\nabla p^{n+1}$

The pressure gradient is calculated as a forward finite difference.

And loop to the next time step.

You can see my code attached

Thank you in advance !

what about the BC.s for the pressure equation?
are you using staggered or colocated grid arrangement?
have you checked the divergence-free constraint in each cell after the correction?

Then, provide a velocity plot

Hi!
Thank you for your answer.

I think it is a collocated grid I'm using in the program.

The BC for pressure are Neumann type, but actually I think the problem might be here.
Here is the code I use for computing the pressure, I didn't know how to implement the BC.

Code:

        

        do j=2,jmax-1

            do i=2,imax-1

                as=w/(Dy*Dy)

                aw=w/(Dx*Dx)

                ae=w/(Dx*Dx)

                an=w/(Dy*Dy)

                ap=2*((1/(Dx*Dx))+(1/(Dy*Dy)))

                p(i,j)=(1-w)*p(i,j)+(1/ap)*(as*p(i,j-1)+aw*p(i-1,j)+ae*p(i+1,j)+&

                    an*p(i,j+1)-w*RHS(i,j))    

            end do

        end do

        

        !Neumann conditions

        do i=1,imax

            do j=1,jmax

                !Bottom

                p(i,1)=(1-w)*p(i,j)+(1/ap)*(as*p(i,j)+aw*p(i-1,j)+ae*p(i+1,j)+&

                        an*p(i,j+1)-w*RHS(i,j)) 

                !Top

                p(i,jmax)=(1-w)*p(i,j)+(1/ap)*(as*p(i,j-1)+aw*p(i-1,j)+ae*p(i+1,j)+&

                        an*p(i,j)-w*RHS(i,j)) 

                !Left

                p(1,j)=(1-w)*p(i,j)+(1/ap)*(as*p(i,j-1)+aw*p(i,j)+ae*p(i+1,j)+&

                        an*p(i,j+1)-w*RHS(i,j)) 

                !Right

                p(imax,j)=(1-w)*p(i,j)+(1/ap)*(as*p(i,j-1)+aw*p(i-1,j)+ae*p(i,j)+&

                        an*p(i,j+1)-w*RHS(i,j)) 

            end do

        end do

I'm not checking the divergence free constraint, what do you mean by that ?

I suggest a deeper study of the method, the original paper in JCP is a good starting point
http://www.sciencedirect.com/science...21999185901482

The main problem that I see in your code is that you dont show us how you calculate divergence for velocity field.

For collocated method to work, the velocity at face that you use for divergence has two contribution one from the interpolated velocity and one from rhie and chow dissipation term (which is function of pressure and AP).

Also if you did so based on Ap you used then you still have problem because from it you can not derive the pressure equation you use.

Hint: the pressure equation you use you will get when your Ap has only time contribution and nothing else.
But i guess you use different Ap for interpolation.

Last make sure the solver is handling all neumann condition too.

Quote:

Originally Posted by kopaka (Post 579833)

Code:

        

        do j=2,jmax-1

            do i=2,imax-1

                as=w/(Dy*Dy)

                aw=w/(Dx*Dx)

                ae=w/(Dx*Dx)

                an=w/(Dy*Dy)

                ap=2*((1/(Dx*Dx))+(1/(Dy*Dy)))

                p(i,j)=(1-w)*p(i,j)+(1/ap)*(as*p(i,j-1)+aw*p(i-1,j)+ae*p(i+1,j)+&

                    an*p(i,j+1)-w*RHS(i,j))    

            end do

        end do

        

        !Neumann conditions

        do i=1,imax

            do j=1,jmax

                !Bottom

                p(i,1)=(1-w)*p(i,j)+(1/ap)*(as*p(i,j)+aw*p(i-1,j)+ae*p(i+1,j)+&

                        an*p(i,j+1)-w*RHS(i,j)) 

                !Top

                p(i,jmax)=(1-w)*p(i,j)+(1/ap)*(as*p(i,j-1)+aw*p(i-1,j)+ae*p(i+1,j)+&

                        an*p(i,j)-w*RHS(i,j)) 

                !Left

                p(1,j)=(1-w)*p(i,j)+(1/ap)*(as*p(i,j-1)+aw*p(i,j)+ae*p(i+1,j)+&

                        an*p(i,j+1)-w*RHS(i,j)) 

                !Right

                p(imax,j)=(1-w)*p(i,j)+(1/ap)*(as*p(i,j-1)+aw*p(i-1,j)+ae*p(i,j)+&

                        an*p(i,j+1)-w*RHS(i,j)) 

            end do

        end do

I'm not checking the divergence free constraint, what do you mean by that ?

the way in which the code is written is not correct... the Neumann condition must change both the matrix entries and the known term (RHS of the Poisson equation). The first "do" cycles must be rewritten differently for nodes adjacent to the boundaries and the second "do" cycles can be eliminated.

Hi!
Thank you everyone :)

These Neumann BC are really annoying me, I think my approach is wrong.
For the problem to be easier, I use the Jacobi method instead of a SOR.

I use the formulas which I found here : http://demonstrations.wolfram.com/So...ferentMethods/ (I have a uniform grid spacing)

For the no BC cells :

http://demonstrations.wolfram.com/So...CellNone_1.gif

For the BC cells :
http://demonstrations.wolfram.com/So...CellNone_2.gif

My code is the following :

p represent p(k+1) and pold represent p(k)

Code:

    !Solve pressure poisson iteratively

    rP=11

    iterationPP=1

    p(:,:)=0

    

    do while ((iterationPP < 200).and.(rP>0.0000001))

        iterationPP=iterationPP+1

        pold=p

        !Neumann Boudary

        do i=2,imax-1

                !Left

                p(i,1)=0.25*(2*pold(i,2)+pold(i-1,1)+pold(i+1,1)-Dx*Dy*RHS(i,1))

                !Right

                p(i,imax)=0.25*(2*pold(i,imax-1)+pold(i-1,imax)+pold(i+1,imax)-&

                        Dx*Dy*RHS(i,imax))        

        end do

        do j=2,jmax-1

                !Top

                p(1,j)=0.25*(2*pold(2,j)+pold(1,j-1)+pold(1,j+1)-Dy*Dy*RHS(1,j))

                !Bottom

                p(imax,j)=0.25*(2*pold(imax-1,j)+pold(imax,j-1)+pold(imax,j+1)-&

                        Dy*Dy*RHS(imax,j))    

        end do

        

        !Jacobi iteration for the cells in the center

        do j=2,jmax-1

            do i=2,imax-1

                p(i,j)=0.25*(pold(i-1,j)+pold(i+1,j)+pold(i,j-1)+pold(i,j+1)-&

                       Dx*Dy*RHS(i,j))    

            end do

        end do

        

        !Corners

        p(1,1)=0.5*(p(1,2)+p(2,1)) !Bot-Left

        p(1,jmax)=0.5*(p(1,jmax-1)+p(2,jmax)) !Top-Left

        p(imax,1)=0.5*(p(imax-1,1)+p(imax,2)) !Bot-Right

        p(imax,jmax)=0.5*(p(imax,jmax-1)+p(imax-1,jmax)) !Top-Right

        

        

        rP=abs(maxval(pold(:,:)-p(:,:)))



        print *, "rP=", rP

    end do

I think my approach to add the Neumann BC is wrong but I really don't understand why. I maybe someone can help me :)

Thanks

Simon

Simple 1d example for the Poisson equation

((dp/dx)i+1/2 - (dp/dx)i-1/2)/h = ((u*)i+1/2 - (u*)i-1/2)/h

Neumann BC at i-1/2

dp/dx = u* - un+1

the equation becomes

((dp/dx)i+1/2 )/h = ((u*)i+1/2 - (un+1)i-1/2)/h

for pure poisson problem your iterative solver shall not work.

SIMPLE algorithm for lid-driven cavity not converging

Hi,

I have written a code in FORTRAN for solving incompressible steady flow in a lid driven square cavity. I have used a staggered grid formulation; I have used upwind scheme for calculating the convective fluxes at the faces (Fe, Fw, Fn, Fs in the code).

I tried playing around with the under-relaxation factors but could not get a converged solution for v velocity.

I am attaching my code here. Please provide 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(:,:)

  double precision, allocatable :: v(:,:), v_temp(:,:), v_aP(:,:), v_aE(:,:), v_aW(:,:), v_aN(:,:), v_aS(:,:), v_b(:,:)

  double precision, allocatable :: p(:,:), pprime(:,:), aP(:,:), aE(:,:), aW(:,:), aN(:,:), aS(:,:), b(:,:), psi(:,:), &

    & d_u(:,:), d_v(:,:), p_temp(:,:)

  

  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.01

  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.52, alpha_v = 0.5, alpha_p = 0.01

  double precision :: u_residual, v_residual, c_residual, tolerance = 1e-6, 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)) 

    

  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))

    

  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))

  

  !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(/I4,T5,I10,T20,'u-error:'e15.6,T55,'u =',T58,f12.6/)

  75 FORMAT(/I4,T5,I10,T20,'v-error:'e15.6,T55,'u =',T58,f12.6/)

  80 FORMAT(/I4,T5,I10,T20,'p''-error:'e15.6,T55,'u =',T58,f12.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

    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_temp(i,j)/alpha_u)*alpha_u/u_aP(i,j)

        end do

      end do 

      temp = 0

      do j = 1, jmax

        do i = 2, imax 

          temp = temp + abs(u(i,j) - u_temp(i,j))

        end do

      end do

      r = temp/((imax-1)*jmax)

      write(*,70),iter,it,r,u((imax/2),(jmax/2))

      if(r < 1e-8) 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_temp(i,j)/alpha_v)*alpha_v/v_aP(i,j)

        end do

      end do 

      temp = 0

      do j = 2, jmax

        do i = 1, imax 

          temp = temp + abs(v(i,j) - v_temp(i,j))

        end do

      end do

      r = temp/((imax)*(jmax-1))

      write(*,75),iter,it,r,u((imax/2),(jmax/2))

      if(r < 1e-10) 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) = d_u(i+1,j)*del_y

        aW(i,j) = d_u(i,j)*del_y

        aN(i,j) = d_v(i,j+1)*del_x

        aS(i,j) = d_v(i,j)*del_x

        aP(i,j) = aE(i,j) + aW(i,j) + aN(i,j) + aS(i,j)

        b(i,j) = -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 aP, aE, aW, aN, 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

          temp = temp + abs(pprime(i,j) - p_temp(i,j))

        end do

      end do

      r = temp/(imax*jmax)

      write(*,80),iter,it,r,u((imax/2),(jmax/2))

      if(r < 1e-5) 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)

  

  !calculating stream function from velocity field

  do j = 1,jmax

    do i = 1, imax

      psi(i,j) = (u(i,j) + u(i+1,j))*del_x/2 + (v(i,j) + v(i,j+1))*del_y/2

    end do

  end do

  

  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.txt')

  do j = 1, jmax+1

    write(*,'(1000f10.6)'),(v(i,j), i = 1, imax)

    write(1,'(1000f10.6)'),(v(i,j), i = 1, imax)

  end do 

  close(2)

  

  write(*,*),'the stream function at cell centers after',iter,'iterations  are:'

  open(3, FILE = 'stream.txt')

  do j = 1, jmax

    write(*,'(1000f10.6)'),(psi(i,j), i = 1, imax)

    write(1,'(1000f10.6)'),(psi(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

Quote:

Originally Posted by kar1209 (Post 671972)

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(:,:)

  double precision, allocatable :: v(:,:), v_temp(:,:), v_aP(:,:), v_aE(:,:), v_aW(:,:), v_aN(:,:), v_aS(:,:), v_b(:,:)

  double precision, allocatable :: p(:,:), pprime(:,:), aP(:,:), aE(:,:), aW(:,:), aN(:,:), aS(:,:), b(:,:), psi(:,:), &

    & d_u(:,:), d_v(:,:), p_temp(:,:)

  

  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.01

  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.52, alpha_v = 0.5, alpha_p = 0.01

  double precision :: u_residual, v_residual, c_residual, tolerance = 1e-6, 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)) 

    

  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))

    

  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))

  

  !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(/I4,T5,I10,T20,'u-error:'e15.6,T55,'u =',T58,f12.6/)

  75 FORMAT(/I4,T5,I10,T20,'v-error:'e15.6,T55,'u =',T58,f12.6/)

  80 FORMAT(/I4,T5,I10,T20,'p''-error:'e15.6,T55,'u =',T58,f12.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

    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_temp(i,j)/alpha_u)*alpha_u/u_aP(i,j)

        end do

      end do 

      temp = 0

      do j = 1, jmax

        do i = 2, imax 

          temp = temp + abs(u(i,j) - u_temp(i,j))

        end do

      end do

      r = temp/((imax-1)*jmax)

      write(*,70),iter,it,r,u((imax/2),(jmax/2))

      if(r < 1e-8) 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_temp(i,j)/alpha_v)*alpha_v/v_aP(i,j)

        end do

      end do 

      temp = 0

      do j = 2, jmax

        do i = 1, imax 

          temp = temp + abs(v(i,j) - v_temp(i,j))

        end do

      end do

      r = temp/((imax)*(jmax-1))

      write(*,75),iter,it,r,u((imax/2),(jmax/2))

      if(r < 1e-10) 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) = d_u(i+1,j)*del_y

        aW(i,j) = d_u(i,j)*del_y

        aN(i,j) = d_v(i,j+1)*del_x

        aS(i,j) = d_v(i,j)*del_x

        aP(i,j) = aE(i,j) + aW(i,j) + aN(i,j) + aS(i,j)

        b(i,j) = -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 aP, aE, aW, aN, 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

          temp = temp + abs(pprime(i,j) - p_temp(i,j))

        end do

      end do

      r = temp/(imax*jmax)

      write(*,80),iter,it,r,u((imax/2),(jmax/2))

      if(r < 1e-5) 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)

  

  !calculating stream function from velocity field

  do j = 1,jmax

    do i = 1, imax

      psi(i,j) = (u(i,j) + u(i+1,j))*del_x/2 + (v(i,j) + v(i,j+1))*del_y/2

    end do

  end do

  

  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.txt')

  do j = 1, jmax+1

    write(*,'(1000f10.6)'),(v(i,j), i = 1, imax)

    write(1,'(1000f10.6)'),(v(i,j), i = 1, imax)

  end do 

  close(2)

  

  write(*,*),'the stream function at cell centers after',iter,'iterations  are:'

  open(3, FILE = 'stream.txt')

  do j = 1, jmax

    write(*,'(1000f10.6)'),(psi(i,j), i = 1, imax)

    write(1,'(1000f10.6)'),(psi(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

Open a new post and describe better your convergence problem

Hi FMDenaro,

I have posted my question in a new thread.

https://www.cfd-online.com/Forums/ma...tml#post671991

Thanks

It seems you are using the forward differencing scheme to compute the divergence of velocity field?
We should use the CD scheme to compute it.
For the pressure equation, a simple "insulation" BC can be used.