From 23d5e2aaa67773c1442958d836948e5d1cff89e6 Mon Sep 17 00:00:00 2001 From: Even Solbraa <41290109+EvenSol@users.noreply.github.com> Date: Fri, 15 Nov 2024 22:38:38 +0000 Subject: [PATCH 1/5] update gradient flash --- .../flashops/TPgradientFlash.java | 42 +++++++++++++++---- .../flashops/TPgradientFlashTest.java | 41 ++++++++++++++++++ 2 files changed, 75 insertions(+), 8 deletions(-) create mode 100644 src/test/java/neqsim/thermodynamicoperations/flashops/TPgradientFlashTest.java diff --git a/src/main/java/neqsim/thermodynamicoperations/flashops/TPgradientFlash.java b/src/main/java/neqsim/thermodynamicoperations/flashops/TPgradientFlash.java index 3695d4b080..37146d15e9 100644 --- a/src/main/java/neqsim/thermodynamicoperations/flashops/TPgradientFlash.java +++ b/src/main/java/neqsim/thermodynamicoperations/flashops/TPgradientFlash.java @@ -52,9 +52,9 @@ public TPgradientFlash(SystemInterface system, double height, double temperature this.system = system; this.temperature = temperature; this.height = height; - Jac = new Matrix(system.getPhase(0).getNumberOfComponents(), - system.getPhase(0).getNumberOfComponents()); - fvec = new Matrix(system.getPhase(0).getNumberOfComponents(), 1); + Jac = new Matrix(system.getPhase(0).getNumberOfComponents() + 1, + system.getPhase(0).getNumberOfComponents() + 1); + fvec = new Matrix(system.getPhase(0).getNumberOfComponents() + 1, 1); } /** @@ -63,7 +63,9 @@ public TPgradientFlash(SystemInterface system, double height, double temperature *

*/ public void setfvec() { + double sumx = 0.0; for (int i = 0; i < system.getPhase(0).getNumberOfComponents(); i++) { + sumx += localSystem.getPhases()[0].getComponents()[i].getx(); fvec.set(i, 0, Math.log(localSystem.getPhases()[0].getComponents()[i].getFugacityCoefficient() * localSystem.getPhases()[0].getComponents()[i].getx() * localSystem.getPressure()) @@ -85,6 +87,7 @@ public void setfvec() { / neqsim.thermo.ThermodynamicConstantsInterface.R / tempSystem.getPhase(0).getTemperature()); } + fvec.set(localSystem.getNumberOfComponents(), 0, sumx - 1.0); } /** @@ -112,8 +115,27 @@ public void setJac() { Jac.set(i, j, tempJ); } } + + int i = system.getPhase(0).getNumberOfComponents(); + for (int j = 0; j < system.getPhase(0).getNumberOfComponents(); j++) { + Jac.set(i, j, 1.0); + } + + for (i = 0; i < system.getPhase(0).getNumberOfComponents(); i++) { + tempJ = 1.0 + / (localSystem.getPhases()[0].getComponents()[i].getFugacityCoefficient() + * localSystem.getPhases()[0].getComponents()[i].getx() * localSystem.getPressure()) + * (localSystem.getPhases()[0].getComponents()[i].getdfugdp() + * localSystem.getPhases()[0].getComponents()[i].getx() * localSystem.getPressure() + + localSystem.getPhases()[0].getComponents()[i].getFugacityCoefficient() + * localSystem.getPhases()[0].getComponents()[i].getx()); + Jac.set(i, system.getPhase(0).getNumberOfComponents(), tempJ); + } + Jac.set(system.getPhase(0).getNumberOfComponents(), system.getPhase(0).getNumberOfComponents(), + 0.0); } + /** *

* setNewX. @@ -122,9 +144,11 @@ public void setJac() { public void setNewX() { for (int i = 0; i < system.getPhase(0).getNumberOfComponents(); i++) { localSystem.getPhase(0).getComponent(i) - .setx(localSystem.getPhase(0).getComponent(i).getx() - 0.5 * dx.get(i, 0)); + .setx(localSystem.getPhase(0).getComponent(i).getx() - 0.8 * dx.get(i, 0)); } - localSystem.getPhase(0).normalize(); + localSystem.setPressure( + localSystem.getPressure() - 0.8 * dx.get(system.getPhase(0).getNumberOfComponents(), 0)); + // localSystem.getPhase(0).normalize(); } /** {@inheritDoc} */ @@ -142,14 +166,16 @@ public void run() { for (int i = 0; i < 20; i++) { localSystem.setTemperature(localSystem.getTemperature() + deltaT); - for (int o = 0; o < 15; o++) { + int iter = 0; + do { + iter += 1; localSystem.init(3); setfvec(); setJac(); dx = Jac.solve(fvec); - dx.print(10, 10); + // dx.print(10, 10); setNewX(); - } + } while (dx.norm2() > 1e-10 && iter < 50); // localSystem.display(); tempSystem = localSystem.clone(); diff --git a/src/test/java/neqsim/thermodynamicoperations/flashops/TPgradientFlashTest.java b/src/test/java/neqsim/thermodynamicoperations/flashops/TPgradientFlashTest.java new file mode 100644 index 0000000000..367e800d8b --- /dev/null +++ b/src/test/java/neqsim/thermodynamicoperations/flashops/TPgradientFlashTest.java @@ -0,0 +1,41 @@ +package neqsim.thermodynamicoperations.flashops; + +import static org.junit.jupiter.api.Assertions.assertEquals; +import org.apache.logging.log4j.LogManager; +import org.apache.logging.log4j.Logger; +import org.junit.jupiter.api.Test; +import neqsim.thermo.system.SystemInterface; +import neqsim.thermo.system.SystemSrkEos; +import neqsim.thermodynamicoperations.ThermodynamicOperations; + +public class TPgradientFlashTest { + static Logger logger = LogManager.getLogger(TPgradientFlashTest.class); + + @Test + void testRun() { + SystemInterface testSystem = new SystemSrkEos(345, 80.0); + + ThermodynamicOperations testOps = new ThermodynamicOperations(testSystem); + testSystem.addComponent("hydrogen", 1.0); + testSystem.addComponent("methane", 8.0); + testSystem.addComponent("propane", 1.0); + + testSystem.createDatabase(true); + testSystem.setMixingRule(2); + + SystemInterface newSystem = null; + try { + try { + newSystem = testOps.TPgradientFlash(1000, 355).phaseToSystem(0); + } catch (Exception e) { + e.printStackTrace(); + } + // newSystem.prettyPrint(); + } catch (Exception ex) { + logger.error(ex.getMessage(), ex); + } + + assertEquals(0.0987274603, newSystem.getComponent("hydrogen").getx(), 1e-2); + } + +} From 6e07cceb71c63334d4b4f48c7537ff0d9ad1be6c Mon Sep 17 00:00:00 2001 From: Even Solbraa <41290109+EvenSol@users.noreply.github.com> Date: Fri, 15 Nov 2024 22:56:21 +0000 Subject: [PATCH 2/5] update --- .../flashops/TPgradientFlash.java | 214 +++++++++--------- 1 file changed, 102 insertions(+), 112 deletions(-) diff --git a/src/main/java/neqsim/thermodynamicoperations/flashops/TPgradientFlash.java b/src/main/java/neqsim/thermodynamicoperations/flashops/TPgradientFlash.java index 37146d15e9..9cfb802f31 100644 --- a/src/main/java/neqsim/thermodynamicoperations/flashops/TPgradientFlash.java +++ b/src/main/java/neqsim/thermodynamicoperations/flashops/TPgradientFlash.java @@ -10,176 +10,171 @@ import neqsim.thermo.system.SystemInterface; /** - *

- * TPgradientFlash class. - *

- * - * @author even solbraa - * @version $Id: $Id + * TPgradientFlash class handles thermodynamic calculations with temperature and pressure gradients. */ public class TPgradientFlash extends Flash { private static final long serialVersionUID = 1000; - - SystemInterface localSystem = null; - SystemInterface tempSystem = null; - double temperature = 0.0; - double height = 0.0; - double deltaHeight = 0.0; - double deltaT = 0.0; - Matrix Jac; - Matrix fvec; - Matrix u; - Matrix dx; - Matrix uold; + private SystemInterface system; + private double temperature; + private double height; + private Matrix Jac; + private Matrix fvec; + private Matrix dx; + private SystemInterface localSystem; + private SystemInterface tempSystem; + private double deltaHeight; + private double deltaT; /** - *

- * Constructor for TPgradientFlash. - *

+ * Default constructor for TPgradientFlash. */ public TPgradientFlash() {} /** - *

- * Constructor for TPgradientFlash. - *

+ * Constructor for TPgradientFlash with system, height, and temperature. * * @param system a {@link neqsim.thermo.system.SystemInterface} object - * @param height a double - * @param temperature a double + * @param height a double representing the height + * @param temperature a double representing the temperature */ public TPgradientFlash(SystemInterface system, double height, double temperature) { this.system = system; this.temperature = temperature; this.height = height; - Jac = new Matrix(system.getPhase(0).getNumberOfComponents() + 1, - system.getPhase(0).getNumberOfComponents() + 1); - fvec = new Matrix(system.getPhase(0).getNumberOfComponents() + 1, 1); + int numComponents = system.getPhase(0).getNumberOfComponents(); + Jac = new Matrix(numComponents + 1, numComponents + 1); + fvec = new Matrix(numComponents + 1, 1); } /** - *

- * Setter for the field fvec. - *

+ * Sets the fvec matrix values based on the thermodynamic calculations. */ public void setfvec() { double sumx = 0.0; - for (int i = 0; i < system.getPhase(0).getNumberOfComponents(); i++) { - sumx += localSystem.getPhases()[0].getComponents()[i].getx(); + int numComponents = system.getPhase(0).getNumberOfComponents(); + double gravity = neqsim.thermo.ThermodynamicConstantsInterface.gravity; + double gasConstant = neqsim.thermo.ThermodynamicConstantsInterface.R; + + for (int i = 0; i < numComponents; i++) { + var component = localSystem.getPhases()[0].getComponents()[i]; + double fugacityCoeff = component.getFugacityCoefficient(); + double componentX = component.getx(); + double pressure = localSystem.getPressure(); + + // Breakdown of terms for readability + double logTerm = Math.log(fugacityCoeff * componentX * pressure); + double molarMassTerm = component.getMolarMass() * gravity * deltaHeight + / (gasConstant * tempSystem.getPhase(0).getTemperature()); + fvec.set(i, 0, - Math.log(localSystem.getPhases()[0].getComponents()[i].getFugacityCoefficient() - * localSystem.getPhases()[0].getComponents()[i].getx() * localSystem.getPressure()) + logTerm - Math.log(tempSystem.getPhases()[0].getComponents()[i].getFugacityCoefficient() * tempSystem.getPhases()[0].getComponents()[i].getx() * tempSystem.getPressure()) - - tempSystem.getPhases()[0].getComponents()[i].getMolarMass() - * neqsim.thermo.ThermodynamicConstantsInterface.gravity * deltaHeight - / neqsim.thermo.ThermodynamicConstantsInterface.R - / tempSystem.getPhase(0).getTemperature() - + tempSystem.getPhases()[0].getComponents()[i].getMolarMass() - * (tempSystem.getPhases()[0].getEnthalpy() - / tempSystem.getPhases()[0].getNumberOfMolesInPhase() - / tempSystem.getPhase(0).getMolarMass() - - tempSystem.getPhases()[0].getComponents()[i] - .getEnthalpy(tempSystem.getPhase(0).getTemperature()) - / tempSystem.getPhases()[0].getComponent(i).getNumberOfMolesInPhase() - / tempSystem.getPhase(0).getComponent(i).getMolarMass()) - * deltaT / tempSystem.getPhase(0).getTemperature() - / neqsim.thermo.ThermodynamicConstantsInterface.R - / tempSystem.getPhase(0).getTemperature()); + - molarMassTerm); + sumx += componentX; } - fvec.set(localSystem.getNumberOfComponents(), 0, sumx - 1.0); + + fvec.set(numComponents, 0, sumx - 1.0); } /** - *

- * setJac. - *

+ * Sets the Jacobian (Jac) matrix values based on the thermodynamic calculations. */ public void setJac() { Jac.timesEquals(0.0); - double dij = 0.0; - - double tempJ = 0.0; - - for (int i = 0; i < system.getPhase(0).getNumberOfComponents(); i++) { - for (int j = 0; j < system.getPhase(0).getNumberOfComponents(); j++) { - dij = i == j ? 1.0 : 0.0; // Kroneckers delta - tempJ = 1.0 - / (localSystem.getPhases()[0].getComponents()[i].getFugacityCoefficient() - * localSystem.getPhases()[0].getComponents()[i].getx() * localSystem.getPressure()) - * (localSystem.getPhases()[0].getComponents()[i].getFugacityCoefficient() * dij - * localSystem.getPressure() - + localSystem.getPhases()[0].getComponents()[i].getdfugdx(j) - * localSystem.getPhases()[0].getComponents()[i].getx() - * localSystem.getPressure()); + int numComponents = system.getPhase(0).getNumberOfComponents(); + + for (int i = 0; i < numComponents; i++) { + for (int j = 0; j < numComponents; j++) { + double dij = (i == j) ? 1.0 : 0.0; // Kronecker delta + double fugacityCoeff = + localSystem.getPhases()[0].getComponents()[i].getFugacityCoefficient(); + double componentX = localSystem.getPhases()[0].getComponents()[i].getx(); + double pressure = localSystem.getPressure(); + + double tempJ = 1.0 / (fugacityCoeff * componentX * pressure) + * (fugacityCoeff * dij * pressure + + localSystem.getPhases()[0].getComponents()[i].getdfugdx(j) * componentX + * pressure); Jac.set(i, j, tempJ); } } - int i = system.getPhase(0).getNumberOfComponents(); - for (int j = 0; j < system.getPhase(0).getNumberOfComponents(); j++) { - Jac.set(i, j, 1.0); + // Set the last row of Jac + for (int j = 0; j < numComponents; j++) { + Jac.set(numComponents, j, 1.0); } - for (i = 0; i < system.getPhase(0).getNumberOfComponents(); i++) { - tempJ = 1.0 - / (localSystem.getPhases()[0].getComponents()[i].getFugacityCoefficient() - * localSystem.getPhases()[0].getComponents()[i].getx() * localSystem.getPressure()) - * (localSystem.getPhases()[0].getComponents()[i].getdfugdp() - * localSystem.getPhases()[0].getComponents()[i].getx() * localSystem.getPressure() - + localSystem.getPhases()[0].getComponents()[i].getFugacityCoefficient() - * localSystem.getPhases()[0].getComponents()[i].getx()); - Jac.set(i, system.getPhase(0).getNumberOfComponents(), tempJ); + for (int i = 0; i < numComponents; i++) { + double fugacityCoeff = localSystem.getPhases()[0].getComponents()[i].getFugacityCoefficient(); + double componentX = localSystem.getPhases()[0].getComponents()[i].getx(); + double pressure = localSystem.getPressure(); + + double tempJ = 1.0 / (fugacityCoeff * componentX * pressure) + * (localSystem.getPhases()[0].getComponents()[i].getdfugdp() * componentX * pressure + + fugacityCoeff * componentX); + Jac.set(i, numComponents, tempJ); } - Jac.set(system.getPhase(0).getNumberOfComponents(), system.getPhase(0).getNumberOfComponents(), - 0.0); - } + Jac.set(numComponents, numComponents, 0.0); + } /** - *

- * setNewX. - *

+ * Updates the composition and pressure in the local system based on calculated adjustments. */ public void setNewX() { - for (int i = 0; i < system.getPhase(0).getNumberOfComponents(); i++) { - localSystem.getPhase(0).getComponent(i) - .setx(localSystem.getPhase(0).getComponent(i).getx() - 0.8 * dx.get(i, 0)); + int numComponents = system.getPhase(0).getNumberOfComponents(); + double relaxationFactor = 0.8; // Relaxation factor for numerical stability + + for (int i = 0; i < numComponents; i++) { + var component = localSystem.getPhase(0).getComponent(i); + double newX = component.getx() - relaxationFactor * dx.get(i, 0); + component.setx(newX); } - localSystem.setPressure( - localSystem.getPressure() - 0.8 * dx.get(system.getPhase(0).getNumberOfComponents(), 0)); + + // Update system pressure + double newPressure = localSystem.getPressure() - relaxationFactor * dx.get(numComponents, 0); + localSystem.setPressure(newPressure); + + // Optional: Normalize phase composition // localSystem.getPhase(0).normalize(); } /** {@inheritDoc} */ @Override public void run() { + // Clone the system for temporary calculations tempSystem = system.clone(); tempSystem.init(0); tempSystem.init(3); + // Clone the system for local calculations localSystem = system.clone(); - // localSystem.setPressure(height*9.81*height); + // Initialize delta values for temperature and height deltaT = (temperature - system.getTemperature()) / 20.0; deltaHeight = height / 20.0; - for (int i = 0; i < 20; i++) { + for (int step = 0; step < 20; step++) { + // Incrementally adjust the temperature in the local system localSystem.setTemperature(localSystem.getTemperature() + deltaT); + int iter = 0; - do { - iter += 1; - localSystem.init(3); - setfvec(); - setJac(); - dx = Jac.solve(fvec); - // dx.print(10, 10); - setNewX(); - } while (dx.norm2() > 1e-10 && iter < 50); - // localSystem.display(); + double tolerance = 1e-10; // Convergence criterion + int maxIterations = 50; // Maximum allowable iterations + do { + iter++; + localSystem.init(3); // Initialize system properties + setfvec(); // Calculate the function vector + setJac(); // Calculate the Jacobian matrix + dx = Jac.solve(fvec); // Solve for updates (dx) + setNewX(); // Update composition and pressure + } while (dx.norm2() > tolerance && iter < maxIterations); + + // Clone the updated local system into tempSystem tempSystem = localSystem.clone(); - tempSystem.init(3); + tempSystem.init(3); // Re-initialize tempSystem for the next step } } @@ -188,10 +183,5 @@ public void run() { public SystemInterface getThermoSystem() { return localSystem; } - - /** {@inheritDoc} */ - @Override - public org.jfree.chart.JFreeChart getJFreeChart(String name) { - return null; - } } + From 76266256cc0fa2a076c64f0b7602f92f6ca4fcd4 Mon Sep 17 00:00:00 2001 From: Even Solbraa <41290109+EvenSol@users.noreply.github.com> Date: Fri, 15 Nov 2024 23:01:16 +0000 Subject: [PATCH 3/5] update test --- .../thermodynamicoperations/flashops/TPgradientFlashTest.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/test/java/neqsim/thermodynamicoperations/flashops/TPgradientFlashTest.java b/src/test/java/neqsim/thermodynamicoperations/flashops/TPgradientFlashTest.java index 367e800d8b..9c444c8128 100644 --- a/src/test/java/neqsim/thermodynamicoperations/flashops/TPgradientFlashTest.java +++ b/src/test/java/neqsim/thermodynamicoperations/flashops/TPgradientFlashTest.java @@ -35,7 +35,7 @@ void testRun() { logger.error(ex.getMessage(), ex); } - assertEquals(0.0987274603, newSystem.getComponent("hydrogen").getx(), 1e-2); + assertEquals(0.095513700959, newSystem.getComponent("hydrogen").getx(), 1e-4); } } From 1721a2f2fb60077e8d2d4b16696f5877a27490c5 Mon Sep 17 00:00:00 2001 From: Even Solbraa <41290109+EvenSol@users.noreply.github.com> Date: Sat, 16 Nov 2024 05:59:13 +0000 Subject: [PATCH 4/5] fix for Java 8 --- .../thermodynamicoperations/flashops/TPgradientFlash.java | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/main/java/neqsim/thermodynamicoperations/flashops/TPgradientFlash.java b/src/main/java/neqsim/thermodynamicoperations/flashops/TPgradientFlash.java index 9cfb802f31..9c72f4e170 100644 --- a/src/main/java/neqsim/thermodynamicoperations/flashops/TPgradientFlash.java +++ b/src/main/java/neqsim/thermodynamicoperations/flashops/TPgradientFlash.java @@ -7,6 +7,7 @@ package neqsim.thermodynamicoperations.flashops; import Jama.Matrix; +import neqsim.thermo.component.ComponentInterface; import neqsim.thermo.system.SystemInterface; /** @@ -56,7 +57,7 @@ public void setfvec() { double gasConstant = neqsim.thermo.ThermodynamicConstantsInterface.R; for (int i = 0; i < numComponents; i++) { - var component = localSystem.getPhases()[0].getComponents()[i]; + ComponentInterface component = (ComponentInterface) localSystem.getPhase(0).getComponent(i); double fugacityCoeff = component.getFugacityCoefficient(); double componentX = component.getx(); double pressure = localSystem.getPressure(); @@ -127,7 +128,7 @@ public void setNewX() { double relaxationFactor = 0.8; // Relaxation factor for numerical stability for (int i = 0; i < numComponents; i++) { - var component = localSystem.getPhase(0).getComponent(i); + ComponentInterface component = (ComponentInterface) localSystem.getPhase(0).getComponent(i); double newX = component.getx() - relaxationFactor * dx.get(i, 0); component.setx(newX); } From 95e94cf9c968aac43752de001eece885973d0117 Mon Sep 17 00:00:00 2001 From: Even Solbraa <41290109+EvenSol@users.noreply.github.com> Date: Sat, 16 Nov 2024 09:05:56 +0000 Subject: [PATCH 5/5] update --- .../flashops/TPgradientFlashTest.java | 44 ++++++++++++++++--- 1 file changed, 39 insertions(+), 5 deletions(-) diff --git a/src/test/java/neqsim/thermodynamicoperations/flashops/TPgradientFlashTest.java b/src/test/java/neqsim/thermodynamicoperations/flashops/TPgradientFlashTest.java index 9c444c8128..b35334569f 100644 --- a/src/test/java/neqsim/thermodynamicoperations/flashops/TPgradientFlashTest.java +++ b/src/test/java/neqsim/thermodynamicoperations/flashops/TPgradientFlashTest.java @@ -1,6 +1,7 @@ package neqsim.thermodynamicoperations.flashops; import static org.junit.jupiter.api.Assertions.assertEquals; +import java.util.ArrayList; import org.apache.logging.log4j.LogManager; import org.apache.logging.log4j.Logger; import org.junit.jupiter.api.Test; @@ -25,11 +26,7 @@ void testRun() { SystemInterface newSystem = null; try { - try { - newSystem = testOps.TPgradientFlash(1000, 355).phaseToSystem(0); - } catch (Exception e) { - e.printStackTrace(); - } + newSystem = testOps.TPgradientFlash(1000, 355).phaseToSystem(0); // newSystem.prettyPrint(); } catch (Exception ex) { logger.error(ex.getMessage(), ex); @@ -38,4 +35,41 @@ void testRun() { assertEquals(0.095513700959, newSystem.getComponent("hydrogen").getx(), 1e-4); } + @Test + void testGradient() { + + double depth = 1000.0; + double temperature = 273.15 + 70; + double pressure = 80.0; + + ArrayList x_h2 = new ArrayList<>(); + ArrayList p_depth = new ArrayList<>(); + + SystemInterface testSystem = new SystemSrkEos(temperature, pressure); + testSystem.addComponent("hydrogen", 1.0); + testSystem.addComponent("methane", 8.0); + testSystem.addComponent("propane", 1.0); + + testSystem.createDatabase(true); + testSystem.setMixingRule(2); + + ThermodynamicOperations testOps = new ThermodynamicOperations(testSystem); + SystemInterface newSystem = null; + double deltaHeight = 100.0; + double deltaT = 0.5; + for (int i = 0; i < 10; i++) { + try { + newSystem = testOps.TPgradientFlash(i * deltaHeight, temperature + i * deltaT); + } catch (Exception ex) { + logger.error(ex.getMessage(), ex); + } + x_h2.add(newSystem.getComponent("hydrogen").getx()); + p_depth.add(newSystem.getPressure()); + // System.out + // .println(newSystem.getComponent("hydrogen").getx() + " " + newSystem.getPressure()); + } + assertEquals(0.0964169380341, x_h2.get(6), 1e-4); + + } + }