Gracias por dedicar tiempo a esta pregunta. Abajo está la ecuación donde necesito su ayuda.
$$\left ( \frac{dh}{dt} \right )^2 + h\frac{d^2h}{dt^2} = F_{u1} - gh $$ $F_{u1}$ y $g$ son constantes.
La anterior ED no lineal de segundo orden describe el siguiente proceso físico:
-
Consideremos un gancho en forma de "T" que tira de un alambre pesado y largo que está inicialmente en el suelo, con una fuerza constante hacia arriba $ F_u \ [N]$ como se muestra en la siguiente imagen.
-
Que la altura del gancho en el momento $t=0$ sea $h=0$ .
-
A medida que se tira del cable, debido al aumento del peso del mismo, la velocidad del gancho disminuiría, y finalmente se detendría en $t_{end}$ cuando se alcanza el equilibrio de la fuerza estática.
-
La variación de la altura del gancho $h(t)$ se desea conocer.
Para determinar $h(t)$ La segunda ley de Newton se utiliza como sigue: $$ 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 + h\frac{d^2h}{dt^2} $$ $$\left ( \frac{dh}{dt} \right )^2 + h\frac{d^2h}{dt^2} = F_{u1} - gh $$ $$ \text{where } \ F_{u1}=\frac{F_u}{\rho A}$$
Como no se dispone de una solución analítica de forma cerrada, se utilizan diferencias finitas como se indica a continuación:
He intentado solucionarlo con el siguiente código. Como la ecuación anterior es un polinomio de segundo grado en $h^{t+\Delta t}$ habrá dos soluciones. Las soluciones serán de signo contrario hasta que se alcance la altura de equilibrio. Así que el código selecciona sólo el valor positivo de la raíz, que tiene un significado físico.
He tomado la altura inicial del cable, h(2), como 0 m. h(1) se trata como un punto fantasma en el tiempo, es decir, el tiempo t = 0 en h(2), t = dt en h(3), t = 2*dt en h(4) y así sucesivamente. Aquí, h(1) está en t = - dt, y por lo tanto h(1) = h(2). Obtengo la solución como se muestra a continuación, donde la línea roja indica la altura del equilibrio de fuerzas estático.
Como se puede ver en el gráfico de la solución, la solución numérica alcanza una altura de equilibrio que es considerablemente mayor que la altura del equilibrio de fuerzas estático. He probado a reducir el paso de tiempo, pero la diferencia sigue siendo la misma.
¿Podría arrojar algo de luz? Gracias.
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'])
