Can not solve system of differential equations with ode45 taking into account continuity equation

I have a schematic of the following situation.

I have described the fluid velocity (c1 and c2) through the 2 pipe diameters, pressure p and volume of air in tank V, all in function of the time, with the following 5 differential equations:

where Q = c1 * S1

Note that the sum of static and dynamic pressure is the total pressure and the total pressure is equal through the whole pipe. This can be described with the following continuity equation.

There are 6 variables that change in time and therefore need to be solved with the differential equations and the continuity equation: p1L, p1R, p, c1, c2 and V.

I have made the following code.

% Functions that calculate friction losses due to pipes
zeta1 = @(c1) (0.11*(abs_rough/d1 + 68./(c1*d1/kin_visc)).^0.25)*L1/d1 * 1.1;
zeta2 = @(c2) (0.11*(abs_rough/d2 + 68./(c2*d2/kin_visc)).^0.25)*L2/d2 * 1.1;

% Continuity equation
syms p1L p1R c1 c2
convgl = p1L + rho .* c1.^2./2 == p1R + rho .* c2.^2./2;
f_equation = matlabFunction(convgl, 'Vars', {p1L, p1R, c1, c2});

% Define differential equations
dPdt = @(t,p) data_debiet(end)*(1-(p_v./p))/(1-(p_v/p_a)).*-p./V_0;
dCdt1 = @(t,c1,p1L) ((p_a - p1L)./rho - g*(H_1 - H_0) - (c1.^2)/2*(1+zeta_1(c1)))/L1;
dCdt2 = @(t,c2,p1R,p) ((p1R - p)./rho - g*(H_2 - H_1) - (c2.^2)/2*(1+zeta_2(c2)))/L2;
dVdt = @(t,V,c1) -c1 * S1;

% Define initial conditions
p0 = p_a;
c0 = 1e-6;
V0 = V_0;

% Solve differential equations
[t, y] = ode45(@(t, y) [dPdt(t, y(3)); dCdt1(t, y(4), y(1)); dCdt2(t, y(5), y(2), y(3)); dVdt(t, y(6), y(4))], tspan1, [p0; p0; p0; c0; c0; V0], options1);

% Extract results from y
p1L = y(:, 1);
p1R = y(:, 2);
p = y(:, 3);
c1 = y(:, 4);
c2 = y(:, 5);
V = y(:, 6);

Assume that all the other variables are known constants. I have the following error message.

Error using odearguments

@(T,Y)[DPDT(T,Y(3));DCDT1(T,Y(4),Y(1));DCDT2(T,Y(5),Y(2),Y(3));DVDT(T,Y(6),Y(4))] returns a vector of length 4, but the length of initial conditions vector is 6. The vector returned by

@(T,Y)[DPDT(T,Y(3));DCDT1(T,Y(4),Y(1));DCDT2(T,Y(5),Y(2),Y(3));DVDT(T,Y(6),Y(4))] and the initial conditions vector must have the same number of elements.

Error in ode45 (line 107)

odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in Model_v4 (line 81)

[t, y] = ode45(@(t, y) [dPdt(t, y(3)); dCdt1(t, y(4), y(1)); dCdt2(t, y(5), y(2), y(3)); dVdt(t, y(6), y(4))], tspan1, [p0; p0; p0; c0; c0; V0], options1);

Answer 1

[Torsten]

Where do you specify the equations for p1L and p1R ?

In your vector of time derivatives

[dPdt(t, y(3)); dCdt1(t, y(4), y(1)); dCdt2(t, y(5), y(2), y(3)); dVdt(t, y(6), y(4))]

I only see equations for P, C1, C2 and V.

Even if you included the Bernoulli equation, one further equation had to be supplied.

[Wout Suykerbuyk]

I have added another continuity equation: c1*S1 = c2*S2.

I also changed the number of initial conditions from 6 to 4, because I have 4 differential equations:

[p0; p0; p0; c0; c0; V0] is changed to [p0; c0; c0; V0].

