Solving non linear ode using ode45

4 次查看（过去 30 天）
Shubham
Shubham 2022-12-16
评论： Sam Chak 2022-12-16
I want to solve a problem of die casting where a liquid alloy contained in a crucible is pressurised by argon into a tube against gravity. The ODE i have derived is P.A- pho.a.y.g= pho.a.y.(y'').
where pho=2375kg/m^3
P=7x10^5 Pa.
A=0.015362 m^2
A/a=581.45
g=9.8m/s^2
y''= 2nd order time derivative of y;
at t=0,y=0,
at t=tf,y=L=1m;
The simulation ends when y=L is reached.
tf=time to fill the tube.
The problem I am facing is during discretization where the variable y is coming in the denominator. so for the first iteration y=0 yeilds an infinitely large no. Please help
  1 个评论
Walter Roberson
Walter Roberson 2022-12-16
I wonder if this could be interpreted in terms of Mass Matrix? https://www.mathworks.com/help/symbolic/sym.massmatrixform.html

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Walter Roberson
Walter Roberson 2022-12-16
P.A- pho.a.y.g= pho.a.y.(y'')
P and A are nonzero finite constants as are everything else except y and y''
y(0)=0
P*A - pho*a*y(0) = pho*a*y(0)*y''(0)
P*A - 0*constant = constant*0*y''(0)
Assuming y''(0) is finite this gives us
P*A - 0 = 0
and the system is inconsistent unless P or A are 0.
However that logic does not work if in fact y''(0) is infinite, in which case the right hand side has 0*infinity in which case it could potentially be the case that even though you have 0*infinity, that potentially the limit of the product is P*A. But for that to be the case, y(0) = 0 would have to be discontinuous with the rest of y(t), which is a situation that the ode routines cannot deal with.
So... you cannot do that. The equations are inconsistent in a way that cannot be worked with.
Sam Chak
Sam Chak 2022-12-16
编辑：Sam Chak 2022-12-16
Ran in:
Hi @Shubham,
I don't know how you derived this, but after rearranging the equation
you can see the state variable is at the denominator.
Thus, the initial value cannot be absolutely zero, or else it triggers a singularity caused by the Division-by-zero event.
If you select a sufficiently small value for , then you can still run the simulation.
rho = 2375; % liquid aluminum density
P = 7e5; % pressure
A = 0.015362; % area 153.62 cm^2
Aoa = 581.45; % ratio A/a
a = A/Aoa;
g = 9.8;
odefcn = @(t, y) [y(2);
(P*A)/(rho*a*y(1)) - g];
tspan = [0 7.15e-4]; % simulation time
y0 = [0.001 0]; % assume that 1 mm thick layer of residue in the tube
[t, y] = ode45(odefcn, tspan, y0);
tubeLg = y(:,1);
tubeLg(end) % length of tube filled
ans = 1.0000
plot(t, y(:,1), 'linewidth', 1.5), grid on
xlabel('Time, (seconds)')
ylabel('Length of Tube filled, (meters)')
  2 个评论
Shubham
Shubham 2022-12-16
Same thing I did, but the results are not practical. Velocity should be decreasing and in the end due to force balane it should be zero.
V v/s t - should be decreasing with concave upward nature.
Sam Chak
Sam Chak 2022-12-16
Ran in:
@Shubham,
If the velocity should be decreasing, then the dynamics of system (ODE) must be intrinsically stable based on logical thinking.
However, this is not case with the given ODE:
.
Due to the physical quantity (length of tube filled with liquid aluminum) of y is non-negative, , and is positive, then first term is adding energy to the system in a reciprocal manner until the term is dominated by the energy dissipation from the effect of gravity.
Based on the above prediction, if you simulate the system long enough, you should observe an arc like this:
rho = 2375; % liquid aluminum density
P = 7e5; % pressure
A = 0.015362; % area 153.62 cm^2
Aoa = 581.45; % ratio A/a
a = A/Aoa;
g = 9.8;
odefcn = @(t, y) [y(2);
(P*A)/(rho*a*y(1)) - g];
tspan = [0 5.49696e2]; % simulation time
y0 = [0.001 0]; % assume that 1 mm thick layer of residue in the tube
[t, y] = ode45(odefcn, tspan, y0);
plot(t, y(:,1), 'linewidth', 1.5), grid on
TLDR: The model must be inaccurate. Can you recheck the derivation steps?

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by