I am following this guide to create a simple cfd solver for incompressible 2D fluid in a square cavity. The following Matlab code sample from the guide creates a Laplacian (coefficient) matrix with zero-derivative (Neumann-type) conditions imposed on all four (left/right/bottom/top) boundaries. L is the laplacian coefficient matrix of the matrix equation L*P=R, where P=P(x,y) is the unknown (solution) 2d pressure function, R(x,y) - source function.
nx = 5;
ny = 5;
Lx = 1;
Ly = 1;
dx = Lx/(nx); %grid x-spacing
dy = Ly/(ny); %grid y-spacing
dxi = 1/dx;
dyi = 1/dy;
D = 2*dxi^2 + 2*dyi^2;
Dx = D - 1/dx^2;
Dy = D - 1/dy^2;
Dxy = D - 1/dx^2 - 1/dy^2;
L = zeros(nx*ny,nx*ny);
% Creat Laplacian operator for solving pressure Poisson equation
for j=1:ny
for i=1:nx
L(i+(j-1)*nx, i+(j-1)*nx)=D;
for ii=i-1:2:i+1
if (ii>0 && ii<=nx) % Interior points
L(i+(j-1)*nx,ii+(j-1)*nx)=-dxi^2;
else % Neuman conditions on boundary
L(i+(j-1)*nx,i +(j-1)*nx)= ...
L(i+(j-1)*nx,i +(j-1)*nx)-dxi^2;
end
end
for jj=j-1:2:j+1
if (jj>0 && jj<=ny) % Interior points
L(i+(j-1)*nx,i+(jj-1)*nx)=-dyi^2;
else % Neuman conditions on boundary
L(i+(j-1)*nx,i +(j-1)*nx)= ...
L(i+(j-1)*nx,i +(j-1)*nx)-dyi^2;
end
end
end
end
L(1,:)=0;
L(1,1)=1; % to avoid singularity
This is the output of the code. All four boundaries are zero-derivative (Neumann-type).
Then I wanted to create the Laplacian matrix with Dirichlet condition (i.e. P(x,y_top)=0) imposed on the top boundary, the other boundaries are with same Neumann conditions. The structure and contents of the (Nx x Ny, Nx x Ny) Laplacian matrix depend on how we arrange the unknowns (solutions) into a (Nx x Ny) vector.
Assuming that in the guide the pressure array elements P(x_i,y_j)=Pij were arranged into a (Nx x Ny) vector in so-called "natural order":
, I have modified the code for the Laplacian matrix creation like this:
% creating 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
When Nx = Ny = 5, the Laplacian matrix looks like this, which I believe is correct:
Then I wrote a complete Matlab program further following the guide:
clear; clc;
% grid parameters
Lx=1; % x length
Ly=1; % y length
nx=5; % number of grid points along x
ny=5; % number of grid points along y
%simulation parameters
dt=0.001; % time step
%t_final=0.1; % simulation run time
t_final=0.2; % simulation run time
% physics parameters
nu=0.001; % kinematic Viscosity
rho=1.0; % density
% 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);
%disp(["dx=",dx]);
%disp(["dy=",dy]);
% Create mesh sizes
dxi=1/dx;
dyi=1/dy;
% Number of timesteps
Nt=t_final/dt;
% preallocate important arrays
p=zeros(imax,jmax); % ok
%p=zeros(nx,ny); % ok
us=zeros(imax+1,jmax+1);
vs=zeros(imax+1,jmax+1);
% R=zeros(imax,1);
R=zeros(nx,1);
u=zeros(imax+1,jmax+1);
v=zeros(imax+1,jmax+1);
L=zeros(nx*ny,nx*ny);
% 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
% L(1,:)=0; L(1,1)=1; % since we have Dirichlet condition,
% we don't need this anymore
disp("L="); disp(num2str(L));
% solver
disp("Starting simulation");
while t <= t_final
% update time
t = t + dt;
disp("t:");
disp(t);
u_top=1;
% u Momentum 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 Momentum 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 the r.h.s (R) of PPEquation using
% the predicted velocities vs and us.
n=0;
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(["Time: ",t]);
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 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^(n+1) & v^(n+1) velocities
for j=jmin:jmax
for i=imin+1:imax
u(i,j)=us(i,j)-dt/rho*(p(i-1,j-1)-p(i-2,j-1))*dxi;
end
end
for j=jmin+1:jmax
for i=imin:imax
v(i,j)=vs(i,j)-dt/rho*(p(i-1,j-1)-p(i-1,j-2))*dyi;
end
end
% 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 p="); disp(num2str(p));
disp("Last u="); disp(num2str(u));
disp("Last v="); disp(num2str(v));
disp("Finita!");
But it leads to NAN result at around 0.006 s time. I am stuck with this situation. Either I can't find a bug, or my take is wrong. Maybe I am missing something?
