Libraries and Models

The main application is rocketMotor which depends on several libraries including the base libraries of OpenFOAM-v2112. The solver is developed based on the native OpenFOAM v2112 solver "reactingMultiphaseEulerFoam". This framework has an additional library that includes the relavant closure models for propellant surface regression and interphase interactions (drag and heat transfer models).

%%{ init: { 'theme': 'neutral', 'fontFamily': 'Exo 2', 'flowchart': { 'curve': 'stepAfter' }, 'themeVariables': { 'primaryColor': '#ffffff', 'edgeLabelBackground': '#ffffff', 'tertiaryColor': '#ffffff', 'lineColor': '#000000' } } }%% flowchart LR %% Left column: base libraries and multiphase group subgraph OpenFOAM_Base_Group [OpenFOAM v2112 Base Libraries] List["libfiniteVolume.so
libOpenFOAM.so
libfluidThermophysicalModels.so
libreactionThermophysicalModels.so
libcompressibleTransportModels.so
libspecie.so
libcompressibleTurbulenceModels.so
libturbulenceModels.so
libthermophysicalProperties.so
libcombustionModels.so
libmeshTools.so
libsurfMesh.so
libdynamicMesh.so"] style List text-align:left end subgraph OpenFOAM_Multiphase_Group [OpenFOAM v2112 Multiphase Library] B["libreactingMultiphaseSystem.so"] end %% Right column: RocPerf framework containing source library and main app subgraph RocPerf_Group [RocPerf Framework] direction TB C["Updates to the base library
(libpropellantRegressionPhaseSystem.so)
Particle Drag Models
Heat Transfer Models
Interface Tracking Model"] D["
rocketMotor
(Main Application)
"] end %% Connections List -->|Mesh, Thermophysical &
Transport properties| C B -->|Phase Models, Phase System &
Phase Pair| C C -->|Multiphase Systems &
Interphase Models| D %% Styling - keep the original color scheme style OpenFOAM_Base_Group fill:#fff8e1,stroke:#ffa000,stroke-width:2px; style OpenFOAM_Multiphase_Group fill:#fff8e1,stroke:#ffa000,stroke-width:2px; style List fill:#ffe0b2,stroke:#ffa000,stroke-width:1px; style B fill:#ffe0b2,stroke:#ffa000,stroke-width:1px; style C fill:#eceff1,stroke:#78909c,stroke-width:2px; style D fill:#c5e1a5,stroke:#7cb342,stroke-width:4px; %% Link Styling linkStyle 0,1,2 stroke:#000000,stroke-width:2px,fill:none,color:#000000
Particle Drag Model ▼
particleDragModel.C
Foam::tmp<Foam::volScalarField> Foam::particleDragModel::Ki() const
{
    const tmp<volScalarField> tmu(pair_.continuous().nu()*pair_.continuous().rho());
    volScalarField muc(tmu());

    if (isInviscid())
    {
        const tmp<volScalarField> tmc(mu());
        muc = tmc();
    }
    return
        0.75*factor_
        *CdRe()
        *muc
        /sqr(pair_.dispersed().d());
}
  • particleDragModel is the base class which is derived by all the drag models for the drag coefficient implementation.
  • The drag parameter $\delta$ is implemented within the function K() present in particleDragModel.C.
    • $$ \delta = \frac{3}{4} \frac{\alpha_p\mu_c}{d_p^2} C_D Re_p $$
  • The dimensionless drag coefficient term, $C_D Re_p$, is computed according to the specific drag model selected at runtime (from phaseProperties file).

Available Drag Models

Drag Models Source File Model/References
Schiller-NaumannSchillerNaumann.C$C_D = \frac{24}{Re_p}\left(1.0 + 0.15(Re_p)^{0.687}\right)$
Clift-GauvinCliftGauvin.C$C_D = \frac{24}{Re_p}\left(1.0 + 0.15(Re_p)^{0.687} + \frac{0.42}{1 + 42500/Re_p^{1.16}}\right)$
Clift-Gauvin (for inviscid simulations)CliftGauvinInviscid.CSame as Clift-Gauvin drag model with viscosity of gas phase computed based on sutherland correlation.
Gidaspow Ergun Wen-YuGidaspowErgunWenYu.CGidaspow (1994)
HendersonHenderson.CHenderson (1976)
Eric LothLoth.CLoth et. al., (2008)
Heat Transfer Model ▼
RanzMarshall.C
Foam::tmp<Foam::volScalarField> 
Foam::sharpInterfaceHeatTransferModels::RanzMarshall::K(const scalar residualAlpha) const
{
    volScalarField Nu(scalar(2) + 0.6*sqrt(pair_.Re())*cbrt(pair_.Pr()));

    return
        6.0
        *pair_.dispersed()*pos(pair_.dispersed() - cutoff)
        *pair_.continuous().kappa()
        *Nu
        /sqr(pair_.dispersed().d());
}
The heat transfer coefficients are calculated within their respective sharpInterfaceHeatTransferModel source files. The implementation of the heat transfer parameter based on the Ranz-Marshall correlation is shown in the code snippet. $$ \mathcal{H} = 6 \text{Nu} \frac{\alpha_p\lambda_c}{d_p^2} $$ According to RanzMarshall correlation, the Nusselt number is computed as, $$ \text{Nu} = 2.0 + 0.6(Re_p)^{1/2}(Pr)^{1/3} $$

Available Heat Transfer Models

Model Name Source File Nusselt Relation ($Nu$)/References
Ranz-Marshall RanzMarshall.C $\text{Nu} = 2 + 0.6 Re_p^{1/2} Pr^{1/3}$
Kavanau-Ranz-Marshall KavanauRanzMarshall.C $\text{Nu} = \frac{\text{Nu}_0}{1 + 3.42\frac{M_p}{Re_pPr}\text{Nu}_0}$ where, $\text{Nu}_0$ is the Nusselt number computed from Ranz-Marshall correlation
DrakeInviscid (for inviscid simulations) DrakeInviscid.C $\text{Nu} = 2 + 0.459 Re_p^{0.55} Pr^{0.33}$ and viscosity of the gas phase is computed from sutherland correlation.
Interface Tracking Model ▼

subCellularInterfaceMotion

A simple surface regression algorithm based on the empirical burning rate law, employing sub-grid level interface tracking. The algorithm identifies the interface location from the propellant phase volume fraction field and calculates rate of regression based on the local pressure (burning rate law) at each time step. The routine is invoked by the regress() function in the subCellularInterfaceMotion.C. This process updates the propellant phase volume fraction and the mass burning rate, which are subsequently integrated as source terms into the governing equations. For further implementation details, users are encouraged to review the source code directly.

entrainedInterfaceMotion

This model facilitates the tracking of dual interfaces, making it particularly effective for simulations where particles are retained within the propellant bed. Extending the logic of the subCellularInterfaceMotion algorithm, this approach independently tracks both the flame surface and the bed surface. The flame surface regresses according to the burning rate law, while the regression of the bed surface is governed by the fraction of particles retained from the flame. Consequently, this model requires the mass retention fraction as an additional input parameter. The algorithm is implemented within the regress() function in the entrainedInterfaceMotion.C. For further information on particle retainment model, the users are referred to Venukumar et. al., (2025)

🚀 CFD solver - rocketMotor

rocketMotor/rocketMotor.C (PIMPLE Loop) 
        // ... previous code ... (For one time step)
        
        while (pimple.loop())
        {
            fluid.solve();
            fluid.correct();
        
            // Find propellant cells to restrict velocity
            #include "propellantCells.H"
        
            // Find particle free cells to restrict Temp
            #include "particleFreeCells.H"
        
            // Solve Conservation Equations
            #include "pU/UEqns.H"
            #include "EEqns.H"
            #include "pU/pEqn.H"
        
            fluid.correctKinematics();
            fluid.correctTurbulence();
        }
        
        // ... remaining code ...

Outer PIMPLE Loop

We enter the main iterative loop for the current time step.

Solve volume fraction equations for propellant, gas and particle phases

At the start of every iteration, volume fraction equations are solved by calling the solve() method on the multiPhaseSystem object (fluid). This triggers a specific sequence:

  1. Propellant Regression:
    First, it invokes the solve() function in PropellantCombustionPhaseSystem.C (or RegressionPhaseSystem). Here, the mass burning rate is computed, and the propellant surface is regressed according to the selected interfaceTrackingModel.
  2. Gas & Particle Phases:
    Next, it solves the volume fractions for the remaining phases. This utilizes the MULES algorithm based on the native OpenFOAM-v2112 implementation found in MultiPhaseSystem.C.

Compute Mass Fluxes

The correct() method in phaseSystem.C computes the mass fluxes for all phases. This update is based on two key inputs:

  • The new volume fraction fields (calculated in the previous step).
  • The velocity field from the previous iteration.

Stability & Robustness

Unique to the rocketMotor solver, we apply explicit constraints to avoid non-physical values in the shared fluid domain:

  • Propellant Cells:
    In cells containing only solid propellant, gas and particle velocities are fixed to zero. This enforces the boundary condition at the regressing surface and prevents non-physical flows inside the solid.
  • Particle-Free Regions:
    Where the particle fraction is negligible ($\alpha_p < 10^{-10}$), the particle temperature is forced to match the gas temperature ($T_p = T_g$). This prevents numerical divergence when heat transfer coefficients become artificially large in cells with no particles.

Constructs Discretized Momentum Equation for gas and particle phase

The base momentum equation (excluding the pressure gradient) is constructed by calling the UEqn() method in the native OpenFOAM-v2112 movingPhaseModel.C.

However, inside UEqn.H, we extend this by calling momentumTransfer() to inject specific source terms sequentially:

  • Regression Source: Added via PropellantRegressionPhaseSystem.C to account for momentum changes due to surface regression.
  • Explicit Drag: Added via MomentumEnergyTransferPhaseSystem.C to include explicit interphase drag forces.

Note: The implicit drag force terms are not added here; they are reserved for the pressure equation (pEqn.H) in the next section.

Energy Correction Loop

The energy equations are constructed via the EEqn() method in the respective phase models:

  • Gas: ThermalEnergyPhaseModel.C
  • Particles: AnisothermalParticlePhaseModel.C

Interphase interactions are added separately via InterphaseHeatTransferPhaseSystem.C. Finally, enthalpy is converted to Temperature using a Newton-Raphson iterator, repeating nEnergyCorrector times.

Pressure-Velocity Coupling

  • Predict Velocity:
    Using the discretized momentum equation matrix, the explicit drag force terms are included by calling phiKdPhis() and KdUByAs() in MomentumEnergyTransferPhaseSystem.C. These predict the velocities and fluxes for both gas and particle phase, which are summed to obtain the total predicted flux.
  • Pressure Poisson Equation:
    The equation for the common pressure $P$ is constructed based on the total predicted flux and the compressibility terms of the gas and particle phases.
  • Correct Velocity:
    Finally, corrections are applied to the predicted velocity and flux fields. In addition, the density is updated here. This process repeats for nCorrectors loops (as set in fvSolution).

Turbulence Correction

Once the mean flow fields ($\alpha_k, \boldsymbol{u}_k, P, T_k$) are updated, the turbulence closure models are updated.

  • Gas Phase Turbulence:
    The turbulence models are corrected to update the eddy viscosity ($\nu_t$) and Reynolds stresses.
  • Kinetic Theory of Granular Flow:
    Simultaneously, the granular temperature ($\Theta$) is solved for the particle phase. This updates the particle pressure, shear stress, and collisional viscosities based on selected KTGF models.

Boundary Conditions

As this library is an extension of the OpenFOAM-v2112 library, it fully supports all standard boundary conditions included in the base distribution. You can use standard conditions (e.g., fixedValue, zeroGradient) interchangeably with the specialized boundary conditions listed below.

Boundary Condition Description Compatible Fields Usage
massLoadingFraction Computes the particle phase volume fraction $\alpha_p$ at the inlet derived from a specified mass loading ratio ($X_p$). $$ \alpha_p = \left[1 + \left(\frac{1 - X_p}{X_p}\right)\left(\frac{\rho_g\boldsymbol{u}_g\cdot\boldsymbol{S}_f}{\rho_p\boldsymbol{u}_p\cdot\boldsymbol{S}_f}\right)\right]^{-1} $$ where $\boldsymbol{S}_f$ is the patch area. $\alpha_p$
(volume fraction) 
inlet
{
    type            massLoadingFraction;
    name            particles; // name of the current phase
    gasPhase        gas; // name of the gas phase
    particlePhase   particles; // name of the particle phase
    massFraction    0.81; // particle phase mass fraction at the inlet
    value           uniform 0;
}
gasFlowRateInletVelocity Computes the gas phase inlet velocity based on a specified mass flow rate ($\dot{m}_g$). The flow rate can be defined as a constant value or a time-dependent function. $$ \boldsymbol{u}_g = \frac{\dot{m}_g}{\alpha_g\rho_g|\boldsymbol{S}_f|}\cdot\hat{\boldsymbol{n}}_f $$ where $\hat{\boldsymbol{n}}_f$ is the patch normal vector. $\boldsymbol{u}_g$
(vector) 
inlet
{
    type           gasFlowRateInletVelocity;
    massFlowRate   constant 0.0001; // flow rate function
    rho            thermo:rho.gas; // density of the gas phase
    value          uniform (0 0 0);
}
gasBurningRateInletVelocity Calculates the gas phase inlet velocity based on the propellant mass burning rate which is computed based on burning rate law $r_b[\text{cm/s}] = a(P[\text{MPa}])^n$, where $P$ is the local pressure in the units of $\text{MPa}$. $$ \boldsymbol{u}_g = \frac{X_g\rho_{prop}r_b}{\alpha_g\rho_g}\cdot\hat{\boldsymbol{n}}_f $$ where $\hat{\boldsymbol{n}}_f$ is the patch normal vector. $\boldsymbol{u}_g$
(Vector) 
inlet
{
    type            gasBurningRateInletVelocity;
    rhoPropellant   1500; // density of the propellant
    a               0.607; // pre-exponential factor
    n               0.798; // pressure exponent
    X               0.5; // gas phase mass fraction
    value           uniform (0 0 0);
}
sameAs A utility condition that copies boundary values from another specified field to the current patch. Useful for coupling fields directly between phases. Any Field 
inlet
{
    type    sameAs;
    field   T.gas; // field to the copied from...
    value   uniform 0;
}
transonic The pressure outlet boundary condition that switches from back pressure condition to extrapolation condition as the flow transitions from subsonic to supersonic flow at the outlet or vice versa from extrapolation condition to back pressure condition. $P$
(scalar) 
outlet
{
    type    transonic;
    backPressure   101325; // back pressure / atmospheric pressure
    gas   gas; // name of the gas phase
    gamma   1.4; // specific heat ratio
    R   287; // gas constant
    supersonic   false; // 'true' if the flow is already supersonic and 
                        //  expect it to change to subsonic
    value   uniform 101325;
}

In addition to the above boundary conditions, there are additional boundary conditions are currently in developmental stage and will be included in the future releases.

Solver Utilities

The following utilities are included to streamline the workflow, from initialization to post-processing.

postProcessRocket

Post-processing

Automatically identifies the domain outlet patch and computes the mass flow rate and thrust evolution for all available time steps.

Usage:
$ postProcessRocket

setRocketInitField

Pre-processing

Sets the initial propellant volume fraction field based on geometric regions.
Note: Requires a setRocketInitFieldDict dictionary to be present in the case directory.

Usage:
$ setRocketInitField

mapVolFields

Grid Analysis

Maps solution fields from a fine grid (source) to the current coarse grid (target) using volume-weighted interpolation. Primarily used for grid independence studies and a posteriori error estimation.
Requirement: The fine mesh must have been generated by refining the coarse mesh using refineMesh.
Note: You must execute this command from within the target case directory (the coarse mesh case).

Usage:
$ mapVolFields <sourceDir> -sourceTime <time>