Solvin Differential Algebraic Equations
4 views (last 30 days)
Show older comments
Please how can I solve a DAE for a gas lift system network. I have 5 differential equations, and 18 algebraic equations. I tried following the codes to create mine but it is not running. Here is a sample:
function out = gas_lift_ftn(t,y)
%Define other parameters%
mu_o = 0.01;
pi = 3.14;
Lw = 1500;
Hw = 1000;
Dw = 0.121;
Lbh = 500;
Hbh = Lbh;
Dbh = Dw;
La = Lw;
Ha = Hw;
Da = 0.189;
Lr = Lbh;
Hr = Lbh;
Dr = Dw;
Aw = (pi/4)*Dw^2;
Abh = Aw;
Ar = Aw;
Aa = (pi/4)*Da;
Va = La*pi/4*(Da^2-Dw^2);
rho_o = 800;
Civ = 0.1e-03; Crh = 10e-3; Cpc = 2e-3;
Pr = 150; Ps = 20; PI = 0.7;
Ta = 28+273; Tw = 32+273; Tr = 30+273; Mw = 20; R = 0.083145;
g = 9.8*10e-5; GOR = 0.1; wgl = 2.217;
out = [ wgl - y(15)
y(15) - y(19) + y(18)
y(20) - y(21)
y(19) - y(22)
y(21) - y(23)
((((Ta*R)/(Va*Mw))+((g*La)/(La*Aa)))*y(1)) - y(6)
(Mw*y(6)) -(Ta*R*y(9))
(((((Tw*R)/Mw))*((y(2))/((Lw*Aw)+(Lbh*Abh)-(y(3)/rho_o)))- 0.5*((y(2)+(y(3))/(Lw*Aw))*g*Hw))-y(7))
y(2)+ y(3)-(rho_o*Lbh*Abh)-(Lw*Aw*y(10))
(y(7)+((g/(Lw*Aw))*(y(3)+y(2)-(rho_o*Lbh*Abh)*Hw))+((128*mu_o*Lw*y(16))/(pi*Dw^4*((y(2) + y(3))*y(7)*Mw*rho_o)/(y(3)*y(7)*Mw + rho_o*R*Tw*y(2)))))-y(8)
(y(8)+ (y(10)*g*Hbh)+((128*mu_o*Lbh*y(20))/(pi*Dbh^4*rho_o)))-y(12)
(Tr*R*y(4))-(Mw*Lr*Ar*y(13))
y(4)+ y(5)-(Lr*Ar*y(11))
(y(13) + y(11)*g*Hr + ((128*mu_o*Lr*y(17))/(pi*Dr^4*((y(4) + y(5))*y(13)*Mw*y(11))/(y(5)*y(13)*Mw + y(11)*R*Tr*y(4)))) - y(14))
(Civ*sqrt(y(9)*(y(6)-y(8))))-y(15)
(Cpc*sqrt(y(10)*(y(7)-y(14))))-y(16)
(Crh*sqrt(y(11)*(y(6)-Ps)))- y(17)
(PI*(Pr -y(12)))- y(20)
GOR* y(20) -y(18)
(y(2)*y(16)) - (y(2)*y(19)) - (y(3)*y(19))
(y(3)*y(16)) - (y(2)*y(21)) - (y(3)*y(21))
(y(4)*y(17)) - (y(4)*y(22)) - (y(4)*y(22))
(y(5)*y(17)) - (y(5)*y(23)) + (y(5)*y(23))];
end
%gaslift systems network
tspan = [0:1:60];
% The script for gas lift network model defining initial conditions and
% boundaries
y0 =[2812.3; 1269.75; 9200.4; 136.9; 2300.2; 160; 80.52; 113.87; 127.9;
340.29; 428.8; 130.54; 30; 50.77; 1.853;63.65; 7.71; 55.93; 13.62; 1.36;
205.86; 11.56; 194.3];
% call
[t,y] = ode15s(@gas_lift_ftn,tspan,y0);
disp('y(1), y(2), y(3), y(4), y(5)')
disp(y(5:1))
4 Comments
Torsten
on 27 Mar 2019
Then please let me know which of your 23 equations read "0 = out(i)" and which read "dy(i)/dt = out(i)".
Accepted Answer
Torsten
on 28 Mar 2019
function main
%gaslift systems network
tspan = [0:1:60];
% The script for gas lift network model defining initial conditions and
% boundaries
y0 =[2812.3; 1269.75; 9200.4; 136.9; 2300.2; 160; 80.52; 113.87; 127.9;
340.29; 428.8; 130.54; 30; 50.77; 1.853;63.65; 7.71; 55.93; 13.62; 1.36;
205.86; 11.56; 194.3];
M = diag([ones(1,5) zeros(1,18)]);
options = odeset('Mass',M);
% call
[t,y] = ode15s(@gas_lift_ftn,tspan,y0,options);
end
function out = gas_lift_ftn(t,y)
%Define other parameters%
mu_o = 0.01;
pi = 3.14;
Lw = 1500;
Hw = 1000;
Dw = 0.121;
Lbh = 500;
Hbh = Lbh;
Dbh = Dw;
La = Lw;
Ha = Hw;
Da = 0.189;
Lr = Lbh;
Hr = Lbh;
Dr = Dw;
Aw = (pi/4)*Dw^2;
Abh = Aw;
Ar = Aw;
Aa = (pi/4)*Da;
Va = La*pi/4*(Da^2-Dw^2);
rho_o = 800;
Civ = 0.1e-03; Crh = 10e-3; Cpc = 2e-3;
Pr = 150; Ps = 20; PI = 0.7;
Ta = 28+273; Tw = 32+273; Tr = 30+273; Mw = 20; R = 0.083145;
g = 9.8*10e-5; GOR = 0.1; wgl = 2.217;
out=zeros(23,1);
out(1) = wgl - y(15);
out(2) = y(15) - y(19) + y(18);
out(3) = y(20) - y(21);
out(4) = y(19) - y(22);
out(5) = y(21) - y(23);
out(6) = ((((Ta*R)/(Va*Mw))+((g*La)/(La*Aa)))*y(1)) - y(6);
out(7) = (Mw*y(6)) -(Ta*R*y(9));
out(8) = (((((Tw*R)/Mw))*((y(2))/((Lw*Aw)+(Lbh*Abh)-(y(3)/rho_o)))- 0.5*((y(2)+(y(3))/(Lw*Aw))*g*Hw))-y(7));
out(9) = y(2)+ y(3)-(rho_o*Lbh*Abh)-(Lw*Aw*y(10));
out(10) = (y(7)+((g/(Lw*Aw))*(y(3)+y(2)-(rho_o*Lbh*Abh)*Hw))+((128*mu_o*Lw*y(16))/(pi*Dw^4*((y(2) + y(3))*y(7)*Mw*rho_o)/(y(3)*y(7)*Mw + rho_o*R*Tw*y(2)))))-y(8);
out(11) = (y(8)+ (y(10)*g*Hbh)+((128*mu_o*Lbh*y(20))/(pi*Dbh^4*rho_o)))-y(12);
out(12) = (Tr*R*y(4))-(Mw*Lr*Ar*y(13));
out(13) = y(4)+ y(5)-(Lr*Ar*y(11));
out(14) = (y(13) + y(11)*g*Hr + ((128*mu_o*Lr*y(17))/(pi*Dr^4*((y(4) + y(5))*y(13)*Mw*y(11))/(y(5)*y(13)*Mw + y(11)*R*Tr*y(4)))) - y(14));
out(15) = (Civ*sqrt(y(9)*(y(6)-y(8))))-y(15);
out(16) = (Cpc*sqrt(y(10)*(y(7)-y(14))))-y(16);
out(17) = (Crh*sqrt(y(11)*(y(6)-Ps)))- y(17);
out(18) = (PI*(Pr -y(12)))- y(20);
out(19) = GOR* y(20) -y(18);
out(20) = (y(2)*y(16)) - (y(2)*y(19)) - (y(3)*y(19));
out(21) = (y(3)*y(16)) - (y(2)*y(21)) - (y(3)*y(21))
out(22) = (y(4)*y(17)) - (y(4)*y(22)) - (y(4)*y(22))
out(23) = (y(5)*y(17)) - (y(5)*y(23)) + (y(5)*y(23))];
end
5 Comments
Jaffar Ali Lone
on 25 Jan 2022
@Torsten How to solve discrete time descriptor system of the following type
Ex(k) =Ax(k) + Bu(k) where E is Singular matrix
For example
E = [1 0 0 0 0 0;0 1 0 0 0 0;0 0 1 0 0 0;0 0 0 1 0 0;0 0 0 0 0 0;0 0 0 0 0 0];
A = [0 0 0 0 0.43478 0;0 -0.16666 0 0 0.000666 0;0 0 0 0 0 0.5;0 0 0 -0.14285 0 0.0005;22.9676 -21.9676 0 -1 0.0025 -0.0015;0 0 0 0 1 1];
C = [22.9676 1 0 0 0.0025 0;0 22.9676 0 1 0 0.0015];
B = [0;0;0;0;1;1];
u = sin(t)
More Answers (0)
See Also
Categories
Find more on Systems of Nonlinear Equations 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!