Class TimeIntegrator

java.lang.Object
neqsim.process.equipment.pipeline.twophasepipe.numerics.TimeIntegrator
All Implemented Interfaces:
Serializable
public class TimeIntegrator extends Object implements Serializable
Time integration methods for the two-fluid transient pipe model.

Provides explicit time stepping algorithms with CFL-based adaptive time step control. The primary method is the classical 4th-order Runge-Kutta (RK4) scheme.

RK4 Algorithm

For dU/dt = R(U):

  • k1 = R(U^n)
  • k2 = R(U^n + 0.5*dt*k1)
  • k3 = R(U^n + 0.5*dt*k2)
  • k4 = R(U^n + dt*k3)
  • U^{n+1} = U^n + dt/6 * (k1 + 2*k2 + 2*k3 + k4)

CFL Condition

The time step is limited by: dt ≤ CFL * dx / (|v| + c) where c is the sound speed and CFL is typically 0.5-0.9 for explicit schemes.

Version:
1.0
Author:
Even Solbraa
See Also:

  • Field Details

    • serialVersionUID

      private static final long serialVersionUID
      See Also:

    • method

      private TimeIntegrator.Method method
      Current integration method.

    • cflNumber

      private double cflNumber
      CFL number for time step control (0 < CFL < 1).

    • minTimeStep

      private double minTimeStep
      Minimum allowed time step (s).

    • maxTimeStep

      private double maxTimeStep
      Maximum allowed time step (s).

    • currentTime

      private double currentTime
      Current simulation time (s).

    • currentDt

      private double currentDt
      Current time step (s).

    • cellSoundSpeeds

      private double[] cellSoundSpeeds
      Cell-averaged sound speeds for IMEX pressure solve. Set by caller before step.

    • cellDensities

      private double[] cellDensities
      Cell-averaged mixture densities for IMEX pressure solve. Set by caller before step.

    • imexDx

      private double imexDx
      Cell size for IMEX pressure solve (m). Set by caller before step.

    • imexOutletPressure

      private double imexOutletPressure
      Pressure boundary (Pa) at outlet for IMEX. Set by caller before step.

    • imexOutletPressureFixed

      private boolean imexOutletPressureFixed
      Whether outlet pressure BC is fixed for IMEX.

  • Constructor Details

    • TimeIntegrator

      public TimeIntegrator()
      Default constructor.

    • TimeIntegrator

      public TimeIntegrator(TimeIntegrator.Method method)
      Constructor with method selection.
      Parameters:
      method - Integration method

  • Method Details

    • step

      public double[][] step(double[][] U, TimeIntegrator.RHSFunction rhs, double dt)
      Advance solution by one time step using selected method.
      Parameters:
      U - Current state [nCells][nVars]
      rhs - Right-hand side function
      dt - Time step
      Returns:
      Updated state at t + dt

    • stepEuler

      public double[][] stepEuler(double[][] U, TimeIntegrator.RHSFunction rhs, double dt)
      Forward Euler method (first-order).

      U^{n+1} = U^n + dt * R(U^n)

      Parameters:
      U - Current state
      rhs - Right-hand side function
      dt - Time step
      Returns:
      Updated state

    • stepRK2

      public double[][] stepRK2(double[][] U, TimeIntegrator.RHSFunction rhs, double dt)
      Second-order Runge-Kutta (Heun's method).
      Parameters:
      U - Current state
      rhs - Right-hand side function
      dt - Time step
      Returns:
      Updated state

    • stepRK4

      public double[][] stepRK4(double[][] U, TimeIntegrator.RHSFunction rhs, double dt)
      Classical 4th-order Runge-Kutta.
      Parameters:
      U - Current state
      rhs - Right-hand side function
      dt - Time step
      Returns:
      Updated state

    • stepSSPRK3

      public double[][] stepSSPRK3(double[][] U, TimeIntegrator.RHSFunction rhs, double dt)
      Strong Stability Preserving RK3 (Shu-Osher).

      Maintains TVD property, good for problems with shocks.

      Parameters:
      U - Current state
      rhs - Right-hand side function
      dt - Time step
      Returns:
      Updated state

    • calcStableTimeStep

      public double calcStableTimeStep(double maxWaveSpeed, double dx)
      Calculate stable time step based on CFL condition.
      Parameters:
      maxWaveSpeed - Maximum wave speed in the domain (|v| + c)
      dx - Minimum cell size
      Returns:
      Stable time step

    • calcTwoFluidTimeStep

      public double calcTwoFluidTimeStep(double[] gasVelocities, double[] liquidVelocities, double[] gasSoundSpeeds, double[] liquidSoundSpeeds, double dx)
      Calculate stable time step for two-fluid system.
      Parameters:
      gasVelocities - Gas velocities at each cell (m/s)
      liquidVelocities - Liquid velocities at each cell (m/s)
      gasSoundSpeeds - Gas sound speeds at each cell (m/s)
      liquidSoundSpeeds - Liquid sound speeds at each cell (m/s)
      dx - Cell size (m)
      Returns:
      Stable time step

    • advanceTime

      public void advanceTime(double dt)
      Advance current time by dt.
      Parameters:
      dt - Time step taken

    • addArrays

      private double[][] addArrays(double[][] A, double[][] B)
      Add two 2D arrays element-wise.
      Parameters:
      A - First array
      B - Second array
      Returns:
      A + B

    • scaleArray

      private double[][] scaleArray(double[][] A, double s)
      Scale a 2D array by a scalar.
      Parameters:
      A - Array
      s - Scalar
      Returns:
      s * A

    • getMethod

      public TimeIntegrator.Method getMethod()
      Get integration method.
      Returns:
      Current method

    • setMethod

      public void setMethod(TimeIntegrator.Method method)
      Set integration method.
      Parameters:
      method - New method

    • getCflNumber

      public double getCflNumber()
      Get CFL number.
      Returns:
      CFL number

    • setCflNumber

      public void setCflNumber(double cflNumber)
      Set CFL number.
      Parameters:
      cflNumber - New CFL number (0 < CFL < 1)

    • getMinTimeStep

      public double getMinTimeStep()
      Get minimum time step.
      Returns:
      Minimum time step (s)

    • setMinTimeStep

      public void setMinTimeStep(double minTimeStep)
      Set minimum time step.
      Parameters:
      minTimeStep - Minimum time step (s)

    • getMaxTimeStep

      public double getMaxTimeStep()
      Get maximum time step.
      Returns:
      Maximum time step (s)

    • setMaxTimeStep

      public void setMaxTimeStep(double maxTimeStep)
      Set maximum time step.
      Parameters:
      maxTimeStep - Maximum time step (s)

    • getCurrentTime

      public double getCurrentTime()
      Get current simulation time.
      Returns:
      Current time (s)

    • setCurrentTime

      public void setCurrentTime(double currentTime)
      Set current simulation time.
      Parameters:
      currentTime - Current time (s)

    • getCurrentDt

      public double getCurrentDt()
      Get current time step.
      Returns:
      Current time step (s)

    • setIMEXProperties

      public void setIMEXProperties(double[] soundSpeeds, double[] densities, double dx, double outletPressure, boolean outletFixed)
      Set cell properties required for the IMEX pressure correction step. Must be called before stepping with IMEX_PRESSURE_CORRECTION.
      Parameters:
      soundSpeeds - sound speed per cell (m/s)
      densities - mixture density per cell (kg/m3)
      dx - cell size (m)
      outletPressure - outlet boundary pressure (Pa)
      outletFixed - true if outlet pressure is a Dirichlet BC

    • stepIMEXPressureCorrection

      public double[][] stepIMEXPressureCorrection(double[][] U, TimeIntegrator.RHSFunction rhs, double dt)
      IMEX (Implicit-Explicit) pressure correction step.

      The algorithm follows a two-stage splitting:

      1. Predictor (explicit): Advance mass and momentum using the explicit RHS (advection + source terms) to obtain intermediate values U*.
      2. Pressure correction (implicit): Solve a tridiagonal Helmholtz equation for the pressure correction dp that enforces mass conservation implicitly. The pressure wave equation dp - (c*dt/dx)^2 * d^2(dp)/dx^2 = RHS removes the acoustic CFL constraint.
      3. Corrector: Update momenta using the pressure correction gradient.

      This allows the convective CFL (based on material velocity |v|) to govern the time step instead of the acoustic CFL (based on |v|+c). For typical gas-liquid flows where c=300 m/s and v=5 m/s, this gives a factor of ~60 speedup.

      Parameters:
      U - Current state [nCells][nVars]
      rhs - Right-hand side function (explicit part)
      dt - Time step
      Returns:
      Updated state at t + dt

    • solveTridiagonal

      private double[] solveTridiagonal(double[] a, double[] bDiag, double[] c, double[] d)
      Solve a tridiagonal system Ax = d using the Thomas algorithm (O(n)).

      The system has the form: a_i * x_{i-1} + b_i * x_i + c_i * x_{i+1} = d_i for i = 0..n-1, where a_0 = 0, c_{n-1} = 0.

      Parameters:
      a - sub-diagonal (a[0] is not used)
      bDiag - main diagonal
      c - super-diagonal (c[n-1] is not used)
      d - right-hand side
      Returns:
      solution vector x

    • calcIMEXTimeStep

      public double calcIMEXTimeStep(double[] gasVelocities, double[] liquidVelocities, double dx)
      Calculate stable time step for IMEX method. The IMEX scheme removes the acoustic CFL constraint, so the time step is limited only by the convective CFL based on material velocities (not sound speed).

      dt_IMEX = CFL * dx / max(|v_G|, |v_L|)

      For typical gas-liquid flows where c=300 m/s and v=5 m/s, this gives a speedup factor of ~60 compared to the standard acoustic CFL.

      Parameters:
      gasVelocities - gas velocities at each cell (m/s)
      liquidVelocities - liquid velocities at each cell (m/s)
      dx - cell size (m)
      Returns:
      stable time step (s), typically 10-100x larger than acoustic CFL

    • reset

      public void reset()
      Reset integrator state.