CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > General Forums > Main CFD Forum

Pressure-projection scheme for a lid-driven cavity

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   January 6, 2016, 15:58
Default Pressure-projection scheme for a lid-driven cavity
  #1
New Member
 
SD
Join Date: Nov 2015
Location: United Kingdom
Posts: 5
Rep Power: 10
kopaka is on a distinguished road
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 ). 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 !
Attached Images
File Type: jpg export.jpg (49.3 KB, 37 views)
Attached Files
File Type: f main.f (4.7 KB, 22 views)
kopaka is offline   Reply With Quote

Old   January 6, 2016, 16:54
Default
  #2
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,760
Rep Power: 71
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
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
FMDenaro is offline   Reply With Quote

Old   January 6, 2016, 19:57
Default
  #3
New Member
 
SD
Join Date: Nov 2015
Location: United Kingdom
Posts: 5
Rep Power: 10
kopaka is on a distinguished road
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 ?
kopaka is offline   Reply With Quote

Old   January 7, 2016, 04:35
Default
  #4
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,760
Rep Power: 71
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
I suggest a deeper study of the method, the original paper in JCP is a good starting point
http://www.sciencedirect.com/science...21999185901482
FMDenaro is offline   Reply With Quote

Old   January 7, 2016, 04:58
Default
  #5
Senior Member
 
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 1,273
Rep Power: 34
arjun will become famous soon enougharjun will become famous soon enough
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.
arjun is offline   Reply With Quote

Old   January 7, 2016, 05:10
Default
  #6
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,760
Rep Power: 71
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by kopaka View Post
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 ?


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.
FMDenaro is offline   Reply With Quote

Old   January 8, 2016, 11:12
Default
  #7
New Member
 
SD
Join Date: Nov 2015
Location: United Kingdom
Posts: 5
Rep Power: 10
kopaka is on a distinguished road
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 :



For the BC cells :



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
kopaka is offline   Reply With Quote

Old   January 8, 2016, 11:46
Default
  #8
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,760
Rep Power: 71
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
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
FMDenaro is offline   Reply With Quote

Old   January 9, 2016, 08:52
Default
  #9
Senior Member
 
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 1,273
Rep Power: 34
arjun will become famous soon enougharjun will become famous soon enough
for pure poisson problem your iterative solver shall not work.
arjun is offline   Reply With Quote

Old   November 18, 2017, 02:45
Default SIMPLE algorithm for lid-driven cavity not converging
  #10
New Member
 
Karthikeyan
Join Date: Nov 2016
Posts: 14
Rep Power: 9
kar1209 is on a distinguished road
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
kar1209 is offline   Reply With Quote

Old   November 18, 2017, 04:45
Default
  #11
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,760
Rep Power: 71
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by kar1209 View Post
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

Open a new post and describe better your convergence problem
FMDenaro is offline   Reply With Quote

Old   November 18, 2017, 08:40
Default
  #12
New Member
 
Karthikeyan
Join Date: Nov 2016
Posts: 14
Rep Power: 9
kar1209 is on a distinguished road
Hi FMDenaro,

I have posted my question in a new thread.

SIMPLE algorithm for lid-driven cavity not converging

Thanks
kar1209 is offline   Reply With Quote

Old   November 18, 2017, 12:10
Default
  #13
Senior Member
 
ztdep's Avatar
 
p ding
Join Date: Mar 2009
Posts: 427
Rep Power: 19
ztdep is on a distinguished road
Send a message via Yahoo to ztdep Send a message via Skype™ to ztdep
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.
ztdep is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Lid driven Cavity suman91 SU2 3 May 24, 2022 06:23
Pressure variation in 2D lid driven cavity Agad15 Main CFD Forum 0 October 9, 2014 05:41
lid-driven cavity in matlab using BiCGStab Don456 Main CFD Forum 1 January 19, 2012 16:00
Lid Driven Cavity velocity profiles new_at_this Main CFD Forum 0 December 22, 2011 17:04
[OpenFOAM] Paraview - Lid Driven Cavity kieranhood ParaView 0 February 13, 2010 17:28


All times are GMT -4. The time now is 08:13.