Skip to content

Commit 19251bb

Browse files
EvenSolEvenSol
authored
start adding more fluid property models (#1075)
* start adding more fluid property models * add java test to sdcode * further implementation * futher updates on Twu model * updated test * added PC, TC and VC
1 parent e9eacdb commit 19251bb

File tree

4 files changed

+237
-6
lines changed

4 files changed

+237
-6
lines changed

‎.devcontainer/devcontainer.json

+2-1
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,8 @@
2424
"ms-vscode.test-adapter-converter",
2525
"shengchen.vscode-checkstyle",
2626
"mechatroner.rainbow-csv",
27-
"redhat.vscode-xml"
27+
"redhat.vscode-xml",
28+
"vscjava.vscode-java-test"
2829
],
2930
"settings": {
3031
"java.configuration.updateBuildConfiguration": "interactive",

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

+177-4
Original file line numberDiff line numberDiff line change
@@ -104,6 +104,16 @@ public double calcCriticalViscosity(double molarMass, double density) {
104104
* 1e-7;
105105
}
106106

107+
@Override
108+
public double calcRacketZ(SystemInterface thermoSystem, double molarMass, double density) {
109+
throw new RuntimeException("calcm() method not defined");
110+
}
111+
112+
@Override
113+
public double calcm(double molarMass, double density) {
114+
throw new RuntimeException("calcm() method not defined");
115+
}
116+
107117
@Override
108118
public boolean isCalcm() {
109119
return calcm;
@@ -154,13 +164,11 @@ public double calcPC(double molarMass, double density) {
154164
public double calcm(double molarMass, double density) {
155165
if (molarMass < 1120) {
156166
TBPfractionCoefs = TBPfractionCoefOil;
157-
return TBPfractionCoefs[2][0] + TBPfractionCoefs[2][1] * molarMass
158-
+ TBPfractionCoefs[2][2] * density + TBPfractionCoefs[2][3] * Math.pow(molarMass, 2.0);
159167
} else {
160168
TBPfractionCoefs = TBPfractionCoefsHeavyOil;
161-
return TBPfractionCoefs[2][0] + TBPfractionCoefs[2][1] * Math.log(molarMass)
162-
+ TBPfractionCoefs[2][2] * density + TBPfractionCoefs[2][3] * Math.sqrt(molarMass);
163169
}
170+
return TBPfractionCoefs[2][0] + TBPfractionCoefs[2][1] * molarMass
171+
+ TBPfractionCoefs[2][2] * density + TBPfractionCoefs[2][3] * Math.pow(molarMass, 2.0);
164172
}
165173

166174
@Override
@@ -312,6 +320,167 @@ public double calcAcentricFactor(double molarMass, double density) {
312320
}
313321
}
314322

323+
/**
324+
* Lee-Kesler property estimation method
325+
*/
326+
public class LeeKesler extends TBPBaseModel {
327+
private static final long serialVersionUID = 1000;
328+
329+
public LeeKesler() {
330+
calcm = false;
331+
}
332+
333+
@Override
334+
public double calcTC(double molarMass, double density) {
335+
double sg = density;
336+
double TB = calcTB(molarMass, density);
337+
double TC =
338+
189.8 + 450.6 * sg + (0.4244 + 0.1174 * sg) * TB + (0.1441 - 1.0069 * sg) * 1e5 / TB;
339+
return TC;
340+
}
341+
342+
@Override
343+
public double calcPC(double molarMass, double density) {
344+
double sg = density;
345+
double TB = calcTB(molarMass, density);
346+
double logpc =
347+
3.3864 - 0.0566 / sg - ((0.43639 + 4.1216 / sg + 0.21343 / sg / sg) * 1e-3 * TB)
348+
+ ((0.47579 + 1.182 / sg + 0.15302 / sg / sg) * 1e-6 * TB * TB)
349+
- ((2.4505 + 9.9099 / sg / sg) * 1e-10 * TB * TB * TB);
350+
double PC = Math.exp(logpc) * 10;
351+
return PC;
352+
}
353+
354+
public double calcAcentricFactor(double molarMass, double density) {
355+
return super.calcAcentricFactorKeslerLee(molarMass, density);
356+
}
357+
}
358+
359+
/**
360+
* Two property estimation method
361+
*/
362+
public class TwuModel extends TBPBaseModel {
363+
private static final long serialVersionUID = 1000;
364+
365+
public TwuModel() {
366+
calcm = false;
367+
}
368+
369+
@Override
370+
public double calcTC(double molarMass, double density) {
371+
double sg = density;
372+
double TB = calcTB(molarMass, density);
373+
double MW = solveMW(TB);
374+
double Tcnalkane = TB * 1.0 / (0.533272 + 0.343831e-3 * TB + 2.526167e-7 * TB * TB
375+
- 1.65848e-10 * TB * TB * TB + 4.60774e24 * Math.pow(TB, -13));
376+
double phi = 1.0 - TB / Tcnalkane;
377+
double SGalkane =
378+
0.843593 - 0.128624 * phi - 3.36159 * Math.pow(phi, 3) - 13749 * Math.pow(phi, 12);
379+
double PCnalkane = Math.pow(0.318317 + 0.099334 * Math.sqrt(phi) + 2.89698 * phi
380+
+ 3.0054 * phi * phi + 8.65163 * Math.pow(phi, 4), 2);
381+
double VCnalkane = Math.pow(
382+
(0.82055 + 0.715468 * phi + 2.21266 * phi * phi * phi + 13411.1 * Math.pow(phi, 14)), -8);
383+
double deltaST = Math.exp(5.0 * (SGalkane - sg)) - 1.0;
384+
double fT = deltaST * (-0.270159 * Math.pow(TB, -0.5)
385+
+ (0.0398285 - 0.706691 * Math.pow(TB, -0.5) * deltaST));
386+
double TC = Tcnalkane * Math.pow(((1 + 2 * fT) / (1 - 2 * fT)), 2);
387+
return TC;
388+
}
389+
390+
public double calculateTfunc(double MW_alkane, double TB) {
391+
double phi = Math.log(MW_alkane);
392+
return Math
393+
.exp(5.1264 + 2.71579 * phi - 0.28659 * phi * phi - 39.8544 / phi - 0.122488 / phi / phi)
394+
- 13.7512 * phi + 19.6197 * phi * phi - TB;
395+
}
396+
397+
public double computeGradient(double MW_alkane, double TB) {
398+
double delta = 1;
399+
double TfuncPlus = calculateTfunc(MW_alkane + delta, TB);
400+
double TfuncMinus = calculateTfunc(MW_alkane - delta, TB);
401+
return (TfuncPlus - TfuncMinus) / (2 * delta);
402+
}
403+
404+
public double solveMW(double TB) {
405+
double MW_alkane = TB / (5.8 - 0.0052 * TB);
406+
double tolerance = 1e-6;
407+
double prevMW_alkane;
408+
double error = 1.0;
409+
int iter = 0;
410+
411+
do {
412+
iter++;
413+
prevMW_alkane = MW_alkane;
414+
double gradient = computeGradient(MW_alkane, TB);
415+
MW_alkane -= 0.5 * calculateTfunc(MW_alkane, TB) / gradient;
416+
error = Math.abs(MW_alkane - prevMW_alkane);
417+
} while (Math.abs(error) > tolerance && iter < 1000 || iter < 3);
418+
419+
return MW_alkane;
420+
}
421+
422+
@Override
423+
public double calcPC(double molarMass, double density) {
424+
double sg = density;
425+
double TB = calcTB(molarMass, density);
426+
double MW = solveMW(TB);
427+
double Tcnalkane = TB * 1.0 / (0.533272 + 0.343831e-3 * TB + 2.526167e-7 * TB * TB
428+
- 1.65848e-10 * TB * TB * TB + 4.60774e24 * Math.pow(TB, -13));
429+
double phi = 1.0 - TB / Tcnalkane;
430+
double SGalkane =
431+
0.843593 - 0.128624 * phi - 3.36159 * Math.pow(phi, 3) - 13749 * Math.pow(phi, 12);
432+
double PCnalkane = Math.pow(0.318317 + 0.099334 * Math.sqrt(phi) + 2.89698 * phi
433+
+ 3.0054 * phi * phi + 8.65163 * Math.pow(phi, 4), 2);
434+
double VCnalkane = Math.pow(
435+
(0.82055 + 0.715468 * phi + 2.21266 * phi * phi * phi + 13411.1 * Math.pow(phi, 14)), -8);
436+
double deltaST = Math.exp(5.0 * (SGalkane - sg)) - 1.0;
437+
double fT = deltaST * (-0.270159 * Math.pow(TB, -0.5)
438+
+ (0.0398285 - 0.706691 * Math.pow(TB, -0.5) * deltaST));
439+
double TC = Tcnalkane * Math.pow(((1 + 2 * fT) / (1 - 2 * fT)), 2);
440+
double deltaSP = Math.exp(0.5 * (SGalkane - sg)) - 1.0;
441+
double deltaSV = Math.exp(4.0 * (SGalkane * SGalkane - sg * sg)) - 1.0;
442+
double fV = deltaSV
443+
* (0.347776 * Math.pow(TB, -0.5) + (-0.182421 + 2.24890 * Math.pow(TB, -0.5)) * deltaSV);
444+
double VC = VCnalkane * Math.pow(((1 + 2 * fV) / (1 - 2 * fV)), 2);
445+
double fP = deltaSP * ((2.53262 - 34.4321 * Math.pow(TB, -0.5) - 0.00230193 * TB)
446+
+ (-11.4277 + 187.934 * Math.pow(TB, -0.5) + 0.00414963 * TB) * deltaSP);
447+
double PC = PCnalkane * (TC / Tcnalkane) * (VCnalkane / VC)
448+
* Math.pow(((1 + 2 * fP) / (1 - 2 * fP)), 2);
449+
return PC * 10.0; // * 10 due to conversion MPa to bar
450+
}
451+
452+
@Override
453+
public double calcCriticalVolume(double molarMass, double density) {
454+
double sg = density;
455+
double TB = calcTB(molarMass, density);
456+
double MW = solveMW(TB);
457+
double Tcnalkane = TB * 1.0 / (0.533272 + 0.343831e-3 * TB + 2.526167e-7 * TB * TB
458+
- 1.65848e-10 * TB * TB * TB + 4.60774e24 * Math.pow(TB, -13));
459+
double phi = 1.0 - TB / Tcnalkane;
460+
double SGalkane =
461+
0.843593 - 0.128624 * phi - 3.36159 * Math.pow(phi, 3) - 13749 * Math.pow(phi, 12);
462+
double PCnalkane = Math.pow(0.318317 + 0.099334 * Math.sqrt(phi) + 2.89698 * phi
463+
+ 3.0054 * phi * phi + 8.65163 * Math.pow(phi, 4), 2);
464+
double VCnalkane = Math.pow(
465+
(0.82055 + 0.715468 * phi + 2.21266 * phi * phi * phi + 13411.1 * Math.pow(phi, 14)), -8);
466+
double deltaST = Math.exp(5.0 * (SGalkane - sg)) - 1.0;
467+
double fT = deltaST * (-0.270159 * Math.pow(TB, -0.5)
468+
+ (0.0398285 - 0.706691 * Math.pow(TB, -0.5) * deltaST));
469+
double TC = Tcnalkane * Math.pow(((1 + 2 * fT) / (1 - 2 * fT)), 2);
470+
double deltaSP = Math.exp(0.5 * (SGalkane - sg)) - 1.0;
471+
double deltaSV = Math.exp(4.0 * (SGalkane * SGalkane - sg * sg)) - 1.0;
472+
double fV = deltaSV
473+
* (0.347776 * Math.pow(TB, -0.5) + (-0.182421 + 2.24890 * Math.pow(TB, -0.5)) * deltaSV);
474+
double VC = VCnalkane * Math.pow(((1 + 2 * fV) / (1 - 2 * fV)), 2);
475+
double fP = deltaSP * ((2.53262 - 34.4321 * Math.pow(TB, -0.5) - 0.00230193 * TB)
476+
+ (-11.4277 + 187.934 * Math.pow(TB, -0.5) + 0.00414963 * TB) * deltaSP);
477+
double PC = PCnalkane * (TC / Tcnalkane) * (VCnalkane / VC)
478+
* Math.pow(((1 + 2 * fP) / (1 - 2 * fP)), 2);
479+
return VC * 1e3; // m3/mol
480+
}
481+
482+
}
483+
315484
/**
316485
* <p>
317486
* getModel.
@@ -334,6 +503,10 @@ public TBPModelInterface getModel(String name) {
334503
return new PedersenTBPModelPRHeavyOil();
335504
} else if (name.equals("RiaziDaubert")) {
336505
return new RiaziDaubert();
506+
} else if (name.equals("Lee-Kesler")) {
507+
return new LeeKesler();
508+
} else if (name.equals("Twu")) {
509+
return new TwuModel();
337510
} else {
338511
// System.out.println("not a valid TBPModelName.................");
339512
return new PedersenTBPModelSRK();

‎src/main/java/neqsim/thermo/system/SystemThermo.java

+3-1
Original file line numberDiff line numberDiff line change
@@ -869,7 +869,9 @@ public void addTBPfraction(String componentName, double numberOfMoles, double mo
869869
molarMass = 1000 * molarMass;
870870
TC = characterization.getTBPModel().calcTC(molarMass, density);
871871
PC = characterization.getTBPModel().calcPC(molarMass, density);
872-
m = characterization.getTBPModel().calcm(molarMass, density);
872+
if (characterization.getTBPModel().isCalcm()) {
873+
m = characterization.getTBPModel().calcm(molarMass, density);
874+
}
873875
acs = characterization.getTBPModel().calcAcentricFactor(molarMass, density);
874876
// TBPfractionCoefs[2][0]+TBPfractionCoefs[2][1]*molarMass+TBPfractionCoefs[2][2]*density+TBPfractionCoefs[2][3]*Math.pow(molarMass,2.0);
875877
TB = characterization.getTBPModel().calcTB(molarMass, density);

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

+55
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,55 @@
1+
package neqsim.thermo.characterization;
2+
3+
import static org.junit.jupiter.api.Assertions.assertEquals;
4+
import org.junit.jupiter.api.Test;
5+
import neqsim.thermo.system.SystemInterface;
6+
import neqsim.thermo.system.SystemPrEos;
7+
import neqsim.thermo.system.SystemSrkEos;
8+
9+
public class TBPfractionModelTest {
10+
@Test
11+
void testTwuModel() {
12+
SystemInterface thermoSystem = new SystemSrkEos(298.0, 10.0);
13+
thermoSystem.getCharacterization().setTBPModel("Twu");
14+
thermoSystem.addTBPfraction("C7", 1.0, 110.0 / 1000.0, 0.73);
15+
assertEquals(536.173400, thermoSystem.getComponent(0).getTC(), 1e-3);
16+
assertEquals(26.52357312690, thermoSystem.getComponent(0).getPC(), 1e-3);
17+
assertEquals(0.56001213933, thermoSystem.getComponent(0).getAcentricFactor(), 1e-3);
18+
assertEquals(437.335493, thermoSystem.getComponent(0).getCriticalVolume(), 1e-3);
19+
}
20+
21+
@Test
22+
void testLeeKeslerModel() {
23+
SystemInterface thermoSystem = new SystemSrkEos(298.0, 10.0);
24+
thermoSystem.getCharacterization().setTBPModel("Lee-Kesler");
25+
thermoSystem.addTBPfraction("C7", 1.0, 110.0 / 1000.0, 0.73);
26+
assertEquals(562.4229803010662, thermoSystem.getComponent(0).getTC(), 1e-3);
27+
assertEquals(28.322987349048354, thermoSystem.getComponent(0).getPC(), 1e-3);
28+
assertEquals(0.3509842412742902, thermoSystem.getComponent(0).getAcentricFactor(), 1e-3);
29+
assertEquals(427.99744457199, thermoSystem.getComponent(0).getCriticalVolume(), 1e-3);
30+
}
31+
32+
@Test
33+
void testPedersenSRKModel() {
34+
SystemInterface thermoSystem = new SystemSrkEos(298.0, 10.0);
35+
thermoSystem.getCharacterization().setTBPModel("PedersenSRK");
36+
thermoSystem.addTBPfraction("C7", 1.0, 110.0 / 1000.0, 0.73);
37+
assertEquals(554.3185637098962, thermoSystem.getComponent(0).getTC(), 1e-3);
38+
assertEquals(26.007549082822628, thermoSystem.getComponent(0).getPC(), 1e-3);
39+
assertEquals(0.508241, thermoSystem.getComponent(0).getAcentricFactor(), 1e-3);
40+
assertEquals(384.6714299777243, thermoSystem.getComponent(0).getCriticalVolume(), 1e-3);
41+
}
42+
43+
@Test
44+
void testPedersenPRModel() {
45+
SystemInterface thermoSystem = new SystemPrEos(298.0, 10.0);
46+
thermoSystem.getCharacterization().setTBPModel("PedersenPR");
47+
thermoSystem.addTBPfraction("C7", 1.0, 110.0 / 1000.0, 0.73);
48+
assertEquals(560.546, thermoSystem.getComponent(0).getTC(), 1e-3);
49+
assertEquals(25.838137535018557, thermoSystem.getComponent(0).getPC(), 1e-3);
50+
assertEquals(0.3838836222383, thermoSystem.getComponent(0).getAcentricFactor(), 1e-3);
51+
assertEquals(405.0890245138075, thermoSystem.getComponent(0).getCriticalVolume(), 1e-3);
52+
}
53+
54+
55+
}

0 commit comments

Comments
 (0)