Trouble with array index in a cfd solver

Hello,

I’m trying to implement a simple cfd solver for the incompressible naiver stokes equation on rectangle domain (Lid-driven cavity)

Boundary conditions are needed for the velocity field and the pressure Poisson equation. The velocity boundary conditions are a bit tricky since some of the velocity components are not defined on the boundary. For example, to specify u = utop at the top of the domain is not straightforward because u is defined dx/2 away from the top boundary. One solution to this problem is to create a fictitious velocity outside the domain such that the velocity on the domain boundary satisfies the boundary condition.

Unfortunately, I made a mistake, which is puzzling me.

Here is the script so far:

clc,clear all,close all
% Number of cells
nx=5;
ny=5;

% lengthof the domain
Lx=1;
Ly=1;

% Material data
rho=1;
nu=1;

% Time 
dt=0.001;
tg=0.01;
kkk=tg/dt;
%CFL=u*dt/dx;

% BC
u_bot=0;
u_top=1;
v_lef=0;
v_rig=0;

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

% Create mesh
x(imin:imax+1)=linspace(0,Lx,nx+1);
y(jmin:jmax+1)=linspace(0,Ly,ny+1);
xm(imin:imax)=0.5*(x(imin:imax)+x(imin+1:imax+1));
ym(jmin:jmax)=0.5*(y(jmin:jmax)+y(jmin+1:jmax+1));
% Create mesh sizes
dx=x(imin+1)-x(imin);
dy=y(jmin+1)-y(jmin);
dxi=1/dx;
dyi=1/dy;

u=zeros(imax,jmax);
v=zeros(imax,jmax);

for ijk=1:1:kkk
    t=ijk*dt;
    % Boundary conditions
    u(:,jmin-1)=2*u_bot-u(:,jmin)
    u(:,jmax+1)=2*u_top-u(:,jmax)
    v(imin-1,:)=2*v_lef-v(imin,:)
    v(imax+1,:)=2*v_rig-v(imax,:)
    % u momentum discretization
    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 discretization
    for j=jmin+1: jmax
        for i=imin : imax
            u_here=0.25*(u(i,j-1)+u(i,j)+u(i+1,j-1)+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

    % Create Laplacian operator for solving pressure Poisson equation
    L=zeros(nx*ny,nx*ny);
    for j=1:ny
        for i=1:nx
            L(i+(j-1)*nx,i+(j-1)*nx)=2*dxi^2+2*dyi^2;
            for ii=i-1:2:i+1
                if(ii>0 && ii<=nx) % Interior point
                    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) % Interiorpoint
                    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
    % Set pressure in first cell (all other pressures w.r.t to this one)
    L(1,:)=0;L(1,1)=1;

    % MATLAB code to compute the right-hand-side R
    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);
        end
    end

    % The pressure is found by solving Lpn+1=R
    pv=LnR;
    % Finally, pv is converted to the mesh representation p(i,j)
    n=0;
    p=zeros(imax,jmax);
    for j=jmin:jmax
        for i=imin:imax
            n=n+1;
            p(i,j)=pv(n);
        end
    end

    % Corrector step
    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
end

Maybe someone has a clue how to solve this

Greetings

Steffen

[Steven Lord]

What is the mistake you've made that's puzzling you?

• Do you receive warning and/or error messages? If so the full and exact text of those messages (all the text displayed in orange and/or red in the Command Window) may be useful in determining what's going on and how to avoid the warning and/or error.
• Does it do something different than what you expected? If so, what did it do and what did you expect it to do?
• Did MATLAB crash? If so please send the crash log file (with a description of what you were running or doing in MATLAB when the crash occured) to Technical Support so we can investigate.

[Steffen B.]

I get an array out of bound error

u =

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

u =

0 0 0 0 0 0 2

0 0 0 0 0 0 2

0 0 0 0 0 0 2

0 0 0 0 0 0 2

0 0 0 0 0 0 2

0 0 0 0 0 0 2

v =

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

v =

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

Index in position 1 exceeds array bounds. Index must not exceed 6.

Error in CFDsolver (line 58)

(nu*(u(i-1,j)-2*u(i,j)+u(i+1,j))*dxi^2 ...

>>