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

2D incompressible cavity flow Matlab code NAN result

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

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   June 6, 2024, 06:48
Default 2D incompressible cavity flow Matlab code NAN result
  #1
New Member
 
Sasha Prival
Join Date: Jun 2024
Posts: 2
Rep Power: 0
SashaPrival is on a distinguished road
Hello. Newbie here. I am following this guide to create my very first simple CFD solvers for incompressible 2D fluid in a lid-driven square cavity. In the guide, the pressure is considered with zero-derivative Neumann BCs imposed on all boundaries. I have managed to write and run a Matlab code following the guide. It runs fine.

Then I wanted to change the top boundary condition for the pressure to 0 Dirichlet like in the Step 11 of famous beginner guide of Lorena Barba. I have modified the Laplacian matrix contents to include top Dirchlet BC. Assigned 0 to the corresponding elements of the source (right-hand side) vector. Forced 0 value to the corresponding elements of the solution (pressure) vector. Here is the complete Matlab code I wrote:
Code:
clear; clc;
% grid parameters
Lx=1; % x length
Ly=1; % y length
nx=41; % number of grid points along x
ny=41; % number of grid points along y

%simulation parameters
dt=0.001; % time step
t_final=5.0; % simulation run time

% Number of timesteps
Nt=t_final/dt;

% physics parameters
rho=1.0; % density 
nu=0.001; % viscosity coefficient (kinematic)

% Index extents
imin=2;
imax=imin+nx-1; % ok
jmin=2;         % ok
jmax=jmin+ny-1;

%create mesh sizes
dx=Lx/(nx);
dy=Ly/(ny);
dxi=1/dx;
dyi=1/dy;

% preallocate important arrays 
p=zeros(imax,jmax); % pressure array

%p=zeros(nx,ny); % ok
us=zeros(imax+1,jmax+1); % u - velocity without pressure
vs=zeros(imax+1,jmax+1); % v - velocity without pressure
R = zeros(nx*ny, 1);  % RHS of Pressure Poisson Equation
u=zeros(imax+1,jmax+1); % corrected u
v=zeros(imax+1,jmax+1); % corrected v
L=zeros(nx*ny,nx*ny); % Laplacian matrix

% initial values
t=0; % initial time
u_bot=0; % initial velocity for bottom wall
u_top=1; % initial velocity for top wall
v_lef=0; % initial velocity for left wall
v_rig=0; % initial velocity for right wall

% Creat Laplacian operator for solving pressure Poisson equation
for j=1:ny %number of block matrices along y
    for i=1:nx %number of block matrices along x
        q = i+(j-1)*nx;
        %if(mod(q,1*nx)==0)
        if (q > (nx-1)*ny) % if top boundary            
            L(q,q)=1; % Zero Dirichlet condition
        else
            L(q,q)=-(2*dxi^2+2*dyi^2);
            for ii=i-1:2:i+1
                if (ii>0 && ii<=nx) % Interior points
                    L(q,ii+(j-1)*nx)=dxi^2;      
                else % Neuman conditions on boundary
                    L(q,q)= L(q,q)+dxi^2;
                end
            end
            for jj=j-1:2:j+1
                if (jj>0 && jj<=ny) % Interior points
                    L(q,i+(jj-1)*nx)=dyi^2;
                else % Neuman conditions on boundary
                    L(q,q)= L(q,q)+dyi^2;
                end
            end
        end
    end
end
%disp("L="); disp(num2str(L));

% solver
disp("Starting simulation ****************************");


