Class GibbsReactor

java.lang.Object
neqsim.util.NamedBaseClass
neqsim.process.SimulationBaseClass
neqsim.process.equipment.ProcessEquipmentBaseClass
neqsim.process.equipment.TwoPortEquipment
neqsim.process.equipment.reactor.GibbsReactor
All Implemented Interfaces:
Serializable, Runnable, ProcessEquipmentInterface, TwoPortInterface, SimulationInterface, NamedInterface
public class GibbsReactor extends TwoPortEquipment
Gibbs reactor for chemical equilibrium calculations using Gibbs free energy minimization.

This reactor computes chemical equilibrium compositions by minimizing the total Gibbs free energy of the system subject to elemental mass balance constraints. The implementation uses the Newton-Raphson method with Lagrange multipliers to solve the constrained optimization problem.

Key Features

  • Supports both isothermal and adiabatic operation modes
  • Handles multi-component systems with element balance constraints
  • Allows marking specific components as inert (excluded from reactions)
  • Provides detailed convergence diagnostics and mass balance verification

Algorithm

The reactor minimizes the objective function: 

G = Σ nᵢ(μᵢ⁰ + RT ln(φᵢyᵢP)) - Σ λⱼ(Σ aᵢⱼnᵢ - bⱼ)

where nᵢ are molar amounts, μᵢ⁰ is standard chemical potential, φᵢ is fugacity coefficient, yᵢ is mole fraction, λⱼ are Lagrange multipliers, aᵢⱼ are stoichiometric coefficients, and bⱼ are element totals.

Usage Example

GibbsReactor reactor = new GibbsReactor("reactor", inletStream);
reactor.setEnergyMode(EnergyMode.ISOTHERMAL);
reactor.setMaxIterations(5000);
reactor.setConvergenceTolerance(1e-6);
reactor.setDampingComposition(0.01);
reactor.run();

