# 2D incompressible cavity flow Matlab code NAN result

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

 June 6, 2024, 06:48 2D incompressible cavity flow Matlab code NAN result #1 New Member   Sasha Prival Join Date: Jun 2024 Posts: 2 Rep Power: 0 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.

 June 9, 2024, 05:49 #2 New Member   Sasha Prival Join Date: Jun 2024 Posts: 2 Rep Power: 0 Is there someone who will give me some advice or help in anyway?

 June 10, 2024, 13:30 #3 Senior Member   Daniel Join Date: Feb 2017 Location: Germany Posts: 162 Rep Power: 9 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

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