# Trouble with array index in a cfd solver

1 vue (au cours des 30 derniers jours)
Steffen B. le 11 Juil 2023
Commenté : Steffen B. le 12 Juil 2023
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
##### 2 commentairesAfficher AucuneMasquer Aucune
Steven Lord le 11 Juil 2023
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. le 12 Juil 2023
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 ...
>>

Connectez-vous pour commenter.

### Catégories

En savoir plus sur Mathematics and Optimization dans Help Center et File Exchange

R2021b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by