Skip to content

Commit 3b7bb38

Browse files
SviatoseEvenSol
Sviatose
and
EvenSol
authored
1078 implement gamma distribution model for plus fraction characterization (#1086)
* started inital implementation * feat: V1 Gamma model * feat: last update * some work for verification of code * added get and set gamma model parameters --------- Co-authored-by: Even Solbraa <[email protected]>
1 parent b08e894 commit 3b7bb38

File tree

6 files changed

+262
-49
lines changed

6 files changed

+262
-49
lines changed

‎src/main/java/neqsim/thermo/characterization/Characterise.java

+1-1
Original file line numberDiff line numberDiff line change
@@ -156,7 +156,7 @@ public void characterisePlusFraction() {
156156
if (plusFractionModel.hasPlusFraction()) {
157157
if (plusFractionModel.getMPlus() > plusFractionModel.getMaxPlusMolarMass()) {
158158
logger.error("plus fraction molar mass too heavy for " + plusFractionModel.getName());
159-
plusFractionModel = plusFractionModelSelector.getModel("heavyOil");
159+
plusFractionModel = plusFractionModelSelector.getModel("Pedersen Heavy Oil");
160160
logger.info("changing to " + plusFractionModel.getName());
161161
}
162162
plusFractionModel.characterizePlusFraction(TBPfractionModel);

‎src/main/java/neqsim/thermo/characterization/LumpingModel.java

+4-3
Original file line numberDiff line numberDiff line change
@@ -306,6 +306,7 @@ public void generateLumpedComposition(Characterise charac) {
306306
system.removeComponent(system.getPhase(0)
307307
.getComponent(charac.getPlusFractionModel().getPlusComponentNumber()).getName());
308308
}
309+
system.init(0);
309310
}
310311

311312
@Override
@@ -398,18 +399,18 @@ public void generateLumpedComposition(Characterise charac) {
398399
pseudoNumber++;
399400
String addName = "C" + Integer.toString(starti) + "-" + Integer.toString(i);
400401
getLumpedComponentNames()[k] = addName;
401-
// System.out.println("adding " + addName);
402+
//System.out.println("adding " + addName);
402403
fractionOfHeavyEnd[k] = zPlus[k] / molFracTot;
403-
404404
system.addTBPfraction(addName, totalNumberOfMoles * zPlus[k], Maverage / zPlus[k],
405405
denstemp1 / denstemp2);
406406
k++;
407-
starti = i + 1;
407+
starti = i;
408408
}
409409
if (charac.getPlusFractionModel().hasPlusFraction()) {
410410
system.removeComponent(system.getPhase(0)
411411
.getComponent(charac.getPlusFractionModel().getPlusComponentNumber()).getName());
412412
}
413+
system.init(0);
413414
}
414415

415416
@Override

‎src/main/java/neqsim/thermo/characterization/PlusFractionModel.java

+181-2
Original file line numberDiff line numberDiff line change
@@ -227,9 +227,8 @@ public void characterizePlusFraction(TBPModelInterface TBPModel) {
227227
z[i] = Math.exp(getCoef(0) + getCoef(1) * i);
228228
M[i] = PVTsimMolarMass[i - 6] / 1000.0;
229229
dens[i] = getCoef(2) + getCoef(3) * Math.log(i);
230-
231-
// System.out.println("z,m,dens " + z[i] + " " + M[i] + " " + dens[i]);
232230
}
231+
// System.out.println("z,m,dens " + z[i] + " " + M[i] + " " + dens[i]);
233232
}
234233

235234
@Override
@@ -287,6 +286,11 @@ public void setCoefs(double[] coefs) {
287286
public void setCoefs(double coef, int i) {
288287
this.coefs[i] = coef;
289288
}
289+
290+
@Override
291+
public void setLastPlusFractionNumber(int fract) {
292+
lastPlusFractionNumber = fract + 1;
293+
}
290294
}
291295

292296
private class PedersenHeavyOilPlusModel extends PedersenPlusModel {
@@ -299,6 +303,179 @@ public PedersenHeavyOilPlusModel() {
299303
}
300304
}
301305