while t <= t_final
    % update time
    t = t + dt;
    clc; disp("t:"); disp(t);


    u_top=1;


    % u velocity predictor
    for j = jmin:jmax
        for i = imin+1:imax
            v_here = 0.25*(v(i-1,j)+v(i-1,j+1)+v(i,j)+v(i,j+1));
            us(i,j)=u(i,j)+dt* ...
                (nu*(u(i-1,j)-2*u(i,j)+u(i+1,j))*dxi^2 ...
                +nu*(u(i,j-1)-2*u(i,j)+u(i,j+1))*dyi^2 ...
                -u(i,j)*(u(i+1,j)-u(i-1,j))*0.5*dxi ...
                -v_here*(u(i,j+1)-u(i,j-1))*0.5*dyi);
        end
     end


    % v velocity predictor
    for j = jmin+1:jmax
        for i = imin:imax
            u_here = 0.25*(u(i,j-1)+u(i+1,j-1)+u(i,j)+u(i+1,j));
            vs(i,j)=v(i,j)+dt* ...
                (nu*(v(i-1,j)-2*v(i,j)+v(i+1,j))*dxi^2 ...
                +nu*(v(i,j-1)-2*v(i,j)+v(i,j+1))*dyi^2 ...
                -u_here*(v(i+1,j)-v(i-1,j))*0.5*dxi...
                -v(i,j)*(v(i,j+1)-v(i,j-1))*0.5*dyi);
        end
    end


    % compute R.H.S. (R) of the PPEquation using vs & us.
    n=0;    
    %R = zeros(nx*ny, 1);  % RHS of PPEquation
    for j=jmin:jmax
        for i=imin:imax
            n=n+1;
            R(n)=(rho/dt)*((us(i+1,j)-us(i,j))*dxi+(vs(i,j+1)-vs(i,j))*dyi);
            % if(mod(n,nx)==0)
            if (n > (nx-1)*ny) % if top boundary
                R(n)=0; % force 0 value
            end
            %disp(["R[",n,"]:",R(n)]);
        end
     end


    % Find pressure solution vector using direct method
    pv=L\R;
    % disp("pv="); disp(pv);


    n=0;    
    p=zeros(imax,jmax);
    % convert pressure vector into mesh array
    for j=jmin:jmax
       for i=imin:imax
           n=n+1;
           % if (mod(n,nx)==0)
           if(n > (nx-1)*ny) % if top boundary
               pv(n)=0; % force Zero Dirichlet
           end
           p(i,j)=pv(n);
       end
    end
    
    %disp("p="); disp(num2str(p));

    % correcting (updating) u & v velocities


    for j=jmin:jmax
        for i=imin+1:imax
            u(i,j)=us(i,j)-(dt/rho)*(p(i,j)-p(i-1,j))*dxi;
        end
    end


    for j=jmin+1:jmax
        for i=imin:imax
            v(i,j)=vs(i,j)-(dt/rho)*(p(i,j)-p(i,j-1))*dyi;
        end
    end

    % set Boundary Conditions
    u(:,jmin-1)=u(:,jmin)-2*(u(:,jmin)-u_bot);
    u(:,jmax+1)=u(:,jmax)-2*(u(:,jmax)-u_top);
    v(imin-1,:)=v(imin,:)-2*(v(imin,:)-v_lef);
    v(imax+1,:)=v(imax,:)-2*(v(imax,:)-v_rig);
    
end;


disp("End of simulation");
disp("Last t="); disp(num2str(t));
disp("Last p="); disp(num2str(p));
disp("Last u="); disp(num2str(u));
disp("Last v="); disp(num2str(v));
disp("Finita!");
But it leads to NAN after some time of running (~3.1 s). I tried to thoroughly scan my code for bugs or other types of errors many times. But I can't find anything. I tried to decrease time step (from 1e-2 to 1e-4). Tried to increase grid resolution (from 11x11 to 51x51). No luck. Maybe there is some fundamental flaw in my realization?
P.S. This is my first post here. Please kindly remind me of any "don'ts" I have made in the post if there are any. Thank you.
SashaPrival is offline   Reply With Quote

Old   June 9, 2024, 05:49
Default
  #2
New Member
 
Sasha Prival
Join Date: Jun 2024
Posts: 2
Rep Power: 0
SashaPrival is on a distinguished road
Is there someone who will give me some advice or help in anyway?
SashaPrival is offline   Reply With Quote

Old   June 10, 2024, 13:30
Default
  #3
Senior Member
 
Daniel
Join Date: Feb 2017
Location: Germany
Posts: 162
Rep Power: 9
zacko is on a distinguished road
Did you check at which loop it crashes? Did you print out your matrices? Maybe check the initialalization first and iterate once and recheck the flow field. Maybe there are any obvious problems already apparent
zacko is offline   Reply With Quote

Reply

Tags
cavity flow, dirichlet pressure, matlab, nan error, neumann bcs

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
Ncrit for a glider Xfoil. How to use it. GPT4 answer AlanMattanó Main CFD Forum 0 April 10, 2023 12:16
code for SIMPLE algorithm - 2D Lid driven cavity flow problem - Collocated grid h_amooie OpenFOAM Programming & Development 1 January 22, 2022 11:33
Matlab code for pipe flow cici Main CFD Forum 72 May 12, 2017 18:05
Need Help Here! nan value error mgkid3310 OpenFOAM Running, Solving & CFD 1 October 6, 2016 00:59
Cavity Flow Code Adeel Khan Main CFD Forum 0 April 20, 2014 03:08


All times are GMT -4. The time now is 16:30.