CFD Online Discussion Forums - Boundary Conditions for FORTRAN Lid Driven Cavity Error

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

I have created a SIMPLE Solver for the the Lid Driven Cavity Problem in MATLAB. It is known working and Validated against the 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 have attached a .zip file with the original Fortran Code, my attempted Lid Driven Cavity Conversion and parameter file.

for U Momentum function

Code:

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

Code:

! 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