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

37 views (last 30 days)
grandhi sai charan on 16 Jun 2020
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?

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:
1. Try to use the t that you are passing as argument in path function. It is not required to create a new variable ‘T’.
2. 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'
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!

grandhi sai charan on 19 Jun 2020
Hi Thanmay das,
Thanks alot for responding to my problem.
U = Umax*((1.5*a)-(a.^3/2));
This is velocity profile of the river or stream which is a parabolic profile so for simplicity I considered this uniform and steady flow and it won't depend on time.
Earlier I tried without using T and also yy but I do get vercat error always so I did this changes and implemented the code.
But thanks after referring to ode45 I tried to use interp1 function and I got result.
here is the code for it.
function sph_path
clear all;
clc; clf;
x0(:)=0;y0(:)=4;Vx0(:)=0;Vy0(:)=0;
IC = [x0,Vx0,y0,Vy0];
tspan = [0 4];
[time, state_values] = ode45(@(time,state_values) path(time,state_values),tspan,IC);
x = state_values(:,1);
Vx1 = state_values(:,2);
y = state_values(:,3);
Vy1 = state_values(:,4);
figure(3);
plot(x,y);
end
function sdot = path(t,s,v1,v2)
T = linspace(0,4,100);
g=9.8;
Cd = 0.44;%drag coefficient
P = 1000; %density of water
P1 = 1600; %density of sphere
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 = P1*V; % 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). * abs (U-VV). * (U-Vx);
v2 = (-K/m).*abs(U-VV).*(U-Vy)+(Fb-m*g)/m;
v1 = interp1(T,v1,t);
v2 = interp1(T,v2,t);
sdot = [ s(2); v1; s(4); v2];
end
I have 2issues here
1) I get the path of the sediment but I can't see the values in my workspace.
2) I am trying to get the relation between the various values of density and the distance travelled by the sediment, is there any idea to do it?
Thanks alot once again.
grandhi sai charan on 22 Jun 2020
Hii,
Thanks alot for the information regarding how to call the values in workspace too.
I went through your suggestions for coding and tried not to use T and used only t in the functions Vx,Vy but I get concatenation error.
Is there any way to solve this issue because using T is not a good way and results are changing alot for various paramaters.