Skip to content

Commit

Permalink
solve speed method for compressor
Browse files Browse the repository at this point in the history
  • Loading branch information
@EvenSol
EvenSol committed Nov 24, 2024
1 parent 6564c83 commit 7acb154
Show file tree
Hide file tree
Showing 2 changed files with 236 additions and 56 deletions.
209 changes: 153 additions & 56 deletions src/main/java/neqsim/process/equipment/compressor/Compressor.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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;
}
}
83 changes: 83 additions & 0 deletions src/test/java/neqsim/process/equipment/compressor/CompressorChartTest.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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());
}

}

0 comments on commit 7acb154

Please sign in to comment.