MATLAB: Solvin Differential Algebraic Equations

differential algebraic equation

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))

Best Answer

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
Related Question