if (reactor.hasConverged()) {
  Stream outlet = reactor.getOutletStream();
  double enthalpyChange = reactor.getEnthalpyOfReactions();
}
Version:
$Id: $Id
Author:
Sviatoslav Eroshkin
See Also:

  • Field Details

    • tempFugacitySystem

      private transient ThreadLocal<SystemInterface> tempFugacitySystem

    • serialVersionUID

      private static final long serialVersionUID
      Serialization version UID.
      See Also:

    • REFERENCE_TEMPERATURE

      private static final double REFERENCE_TEMPERATURE
      Reference temperature in Kelvin for thermodynamic calculations (298.15 K = 25°C).
      See Also:

    • R_KJ

      private static final double R_KJ
      Universal gas constant in kJ/(mol·K).
      See Also:

    • MIN_MOLES

      private static final double MIN_MOLES
      Minimum mole amount to prevent numerical issues with log(0).
      See Also:

    • MIN_JACOBIAN_MOLES

      private static final double MIN_JACOBIAN_MOLES
      Minimum value for Jacobian calculations to avoid division by zero.
      See Also:

    • ELEMENT_ZERO_THRESHOLD

      private static final double ELEMENT_ZERO_THRESHOLD
      Small threshold for detecting zero element coefficients.
      See Also:

    • logger

      private static final org.apache.logging.log4j.Logger logger
      Logger object for class.

    • ION_NAME_PATTERN

      private static final Pattern ION_NAME_PATTERN
      Pattern to identify ionic species (names ending with + or -).

    • energyMode

      private GibbsReactor.EnergyMode energyMode

    • componentAliases

      private static final Map<String,String> componentAliases

    • method

      private String method

    • useAllDatabaseSpecies

      private boolean useAllDatabaseSpecies

    • gibbsDatabase

      private List<GibbsReactor.GibbsComponent> gibbsDatabase

    • componentMap

      private Map<String, GibbsReactor.GibbsComponent> componentMap

    • lambda

      private double[] lambda

    • lagrangeContributions

      private Map<String,Double> lagrangeContributions

    • elementNames

      private String[] elementNames

    • processedComponents

      private List<String> processedComponents

    • objectiveFunctionValues

      private Map<String,Double> objectiveFunctionValues

    • inertComponents

      private final Set<String> inertComponents

    • initialMoles

      private Map<String,Double> initialMoles

    • finalMoles

      private Map<String,Double> finalMoles

    • elementMoleBalanceIn

      private double[] elementMoleBalanceIn

    • elementMoleBalanceOut

      private double[] elementMoleBalanceOut

    • elementMoleBalanceDiff

      private double[] elementMoleBalanceDiff

    • inlet_mole

      private List<Double> inlet_mole

    • outlet_mole

      private List<Double> outlet_mole

    • processedComponentIndexMap

      private Map<String,Integer> processedComponentIndexMap

    • objectiveMinimizationVector

      private double[] objectiveMinimizationVector

    • objectiveMinimizationVectorLabels

      private List<String> objectiveMinimizationVectorLabels

    • jacobianMatrix

      private double[][] jacobianMatrix

    • jacobianInverse

      private double[][] jacobianInverse

    • jacobianRowLabels

      private List<String> jacobianRowLabels

    • jacobianColLabels

      private List<String> jacobianColLabels

    • variableComponents

      private List<String> variableComponents

    • maxIterations

      private int maxIterations

    • convergenceTolerance

      private double convergenceTolerance

    • dampingComposition

      private double dampingComposition

    • actualIterations

      private int actualIterations

    • converged

      private boolean converged

    • finalConvergenceError

      private double finalConvergenceError

    • system

      private SystemInterface system

    • inletEnthalpy

      private double inletEnthalpy

    • outletEnthalpy

      private double outletEnthalpy

    • dT

      private double dT

    • tempUpdateIter

      private int tempUpdateIter

    • enthalpyOld

      double enthalpyOld

    • enthalpyOfReactions

      private double enthalpyOfReactions

    • temperatureChange

      private double temperatureChange

    • deltaNorm

      private double deltaNorm

    • GOLD

      private double GOLD

    • G

      private double G

    • dG

      private double dG

    • formulaToComponent

      private static final Map<String,String> formulaToComponent

  • Constructor Details

    • GibbsReactor

      public GibbsReactor(String name)
      Constructor for GibbsReactor.
      Parameters:
      name - Name of GibbsReactor

    • GibbsReactor

      public GibbsReactor(String name, StreamInterface stream)
      Constructor for GibbsReactor.
      Parameters:
      name - Name of GibbsReactor
      stream - Stream to set as inlet Stream. A clone of stream is set as outlet stream.

  • Method Details

    • getMassBalanceError

      public double getMassBalanceError()
      Get the absolute mass balance error (difference between inlet and outlet) in kg/sec.
      Returns:
      absolute difference in total mass flow rate (kg/sec)

    • getMassBalanceConverged

      public boolean getMassBalanceConverged()
      Returns true if the absolute mass balance error is less than 1e-3 kg/sec.
      Returns:
      true if mass balance is converged, false otherwise

    • ensureTempFugacitySystemInitialized

      private void ensureTempFugacitySystemInitialized()
      Ensures tempFugacitySystem is initialized (handles deserialization case).

    • getEnthalpyOfReactions

      public double getEnthalpyOfReactions()
      Get the cumulative enthalpy of reaction (sum of dH for all iterations).
      Returns:
      enthalpyOfReactions in kJ

    • getTemperatureChange

      public double getTemperatureChange()
      Get the cumulative temperature change during the reaction (sum of dT for all iterations).
      Returns:
      temperatureChange in K

    • getPower

      public double getPower()
      Get the reactor power (default: W, negative enthalpyOfReactions*1000).
      Returns:
      Power in Watts (W)

    • getPower

      public double getPower(String unit)
      Get the reactor power in the specified unit ("W", "kW", or "MW").
      Parameters:
      unit - Power unit: "W", "kW", or "MW" (case-insensitive, default is W)
      Returns:
      Power in the specified unit

    • calculateMixtureEnthalpy

      public double calculateMixtureEnthalpy(List<String> componentNames, List<Double> n, double T, Map<String, GibbsReactor.GibbsComponent> componentMap)
      Calculate the total enthalpy of a mixture: sum_i n_i * enthalpy_i(T).
      Parameters:
      componentNames - List of component names (order matches n_i)
      n - List of moles for each component
      T - Temperature in K
      componentMap - Map from component name (lowercase) to GibbsComponent
      Returns:
      Total enthalpy (kJ)

    • calculateMixtureGibbsEnergy

      public double calculateMixtureGibbsEnergy(List<String> componentNames, List<Double> n, Map<String, GibbsReactor.GibbsComponent> componentMap, double T)
      Calculate the total Gibbs energy of a mixture: sum_i n_i * gibbs_i(T).
      Parameters:
      componentNames - List of component names (order matches n_i)
      n - List of moles for each component
      componentMap - Map from component name (lowercase) to GibbsComponent
      T - Temperature in K
      Returns:
      Total Gibbs energy (kJ)

    • getMixtureEnthalpy

      public double getMixtureEnthalpy()
      Get the last calculated mixture enthalpy (kJ).
      Returns:
      a double

    • getMixtureGibbsEnergy

      public double getMixtureGibbsEnergy()
      Get the last calculated mixture Gibbs energy (kJ).
      Returns:
      a double

    • calculateMixtureEnthalpyStandard

      public double calculateMixtureEnthalpyStandard(List<String> componentNames, List<Double> n, Map<String, GibbsReactor.GibbsComponent> componentMap)
      Calculate the total standard enthalpy of a mixture: sum_i n_i * enthalpy_i(T).
      Parameters:
      componentNames - List of component names (order matches n_i)
      n - List of moles for each component
      componentMap - Map from component name (lowercase) to GibbsComponent
      Returns:
      Total enthalpy (kJ)

    • calculateMixtureEnthalpy

      public double calculateMixtureEnthalpy(List<String> componentNames, List<Double> n, Map<String, GibbsReactor.GibbsComponent> componentMap, double T)
      Calculate the total standard enthalpy of a mixture: sum_i n_i * enthalpy_i(T).
      Parameters:
      componentNames - List of component names (order matches n_i)
      n - List of moles for each component
      componentMap - Map from component name (lowercase) to GibbsComponent
      T - a double
      Returns:
      Total enthalpy (kJ)

    • setEnergyMode

      public void setEnergyMode(GibbsReactor.EnergyMode mode)
      Set the energy mode of the reactor (isothermal or adiabatic).
      Parameters:
      mode - EnergyMode.ISOTHERMAL or EnergyMode.ADIABATIC

    • setEnergyMode

      public void setEnergyMode(String mode)
      Set the energy mode of the reactor using a string (case-insensitive). Accepts "adiabatic" or "isothermal" (case-insensitive).
      Parameters:
      mode - String representing the energy mode
      Throws:
      IllegalArgumentException - if the mode is not recognized

    • getEnergyMode

      public GibbsReactor.EnergyMode getEnergyMode()
      Gets the current energy mode of the reactor.
      Returns:
      the current EnergyMode (ISOTHERMAL or ADIABATIC)

    • setComponentAsInert

      public void setComponentAsInert(String componentName)
      Mark a component as inert by name.
      Parameters:
      componentName - component name as in the thermo system

    • setComponentAsInert

      public void setComponentAsInert(int index)
      Mark a component as inert by index in inlet system.
      Parameters:
      index - component index in inlet system

    • isComponentInert

      public boolean isComponentInert(String componentName)
      Check if a component is inert.
      Parameters:
      componentName - the name of the component to check
      Returns:
      true if the component is inert, false otherwise

    • isIonicComponent

      private boolean isIonicComponent(String moleculeName)

    • loadGibbsDatabase

      private void loadGibbsDatabase()
      Load the Gibbs reaction database from resources.

    • setUseAllDatabaseSpecies

      public void setUseAllDatabaseSpecies(boolean useAllDatabaseSpecies)
      Set whether to use all database species or only species in the system.
      Parameters:
      useAllDatabaseSpecies - true to use all database species, false to use only system species

    • getUseAllDatabaseSpecies

      public boolean getUseAllDatabaseSpecies()
      Get whether using all database species or only species in the system.
      Returns:
      true if using all database species, false if using only system species

    • getMethod

      public String getMethod()
      Get the method used for Gibbs minimization.
      Returns:
      The method name

    • setMethod

      public void setMethod(String method)
      Set the method used for Gibbs minimization.
      Parameters:
      method - The method name

    • getLagrangeContributions

      public Map<String,Double> getLagrangeContributions()
      Get the Lagrange contributions (legacy method).
      Returns:
      Map of Lagrange contributions

    • run

      public void run(UUID id)

      In this method all thermodynamic and unit operations will be calculated in a steady state calculation.

      Parameters:
      id - UUID

    • performGibbsMinimization

      private void performGibbsMinimization(SystemInterface system)
      Perform Gibbs free energy minimization.
      Parameters:
      system - The thermodynamic system

    • calculateElementMoleBalance

      private void calculateElementMoleBalance(SystemInterface system, double[] elementBalance, boolean isInput)
      Calculate element mole balance for a system.
      Parameters:
      system - The thermodynamic system
      elementBalance - Array to store the element balance
      isInput - true if this is input balance, false if output balance

    • calculateObjectiveFunctionValues

      private void calculateObjectiveFunctionValues(SystemInterface system)
      Calculate objective function values for each component.
      Parameters:
      system - The thermodynamic system

    • setLagrangeMultiplier

      public void setLagrangeMultiplier(int index, double value)
      Set a Lagrange multiplier value.
      Parameters:
      index - The element index (0=O, 1=N, 2=C, 3=H, 4=S, 5=Ar)
      value - The Lagrange multiplier value

    • getLagrangianMultipliers

      public double[] getLagrangianMultipliers()
      Get the Lagrange multiplier values.
      Returns:
      Array of Lagrange multiplier values [O, N, C, H, S, Ar]

    • getElementNames

      public String[] getElementNames()
      Get the element names.
      Returns:
      Array of element names

    • getLagrangeMultiplierContributions

      public Map<String, Map<String,Double>> getLagrangeMultiplierContributions()
      Get the Lagrange multiplier contributions for each component.
      Returns:
      Map with component names and their Lagrange multiplier contributions

    • getObjectiveFunctionValues

      public Map<String,Double> getObjectiveFunctionValues()
      Get the objective function values for each component.
      Returns:
      Map with component names and their objective function values

    • getElementMoleBalanceIn

      public double[] getElementMoleBalanceIn()
      Get the element mole balance for input streams.
      Returns:
      Array of total moles for each element in input [O, N, C, H, S, Ar]

    • getElementMoleBalanceOut

      public double[] getElementMoleBalanceOut()
      Get the element mole balance for output streams.
      Returns:
      Array of total moles for each element in output [O, N, C, H, S, Ar]

    • getElementMoleBalanceDiff

      public double[] getElementMoleBalanceDiff()
      Get the element mole balance difference (in - out).
      Returns:
      Array of mole differences for each element [O, N, C, H, S, Ar]

    • getDetailedMoleBalance

      public Map<String, Map<String,Double>> getDetailedMoleBalance()
      Get detailed mole balance information for each component.
      Returns:
      Map with component names and their element contributions to mole balance

    • calculateObjectiveMinimizationVector

      private void calculateObjectiveMinimizationVector()
      Calculate the objective minimization vector. This vector contains the F values for each component and the mass balance constraints. The system is in equilibrium when this vector is zero. Only includes elements that are actually present in the system.

    • getObjectiveMinimizationVector

      public double[] getObjectiveMinimizationVector()
      Get the objective minimization vector.
      Returns:
      The objective minimization vector

    • getObjectiveVectorForVariables

      private double[] getObjectiveVectorForVariables()
      Build objective vector matching the optimization variables (variableComponents) and active element balances. This excludes inert components from the variable part but keeps balances.
      Returns:
      the objective vector for optimization variables

    • getObjectiveMinimizationVectorLabels

      public List<String> getObjectiveMinimizationVectorLabels()
      Get the labels for the objective minimization vector.
      Returns:
      List of labels for each element in the objective minimization vector

    • calculateJacobian

      private void calculateJacobian()
      Calculate the Jacobian matrix for the Newton-Raphson method. The Jacobian represents the derivatives of the objective function with respect to the variables. Only includes elements that are actually present in the system to avoid singular matrices.

    • getJacobianMatrix

      public double[][] getJacobianMatrix()
      Get the Jacobian matrix.
      Returns:
      The Jacobian matrix

    • getJacobianRowLabels

      public List<String> getJacobianRowLabels()
      Get the Jacobian row labels.
      Returns:
      List of row labels for the Jacobian matrix

    • getJacobianColLabels

      public List<String> getJacobianColLabels()
      Get the Jacobian column labels.
      Returns:
      List of column labels for the Jacobian matrix

    • getJacobianInverse

      public double[][] getJacobianInverse()
      Get the Jacobian inverse matrix.
      Returns:
      The Jacobian inverse matrix, or null if it couldn't be calculated

    • calculateJacobianInverse

      private double[][] calculateJacobianInverse()
      Calculate the inverse of the Jacobian matrix using JAMA Matrix library.
      Returns:
      Inverse matrix, or null if matrix is singular

    • enforceMinimumConcentrations

      private void enforceMinimumConcentrations(SystemInterface system)
      Enforce minimum concentration threshold to prevent numerical issues. If any component has moles less than MIN_MOLES, it will be set to MIN_MOLES.
      Parameters:
      system - The thermodynamic system to check and modify

    • getInletMole

      public List<Double> getInletMole()
      Get the inlet mole list.
      Returns:
      List of inlet moles for each component

    • getOutletMole

      public List<Double> getOutletMole()
      Get the outlet mole list.
      Returns:
      List of outlet moles for each component (with 0 values replaced by 1E-6)

    • getInletMoles

      public List<Double> getInletMoles()
      Get the inlet mole list for debugging.
      Returns:
      List of inlet moles

    • getOutletMoles

      public List<Double> getOutletMoles()
      Get the outlet mole list for debugging.
      Returns:
      List of outlet moles

    • printDatabaseComponents

      public void printDatabaseComponents()
      Print loaded database components for debugging.

    • findActiveElements

      private List<Integer> findActiveElements()
      Find which elements are actually present in the system (have non-zero coefficients).
      Returns:
      List of indices of active elements

    • getActiveElementIndices

      private List<Integer> getActiveElementIndices()
      Get the indices of active elements (elements that have non-zero coefficients in any component).
      Returns:
      List of active element indices

    • verifyJacobianInverse

      public boolean verifyJacobianInverse()
      Verify that the Jacobian inverse is correct by multiplying J * J^-1. Should return the identity matrix if the inverse is correct.
      Returns:
      True if the inverse is correct (within tolerance)

    • performNewtonRaphsonIteration

      public double[] performNewtonRaphsonIteration()
      Perform one Newton-Raphson iteration step to calculate the delta vector (dX). Uses the formula: dX = -J^(-1) * F where J is the Jacobian matrix and F is the objective function vector.
      Returns:
      The delta vector (dX) for updating variables, or null if calculation fails

    • performIterationUpdate

      public boolean performIterationUpdate(double[] deltaX, double alphaComposition)
      Perform a Newton-Raphson iteration update. Updates outlet compositions with damping factor and Lagrange multipliers directly.
      Parameters:
      deltaX - The delta vector from Newton-Raphson iteration
      alphaComposition - Damping factor for composition updates (e.g., 0.0001)
      Returns:
      True if update was successful, false otherwise

    • updateSystemWithNewCompositions

      private boolean updateSystemWithNewCompositions()
      Update the thermodynamic system with the new outlet compositions.
      Returns:
      True if update was successful, false otherwise

    • getFugacityCoefficient

      public double[] getFugacityCoefficient(Object phaseNameOrIndex)
      Get the fugacity coefficient array for all components in a specified phase using the current outlet composition. Uses direct phase composition assignment for efficiency.
      Parameters:
      phaseNameOrIndex - Name or index of the phase (e.g., "gas", "oil", "aqueous", or 0/1/2)
      Returns:
      Fugacity coefficient (phi) array for all components in the specified phase, or Double.NaN if not found

    • setMaxIterations

      public void setMaxIterations(int maxIterations)
      Set maximum number of Newton-Raphson iterations.
      Parameters:
      maxIterations - Maximum number of iterations

    • getMaxIterations

      public int getMaxIterations()
      Get maximum number of Newton-Raphson iterations.
      Returns:
      Maximum number of iterations

    • setConvergenceTolerance

      public void setConvergenceTolerance(double convergenceTolerance)
      Set convergence tolerance for Newton-Raphson iterations.
      Parameters:
      convergenceTolerance - Convergence tolerance

    • getConvergenceTolerance

      public double getConvergenceTolerance()
      Get convergence tolerance for Newton-Raphson iterations.
      Returns:
      Convergence tolerance

    • setDampingComposition

      public void setDampingComposition(double dampingComposition)
      Set damping factor for composition updates in Newton-Raphson iterations.
      Parameters:
      dampingComposition - Damping factor (typically 0.0001 to 0.01)

    • getDampingComposition

      public double getDampingComposition()
      Get damping factor for composition updates.
      Returns:
      Damping factor for composition updates

    • getActualIterations

      public int getActualIterations()
      Get actual number of iterations performed.
      Returns:
      Actual iterations performed

    • hasConverged

      public boolean hasConverged()
      Check if iterations converged.
      Returns:
      true if converged, false otherwise

    • getFinalConvergenceError

      public double getFinalConvergenceError()
      Get final convergence error.
      Returns:
      Final convergence error

    • solveGibbsEquilibrium

      public boolean solveGibbsEquilibrium(double alphaComposition)
      Solve Gibbs equilibrium using Newton-Raphson iterations with specified step size.
      Parameters:
      alphaComposition - Step size for composition updates
      Returns:
      true if converged, false otherwise

    • solveGibbsEquilibrium

      public boolean solveGibbsEquilibrium()
      Solve Gibbs equilibrium using Newton-Raphson iterations with default damping factor.
      Returns:
      true if converged, false otherwise