Skip to content

Commit

Permalink
PFCT-bug-fix (#1271)
Browse files Browse the repository at this point in the history 
* Fixed PFCT viscosity function so that it can handle higher mole fractions of hydrogen. Also changed test-class so that it is easier to test the viscosity functions.

* PFCT viscosity bug fix
  • Loading branch information
@akselrs
akselrs authored Feb 7, 2025
1 parent 172440b commit 0a560d0
Show file tree
Hide file tree
Showing 2 changed files with 34 additions and 49 deletions.
72 changes: 28 additions & 44 deletions ...lproperties/methods/commonphasephysicalproperties/viscosity/PFCTViscosityMethodMod86.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,12 +20,12 @@ public class PFCTViscosityMethodMod86 extends Viscosity {
// SystemInterface referenceSystem = new SystemBWRSEos(273.15,
// ThermodynamicConstantsInterface.referencePressure);
SystemInterface referenceSystem =
new SystemSrkEos(273.0, ThermodynamicConstantsInterface.referencePressure);
new SystemSrkEos(273.15, ThermodynamicConstantsInterface.referencePressure);

// todo: is this parameter required?
int phaseTypeNumb = 1;
double[] GVcoef = {-2.090975e5, 2.647269e5, -1.472818e5, 4.716740e4, -9.491872e3, 1.219979e3,
-9.627993e1, 4.274152, -8.141531e-2};
-9.627993e1, 4.274152, -8.141531e-2};
double visRefA = 1.696985927;
double visRefB = -0.133372346;
double visRefC = 1.4;
Expand All @@ -34,7 +34,7 @@ public class PFCTViscosityMethodMod86 extends Viscosity {
double visRefG = 0.0;

double[] viscRefJ = {-1.035060586e1, 1.7571599671e1, -3.0193918656e3, 1.8873011594e2,
4.2903609488e-2, 1.4529023444e2, 6.1276818706e3};
4.2903609488e-2, 1.4529023444e2, 6.1276818706e3};
double[] viscRefK = {-9.74602, 18.0834, -4126.66, 44.6055, 0.976544, 81.8134, 15649.9};

/**
Expand Down Expand Up @@ -93,37 +93,29 @@ public double calcViscosity() {
}
PCmix = 8.0 * tempPC1 / (tempPC2 * tempPC2);
double TCmix = tempTC1 / tempTC2;
double Mmix =
(Mmtemp + 1.304e-4 * (Math.pow(Mwtemp / Mmtemp, 2.303) - Math.pow(Mmtemp, 2.303))) * 1e3; // phase.getPhase().getMolarMass();
double Mmix = (Mmtemp + 1.304e-4 * (Math.pow(Mwtemp / Mmtemp, 2.303) - Math.pow(Mmtemp, 2.303))) * 1e3; // phase.getPhase().getMolarMass();

referenceSystem.setTemperature(phase.getPhase().getTemperature()
* referenceSystem.getPhase(0).getComponent(0).getTC() / TCmix);
referenceSystem.setPressure(phase.getPhase().getPressure()
* referenceSystem.getPhase(0).getComponent(0).getPC() / PCmix);
referenceSystem.setTemperature(phase.getPhase().getTemperature());
referenceSystem.setPressure(phase.getPhase().getPressure());
referenceSystem.init(1);
// todo: mixing phasetype and phase index?
// double molDens = 1.0 /
// referenceSystem.getPhase(phaseTypeNumb).getMolarVolume() * 100.0;

double molDens = 1.0 / referenceSystem.getLowestGibbsEnergyPhase().getMolarVolume() * 100.0;
double critMolDens = 10.15; // 1.0/referenceSystem.getPhase(0).getComponent(0).getCriticalVolume();
double redDens = molDens / critMolDens;

double alfaMix = 1.0 + 7.378e-3 * Math.pow(redDens, 1.847) * Math.pow(Mmix, 0.5173);
double alfa0 = 1.0 + 7.378e-3 * Math.pow(redDens, 1.847)
* Math.pow(referenceSystem.getMolarMass() * 1.0e3, 0.5173);
// alfa0 = 1.0 + 8.374e-4 * Math.pow(redDens, 4.265);
// System.out.println("func " + 7.475e-5*Math.pow(16.043, 0.8579));

double T0 = phase.getPhase().getTemperature()
* referenceSystem.getPhase(0).getComponent(0).getTC() / TCmix * alfa0 / alfaMix;
double P0 = phase.getPhase().getPressure() * referenceSystem.getPhase(0).getComponent(0).getPC()
/ PCmix * alfa0 / alfaMix;

double refVisosity = getRefComponentViscosity(T0, P0);
// System.out.println("m/mix " + Mmix/M0);
// System.out.println("a/amix " + alfaMix/alfa0);

double viscosity = refVisosity * Math.pow(TCmix / Tc0, -1.0 / 6.0)
* Math.pow(PCmix / Pc0, 2.0 / 3.0) * Math.pow(Mmix / M0, 0.5) * alfaMix / alfa0;
// System.out.println("viscosity " + refVisosity);
return viscosity;
}

Expand All @@ -138,32 +130,19 @@ public double calcViscosity() {
*/
public double getRefComponentViscosity(double temp, double pres) {
referenceSystem.setTemperature(temp);
// System.out.println("ref temp " + temp);
referenceSystem.setPressure(pres);
// System.out.println("ref pres " + pres);
referenceSystem.init(1);
// referenceSystem.display();
// todo: mixing phasetype and phase index?
// double molDens = 1.0 /
// referenceSystem.getPhase(phaseTypeNumb).getMolarVolume() * 100.0;
double molDens = 1.0 / referenceSystem.getLowestGibbsEnergyPhase().getMolarVolume() * 100.0; // mol/dm^3
// System.out.println("mol dens " + molDens);
double critMolDens = 10.15; // 1.0/referenceSystem.getPhase(0).getComponent(0).getCriticalVolume();
double molDens = referenceSystem.getLowestGibbsEnergyPhase().getDensity() * 1e-3; //[kg/L]
double critMolDens = 162.66e-3; //[kg/L] (Source: NIST)
double redMolDens = (molDens - critMolDens) / critMolDens;
// System.out.println("gv1 " +GVcoef[0]);

molDens = referenceSystem.getLowestGibbsEnergyPhase().getDensity() * 1e-3;

double viscRefO = GVcoef[0] * Math.pow(temp, -1.0) + GVcoef[1] * Math.pow(temp, -2.0 / 3.0)
+ GVcoef[2] * Math.pow(temp, -1.0 / 3.0) + GVcoef[3] + GVcoef[4] * Math.pow(temp, 1.0 / 3.0)
+ GVcoef[5] * Math.pow(temp, 2.0 / 3.0) + GVcoef[6] * temp
+ GVcoef[7] * Math.pow(temp, 4.0 / 3.0) + GVcoef[8] * Math.pow(temp, 5.0 / 3.0);

// System.out.println("ref visc0 " + viscRefO);
double viscRef1 =
(visRefA + visRefB * Math.pow(visRefC - Math.log(temp / visRefF), 2.0)) * molDens;
// System.out.println("ref visc1 " + viscRef1);

double viscRefO = 0.0;
for (int i = 0; i < GVcoef.length; i++) {
viscRefO += GVcoef[i] * Math.pow(temp, ((i + 1) - 4) / 3.0);
}


//Calculating the reference viscosity contributions:
double temp1 = Math.pow(molDens, 0.1) * (viscRefJ[1] + viscRefJ[2] / Math.pow(temp, 3.0 / 2.0));
double temp2 = redMolDens * Math.pow(molDens, 0.5)
* (viscRefJ[4] + viscRefJ[5] / temp + viscRefJ[6] / Math.pow(temp, 2.0));
Expand All @@ -172,21 +151,26 @@ public double getRefComponentViscosity(double temp, double pres) {
double dTfreeze = temp - 90.69;
double HTAN =
(Math.exp(dTfreeze) - Math.exp(-dTfreeze)) / (Math.exp(dTfreeze) + Math.exp(-dTfreeze));
visRefE = (HTAN + 1.0) / 2.0;

//This compensates for that the HTAN function is not defined for dTfreeze > 709.0:
if (dTfreeze > 709.0) {
double visRefE = 1.0;
double visRefG = 0.0;
} else {
double visRefE = (HTAN + 1.0) / 2.0;
double visRefG = (1.0 - HTAN) / 2.0;
}

double viscRef1 = (visRefA + visRefB * Math.pow(visRefC - Math.log(temp / visRefF), 2.0)) * molDens;
double viscRef2 = visRefE * Math.exp(viscRefJ[0] + viscRefJ[3] / temp) * (temp3 - 1.0);

double temp4 = Math.pow(molDens, 0.1) * (viscRefK[1] + viscRefK[2] / Math.pow(temp, 3.0 / 2.0));
double temp5 = redMolDens * Math.pow(molDens, 0.5)
* (viscRefK[4] + viscRefK[5] / temp + viscRefK[6] / Math.pow(temp, 2.0));
double temp6 = Math.exp(temp4 + temp5);
visRefG = (1.0 - HTAN) / 2.0;
double viscRef3 = visRefG * Math.exp(viscRefK[0] + viscRefK[3] / temp) * (temp6 - 1.0);

// System.out.println("ref visc2 " + viscRef2);
// System.out.println("ref visc3 " + viscRef3);
double refVisc = (viscRefO + viscRef1 + viscRef2 + viscRef3) / 1.0e7;
// System.out.println("ref visc " + refVisc);
return refVisc;
}
}
11 changes: 6 additions & 5 deletions ...alproperties/methods/commonphasephysicalproperties/viscosity/PFCTViscosityMethodTest.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,14 +9,15 @@ public class PFCTViscosityMethodTest extends neqsim.NeqSimTest {

@Test
void testCalcViscosityHydrogen() {
testSystem = new neqsim.thermo.system.SystemSrkEos(273.15 + 25.0, 42.0);
testSystem.addComponent("hydrogen", 0.5);
double T = 298.15;
double P = 42.00; //Pressure in bar(a)
testSystem = new neqsim.thermo.system.SystemSrkEos(T, P);
testSystem.addComponent("hydrogen", 1.0);
testSystem.setMixingRule("classic");
ThermodynamicOperations testOps = new ThermodynamicOperations(testSystem);
testOps.TPflash();
testSystem.getPhase("gas").getPhysicalProperties().setViscosityModel("PFCT");
testSystem.initProperties();
double expected = 7.8e-6;
double actual = testSystem.getPhase("gas").getViscosity("kg/msec");
assertEquals(expected, actual, 1e-6);
assertEquals(7.8e-6, testSystem.getPhase(0).getPhysicalProperties().getViscosity(), 1e-6);
}
}

0 comments on commit 0a560d0

Please sign in to comment.