From 7acb15426e38c0fc2afde99a94b161b1437d2fdf Mon Sep 17 00:00:00 2001 From: Even Solbraa <41290109+EvenSol@users.noreply.github.com> Date: Sun, 24 Nov 2024 14:54:35 +0000 Subject: [PATCH 1/3] solve speed method for compressor --- .../equipment/compressor/Compressor.java | 209 +++++++++++++----- .../compressor/CompressorChartTest.java | 83 +++++++ 2 files changed, 236 insertions(+), 56 deletions(-) diff --git a/src/main/java/neqsim/process/equipment/compressor/Compressor.java b/src/main/java/neqsim/process/equipment/compressor/Compressor.java index 45e1415558..ed25de3c81 100644 --- a/src/main/java/neqsim/process/equipment/compressor/Compressor.java +++ b/src/main/java/neqsim/process/equipment/compressor/Compressor.java @@ -46,6 +46,7 @@ public class Compressor extends TwoPortEquipment implements CompressorInterface private CompressorPropertyProfile propertyProfile = new CompressorPropertyProfile(); public double dH = 0.0; public double inletEnthalpy = 0; + private boolean solveSpeed = false; public double pressure = 0.0; private double speed = 3000; private double maxspeed = 30000; @@ -406,67 +407,155 @@ public void run(UUID id) { } } if (compressorChart.isUseCompressorChart()) { - do { - double actualFlowRate = thermoSystem.getFlowRate("m3/hr"); - double z_inlet = thermoSystem.getZ(); - double MW = thermoSystem.getMolarMass(); + if (solveSpeed) { + double targetPressure = getOutletPressure(); // Desired outlet pressure + double tolerance = 1e-3; // Tolerance for pressure difference + double minSpeed = getMinimumSpeed(); // Minimum speed for the compressor + double maxSpeed = getMaximumSpeed(); // Maximum speed for the compressor + double currentSpeed = getSpeed(); // Initial guess for speed + double maxIterations = 100; // Maximum number of iterations + double deltaSpeed = 10.0; // Small increment for numerical derivative + int iteration = 0; + + while (iteration < maxIterations) { + // Calculate the pressure at the current speed + double actualFlowRate = thermoSystem.getFlowRate("m3/hr"); + double z_inlet = thermoSystem.getZ(); + double MW = thermoSystem.getMolarMass(); + if (getCompressorChart().useRealKappa()) { + kappa = thermoSystem.getGamma(); + } else { + kappa = thermoSystem.getGamma2(); + } - if (getCompressorChart().useRealKappa()) { - kappa = thermoSystem.getGamma(); - } else { - kappa = thermoSystem.getGamma2(); - } - if (useGERG2008 && inStream.getThermoSystem().getNumberOfPhases() == 1) { - double[] gergProps; - gergProps = getThermoSystem().getPhase(0).getProperties_GERG2008(); - actualFlowRate *= gergProps[1] / z_inlet; - kappa = gergProps[14]; - z_inlet = gergProps[1]; - } + if (useGERG2008 && inStream.getThermoSystem().getNumberOfPhases() == 1) { + double[] gergProps = getThermoSystem().getPhase(0).getProperties_GERG2008(); + actualFlowRate *= gergProps[1] / z_inlet; + kappa = gergProps[14]; + z_inlet = gergProps[1]; + } - double polytropEff = - getCompressorChart().getPolytropicEfficiency(actualFlowRate, getSpeed()); - setPolytropicEfficiency(polytropEff / 100.0); - polytropicHead = getCompressorChart().getPolytropicHead(actualFlowRate, getSpeed()); - double temperature_inlet = thermoSystem.getTemperature(); - double n = 1.0 / (1.0 - (kappa - 1.0) / kappa * 1.0 / (polytropEff / 100.0)); - polytropicExponent = n; - if (getCompressorChart().getHeadUnit().equals("meter")) { - polytropicFluidHead = polytropicHead / 1000.0 * 9.81; - polytropicHeadMeter = polytropicHead; - } else { - polytropicFluidHead = polytropicHead; - polytropicHeadMeter = polytropicHead * 1000.0 / 9.81; - } - double pressureRatio = Math.pow((polytropicFluidHead * 1000.0 + (n / (n - 1.0) * z_inlet - * ThermodynamicConstantsInterface.R * (temperature_inlet) / MW)) - / (n / (n - 1.0) * z_inlet * ThermodynamicConstantsInterface.R * (temperature_inlet) - / MW), - n / (n - 1.0)); - setOutletPressure(thermoSystem.getPressure() * pressureRatio); - if (getAntiSurge().isActive()) { - logger.info("surge flow " - + getCompressorChart().getSurgeCurve().getSurgeFlow(polytropicHead) + " m3/hr"); - surgeCheck = isSurge(polytropicHead, actualFlowRate); - } - if (getCompressorChart().getStoneWallCurve().isActive()) { - // logger.info("stone wall? " + isStoneWall(polytropicHead, - // thermoSystem.getFlowRate("m3/hr"))); + double polytropEff = + getCompressorChart().getPolytropicEfficiency(actualFlowRate, currentSpeed); + setPolytropicEfficiency(polytropEff / 100.0); + polytropicHead = getCompressorChart().getPolytropicHead(actualFlowRate, currentSpeed); + double temperature_inlet = thermoSystem.getTemperature(); + double n = 1.0 / (1.0 - (kappa - 1.0) / kappa * 1.0 / (polytropEff / 100.0)); + double polytropicFluidHead = + (getCompressorChart().getHeadUnit().equals("meter")) ? polytropicHead / 1000.0 * 9.81 + : polytropicHead; + double pressureRatio = Math.pow((polytropicFluidHead * 1000.0 + (n / (n - 1.0) * z_inlet + * ThermodynamicConstantsInterface.R * (temperature_inlet) / MW)) + / (n / (n - 1.0) * z_inlet * ThermodynamicConstantsInterface.R * (temperature_inlet) + / MW), + n / (n - 1.0)); + double currentPressure = thermoSystem.getPressure() * pressureRatio; + + // Calculate the derivative of pressure with respect to speed + double polytropEffDelta = getCompressorChart().getPolytropicEfficiency(actualFlowRate, + currentSpeed + deltaSpeed); + double polytropicHeadDelta = + getCompressorChart().getPolytropicHead(actualFlowRate, currentSpeed + deltaSpeed); + double nDelta = 1.0 / (1.0 - (kappa - 1.0) / kappa * 1.0 / (polytropEffDelta / 100.0)); + double polytropicFluidHeadDelta = (getCompressorChart().getHeadUnit().equals("meter")) + ? polytropicHeadDelta / 1000.0 * 9.81 + : polytropicHeadDelta; + double pressureRatioDelta = + Math.pow((polytropicFluidHeadDelta * 1000.0 + (nDelta / (nDelta - 1.0) * z_inlet + * ThermodynamicConstantsInterface.R * (temperature_inlet) / MW)) + / (nDelta / (nDelta - 1.0) * z_inlet * ThermodynamicConstantsInterface.R + * (temperature_inlet) / MW), + nDelta / (nDelta - 1.0)); + double pressureDelta = thermoSystem.getPressure() * pressureRatioDelta; + + double dPressure_dSpeed = (pressureDelta - currentPressure) / deltaSpeed; + + // Update speed using Newton-Raphson method + double speedUpdate = (targetPressure - currentPressure) / dPressure_dSpeed; + currentSpeed += 0.8 * speedUpdate; + + // Check if speed is within bounds + if (currentSpeed < minSpeed || currentSpeed > maxSpeed) { + throw new IllegalArgumentException( + "Speed out of bounds during Newton-Raphson iteration."); + } + + // Check for convergence + if (Math.abs(currentPressure - targetPressure) <= tolerance) { + setSpeed(currentSpeed); // Update the final speed + break; + } + + iteration++; } - if (surgeCheck && getAntiSurge().isActive()) { - thermoSystem.setTotalFlowRate( - getAntiSurge().getSurgeControlFactor() - * getCompressorChart().getSurgeCurve().getSurgeFlow(polytropicFluidHead), - "Am3/hr"); - thermoSystem.init(3); - fractionAntiSurge = thermoSystem.getTotalNumberOfMoles() / orginalMolarFLow - 1.0; - getAntiSurge().setCurrentSurgeFraction(fractionAntiSurge); + + if (iteration == maxIterations) { + throw new RuntimeException( + "Newton-Raphson method did not converge within the maximum iterations."); } + } else { + do { + double actualFlowRate = thermoSystem.getFlowRate("m3/hr"); + double z_inlet = thermoSystem.getZ(); + double MW = thermoSystem.getMolarMass(); - powerSet = true; - dH = polytropicFluidHead * 1000.0 * thermoSystem.getMolarMass() / getPolytropicEfficiency() - * thermoSystem.getTotalNumberOfMoles(); - } while (surgeCheck && getAntiSurge().isActive()); + if (getCompressorChart().useRealKappa()) { + kappa = thermoSystem.getGamma(); + } else { + kappa = thermoSystem.getGamma2(); + } + if (useGERG2008 && inStream.getThermoSystem().getNumberOfPhases() == 1) { + double[] gergProps; + gergProps = getThermoSystem().getPhase(0).getProperties_GERG2008(); + actualFlowRate *= gergProps[1] / z_inlet; + kappa = gergProps[14]; + z_inlet = gergProps[1]; + } + + double polytropEff = + getCompressorChart().getPolytropicEfficiency(actualFlowRate, getSpeed()); + setPolytropicEfficiency(polytropEff / 100.0); + polytropicHead = getCompressorChart().getPolytropicHead(actualFlowRate, getSpeed()); + double temperature_inlet = thermoSystem.getTemperature(); + double n = 1.0 / (1.0 - (kappa - 1.0) / kappa * 1.0 / (polytropEff / 100.0)); + polytropicExponent = n; + if (getCompressorChart().getHeadUnit().equals("meter")) { + polytropicFluidHead = polytropicHead / 1000.0 * 9.81; + polytropicHeadMeter = polytropicHead; + } else { + polytropicFluidHead = polytropicHead; + polytropicHeadMeter = polytropicHead * 1000.0 / 9.81; + } + double pressureRatio = Math.pow((polytropicFluidHead * 1000.0 + (n / (n - 1.0) * z_inlet + * ThermodynamicConstantsInterface.R * (temperature_inlet) / MW)) + / (n / (n - 1.0) * z_inlet * ThermodynamicConstantsInterface.R * (temperature_inlet) + / MW), + n / (n - 1.0)); + setOutletPressure(thermoSystem.getPressure() * pressureRatio); + if (getAntiSurge().isActive()) { + logger.info("surge flow " + + getCompressorChart().getSurgeCurve().getSurgeFlow(polytropicHead) + " m3/hr"); + surgeCheck = isSurge(polytropicHead, actualFlowRate); + } + if (getCompressorChart().getStoneWallCurve().isActive()) { + // logger.info("stone wall? " + isStoneWall(polytropicHead, + // thermoSystem.getFlowRate("m3/hr"))); + } + if (surgeCheck && getAntiSurge().isActive()) { + thermoSystem.setTotalFlowRate( + getAntiSurge().getSurgeControlFactor() + * getCompressorChart().getSurgeCurve().getSurgeFlow(polytropicFluidHead), + "Am3/hr"); + thermoSystem.init(3); + fractionAntiSurge = thermoSystem.getTotalNumberOfMoles() / orginalMolarFLow - 1.0; + getAntiSurge().setCurrentSurgeFraction(fractionAntiSurge); + } + + powerSet = true; + dH = polytropicFluidHead * 1000.0 * thermoSystem.getMolarMass() + / getPolytropicEfficiency() * thermoSystem.getTotalNumberOfMoles(); + } while (surgeCheck && getAntiSurge().isActive()); + } } if (usePolytropicCalc) { @@ -1484,4 +1573,12 @@ public void setCompressorChartType(String type) { compressorChart = new CompressorChart(); } } + + public boolean isSolveSpeed() { + return solveSpeed; + } + + public void setSolveSpeed(boolean solveSpeed) { + this.solveSpeed = solveSpeed; + } } diff --git a/src/test/java/neqsim/process/equipment/compressor/CompressorChartTest.java b/src/test/java/neqsim/process/equipment/compressor/CompressorChartTest.java index 553252a691..1fe0f84e67 100644 --- a/src/test/java/neqsim/process/equipment/compressor/CompressorChartTest.java +++ b/src/test/java/neqsim/process/equipment/compressor/CompressorChartTest.java @@ -279,4 +279,87 @@ public void runCurveTest() { Assertions.assertEquals(158.7732888, comp1.getOutletPressure(), 1e-3); } + @Test + public void runCurveTest2() { + SystemInterface testFluid = new SystemPrEos(273.15 + 29.96, 75.73); + + testFluid.addComponent("nitrogen", 0.7823); + testFluid.addComponent("CO2", 1.245); + testFluid.addComponent("methane", 88.4681); + testFluid.addComponent("ethane", 4.7652); + testFluid.addComponent("propane", 2.3669); + testFluid.addComponent("i-butane", 0.3848); + testFluid.addComponent("n-butane", 0.873); + testFluid.setMixingRule("classic"); + + testFluid.setTemperature(29.96, "C"); + testFluid.setPressure(75.73, "bara"); + testFluid.setTotalFlowRate(559401.4, "kg/hr"); + + Stream stream_1 = new Stream("Stream1", testFluid); + stream_1.run(); + + Compressor comp1 = new Compressor("compressor 1", stream_1); + comp1.setCompressorChartType("interpolate and extrapolate"); + comp1.setUsePolytropicCalc(true); + comp1.setPolytropicEfficiency(0.85); + comp1.setSpeed(9000); + double[] chartConditions = new double[] {0.3, 1.0, 1.0, 1.0}; + + double[] speed = new double[] {7000, 7500, 8000, 8500, 9000, 9500, 9659, 10000, 10500}; + + double[][] flow = new double[][] {{4512.7, 5120.8, 5760.9, 6401, 6868.27}, + {4862.47, 5486.57, 6172.39, 6858.21, 7550.89}, + {5237.84, 5852.34, 6583.88, 7315.43, 8046.97, 8266.43}, + {5642.94, 6218.11, 6995.38, 7772.64, 8549.9, 9000.72}, + {6221.77, 6583.88, 7406.87, 8229.85, 9052.84, 9768.84}, + {6888.85, 6949.65, 7818.36, 8687.07, 9555.77, 10424.5, 10546.1}, + {7109.83, 7948.87, 8832.08, 9715.29, 10598.5, 10801.6}, + {7598.9, 8229.85, 9144.28, 10058.7, 10973.1, 11338.9}, + {8334.1, 8641.35, 9601.5, 10561.6, 11521.8, 11963.5}}; + + double[][] head = new double[][] {{61.885, 59.639, 56.433, 52.481, 49.132}, + {71.416, 69.079, 65.589, 61.216, 55.858}, {81.621, 79.311, 75.545, 70.727, 64.867, 62.879,}, + {92.493, 90.312, 86.3, 81.079, 74.658, 70.216}, + {103.512, 102.073, 97.83, 92.254, 85.292, 77.638}, + {114.891, 114.632, 110.169, 104.221, 96.727, 87.002, 85.262}, + {118.595, 114.252, 108.203, 100.55, 90.532, 87.54}, + {126.747, 123.376, 117.113, 109.056, 98.369, 92.632}, + {139.082, 137.398, 130.867, 122.264, 110.548, 103.247},}; + double[][] polyEff = new double[][] {{78.3, 78.2, 77.2, 75.4, 73.4}, + + {78.3, 78.3, 77.5, 75.8, 73}, {78.2, 78.4, 77.7, 76.1, 73.5, 72.5}, + {78.2, 78.4, 77.9, 76.4, 74, 71.9}, {78.3, 78.4, 78, 76.7, 74.5, 71.2}, + {78.3, 78.4, 78.1, 77, 74.9, 71.3, 70.5}, {78.4, 78.1, 77.1, 75, 71.4, 70.2}, + {78.3, 78.2, 77.2, 75.2, 71.7, 69.5}, {78.2, 78.2, 77.3, 75.5, 72.2, 69.6}}; + + comp1.getCompressorChart().setCurves(chartConditions, speed, flow, head, polyEff); + comp1.getCompressorChart().setHeadUnit("kJ/kg"); + + double[] surgeflow = + new double[] {4512.7, 4862.47, 5237.84, 5642.94, 6221.77, 6888.85, 7109.83, 7598.9}; + double[] surgehead = + new double[] {61.885, 71.416, 81.621, 92.493, 103.512, 114.891, 118.595, 126.747}; + comp1.getCompressorChart().getSurgeCurve().setCurve(chartConditions, surgeflow, surgehead); + // comp1.getAntiSurge().setActive(true); + comp1.getAntiSurge().setSurgeControlFactor(1.0); + comp1.getCompressorChart().setUseCompressorChart(true); + comp1.setOutletPressure(220.0, "bara"); + comp1.setSolveSpeed(true); + comp1.run(); + + org.apache.logging.log4j.LogManager.getLogger(CompressorChartTest.class) + .debug("feed flow " + (comp1.getInletStream().getFlowRate("m3/hr"))); + org.apache.logging.log4j.LogManager.getLogger(CompressorChartTest.class) + .debug("out pressure " + (comp1.getOutletStream().getPressure("bara"))); + org.apache.logging.log4j.LogManager.getLogger(CompressorChartTest.class) + .debug("power " + comp1.getPower("MW")); + org.apache.logging.log4j.LogManager.getLogger(CompressorChartTest.class) + .debug("polytropic head " + comp1.getPolytropicFluidHead()); + org.apache.logging.log4j.LogManager.getLogger(CompressorChartTest.class) + .debug("polytropic efficiency " + comp1.getPolytropicEfficiency()); + org.apache.logging.log4j.LogManager.getLogger(CompressorChartTest.class) + .debug("speed " + comp1.getSpeed()); + } + } From f217ee87042ce2355a39566a79f848a3b43c0843 Mon Sep 17 00:00:00 2001 From: Even Solbraa <41290109+EvenSol@users.noreply.github.com> Date: Sun, 24 Nov 2024 16:53:44 +0000 Subject: [PATCH 2/3] update --- .../equipment/compressor/Compressor.java | 2 +- .../compressor/CompressorChartTest.java | 30 ++++++++++++++++--- 2 files changed, 27 insertions(+), 5 deletions(-) diff --git a/src/main/java/neqsim/process/equipment/compressor/Compressor.java b/src/main/java/neqsim/process/equipment/compressor/Compressor.java index ed25de3c81..5369d1ec20 100644 --- a/src/main/java/neqsim/process/equipment/compressor/Compressor.java +++ b/src/main/java/neqsim/process/equipment/compressor/Compressor.java @@ -414,7 +414,7 @@ public void run(UUID id) { double maxSpeed = getMaximumSpeed(); // Maximum speed for the compressor double currentSpeed = getSpeed(); // Initial guess for speed double maxIterations = 100; // Maximum number of iterations - double deltaSpeed = 10.0; // Small increment for numerical derivative + double deltaSpeed = 100.0; // Small increment for numerical derivative int iteration = 0; while (iteration < maxIterations) { diff --git a/src/test/java/neqsim/process/equipment/compressor/CompressorChartTest.java b/src/test/java/neqsim/process/equipment/compressor/CompressorChartTest.java index 1fe0f84e67..2f3f17cd6c 100644 --- a/src/test/java/neqsim/process/equipment/compressor/CompressorChartTest.java +++ b/src/test/java/neqsim/process/equipment/compressor/CompressorChartTest.java @@ -292,11 +292,10 @@ public void runCurveTest2() { testFluid.addComponent("n-butane", 0.873); testFluid.setMixingRule("classic"); - testFluid.setTemperature(29.96, "C"); - testFluid.setPressure(75.73, "bara"); - testFluid.setTotalFlowRate(559401.4, "kg/hr"); - Stream stream_1 = new Stream("Stream1", testFluid); + stream_1.setTemperature(29.96, "C"); + stream_1.setPressure(75.73, "bara"); + stream_1.setFlowRate(559401.4, "kg/hr"); stream_1.run(); Compressor comp1 = new Compressor("compressor 1", stream_1); @@ -360,6 +359,29 @@ public void runCurveTest2() { .debug("polytropic efficiency " + comp1.getPolytropicEfficiency()); org.apache.logging.log4j.LogManager.getLogger(CompressorChartTest.class) .debug("speed " + comp1.getSpeed()); + + stream_1.setFlowRate(309401.4, "kg/hr"); + stream_1.run(); + comp1.setOutletPressure(170.0, "bara"); + comp1.run(); + org.apache.logging.log4j.LogManager.getLogger(CompressorChartTest.class) + .debug("feed flow " + (comp1.getInletStream().getFlowRate("m3/hr"))); + org.apache.logging.log4j.LogManager.getLogger(CompressorChartTest.class) + .debug("out pressure " + (comp1.getOutletStream().getPressure("bara"))); + org.apache.logging.log4j.LogManager.getLogger(CompressorChartTest.class) + .debug("power " + comp1.getPower("MW")); + org.apache.logging.log4j.LogManager.getLogger(CompressorChartTest.class) + .debug("polytropic head " + comp1.getPolytropicFluidHead()); + org.apache.logging.log4j.LogManager.getLogger(CompressorChartTest.class) + .debug("polytropic efficiency " + comp1.getPolytropicEfficiency()); + org.apache.logging.log4j.LogManager.getLogger(CompressorChartTest.class) + .debug("speed " + comp1.getSpeed()); + org.apache.logging.log4j.LogManager.getLogger(CompressorChartTest.class) + .debug("dist to surge " + comp1.getDistanceToSurge()); + org.apache.logging.log4j.LogManager.getLogger(CompressorChartTest.class) + .debug("surge flow rate margin " + comp1.getSurgeFlowRateMargin()); + org.apache.logging.log4j.LogManager.getLogger(CompressorChartTest.class) + .debug("surge flow rate " + comp1.getSurgeFlowRate()); } } From 6f5996e8d9d6a6281388fc1d59ec0ba4e3c0a1f4 Mon Sep 17 00:00:00 2001 From: Even Solbraa <41290109+EvenSol@users.noreply.github.com> Date: Sun, 24 Nov 2024 17:20:28 +0000 Subject: [PATCH 3/3] update --- .../processmodel/CompressorModule.java | 33 +++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/src/test/java/neqsim/process/processmodel/CompressorModule.java b/src/test/java/neqsim/process/processmodel/CompressorModule.java index a8d3260903..f83bb6e284 100644 --- a/src/test/java/neqsim/process/processmodel/CompressorModule.java +++ b/src/test/java/neqsim/process/processmodel/CompressorModule.java @@ -1,6 +1,7 @@ package neqsim.process.processmodel; import static org.junit.jupiter.api.Assertions.assertEquals; +import static org.junit.jupiter.api.Assertions.assertFalse; import static org.junit.jupiter.api.Assertions.assertTrue; import org.junit.jupiter.api.Test; import neqsim.process.equipment.compressor.Compressor; @@ -154,7 +155,39 @@ public void testProcess() { feedStream.setFlowRate(204094, "kg/hr"); operations.run(); + assertTrue(seccondStageCompressor.isSurge(seccondStageCompressor.getPolytropicFluidHead(), seccondStageCompressor.getInletStream().getFlowRate("m3/hr"))); + + + double pressurespeedclac = seccondStageCompressor.getOutletPressure(); + double speedcomp = seccondStageCompressor.getSpeed(); + + seccondStageCompressor.setSolveSpeed(true); + operations.run(); + + assertEquals(pressurespeedclac, seccondStageCompressor.getOutletStream().getPressure("bara"), + 0.5); + assertEquals(speedcomp, seccondStageCompressor.getSpeed(), 0.5); + assertEquals(259.8380255, seccondStageCompressor.getInletStream().getFlowRate("m3/hr"), 0.2); + + feedStream.setFlowRate(304094, "kg/hr"); + operations.run(); + + assertEquals(pressurespeedclac, seccondStageCompressor.getOutletStream().getPressure("bara"), + 0.5); + assertEquals(3526, seccondStageCompressor.getSpeed(), 10); + assertEquals(386.5, seccondStageCompressor.getInletStream().getFlowRate("m3/hr"), 0.2); + assertTrue(seccondStageCompressor.isSurge(seccondStageCompressor.getPolytropicFluidHead(), + seccondStageCompressor.getInletStream().getFlowRate("m3/hr"))); + + feedStream.setFlowRate(704094, "kg/hr"); + operations.run(); + + assertEquals(pressurespeedclac, seccondStageCompressor.getOutletStream().getPressure("bara"), + 0.5); + assertEquals(4177, seccondStageCompressor.getSpeed(), 10); + assertFalse(seccondStageCompressor.isSurge(seccondStageCompressor.getPolytropicFluidHead(), + seccondStageCompressor.getInletStream().getFlowRate("m3/hr"))); } }