Trouble with array index in a cfd solver

2 次查看（过去 30 天）
Steffen B.
Steffen B. 2023-7-11
评论： Steffen B. 2023-7-12
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 个评论
Steven Lord
Steven Lord 2023-7-11
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.
Steffen B. 2023-7-12
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 ...
>>

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Matrix Indexing 的更多信息

标签

产品


版本

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by