Ordinary differential equation - pratical problem
1 回表示 (過去 30 日間)
古いコメントを表示
Hello,
I have a single equation of the form: u'=a*u^2+b*u^1.5+c*u+d. I am trying to solve it in Matlab but it does not find the solution. The problem is a body falling down. The body is under the effect of four forces: 1) weight (m*g) 2) flow in opposite direction of the weight (0.5*da*A*U0^2) 3) drag force (0.5*da*C_d*A*u^2), which would change verse depending on flow direction 4) acceleration (m*u') u(t=0)=0;
The C_d changes with velocity (relationship from the sphere in free flow C_d = 21.12/Reynolds+6.3/sqrt(Reynolds)+.25; % drag coefficient)
I coded this problem, but I can't get it to give a result, unless I fix C_d to a value. Would you be able to give a help, please?
da = 1.161; %density of air
g = 9.81; % Gravitational acceleration
d =1E-6*100; % characteristic length of the body
m = 1000*.957251*d.^3.09275; %Mass of the the body
A = 3.3108*d.^2.21672; %Surface area of the body
u_flow = .5; % Velocity air flow
syms u(t) c_D Re_ver
Re_ver = da*u(t)*d/1.8E-5; % Reynolds number
c_D = 21.12/Re_ver+6.3/sqrt(Re_ver)+.25; % drag coefficient;
a = m;
b = -.5*da*A*c_D;
c = m*g-.5*da*u_flow^2*A;
if c<0
c = -c;
end
u(t) = dsolve(a*diff(u,t) == b*u^2+c,u(0)==0);
t = linspace(0,10,500);
figure, plot(t,double(u(t)))
I am using Matlab R2014a
Thanks Antonio
8 件のコメント
Karan Gill
2016 年 7 月 13 日
In the ODE a*u' = b*u^2+c , b has a term 1/sqrt(u). My guess is that this term complicates the calculation. The ODE is of the form:
C*u' = u^2(C + C/u + C/sqrt(C*u)) + C
|----'b' term-------|
where C represents the constants. I'm afraid I don't see how b = constant/u. Could you point out where I misunderstand?
Walter Roberson
2016 年 7 月 13 日
At your initial condition, u(0) = 0, your Reynold's number is 0. That is giving you a division by 0 for your drag coefficient, which is giving a singularity right from the start, which is preventing any resolution of the system.
I am not familiar with fluid dynamics, but when I look at various material it appears to me that the equation you are using might perhaps not apply at very low velocities.
回答 (0 件)
参考
製品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!