How do I solve Differential Algebraic Equations?

Question

ojonugwa adukwu 2019-10-28

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/487859-how-do-i-solve-differential-algebraic-equations

评论： ojonugwa adukwu 2019-10-31

How do I solve my DAE systems in gas lift system? My codes are giving results that are not acceptable. For example, my densities sometimes give negative values.

Thanks so much.

The code for my function is:

mga = x(1); %mass of gas in the annulus

mgt = x(2); %mass of gas in the well tubing

mot = x(3); %mass of oil in the well tubing

Pa = x(4); % pressure of gas in the annulus

rho_a = x(5); % density of gas in the annulus

Pwh = x(6); % Wellhead pressure

rho_w = x(7); % density of mixture in the well

wpc = x(8);

Pw = x(9); % Well pressure

Pbh = x(10); % bottomhole pressure

wiv = x(11);

wpg = x(12);

wpo = x(13);

wro = x(14);

wrg = x(15);

Lbh = 500;

La = 1500;

Hw=1000;

Dw=0.121;

Da=0.189;

Aw = 0.25*pi*Dw^2;

Aa=0.25*pi*(Da^2-Dw^2);

Va=Aa*La;

Abh = Aw;

uk=01; us=0.4;

rho_o = 800;

Civ = 1*1e-4; Cpc = 2e-3;Crh= 10^-2; Cr = 2.623*10^-4;

PI=0.7;

Pr = 150; Ps= 20; tf=300;

Ta = 301; Tw = 305; Mw = 0.02; R = 8.314;

g = 9.8;wgl = us*1; GOR = 0.1;

x(1,1) = wgl - wiv;

x(2,1) = wiv - wpg+ wrg;

x(3,1) = wro - wpo;

x(4,1) = -Pa*1e5+((Ta*R/(Va*Mw)+g*La/Va)*mga*1e3) ;

x(5,1) = -rho_a*10^3 +((Mw/(Ta*R))*Pa);

x(6,1) = -Pwh*1e5 + ((Tw*R/Mw)*(mgt*1e3/(Lw*Aw+Lbh*Abh-(mot*1e3/rho_o))))-(0.5*((mgt*1e3 + mot*1e3)*g*Hw/(Lw*Aw)));%

x(7,1) = -rho_w*10^3 +((mgt*1e3+ mot*1e3)*Pwh*Mw*rho_o)/(mot*1e3*Pwh*Mw +rho_o*R*Tw*mgt*1e3);

x(8,1) = -wpc +Cpc*sqrt(rho_w*(max(0.001,(Pwh - Ps*1e5))));

x(9,1) = -Pw*1e5 + Pwh +((g*Hw/(Lw*Aw))*max(0.001,(mot*1e3 + mgt*1e3 -(rho_o*Lbh*Abh))))+(128*0.003*Lw*wpc/(pi*Dw^4*((mgt*1e3+ mot*1e3)*Pwh*Mw*rho_o)/(mot*1e3*Pwh*Mw +rho_o*R*Tw*mgt*1e3)));

x(10,1) = -Pbh*1e5 + Pw + (rho_o*g*Lbh)+ (128*0.003*wro*Lbh/(3.14*Dw^4*rho_o));

x(11,1) = -wiv +Civ*sqrt(rho_a*(max (0.001,(Pa-Pw))));

x(12,1) = -wpg + ((mgt*1e3/max(0.001,(mgt*1e3+mot*1e3)))*wpc);

x(13,1) = -wpo + ((mot*1e3/max(0.001,(mgt*1e3+mot*1e3)))*wpc);

x(14,1) = -wro + 0.7*10^-5*(Pr*1e5-Pbh);

x(15,1) = -wrg + (GOR*wro);

The corresponding code for the script is

tspan = 0:36000;

x0 =[1.3268; 0.8093; 6.2694; 74.70300; 59.7027;

47.57000; 240.7129; 51.5223 ; 63.48600; 102.69000;

0.8184; 5.8905; 45.6318; 33.1200; 3.3120];

M = diag([ones(1,3) zeros(1,12)]);

options = odeset('Mass',M,'RelTol',1e-3,'AbsTol',1e-5);

[t,x] = ode23t(@gasLiftFnReal2,tspan,x0,options);

subplot(2,2,1)

plot(t/3600, 1000*x(:,5), '--')

xlabel('tspan')

ylabel('density of gas in annulus')

subplot(2,2,2)

plot(t/3600,1000*x(:,5), '-')

xlabel('tspan')

ylabel('density of gas in annulus')

subplot(2,2,3)

plot(t/3600, 1000*x(:,7), '-')

xlabel('tspan')

ylabel('density of mixtures in tubing')

subplot(2,2,4)

plot(t/3600, 1000*x(:,7), '-')

xlabel('tspan')

ylabel('density of mixtures in tubing')

% title

title ('Gas-lift Network')

6 个评论
显示 4更早的评论隐藏 4更早的评论

ojonugwa adukwu 2019-10-28

Thanks Fabio.

Neither Text, code, insert nor help is active when i paste the code. I have added the Lw. I have also added the function name before the function so it can be copied and pasted to MATLAB script.

The code is intended to simulate in open loop the nonlinear gas lift system with the details provided. I have gone through everything but it appears I still can't see where the problem is.

The system is starting close to it's steady state and I expected that it settles not too far away. In particular, densities, flowrates and masses can not be negative but these are what I am seeing.

