How do I solve Differential Algebraic Equations?

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

[Fabio Freschi]

Difficult to say

1) the code is not properly formatted: use the code button

2) Lw is not defined, as well as the function gasLiftFnReal2

please provide a minimum working example to help people in the community helping you

[ojonugwa adukwu]

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]

One ODE option is to set a specific component to be nonnegative. From the documentation:

'NonNegative' — Nonnegative solution components

[] (default) | scalar | vector

Nonnegative solution components, specified as the comma-separated pair consisting of 'NonNegative' and a scalar or vector. The scalar or vector selects which solution components must be nonnegative.

Note

NonNegative is not available for ode23s or ode15i. Additionally, for ode15s, ode23t, and ode23tb it is not available for problems where there is a mass matrix.

Example: opts = odeset('NonNegative',1) specifies that the first solution component must be nonnegative.

You may have a look at this option to see if it fixes your problem

[ojonugwa adukwu]

Many thanks Fabio.

I think they should naturally produce nonnegative values if the code is properly debugged, considering order of various lines, intial conditions, commas, etc. ODE 23t which i am using can solve DAE and it won's accept 'nonNegativity'.

I am still researching to see where the problem is.

Thanks so much.

[Fabio Freschi]

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]

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