Solvin Differential Algebraic Equations

Question

Patience Shamaki 2019-3-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/452777-solvin-differential-algebraic-equations

编辑： Torsten 2022-1-25

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 个评论
显示 2更早的评论隐藏 2更早的评论

Torsten 2019-3-27

Then please let me know which of your 23 equations read "0 = out(i)" and which read "dy(i)/dt = out(i)".

Patience Shamaki 2019-3-27

Differential equations: (dx(1)/dt = wgl - wiv; dx(2)/dt = wiv - wpg+ wrg .....)

dy(1)/dt= wgl - y(15);

dy(2)/dt= y(15) - y(19) + y(18);

dy(3)/dt= y(20) - y(21);

dy(4)/dt= y(19) - y(22);

dy(5)/dt= y(21) - y(23);

0= ((((Ta*R)/(Va*Mw))+((g*La)/(La*Aa)))*y(1)) - y(6);

0= (Mw*y(6)) -(Ta*R*y(9));

0=(((((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));

0= y(2)+ y(3)-(rho_o*Lbh*Abh)-(Lw*Aw*y(10));

0= (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);

0= (y(8)+ (y(10)*g*Hbh)+((128*mu_o*Lbh*y(20))/(pi*Dbh^4*rho_o)))-y(12);

0= (Tr*R*y(4))-(Mw*Lr*Ar*y(13));

0= y(4)+ y(5)-(Lr*Ar*y(11));

0= (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));

0= (Civ*sqrt(y(9)*(y(6)-y(8))))-y(15);

0= (Cpc*sqrt(y(10)*(y(7)-y(14))))-y(16);

0= (Crh*sqrt(y(11)*(y(6)-Ps)))- y(17);

0= (PI*(Pr -y(12)))- y(20);

0= GOR* y(20) -y(18);

0= (y(2)*y(16)) - (y(2)*y(19)) - (y(3)*y(19));

0= (y(3)*y(16)) - (y(2)*y(21)) - (y(3)*y(21));

0= (y(4)*y(17)) - (y(4)*y(22)) - (y(4)*y(22));

0=(y(5)*y(17)) - (y(5)*y(23)) + (y(5)*y(23));

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2019-3-28

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/452777-solvin-differential-algebraic-equations#answer_367853

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