So now I have 4 diff eq and 2 extra equations for 6 unknown variables but it still doesn't work. Do I need to 'connect' the differential equations with the extra equations in some way?

This is the new error:

Index exceeds the number of array elements. Index must not exceed 4.

Error in Model_v4>@(t,y)[dPdt(t,y(3));dCdt1(t,y(4),y(1));dCdt2(t,y(5),y(2),y(3));dVdt(t,y(6),y(4))] (line 86)

[t, y] = ode45(@(t, y) [dPdt(t, y(3)); dCdt1(t, y(4), y(1)); dCdt2(t, y(5), y(2), y(3)); dVdt(t, y(6), y(4))], tspan1, [p0; c0; c0; V0], options1);

Error in odearguments (line 92)

f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 107)

odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in Model_v4 (line 86)

[t, y] = ode45(@(t, y) [dPdt(t, y(3)); dCdt1(t, y(4), y(1)); dCdt2(t, y(5), y(2), y(3)); dVdt(t, y(6), y(4))], tspan1, [p0; c0; c0; V0], options1);

[Torsten]

Please include your modified code.

[Wout Suykerbuyk]

Meanwhile I changed something else: the two continuity equations are included in the ode45 solver. This is my complete code except for the given constants.

% Formulas for pressure losses
zeta1 = @(c1) (0.11*(abs_rough1/d1 + 68./(c1*d1/kin_visc)).^0.25)*L1/d1 * 1.1;
zeta2 = @(c2) (0.11*(abs_rough2/d2 + 68./(c2*d2/kin_visc)).^0.25)*L2/d2 * 1.1;
% options for ODE solvers
const.p_a = p_a;
const.rho = rho;
const.g = g;
const.H_0 = H_0;
options1 = odeset('Events', @(t, y) event_function1_m4(t, y, const)); % stop solver when pressure difference > head pressure
options2 = odeset('Events', @event_function2_m4); % stop solver if tank is completely filled (V=0)
%% Solution part 1: creating underpressure (underpressure < p_head)
% In this part, there will be no liquid suction so dCdt1, dCdt2 and the two continuity equations are temporary equal to zero.
tspan1 = [0 50]; % normally head pressure is reached earlier
% Define differential equations
dPdt = @(t,p) data_debiet(end)*(1-(p_v./p))/(1-(p_v/p_a)).*-p./V_0;
% Define initial conditions
p0 = p_a;
c0 = 1e-6;
V0 = V_0;
% Solve differential equations until boundary value
systeem = @(t,y) [dPdt(t,y(1))];
[t, y] = ode45(systeem, tspan1, p0, options1);
% Extract solutions
% c1, c2 and V doesn't change during part 1 because there is no fluid flow
t_1 = t;
p_1 = y(:, 1);
p1L_1 = p_1;
p1R_1 = p_1;
c1_1 = repmat(c0, length(p_1),1);
c2_1 = repmat(c0, length(p_1),1);
V_1 = repmat(V0, length(p_1),1);
%% Solution part 2: suction of liquid (underpressure > p_head)
tspan2 = [t(end) 50];
% Define new differential equations
dCdt1 = @(t,c1,p1L) ((p_a - p1L)./rho - g*(H_1 - H_0) - (c1.^2)/2*(1+zeta1(c1)))/L1;
dCdt2 = @(t,c2,p1R,p) ((p1R - p)./rho - g*(H_2 - H_1) - (c2.^2)/2*(1+zeta2(c2)))/L2;
dVdt = @(t,V,c1) -c1 * S1;
% Define continuity equations
convgl1 = @(p1L,c1,p1R,c2) p1L + rho .* c1.^2./2 - p1R - rho .* c2.^2./2;
convgl2 = @(c1,c2) c1 .* S1 - c2 .* S2;
% Define new initiial condition for p0
p0 = p_1(end);
% Combine differential and continuity equations in one system
systeem = @(t,y) [dPdt(t,y(1)); dCdt1(t,y(4),y(2)); dCdt2(t,y(5),y(3),y(1)); dVdt(t,y(6),y(4)); convgl1(y(2),y(4),y(3),y(5)); convgl2(y(4),y(5))];
% Solve differential equations starting from boundary value
[t, y] = ode45(systeem, tspan2, [p0; p0; p0; c0; c0; V0], options2);
% Extract solutions
t_2 = t;
p_2 = y(:, 1);
p1L_2 = y(:, 2);
p1R_2 = y(:, 3);
c1_2 = y(:, 4);
c2_2 = y(:, 5);
V_2 = y(:, 6);
%% Combine results of part 1 and 2
% note: last value of part 1 = first value of part 2 → first value overlaps
t = [t_1; t_2(2:end)];
p1L = [p1L_1; p1L_2(2:end)];
p1R = [p1R_1; p1R_2(2:end)];
p = [p_1; p_2(2:end)]; % pressure in tank
c1 = [c1_1; c1_2(2:end)]; % velocity through hose
c2 = [c2_1; c2_2(2:end)]; % velocity through riser
V = [V_1; V_2(2:end)];
% Calculate flow rates with Q = c*S, they should be equal to each other
Q1 = c1 * S1; % flow rate through hose
Q2 = c2 * S2; % flow rate through riser

Note that tspan2 is limited to 50 seconds while normally it would take longer. If I should make tspan2 longer, Matlab would give an internal error of the event function (which is set to V=0) because it would take to many steps so reach V = 0. As you can see in the results for tspan2 = 50s, flow rate Q through the hose remains zero. Normally, the tank should be filled within 200s.

Now I made two functions of the continuity equations to solve them within ode45 but I don't know if it is allowed to do that, knowing that the continuity equations are

Results: (ignore graph at right bottom)

[Torsten]

The equation

c1 * S1 - c2 * S2 = 0

does not help.

On the contrary: Most probably, it will contradict the results of the two differential equations where you determine c1 and c2.

You have the Bernoulli equation in which p_1,L and p_1,R appear. But you need a second equation that relates p_1,L and p_1,R to close the system.

There are empirical correlations for pressure losses over regions where the tube radius changes. You should look them up in an engineering handbook (e.g. VDI Wärmeatlas).

What I don't understand is why you use time-dependent equations (p,c1,c2,V) together with stationary ones (Bernoulli equation, static pressure loss over a tube expansion). Isn't that inconsistent ?

By the way:

If you want to include algebraic equations in the system for ODE45 (like equations 5 and 6 of your system), you must modify the mass matrix M.

Usually, ODE45 solves a system of the form

M*y' = f(t,y)

where M is the identity matrix (pure system of ordinary differential equations).

If equations 5 and 6 are algebraic equations like in your case, you have to define

M = eye(6)
M(5,5) = 0;
M(6,6) = 0;

and pass M to ODE45 via the mass matrix option:

options = odeset('Mass',M)

[Wout Suykerbuyk]

I found the following equation:

p1L - p1R = rho*c1^2/2 * (1-S1/S2)^2

I filled it in in dCdt1 so p1R is gone.

But the formula above is actually an indirect formula to Bernoulli's principle which I already applied.

I really don't know how to fix this.

I know the pressure at the beginning, which is the atmospheric pressure p_a, and the pressure at the end, which is the tank pressure p(t). For my goal I actually don't need c1(t) or c2(t). I only need the flow rate Q(t) to determine how fast the tank is filled. But I think it's not possible to determine Q(t) without determining c1(t) and c2(t).

[Torsten]

https://www.youtube.com/watch?v=eZNRsRjX7rs

[Wout Suykerbuyk]

I tried something new: I included the enlargement of the pipe in the first section. So the enlargement is now a local friction loss of the first section. As a result, p1L and p1R are replaced by just p1.

So I have 4 differential equations (dPdt, dCdt1, dCdt2 and dVdt), 1 equation (c1*S1 = c2*S2) and 5 variables (p, p1, c1, c2 and V). I implemented the equation in ode45 by modifying it to

d(lambda)/dt = c1 * S1 - c2 * S2

where d(lambda)/dt = 0

and initial condition lambda0 = c0*S1 - c0*S2

I have tried it in this way because I don't understand how to implement the mass matrix M.

Here is my modified code:

% Formulas for pressure losses
zeta1 = @(f1,K_exp) f1 * L1 / d1 + K_exp;
zeta2 = @(f2,K_bocht) f2 * L2 / d2 + 2 * K_bocht;
f1 = @(Re1) 0.11 * (abs_rough1/d1 + 68/Re1)^0.25;
f2 = @(Re2) 0.11 * (abs_rough2/d2 + 68/Re2)^0.25;
Re1 = @(c1) c1 * d1 / kin_visc;
Re2 = @(c2) c2 * d2 / kin_visc;
K_exp = @(f1) (1 + 0.8 * f1) * (1 - (d1/d2)^2)^2;      % plotse verbreding
K_bocht = 0.45;                                        % R/D = 1.5 van 90°
% Continuity equation
% c2 = @(c1) c1 * S1 / S2;
% options for ODE solvers
const.p_a = p_a;
const.rho = rho;
const.g = g;
const.H_0 = H_0;
const.H_1 = H_1;
options1 = odeset('Events', @(t, y) event_function1_m4(t, y, const));
options2 = odeset('Events', @event_function2_m4);
%% Solution part 1: onderdruk creëren (onderdruk < p_opvoer)
tspan1 = [0 50];
% Define differential equations
dPdt = @(t,p) data_debiet(end)*(1-(p_v./p))/(1-(p_v/p_a)).*-p./V_0;
% Define initial conditions
p0 = p_a;
c0 = 1e-6;
V0 = V_0;
lambda0 = c0*S1 - c0*S2;
% Solve differential equations until head pressure is reached
systeem = @(t,y) [dPdt(t,y(1))];
[t, y] = ode45(systeem, tspan1, p0, options1);
% Extract the solutions
t_1 = t;
p_1 = y(:, 1);
p1_1 = y(:, 1);
c1_1 = repmat(c0, length(p_1),1);
c2_1 = repmat(c0, length(p_1),1);
V_1 = repmat(V0, length(p_1),1);
%% Solution part 2: vloeistof aanzuigen (onderdruk > p_opvoer)
tspan2 = [t(end) 200];
% Define new differential equations
dCdt1 = @(t,c1,p1) ((p_a - p1) - rho * g * (H_1 - H_0) - rho * c1^2 / 2 * (1 + zeta1(f1(Re1(c1)),K_exp(f1(f1(Re1(c1))))))) / L1;
dCdt2 = @(t,c2,p1,p) ((p1 - p) - rho * g * (H_2 - H_1) - rho * c2^2 / 2 * (1 + zeta2(f2(Re2(c2)),K_bocht))) / L2;
dVdt = @(t,V,c1) -c1 * S1;
% Define new initial conditions for p0
p0 = p_1(end);
% Solve differential equations after head pressure is reached
systeem = @(t,y) [dPdt(t,y(1)); dCdt1(t,y(3),y(2)); dCdt2(t,y(4),y(2),y(1)); dVdt(t,y(5),y(3)); y(3)*S1 - y(4)*S2];
[t, y] = ode45(systeem, tspan2, [p0; c0; c0; V0; lambda0], options2); % beginvoorwaarden ordenen volgens volgorde in y
% Extract solutions
t_2 = t;
p_2 = y(:, 1);
p1_2 = y(:, 2);
c1_2 = y(:, 3);
c2_2 = y(:, 4);
V_2 = y(:, 5);
%% combine results
t = [t_1; t_2(2:end)];
p = [p_1; p_2(2:end)]; % druk in tank
p1 = [p1_1; p1_2(2:end)];
c1 = [c1_1; c1_2(2:end)]; % snelheid door de slang
c2 = [c2_1; c2_2(2:end)]; % snelheid door stijgbuis
V = [V_1; V_2(2:end)];

Q1 = c1 * S1; % debiet door slang
Q2 = c2 * S2; % debiet door stijgbuis

Results are still not as expected. Q(t) should increase and V(t) should decrease gradually after the head pressure is reached but instead Q(t) doesn't change and V(t) is zero immediately.

[Torsten]

As you can see in the video, a simplifying assumption is that there is a pressure jump from p_inlet to p_outlet at the position of pipe enlargement. Friction losses in the pipe are neglected, but you can include them in the Bernoulli equation together with the influence of gravity because of height differences. I think this stationary approach to compute the volume flow rate in steady-state would be a good starting point.

Concerning your code:

As I already wrote, c1 and c2 are determined by the differential equations dCdt1 and dCdt2. The continuity equation c1*S1 - c2*S2 cannot be used to determine a variable different from c1 or c2. Thus either you remove one of the differential equations and replace it by the continuity equation or you will get problems with your system. You definitly need an equation that computes p1.

[Wout Suykerbuyk]

I have tried a new way.

The following differential equation determines the fluid velocity through the hose with diameter d1 in function of time.

All parameters in the above differential equation dc1/dt are known, except p1(t). To determine this parameter in relation to time, I use Bernoulli's equation between p1 and p2.

The two unknowns in this equation are p2(t) and c2(t).

• the pressure in the tank p2(t) is calculated with the differential equation dp/dt. The solutions for dp/dt is an array p2(t) with values in function of time.
• c2(t) can be calculated with the continuity equation c1(t) * S1 = c2(t) * S2

Now, the formulas for p1(t), zeta1, zeta2 and the continuity equation are defined as functions in my script. Maybe that's why it does not work. Could it be solved by include them in the ode45 or some other way?

% Formulas for pressure losses
zeta1 = @(f1,K_exp) f1 * L1 / d1 + K_exp;
zeta2 = @(f2,K_bocht) f2 * L2 / d2 + 2 * K_bocht;
f1 = @(Re1) 0.11 * (abs_rough1/d1 + 68/Re1)^0.25;
f2 = @(Re2) 0.11 * (abs_rough2/d2 + 68/Re2)^0.25;
Re1 = @(c1) c1 * d1 / kin_visc;
Re2 = @(c2) c2 * d2 / kin_visc;
K_exp = @(f1) (1 + 0.8 * f1) * (1 - (d1/d2)^2)^2;      % sudden widening
K_bocht = 0.45;                                        % R/D = 1.5 van 90°
% Continuity equation
% c2 = @(c1) c1 * S1 / S2;
% Formula for p1
p1 = @(p,c1) p + rho * g * (H_2 - H_1) + rho * c2(c1)^2 / 2 * (1 + zeta2(f2(Re2(c2(c1))),K_bocht));
% options for ODE solvers
const.p_a = p_a;
const.rho = rho;
const.g = g;
const.H_0 = H_0;
const.H_1 = H_1;
options1 = odeset('Events', @(t, y) event_function1_m4(t, y, const));
options2 = odeset('Events', @event_function2_m4);
%% Solution part 1
tspan1 = [0 50];
% Define differential equations
dPdt = @(t,p) data_debiet(end)*(1-(p_v./p))/(1-(p_v/p_a)).*-p./V_0;
% Define initial conditions
p0 = p_a;
c0 = 1e-6;
V0 = V_0;
% Solve differential equations until head pressure is reached
systeem = @(t,y) [dPdt(t,y(1))];
[t, y] = ode45(systeem, tspan1, p0, options1);
% Extract the solutions
t_1 = t;
p_1 = y(:, 1);
p1_1 = y(:, 1);
c1_1 = repmat(c0, length(p_1),1);
c2_1 = repmat(c0, length(p_1),1);
V_1 = repmat(V0, length(p_1),1);
%% Solution part 2
tspan2 = [t(end) 200];
% Define new differential equations
dCdt1 = @(t,c1,p) ((p_a - p1(p,c2(c1))) - rho * g * (H_1 - H_0) - rho * c1^2 / 2 * (1 + zeta1(f1(Re1(c1)),K_exp(f1(Re1(c1)))))) / L1;
dVdt = @(t,V,c1) -c1 * S1;
% Define new initial conditions for p0
p0 = p_1(end);
% Solve differential equations after head pressure is reached
systeem = @(t,y) [dPdt(t,y(1)); dCdt1(t,y(2),y(1)); dVdt(t,y(3),y(2))];
[t, y] = ode45(systeem, tspan2, [p0; c0; V0], options2);
% Extract solutions
t_2 = t;
p_2 = y(:, 1);
c1_2 = y(:, 2);
V_2 = y(:, 3);
%% combine results
t = [t_1; t_2(2:end)];
p = [p_1; p_2(2:end)]; % druk in tank
c1 = [c1_1; c1_2(2:end)]; % snelheid door de slang
V = [V_1; V_2(2:end)];

[Torsten]

This cannot work since p0 and p2 determine p1, but p0 is not used in your new system as far as I can see (at least in the description of your new way).

[Wout Suykerbuyk]

My fault... p0 (= p_a) is the atmospheric pressure and is defined in the upper part where I put all the constants (which I didn't show in the above script).

% Given constants

V_0 = 12; % [m³]

p_a = 101325; % [Pa]

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

H_0 = 1; % m

H_1 = 4; % m

H_2 = 6; % m

L1 = 10; % m

L2 = 5; % m

d1 = 0.10; % m

d2 = 0.15; % m

abs_rough1 = 0.0002; % m

abs_rough2 = 0.0002; % m

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

rho = 1000; % kg/m³

S1 = pi*d1^2/4; % m²

S2 = pi*d2^2/4; % m²

p_v = 20000; % Pa

[Torsten]

y0 = [P20,c10];   % Initial values for P2 and c1
[T,Y] = ode45(@fun,[0 10],y0)

function dy = fun(t,y)
  % Given constants
  V_0 = 12; % [m³]
  p_a = 101325; % [Pa]
  g = 9.81; % N/kg of m/s²
  H_0 = 1; % m
  H_1 = 4; % m
  H_2 = 6; % m
  L1 = 10; % m
  L2 = 5; % m
  d1 = 0.10; % m
  d2 = 0.15; % m
  abs_rough1 = 0.0002; % m
  abs_rough2 = 0.0002; % m
  kin_visc = 10^-6; % m²/s
  rho = 1000; % kg/m³
  S1 = pi*d1^2/4; % m²
  S2 = pi*d2^2/4; % m²
  p_v = 20000; % Pa
  zeta1 = @(f1,K_exp) f1 * L1 / d1 + K_exp;
  zeta2 = @(f2,K_bocht) f2 * L2 / d2 + 2 * K_bocht;
  f1 = @(Re1) 0.11 * (abs_rough1/d1 + 68/Re1)^0.25;
  f2 = @(Re2) 0.11 * (abs_rough2/d2 + 68/Re2)^0.25;
  Re1 = @(c1) c1 * d1 / kin_visc;
  Re2 = @(c2) c2 * d2 / kin_visc;
  K_exp = @(f1) (1 + 0.8 * f1) * (1 - (d1/d2)^2)^2;      % sudden widening
  K_bocht = 0.45;                                        % R/D = 1.5 van 90°
  P2 = y(1);
  c1 = y(2);
  c2 = c1*S1/S2;
  P1 = P2 + rho * g * (H_2 - H_1) + rho * c2^2 / 2 * (1 + zeta2(f2(Re2(c2)),K_bocht));
  dP2 = data_debiet(end)*(1-(p_v./P2))/(1-(p_v/p_a)).*-P2./V_0;
  dc1 = ((p_a - P1) - rho * g * (H_1 - H_0) - rho * c1^2 / 2 * (1 + zeta1(f1(Re1(c1)),K_exp(f1(Re1(c1)))))) / L1;
  dy = [dP2;dc1];
end

[Wout Suykerbuyk]

Hey Torsten,

I wanted to let you know my problem is solved. One of the event functions was incorrect defined.

Thanks for your help!