- 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
15 次查看(过去 30 天)
显示 更早的评论
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 个评论
回答(2 个)
Tanmay Das
2020-6-18
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 个评论
grandhi sai charan
2020-6-19
4 个评论
Tanmay Das
2020-6-22
Hi,
As I said, every function has its own workspace and it gets flushed once the control goes out of the function. So one way to retain variables can be to return it. Previously, sph_path was not returning anything. After changing the first line, the function now returns x, Vx, y and Vy. Thats why the variables are now visible in the workspace after the function is called.
另请参阅
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!