Solving non linear ode using ode45

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

[Walter Roberson]

I wonder if this could be interpreted in terms of Mass Matrix? https://www.mathworks.com/help/symbolic/sym.massmatrixform.html

Answer 1

[Walter Roberson]

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.

Answer 2

[Sam Chak]

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)')

[Shubham]

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]

@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?