- Try to use the t that you are passing as argument in path function. It is not required to create a new variable ‘T’.
- Try to use the argument s to find the first order derivatives with respect to s in path function. It is not required to create a new variable 'yy'

# simulation for the path of a particle(sphere) in a water flow

37 views (last 30 days)

Show older comments

My project is to simulate the path of a small particle(sediment) of radius 1mm and weight is 1mg in a steady flow(ex:river) which is released in river with is flowing with a speed of 5m/sec so I derived the velocities for the semdiment in X&Y-directions and get the below equations and want use this values in simulation of this sediment when released in water in which depth of water is 4m and intial velocities are zero.

I am trying to solve the ode's in matlab with ode45.

My modelling equations are m* x''(t) = k*(u-v)*(u-vx),

m*y''(t) = -k*(u-v)*(u-vy)+g(Fb-m).

the flow is steady flow and no slip boundary condition.

k = cd*p*A*0.5

u - velocity profile u = {2(y/h) - (y/h)^2)} * Umax

Umax =5m/sec

v=sqrt(vx^2+vy^2).

vx - velocity particle in X-direction - (u^2*k*t)/(u*k*t-m)

vy - velocity of particle in Y-direction - (Vt*tanhZ)

Vt - terminal velocity - (m*g - Fb)/k

Fb - bouyancy force

Z= Vt*k*t/m.

function sph_path

clear all;

clc;

x0=0;y0=4;Vx0=0;Vy0=0;

IC = [x0,Vx0,y0,Vy0];

tspan = linspace(0,4,100);

[time, state_values] = ode45(@path,tspan,IC);

x = state_values(:,1);

Vx1 = state_values(:,2);

y = state_values(:,3);

Vy1 = state_values(:,4);

figure(3);

plot(x,y);

figure(4);

plot(time,x);

axis([0 10 0 8]);

end

function sdot = path(t,s)

T = linspace(0,4,100);

g=9.8;

Cd = 0.44;%drag coefficient

P = 1000; %density of water

r = 10^-3; %radius od sphere

A = 4*pi*r*r; % radius of sphere

V = (4*pi*r^3)/3; % volume of sphere

Fb = P*V*g; % buoyancy force

K = Cd*P*A/2; %constant

m = 10*10^-5; % mass of sphere

h=4; % depth of the river

yy = linspace(0,4,100);

a = yy/h;

Umax = 5;

U = Umax*((1.5*a)-(a.^3/2)); %velocity of stream

Vt = (m*g - Fb)/K;

Z = Vt*K.*T/m;

Vy = -Vt.*tanh(Z); %velocity of the particle in Y-direction

Vx = (K.*T.*U.*U)./((K.*T.*U)-m);% velocity of particle X-direction

VV = sqrt(Vx.^2+Vy.^2);

hold on

plot(T,Vy,'r-');

plot(T,-Vt,'bo');

hold off

figure(2)

plot(T,Vx);

v1 = (K/m).*(U-VV).*(U-Vx);

v2 = (-K/m).*(U-VV).*(U-Vy)+(Fb-m*g)/m;

v3 = v1'; v4 = v2';

sdot = [ s(2); v3; s(4); v4];

end

I some how tried to improve my code but still I have this error.

Error using odearguments (line 92)

PATH returns a vector of length 202, but the length of initial conditions vector is 4. The vector returned by PATH

and the initial conditions vector must have the same number of elements.

How to solve this issue?

##### 0 Comments

### Answers (2)

Tanmay Das
on 18 Jun 2020

Hi,

I have the understanding that you have written the equations correctly in your code and the equations are differentiable at every point between the time span. For e.g., I found that the following line:

U = Umax*((1.5*a)-(a.^3/2));

does not matches with the equations that you have provided above. Although it has nothing to do with the error that you got, but you may try to write the correct equations in order to get correct result.

However, the error is because ode45 considers the dimension of IC and expects that the function 'path' will return an array of same dimension. ode45 has the following syntax:

ode45(odefun,tspan,y0,options) % where y0 is an array of initial conditions and tspan is the time span

In the above case, ode45 passes tspan (or t) and s into the path function that you have created where s is an array of variables that is required to be computed along with their derivatives. In the above case its x, y, Vx, Vy. Now you should map corresponding derivatives with respect to the variables present in s to sdot.

You may refer to ode45 to get a detailed overview of syntax and working of ode45 solver along with examples.

After going through the documentation, you may note the following points:

The below code block might be useful for understanding:

sdot(1,1) = s(2);

sdot(2,1) = (K/m)*(U-VV)*(U-Vx);

sdot(3,1) = s(4);

sdot(4,1) = (-K/m)*(U-VV)*(U-Vy)+(Fb-m*g)/m;

Hope that helps!

##### 0 Comments

### See Also

### Products

### Community Treasure Hunt

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

Start Hunting!