# function call and getting array error

4 views (last 30 days)

Show older comments

Thin Rupar Win
on 9 Nov 2021

Commented: Steven Lord
on 10 Jan 2024

Dear Sir or Madam,

I recently learn matlab for my education purpose. I would like to know how to solve the Index error in my code as I call function in for loop and getting error. I have already check array boundary and can't find the answer for the case. I deeply request you how to solve this case.

clear

% nx and ny must be the same size with t, n_part of the DEM loop

nx=4;ny=9;

% f_ equation for collisions

f=zeros(nx,ny,9);

feq=zeros(nx,ny,9);feq_f=zeros(nx,ny,9);

% velocities for x and y direction

u=zeros(nx,ny);v=zeros(nx,ny);

% velocities on solid particle

% for x direction

U_s_x=zeros(nx,ny,9);U_px=[0,0];Omega_px=[0,0];Xpx=[0,0];

% for y direction

U_s_y=zeros(nx,ny,9);U_py=[0,0];Omega_py=[0,0];Xpy=[0,0];

% position of particle and omega_s

omega_s=zeros(nx,ny,9);

% force term Ff,Tf on solid particle

Ffx=zeros(nx,ny,9);Ffy=zeros(nx,ny,9);

Tfx=zeros(nx,ny,9);Tfy=zeros(nx,ny,9);

rho=ones(nx,ny);x=zeros(nx);y=zeros(ny);

w(9)=zeros;

% initialization of the system

w=[1/9 1/9 1/9 1/9 1/36 1/36 1/36 1/36 4/9];

cx=[1 0 -1 0 1 -1 -1 1 0];

cy=[0 1 0 -1 1 1 -1 -1 0];

c2=1./3.;

% derivative rate

dx=1.0;dy=1.0;

x1=(nx-1)/(ny-1);y1=1.0;

dx=x1/(nx-1);

dy=y1/(ny-1);

x=(0:dx:x1);

y=(0:dy:y1);

u0=0.2;

% relaxiation time

alpha=0.02;

Re=u0*(ny-1)/alpha;

% omega = 1/tau

omega=1./(3.*alpha+0.5);

% toleriant and error

%count=0;tol=1.0e-4;error=10.;erso=0.0;

% setting velocity

for j=2:ny-1

u(1,j)=u0;

end

t_lbm=50;

for t=1:t_lbm

[f,Ffx,Ffy,Tfx,Tfy]=up_collition(nx,ny,u,v,cx,cy,U_s_x,U_px,Omega_px,Xpx,U_s_y,U_py,Omega_py,Xpy,feq,rho,w,feq_f,f,omega_s,Ffx,Ffy,Tfx,Tfy);

end

function[f,Ffx,Ffy,Tfx,Tfy]=up_collition(nx,ny,u,v,cx,cy,U_s_x,U_px,Omega_px,Xpx,U_s_y,U_py,Omega_py,Xpy,feq,rho,w,feq_f,f,omega_s,Ffx,Ffy,Tfx,Tfy)

% the weighting function of solid particles k Bk,

% the local solid ratio e_k divided

e_k=[1/8 1/8 1/8 1/8 1/8 1/8 1/8 1/8 0];

e_total=1;

for j=1:ny

for i=1:nx

t1=u(i,j)*u(i,j)+v(i,j)*v(i,j);

for k=1:9

x=zeros(length(i));% later add with dx term, let dt =1;

t2=u(i,j)*cx(k)+v(i,j)*cy(k);

Bk(k)=e_k(k).*(0.56-0.5)/(1-e_total+(0.56-0.5));

% velocity on solid particle

% the equaiton is Us=UP+omega_px *(x +0.5 *c(k)*dt)-Xp

U_s_x(i,j,k)=U_px(i,j)+Omega_px(i,j)*(x(i)+0.5*(cx(k)))-Xpx(i,j);

U_s_y(i,j,k)=U_py(i,j)+Omega_py(i,j)*(x(i)+0.5*(cy(k)))-Xpy(i,j);

t22=U_s_x(i,j,k)*cx(k)+U_s_y(i,j,k)*cy(k);

t11=U_s_x(i,j,k)*U_s_x(i,j,k)+U_s_y(i,j,k)*U_s_y(i,j,k);

feq(i,j,k)=rho(i,j)*w(k)*(1.0+3.0*t2+4.5*t2*t2-1.5*t1);

feq_f(i,j,k)=rho(i,j)*w(k)*(1.0+3.0*t22+4.5*t22*t22-1.5*t11);

f(i,j,k)=f(i,j,k)-(0.9091*(f(i,j,k)-feq(i,j,k)))+(Bk(k)*(feq_f(i,j,k)-feq(i,j,k)));

% additional collision terms

omega_s(i,j,k)=feq_f(i,j,k)-f(i,j,k)+((f(i,j,k)-feq(i,j,k)).*0.42);

% force term and torque term for DEM

Ffx(i,j,k)=sum(Bk(k)).*sum(omega_s(i,j,k))*cx(k);

Ffy(i,j,k)=sum(Bk(k)).*sum(omega_s(i,j,k))*cy(k);

Tfx(i,j,k)=sum((x(i)-Xpx)*Bk(k)).*sum(omega_s(i,j,k)).*cx(k);

Tfy(i,j,k)=sum((x(i)-Xpy)*Bk(k)).*sum(omega_s(i,j,k)).*cy(k);

end

end

end

end

##### 4 Comments

Rik
on 9 Nov 2021

### Accepted Answer

Sudharsana Iyengar
on 9 Nov 2021

U_s_x is a 4x9x9 matrix but U_px is 1x2 matrix. So there is no U_px(3,4) did you check on that.

U_s_x(i,j,k)=U_px(i,j)+Omega_px(i,j)*(x(i)+0.5*(cx(k)))-Xpx(i,j);

U_s_y(i,j,k)=U_py(i,j)+Omega_py(i,j)*(x(i)+0.5*(cy(k)))-Xpy(i,j);

##### 3 Comments

Steven Lord
on 10 Jan 2024

### More Answers (0)

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!