Mathematics Asked by curious 2 learn on December 5, 2020
Thank you for taking time to look at this question. Below is the equation where I need your help.
$$left ( frac{dh}{dt} right )^2 + hfrac{d^2h}{dt^2} = F_{u1} – gh $$
$F_{u1}$ and $g$ are constants.
The above second-order non-linear DE describes the following physical process:
Consider a "T" shaped hook which pulls a lengthy heavy wire that is initially on the ground, with a constant upward force $ F_u [N]$, as shown in the below image.
Let the height of the hook at time $t=0$ be $h=0$.
As the wire is being pulled, due to increasing weight of the wire, the velocity of the hook would slow down, and it would finally come to rest at $t_{end}$ when static force balance is achieved.
The variation of the hook’s height $h(t)$ is desired to be known.
To determine $h(t)$, Newton’s second law is used as follows:
$$ F_{net} = frac{d (mv)}{dt}, text{where} v=frac{dh}{dt} $$
$$ F_u – F_{wire} = frac{d (rho A h frac{dh}{dt})}{dt}$$
$$ F_u – rho A h g = rho A frac{d (h frac{dh}{dt})}{dt}$$
$$ frac{F_u}{rho A } – h g = left ( frac{dh}{dt} right )^2 + hfrac{d^2h}{dt^2} $$
$$left ( frac{dh}{dt} right )^2 + hfrac{d^2h}{dt^2} = F_{u1} – gh $$
$$ text{where } F_{u1}=frac{F_u}{rho A}$$
As closed-form analytical solution is not available, finite difference is used as follows:
I tried to solve it using the below code. As the above equation is second-degree polynomial in $h^{t+Delta t}$, there will be two solutions. The solutions will be of opposite sign until the equilibrium height is reached. So the code selects only the positive value of the root, which has physical meaning.
I took the initial height of the wire, h(2), as 0 m. h(1) is treated as a ghost point in time, i.e. Time t = 0 at h(2), t = dt at h(3), t = 2*dt at h(4) and so on. Here, h(1) is at t = – dt, and hence h(1) = h(2). I get the solution as shown below, where the red line indicates the height from static force balance.
As you could see from the solution plot, the numerical solution reaches an equilibrium height that is considerably higher than the height from static force balance. I tried reducing the time step, but the difference remains the same.
Could you please shed some light ? Thank you.
clc;
clear all;
rho =1000;
radius =3;
area =pi*(radius/1000)^2;
g = 9.81;
F_u=1;
F_u1 = F_u/(rho*area);
%static force balance
h_static = F_u/(rho*area*g);
dt = 0.00001;
h(1) = 0;
h(2) = 0;
%Model equation Ax^2 + Bx + C = 0
dt_steps = 8/dt;
for i =2:1:dt_steps
A=1;
B=h(i);
C = (h(i))^2-(h(i)*h(i-1))+((F_u1-(g*h(i)))*dt^2);
polynomial_roots=roots([A -B -C]);
if (length(polynomial_roots(polynomial_roots>0))==1)
h(i+1)=polynomial_roots(polynomial_roots>0);
else
break
end
end
scatter(((1:length(h))-2)*dt,h);
hold on;
plot([0 (length(h)-2)*dt],[h_static h_static],'-r');
xlim([0 (length(h)-2)*dt]);
xlabel('Time [s]')
ylabel('Height [m]')
title(['dt = ' num2str(dt),' s'])
You seem to think in terms of a first order DE, where indeed one would expect a limit at the root $A/B$ of the right side. However, you have a second order DE that has a potential term and is conservative. More specifically, multiply with $h'$ and integrate to get $$ frac12h'(t)^2=Aln h(t)-Bh(t)+C $$ which has the form $frac12h'(t)^2+V(h(t))=C$ with a convex function $V$. The expected behavior is an oscillation in the segment of the level set ${h:V(h)le C}$ that contains the initial point, here there is only one segment.
You will have to contemplate again on the physics of your model, you need either a friction term, some other kind of energy dissipation, or a good reason why the behavior can be modeled with a first order DE.
Answered by Lutz Lehmann on December 5, 2020
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP