Problem with solving two differential equations depending on each other
11 views (last 30 days)
Show older comments
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)
0 Comments
Answers (1)
Torsten
on 14 Mar 2024
Edited: Torsten
on 14 Mar 2024
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 ?
4 Comments
Torsten
on 15 Mar 2024
Edited: Torsten
on 15 Mar 2024
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.
See Also
Categories
Find more on Sources in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!