Problem with solving two differential equations depending on each other

I have a case where I want to simulate the transport of a certain liquid from a silo to a tank. This is done in two parts.

Part 1: Valve 4 is closed and valve 2 is open. The water ring vacuum pump (WRVP) will suck the air out of the tank, causing the pressure in the tank to decrease.

Part 2: Valve 4 is open, valve 2 is closed and the WRVP is turned off. The present underpressure, created during step 1, will cause the liquid in the silo to suck up in the tank. During the process, the underpressure will decrease and the flow rate through the pipeline (5) will decrease as well. At a certain point, the flow rate is zero.

Ignore valve 7.

My question is about part 2. The next differential equation describes the flow rate through the pipeline (5).

Zeta describes the friction losses in the pipeline as follows:

The differential equation is depanded on the variables V(t), W(t) and zeta(W). V(t) is the volume of air and W(t) is the fluid velocity in the pipeline at moment t. The following differential equation dVdt and "normal" equation W(t) make the link with dWdt.

All the following variables are constants: L, rho, p_a, p_0, V_0, g, H_0, d, abs_rough (=triangle), kin_visc (=v), S=pi*d^2/4. So all the variables are constants, except V, W, Q and zeta. Now I want to plot Q(t). My code:

% Given constants

g = 9.81; % N/kg of m/s²

H_0 = 3; % m

p_a = 101325; % Pa

L = 100; % m

d = 0.2; % m

V_0 = 12; % [m³]

p1_end = 1.9630e+04; % [Pa]

abs_rough = 0.0002; % m

kin_visc = 10^-6; % m²/s

rho = 1000; % kg/m³

S = pi*d^2/4; % m²

% Constants dependent on W

delta = abs_rough/d;

Re = @(W) W*d/kin_visc;

lambda = @(Re) 0.11*(delta + 68./Re).^0.25;

zeta = @(lambda) lambda*L/d * 1.1; % assume 10% more to cover the local losses

% Define differential equations

dWdt = @(t,W,V) (((p_a - p1_end*V_0./V)./rho) - g*H_0 - (W.^2)/2*(1+zeta(lambda(Re(W)))))/L;

dVdt = @(t,V,W) -W * S;

% Define initial conditions

W0 = 0;

V0 = V_0;

% Define time interval

tspan = [0 100]; % s

% Solve differential equations

[t, y] = ode45(@(t, y) [dWdt(t, y(1),y(2)); dVdt(t, y(2),y(1))], tspan, [W0; V0]);

% Extract de solution for W(t) and V(t)

W = y(:, 1);

V = y(:, 2);

% Calculate Q(t) from W(t) = Q(t)/S

Q = W * S;

I have no errors but the variables W and V are not as expected:

W = [0 NaN NaN ... NaN] - 41x1 (first value is initial condition)

V = [12 NaN NaN ... NaN] - 41x1 (first value is initial condition)

The solution should approximately be like this: (in my case just one curve)

Answer 1

[Torsten]

Since W0 = 0, you get Re = 0 and thus lambda = something undefined (you divide by Re !).

What does the c stand for in dm^3/c for the unit for Q ?

[Wout Suykerbuyk]

The graph comes from a paper. I think it's a typo and has to be dm³/s.

I changed W0 = 0 to W0 = 1e-6 and it worked!

Meanwhile I have a new challenge...

Until now, the first and second step where seperated: the first step is solved after which the second step is solved. The initial pressure of step two is the final pressure of step 1.

Now I want to start both steps at the same time: the WRVP sucks air (creates underpressure) and, when the pressure is low enough, the liquid starts to transport from the silo to the tank while the WRVP keeps running.

The differential equation for the first step, which I didn't show you yet, is

where V0 = volume of the tank, k = leakage and can be assumed as k = 0, p_a = atmospheric pressure, p(t) = pressure in tank and G = f(p) - characteristic of WRVP. The function G was determined by applying curve fitting on some data points from a datasheet.

Also good to know about part 2: p(t) = p_0*V_0/V(t)

where p_0 = pressure at the end of step 1 (= p1_end).

How can I do this and plot p(t), Q(t), V(t)?

[Torsten]

Use an event function that stops the first step as soon as the pressure is "low enough".

Then restart the integration with the new equation and the results from the first step as new initial conditions.

Look at

https://uk.mathworks.com/help/matlab/math/ode-event-location.html

for an example code.

[Wout Suykerbuyk]

Thanks for your answer, Torsten. I think my explanation wasn't clear enough. The WRVP must continue to run (i.e. continue to suck air from the tank) while the liquid is drawn from the tank due to the negative pressure. From your suggestion I deduce that you want to shut down the WRVP (because you want to stop step 1 as soon as the pressure is low enough).

In the meantime I have made a script myself. dPdt from step 1 and dWdt and dVdt from step 2 are solved together. To be able to suck in the liquid, the negative pressure must be lower than the head pressure. That's why I programmed that the two differential equations associated with transporting the fluid (dWdt and dVdt) should only be active if p <= p_boundary = p_a - rho*g*H_0. With this I simulate that valve 4 only opens when p_limit is reached. That's my intention anyway.

For p(t) and Q(t) I get these results:

However, V(t) remains equal to the initial volume all the time. This means that no liquid is transported and the tank therefore remains empty. However, I think I made a mistake in the script.

% Functions

delta = abs_rough/d;

zeta = @(lambda) lambda*L/d * 1.1; % 10% meer voor de lokale drukverliezen (nog bep)

lambda = @(Re) 0.11*(delta + 68./Re).^0.25;

Re = @(W) W*d/kin_visc;

p = @(V) p_0*V_0./V;

tspan = [0 100];

p_boundary = p_a - rho*g*H_0;

p = p_a; % initial pressure

% Define differential equations

dPdt = @(t,p) G_m*(1-(p_v./p))/(1-(p_v/p_a)).*-p./V_0;

dWdt = @(t,W,p) (((p_a - p)./rho) - g*H_0 - (W.^2)/2*(1+zeta(lambda(Re(W)))))/L*(p<=p_boundary);

dVdt = @(t,V,W) -W * S * (p <= p_grens); % only active when p <= p_boundary

% Define initial conditions

p0 = p_a;

W0 = 1e-6;

V0 = V_0;

% Solve differential equations

[t, y] = ode45(@(t, y) [dPdt(t, y(1)); dWdt(t, y(2), y(1)); dVdt(t, y(3), y(2))], tspan, [p0; W0; V0]);

% Extract p(t), W(t) and V(t)

p = y(:, 1);

W = y(:, 2);

V = y(:, 3);

Q = W * S; % m³/s

Q1 = Q*1000; % dm³/s

[Torsten]

Each time the equations change due to an external event, you should stop the solver using the "event" facility of the ODE integrators, redefine your equations (e.g. add or modify terms) and restart from the solutions obtained so far. This does not mean that you cancel any of the equations, but that you get the chance to modify them.

Discontinuities in the equations resulting from terms like

(W.^2)/2*(1+zeta(lambda(Re(W)))))/L*(p<=p_boundary);
-W * S * (p <= p_grens)

should be avoided.

Integrate up to p_boundary without the term (W.^2)/2*(1+zeta(lambda(Re(W)))))/L, interrupt integration by the event function, modify your function handle by adding the term (W.^2)/2*(1+zeta(lambda(Re(W)))))/L to your equation and continue integration starting with the values obtained so far for the variables.