Filter löschen
Filter löschen

How do I solve Differential Algebraic Equations?

4 Ansichten (letzte 30 Tage)
ojonugwa adukwu
ojonugwa adukwu am 28 Okt. 2019
Kommentiert: ojonugwa adukwu am 31 Okt. 2019
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 Kommentare
Fabio Freschi
Fabio Freschi am 31 Okt. 2019
Bearbeitet: Fabio Freschi am 31 Okt. 2019
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
You maybe be interested in using other solvers than ode23t.
ojonugwa adukwu
ojonugwa adukwu am 31 Okt. 2019
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

Melden Sie sich an, um zu kommentieren.

Antworten (0)

Kategorien

Mehr zu Programming finden Sie in Help Center und File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by