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);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); 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); ode45, fluid dynamics MATLAB Answers — New Questions