The first three variable are differential while the last 12 are algebraic. All variable are decribed briefly.

Thanks so much

The Function is;

function x = gasLiftFnReal2(t,x)

mga = x(1); %mass of gas in the annulus

mgt = x(2); %mass of gas in the well tubing

mot = x(3); %mass of oil in the well tubing

Pa = x(4); % pressure of gas in the annulus

rho_a = x(5); % density of gas in the annulus

Pwh = x(6); % Wellhead pressure

rho_w = x(7); % density of mixture in the well

wpc = x(8); %flow rate through choke

Pw = x(9); % Well pressure

Pbh = x(10); % bottomhole pressure

wiv = x(11); flow thyrough

wpg = x(12); produced gas flow rate

wpo = x(13); produced oil flow rate

wro = x(14); flow rate of oil from the reservoir into bottomhole

wrg = x(15); flow rate of gas from the reservoir into bottomhole

Lbh = 500;

Lw=Lbh;

La = 1500;

Hw=1000;

Dw=0.121;

Da=0.189;

Aw = 0.25*pi*Dw^2;

Aa=0.25*pi*(Da^2-Dw^2);

Va=Aa*La;

Abh = Aw;

rho_o = 800;

Civ = 1*1e-4; Cpc = 2e-3;Crh= 10^-2; Cr = 2.623*10^-4;

PI=0.7;

Pr = 150; Ps= 20; tf=300;

Ta = 301; Tw = 305; Mw = 0.02; R = 8.314;

g = 9.8;wgl = 1; GOR = 0.1;

x(1,1) = wgl - wiv;

x(2,1) = wiv - wpg+ wrg;

x(3,1) = wro - wpo;

x(4,1) = -Pa*1e5+((Ta*R/(Va*Mw)+g*La/Va)*mga*1e3) ;

x(5,1) = -rho_a*10^3 +((Mw/(Ta*R))*Pa);

x(6,1) = -Pwh*1e5 + ((Tw*R/Mw)*(mgt*1e3/(Lw*Aw+Lbh*Abh-(mot*1e3/rho_o))))-(0.5*((mgt*1e3 + mot*1e3)*g*Hw/(Lw*Aw)));

x(7,1) = -rho_w*10^3 +((mgt*1e3+ mot*1e3)*Pwh*Mw*rho_o)/(mot*1e3*Pwh*Mw +rho_o*R*Tw*mgt*1e3);

x(8,1) = -wpc +Cpc*sqrt(rho_w*(max(0.001,(Pwh - Ps*1e5))));

x(9,1) = -Pw*1e5 + Pwh +((g*Hw/(Lw*Aw))*max(0.001,(mot*1e3 + mgt*1e3 -(rho_o*Lbh*Abh))))+(128*0.003*Lw*wpc/(pi*Dw^4*((mgt*1e3+ mot*1e3)*Pwh*Mw*rho_o)/(mot*1e3*Pwh*Mw +rho_o*R*Tw*mgt*1e3)));

x(10,1) = -Pbh*1e5 + Pw + (rho_o*g*Lbh)+ (128*0.003*wro*Lbh/(3.14*Dw^4*rho_o));

x(11,1) = -wiv +Civ*sqrt(rho_a*(max (0.001,(Pa-Pw))));

x(12,1) = -wpg + ((mgt*1e3/max(0.001,(mgt*1e3+mot*1e3)))*wpc);

x(13,1) = -wpo + ((mot*1e3/max(0.001,(mgt*1e3+mot*1e3)))*wpc);

x(14,1) = -wro + 0.7*10^-5*(Pr*1e5-Pbh);

x(15,1) = -wrg + (GOR*wro);

The Script is;

tspan = 0:36000;

x0 =[1.3268; 0.8093; 6.2694; 74.70300; 59.7027;

47.57000; 240.7129; 51.5223 ; 63.48600; 102.69000;

0.8184; 5.8905; 45.6318; 33.1200; 3.3120];

M = diag([ones(1,3) zeros(1,12)]);

options = odeset('Mass',M,'RelTol',1e-3,'AbsTol',1e-5);

[t,x] = ode23t(@gasLiftFnReal2,tspan,x0,options);

subplot(2,2,1)

plot(t/3600, 1000*x(:,5), '--')

xlabel('tspan')

ylabel('density of gas in annulus')

subplot(2,2,2)

plot(t/3600,1000*x(:,5), '-')

xlabel('tspan')

ylabel('density of gas in annulus')

subplot(2,2,3)

plot(t/3600, 1000*x(:,7), '-')

xlabel('tspan')

ylabel('density of mixtures in tubing')

subplot(2,2,4)

plot(t/3600, 1000*x(:,7), '-')

xlabel('tspan')

ylabel('density of mixtures in tubing')

Fabio Freschi 2019-10-31

编辑：Fabio Freschi 2019-10-31

Do not undersetimate the power of numerical approximation: you could have negative values because of the local and global discretization error of the ODE: look the example here

https://it.mathworks.com/help/matlab/math/nonnegative-ode-solution.html

You maybe be interested in using other solvers than ode23t.

ojonugwa adukwu 2019-10-31

Many thanks Fabio. It is indeed appreciated.

I observed that this nonnegativity options is not available to DAE systems which I am working with. I have been able to figure out what happened now though. Everything works well now. The problems were with some values provided, proper use of commas and brackets.

Thanks so much bro

请先登录，再进行评论。

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