306+
class WhitsonGammaModel extends PedersenPlusModel {
307+
308+
public double[] zValues;
309+
public double[] molarMasses;
310+
public double[] densities;
311+
public double eta = 90; // minimum molecular weight in C7+
312+
public String model = "Whitson";
313+
public double alfa = 1.0;
314+
public double betta = Double.NaN;
315+
316+
public WhitsonGammaModel() {
317+
name = "Whitson Gamma";
318+
}
319+
320+
public void setCalculationModel(String model) {
321+
this.model = model;
322+
}
323+
324+
public void characterizePlusFractionWhitsonGamma() {
325+
326+
}
327+
328+
public double gamma(double X) {
329+
double[] dataB = {-0.577191652, 0.988205891, -0.897056937, 0.918206857, -0.756704078,
330+
0.482199394, -0.193527818, 0.035868343};
331+
double const_ = 1.0;
332+
double XX = X;
333+
if (X < 1.0) {
334+
XX = X + 1.0;
335+
}
336+
while (XX >= 2.0) {
337+
XX -= 1.0;
338+
}
339+
const_ = XX * const_;
340+
XX -= 1.0;
341+
double Y = 1.0;
342+
for (int i = 1; i <= 8; i++) {
343+
Y += dataB[i - 1] * XX * i;
344+
}
345+
double GAMMA = const_ * Y;
346+
if (X < 1.0) {
347+
GAMMA /= X;
348+
}
349+
return GAMMA;
350+
}
351+
352+
public double[] P0P1(double MWB) {
353+
double P0 = 0.0;
354+
double P1 = 0.0;
355+
if (MWB == eta) {
356+
return new double[] {P0, P1};
357+
}
358+
double Y = (MWB - eta) / betta;
359+
double Q = Math.exp(-Y) * Math.pow(Y, alfa) / gamma(alfa);
360+
double TERM = 1.0 / alfa;
361+
double S = TERM;
362+
for (int j = 1; j <= 10000; j++) {
363+
TERM *= Y / (alfa + j);
364+
S += TERM;
365+
if (Math.abs(TERM) <= 1e-8) {
366+
P0 = Q * S;
367+
P1 = Q * (S - 1.0 / alfa);
368+
break;
369+
}
370+
371+
}
372+
return new double[] {P0, P1};
373+
}
374+
375+
public void densityUOP() {
376+
// Calculates density of the C+ function with Watson or Universal Oil Products
377+
// characterization factor
378+
// Experience has shown, that the method is not very accurate for C20+
379+
double Kw = 4.5579 * Math.pow(MPlus * 1000, 0.15178) * Math.pow(densPlus, -1.18241);
380+
for (int i = firstPlusFractionNumber; i < lastPlusFractionNumber; i++) {
381+
densities[i - 1] = 6.0108 * Math.pow(molarMasses[i - 1], 0.17947) * Math.pow(Kw, -1.18241);
382+
}
383+
}
384+
385+
386+
387+
@Override
388+
public void characterizePlusFraction(TBPModelInterface TBPModel) {
389+
system.init(0);
390+
double MWBU = Double.NaN;
391+
double MWBL = Double.NaN;
392+
double sumZ = 0.0;
393+
394+
betta = (MPlus * 1000 - eta) / alfa;
395+
// Implement the Gamma distribution for the plus fraction
396+
zValues = new double[lastPlusFractionNumber];
397+
molarMasses = new double[lastPlusFractionNumber];
398+
densities = new double[lastPlusFractionNumber];
399+
400+
if (model.equals("Whitson")) {
401+
for (int i = firstPlusFractionNumber; i < lastPlusFractionNumber; i++) {
402+
if (i == 1) {
403+
MWBU = eta;
404+
}
405+
MWBL = MWBU;
406+
MWBU = MWBL + 14;
407+
if (i == lastPlusFractionNumber) {
408+
MWBU = 10000.0;
409+
}
410+
double[] P0LP1L = P0P1(MWBL);
411+
double P0L = P0LP1L[0];
412+
double P1L = P0LP1L[1];
413+
414+
double[] P0UP1U = P0P1(MWBU);
415+
double P0U = P0UP1U[0];
416+
double P1U = P0UP1U[1];
417+
418+
double Z = P0U - P0L;
419+
if (Z < 1E-15) {
420+
Z = 1E-15;
421+
}
422+
zValues[i] = Z * zPlus;
423+
double MWAV = eta + alfa * betta * (P1U - P1L) / (P0U - P0L);
424+
molarMasses[i] = MWAV / 1000;
425+
sumZ = sumZ + zValues[i];
426+
}
427+
densityUOP();
428+
429+
430+
}
431+
432+
// Normalize z values to ensure sumZ equals zPlus
433+
for (int i = firstPlusFractionNumber; i < lastPlusFractionNumber; i++) {
434+
zValues[i] *= zPlus / sumZ;
435+
}
436+
}
437+
438+
public void setGammaParameters(double shape, double minMW) {
439+
this.alfa = shape;
440+
this.eta = minMW;
441+
}
442+
443+
public double[] getGammaParameters() {
444+
return new double[] {this.alfa, this.eta};
445+
}
446+
447+
@Override
448+
public double[] getCoefs() {
449+
return new double[] {alfa, eta};
450+
}
451+
452+
@Override
453+
public double getCoef(int i) {
454+
if (i == 0)
455+
return alfa;
456+
if (i == 1)
457+
return eta;
458+
return 0;
459+
}
460+
461+
@Override
462+
public double[] getZ() {
463+
return zValues;
464+
}
465+
466+
@Override
467+
public double[] getM() {
468+
return molarMasses;
469+
}
470+
471+
@Override
472+
public double[] getDens() {
473+
return densities;
474+
}
475+
476+
477+
}
478+
302479
/**
303480
* <p>
304481
* getModel.
@@ -312,6 +489,8 @@ public PlusFractionModelInterface getModel(String name) {
312489
return new PedersenPlusModel();
313490
} else if (name.equals("Pedersen Heavy Oil")) {
314491
return new PedersenHeavyOilPlusModel();
492+
} else if (name.equals("Whitson Gamma Model")) {
493+
return new WhitsonGammaModel();
315494
} else {
316495
return new PedersenPlusModel();
317496
}

‎src/main/java/neqsim/thermo/characterization/PlusFractionModelInterface.java

+2
Original file line numberDiff line numberDiff line change
@@ -72,6 +72,8 @@ public interface PlusFractionModelInterface extends java.io.Serializable {
7272
*/
7373
public int getLastPlusFractionNumber();
7474

