[BaklavaBaby29]

I am studying Numerical Methods for incompressible flows. Part of the tasks is to model the lid-driven cavity problem in 2D using the SIMPLE method.

I have been provided with Fortran code that is working and solved. However, is for a channel flow problem. I was informed that by modifying the boundary conditions I can convert it from channel flow to cavity flow. (by manipulating which walls as velocity etc).

I have created a SIMPLE Solver for the Lid Driven Cavity Problem in MATLAB. It is known to work and Validated against Ghia et al.(1982). However when I modify the Fortran code to what I believe to be the correct boundary conditions my results.dat file contains NaN for all values of U and V for all X and Y locations

I believe by modifying the boundary conditions in UMomentum,Vmomentum, Pressure Correction and Scalar I can convert it to a lid driven cavity. as south is a known wall in the channel flow I based my edits on this fact when converting the east and west inlets to walls as such I used these as the correct boundary conditions.

for U Momentum function

! set boundary conditions for coefficients in the equations
! south
    aw(:,1)=0.
    ae(:,1)=0.
    as(:,1)=0.
    an(:,1)=0.
    ap(:,1)=1.
    su(:,1)=0.
! west
    aw(1,:)=0.
    ae(1,:)=0.
    as(1,:)=0.
    an(1,:)=0.
    ap(1,:)=1.
    su(1,:)=0.
! north
    aw(:,iNyUNodes)=0.
    ae(:,iNyUNodes)=0.
    as(:,iNyUNodes)=0.
    an(:,iNyUNodes)=0.
    ap(:,iNyUNodes)=1.
    su(:,iNyUNodes)=1.
! east
    aw(iNxUNodes,:)=0.
    ae(iNxUNodes,:)=0.
    as(iNxUNodes,:)=0.
    an(iNxUNodes,:)=0.
    ap(iNxUNodes,:)=1.
    su(iNxUNodes,:)=0

.

as there is no momentum in the Y direction I believe these to be the correct Boundary Conditions for VMomentum Function.

! south
    aw(:,1)=0.
    ae(:,1)=0.
    as(:,1)=0.
    an(:,1)=0.
    ap(:,1)=1.
    su(:,1)=0.
! west
    aw(1,:)=0.
    ae(1,:)=0.
    as(1,:)=0.
    an(1,:)=0.
    ap(1,:)=1.
    su(1,:)=0.
! north
    aw(:,iNyVNodes)=0.
    ae(:,iNyVNodes)=0.
    as(:,iNyVNodes)=0.
    an(:,iNyVNodes)=0.
    ap(:,iNyVNodes)=1.
    su(:,iNyVNodes)=0.
! east
    aw(iNxVNodes,:)=0.
    ae(iNxVNodes,:)=0.
    as(iNxVNodes,:)=0.
    an(iNxVNodes,:)=0.
    ap(iNxVNodes,:)=1.
    su(iNxVNodes,:)=0

.


I am unsure how to manipulate the Pressure Corrections and Scalar Function's Boundary Conditions.

and the output of my code is either incorrect values or NaN with different attempts of the boundary conditions

Apologies for this being long, I wanted to be as thorough as possible as I cannot understand what is incorrect.

(it is known that the Fortran solver is validated as accurate, and my MATLAB solver) however these results do not matchup after my boundary edits.

Thanks in advance
EDIT: code is attached

[naffrancois]

If not already done, the first step is to compile your code with debugging flags. Deactivate all optimizations and add checks for array indices, overflows etc. I give you some of them, possibly redundant, for gfortran:

-O0
-fbounds-check
-fbacktrace
-fcheck=all
-ffpe-trap=invalid,zero,overflow,underflow
-Wuninitialized

fbacktrace should tell you where things start to go wrong.

edit: sorry I read again your post, this is a matlab code. I can not help much here, but there should be some specialized tools as well to track down errors in the code

[BaklavaBaby29]

Originally Posted by naffrancois

If not already done, the first step is to compile your code with debugging flags. Deactivate all optimizations and add checks for array indices, overflows etc. I give you some of them, possibly redundant, for gfortran:

-O0
-fbounds-check
-fbacktrace
-fcheck=all
-ffpe-trap=invalid,zero,overflow,underflow
-Wuninitialized

fbacktrace should tell you where things start to go wrong.

edit: sorry I read again your post, this is a matlab code. I can not help much here, but there should be some specialized tools as well to track down errors in the code

Thanks for the response! It definitely is fotran, I can provide the .f90 when I get home in 10 mins

[arjun]

Originally Posted by BaklavaBaby29

Thanks for the response! It definitely is fotran, I can provide the .f90 when I get home in 10 mins

use very small under-relaxation factors first. This shall get you results. Then increase urfs slowly and see when it diverges.


It is normal in SIMPLE method to have divergence if the code is a simple implementation.

[BaklavaBaby29]

Originally Posted by naffrancois

If not already done, the first step is to compile your code with debugging flags. Deactivate all optimizations and add checks for array indices, overflows etc. I give you some of them, possibly redundant, for gfortran:

-O0
-fbounds-check
-fbacktrace
-fcheck=all
-ffpe-trap=invalid,zero,overflow,underflow
-Wuninitialized

fbacktrace should tell you where things start to go wrong.

edit: sorry I read again your post, this is a matlab code. I can not help much here, but there should be some specialized tools as well to track down errors in the code

it is in the edit

[BaklavaBaby29]

Originally Posted by arjun

use very small under-relaxation factors first. This shall get you results. Then increase urfs slowly and see when it diverges.


It is normal in SIMPLE method to have divergence if the code is a simple implementation.

Code is Attached

[naffrancois]

It is good to post your code, but it is better if you try yourself to apply the advice people give you . Tell us what comes out of your trials and we may guide you.

[arjun]

Originally Posted by BaklavaBaby29

Code is Attached

I am not going to look into it for the reasons that last 18 years SIMPLE is what i have mostly done including this weekend. I am tired of looking into SIMPLE code.


Based on my experience code can be very unstable when coded without safe guards. You may not be doing anything wrong.

[FMDenaro]

My suggestion is to debug carefully your code when you add the two walls but start a very low Reynolds number to check that you are able to get convergence towards a symmetric flow field.

Only after you can increase the Re number and compare to the results in literature. But that will require you are able to increas the grid resolution.