75+
public void setLastPlusFractionNumber(int fract);
76+
7577
/**
7678
* <p>
7779
* getPlusComponentNumber.

‎src/test/java/neqsim/thermo/characterization/CharacteriseTest.java

+31-2
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@
44
import org.junit.jupiter.api.Test;
55
import neqsim.thermo.system.SystemInterface;
66
import neqsim.thermo.system.SystemSrkEos;
7+
import neqsim.thermodynamicOperations.ThermodynamicOperations;
78

89
public class CharacteriseTest extends neqsim.NeqSimTest {
910
static SystemInterface thermoSystem = null;
@@ -33,9 +34,37 @@ void testCharacterisePlusFraction() {
3334
thermoSystem.addTBPfraction("C8", 1.0, 120.0 / 1000.0, 0.76);
3435
thermoSystem.addTBPfraction("C9", 1.0, 140.0 / 1000.0, 0.79);
3536
thermoSystem.addPlusFraction("C10", 11.0, 290.0 / 1000.0, 0.82);
36-
thermoSystem.getCharacterization().getLumpingModel().setNumberOfPseudoComponents(12);
37+
//thermoSystem.getCharacterization().getLumpingModel().setNumberOfLumpedComponents(6);
38+
//thermoSystem.getCharacterization().setLumpingModel("PVTlumpingModel");
39+
thermoSystem.getCharacterization().setLumpingModel("no lumping");
40+
thermoSystem.getCharacterization().characterisePlusFraction();
41+
// assertEquals(15, thermoSystem.getNumberOfComponents());
42+
thermoSystem.prettyPrint();
43+
}
44+
45+
46+
@Test
47+
void testCharacterisePlusFractionGAMMA() {
48+
thermoSystem = new SystemSrkEos(298.0, 10.0);
49+
thermoSystem.addComponent("methane", 51.0);
50+
thermoSystem.addPlusFraction("C10", 11.0, 290.0 / 1000.0, 0.82);
51+
thermoSystem.getCharacterization().setPlusFractionModel("Whitson Gamma Model");
3752
thermoSystem.getCharacterization().setLumpingModel("PVTlumpingModel");
53+
thermoSystem.getCharacterization().getLumpingModel().setNumberOfLumpedComponents(15);
3854
thermoSystem.getCharacterization().characterisePlusFraction();
39-
assertEquals(15, thermoSystem.getNumberOfComponents());
55+
// logger.info("number of components " + thermoSystem.getNumberOfComponents());
56+
//assertEquals(86, thermoSystem.getNumberOfComponents());
57+
//System.out.println(thermoSystem.getComponent("C1-2_PC").getz());
58+
thermoSystem.prettyPrint();
59+
60+
thermoSystem.setPressure(1, "bara");
61+
62+
63+
ThermodynamicOperations thermoOps = new ThermodynamicOperations(thermoSystem);
64+
thermoOps.TPflash();
65+
66+
thermoSystem.initProperties();
67+
thermoSystem.prettyPrint();
4068
}
69+
4170
}

0 commit comments

Comments
 (0)