Four-field solver

The Solvers.FourField module in OpenSTREAM provides tools and classes for simulating annular two-phase flows using a four-field approach. In this method, the liquid phase is represented by three distinct fields (base film, disturbance waves, and droplets) each governed by its own set of conservation equations for mass, momentum, and energy. The vapor phase is solved using the Solvers.Mixture module. This framework enables detailed modeling of liquid film structures, interfacial interactions, slip velocities, and hydrodynamic non-equilibrium effects, making it particularly suitable for annular boiling two-phase flow regimes where field separation and wave dynamics play a significant role.

This module includes:

Solvers.FourField.FourFieldSolver.solve(ffSolver)

SOLVE Executes the four-field solver for steady-state and transient simulations.

Runs the full solution process for the Solvers.FourField.FourFieldSolver object. It manages initialization, time stepping, axial sweeps, inner iterations, convergence checks, and logging.

Workflow:

  • Enables logging and opens persistent log file

  • Validates solver state before execution

  • Runs steady-state initialization (solveINIT = true)

  • If initial step converges, proceeds with transient simulation

  • Handles exceptions and ensures proper log closure

Notes:

  • Uses internal solver function to handle both steady-state and transient modes

  • Applies relaxation factors for phase flows, velocities, and enthalpies

  • Supports thermal non-equilibrium modeling

  • Logs progress and outputs to session directory

class Solvers.FourField.FourFieldSolver.FourFieldSolver(inputSet, mixSolver)

Bases: Solvers.ThreeField.ThreeFieldSolver

FOURFIELDSOLVER Solver for initializing, solving, and visualizing four-field flow.

The FourFieldSolver class extends Solvers.ThreeField.ThreeFieldSolver to handle thermal-hydraulic simulations using a four-field approach.

Responsibilities:

Key Components:

  • Liquid film (base and wave) and drop field construction

  • Field mass/momentum/energy transport modeling

  • Pressure drop gradient obtained from mixture solver

  • Vapor solution obtained from mixture solver

  • Visualization of results across time and axial domains

Constructor Summary
FourFieldSolver(inputSet, mixSolver)

FOURFIELDSOLVER Constructor for the FourFieldSolver class.

Creates a new instance of the FourFieldSolver class using the provided Inputs.InputSet. This constructor initializes the solver by calling the Solvers.ThreeField.ThreeFieldSolver constructor and then sets up all necessary solver parameters and internal data structures via Solvers.FourField.FourFieldSolver.initializeSolver

Input:

  • inputSet — An Inputs.InputSet object containing model configuration, geometry, boundary conditions, and solver options.

Notes:

  • If mixSolver is unsolved, it will be solved automatically.

  • Initialization includes film, wave, base, and droplet objects.

Method Summary
EQUIL(ffSolver, flm, drp, mix, zIdx)

EQUIL Find entrained ratio at film/drop equilibrium state [-]

Iteratively computes the entrained liquid ratio at the onset of annular flow by balancing droplet deposition and film entrainment rates.

Inputs:

Algorithm:

  • Starts with an initial guess for droplet mass flow (50% of liquid mass)

  • Updates droplet and film flow rates iteratively

  • Computes difference between deposition and entrainment rates

  • Stops when error < 1e-4 kg/s/m or after 100 iterations

Notes:

  • Uses interpolation for improved convergence

  • Falls back to ad-hoc update if interpolation fails

initializeSolver(ffSolver)

INITIALIZESOLVER Initializes solver parameters and constructs film, wave, base, and drop objects.

Sets up the solver’s internal state using the provided Inputs.InputSet and mixture solver solutions. Includes initialization of data structures for transient and steady-state simulations.

Operations performed:

  • Initializes film, wave, and base arrays for each time step

  • Initializes fluid property objects

  • Sets up mixture object

  • Assigns initial conditions for field mass flow rates, velocities, enthalpies

Notes:

  • Assumes uniform axial discretization.

  • Wall heat flux computed from interpolated boundary conditions.

plott(ffSolver, zIdx, opt)

PLOTT Plots temporal distributions of four-field parameters at a given axial location.

Generates time-series plots of selected four-field parameters at a specified axial index. The method supports multiple display modes, wall selections, and plot arrangements.

Inputs:

  • ffSolver — Solvers.FourField.FourFieldSolver object containing simulation data

  • zIdx — Axial index (scalar, positive integer)

  • opt.display — Parameters to display (e.g., ‘HFLUX’, ‘W’, ‘U’, etc.)

  • opt.solveMode — Solve mode: ‘REAL’ or ‘NULL’

  • opt.wall — Wall index(es) to plot

  • opt.tIdx — Time indices to include in the plot

  • opt.reverseTime — Logical flag to reverse time axis

  • opt.unitTemp — Temperature unit: ‘K’ or ‘C’

  • opt.arrangement — Plot arrangement: ‘flow’, ‘vertical’, or ‘horizontal’

  • opt.resize — Resize factor for plot scaling

Notes:

  • If opt.display is set to ‘ALL’, all supported parameters are plotted.

  • The function validates that at least two time indices are provided.

  • Temperature unit conversion is applied if ‘C’ is selected.

  • Each wall is plotted in a separate tile with appropriate legends and axis scaling.

plotz(ffSolver, tIdx, opts)

PLOTZ Plots spatial (axial) distributions of four-field parameters.

Generates axial plots of selected four-field parameters at specified time indices. The function supports multiple display modes, wall selections, and optional animation over time.

Inputs:

  • ffSolver — Solvers.FourField.FourFieldSolver object containing simulation data

  • tIdx — Time index or indices (vector of positive integers)

  • opts.display — Parameters to display (e.g., ‘HFLUX’, ‘W’, ‘U’, etc.)

  • opts.solveMode — Solve mode: ‘REAL’ or ‘NULL’

  • opts.wall — Wall index(es) to plot

  • opts.zIdx — Axial indices to include in the plot

  • opts.obstructions — Logical flag to display obstruction locations

  • opts.annular — Logical flag to plot annular two-phase flow parameters only from the onset of annular flow

  • opts.unitTemp — Temperature unit: ‘K’ or ‘C’

  • opts.arrangement — Plot arrangement: ‘flow’, ‘vertical’, or ‘horizontal’

  • opts.resize — Resize factor for plot scaling

Notes:

  • If multiple time indices are provided, the function creates an animated plot.

  • At least two axial indices must be specified to enable spatial plotting.

  • Temperature unit conversion is applied if ‘C’ is selected.

  • Obstruction locations are plotted if opts.obstructions is true.

  • Each wall is plotted in a separate tile with appropriate legends and axis scaling.

plotzt(ffSolver, opt)

PLOTZT Plots 2D time/elevation distributions of four-field parameters.

Generates surface plots of selected four-field related parameters over time and axial (elevation) positions. It supports plotting for base and wave fields, and can handle both real and null (initial) transient data.

Inputs:

  • ffSolver — Solvers.FourField.FourFieldSolver object containing simulation data

  • opt.display — Parameter(s) to display (e.g., ‘HFLUX’, ‘U’, etc.)

  • opt.label — Corresponding labels for display parameters

  • opt.unit — Units for each parameter

  • opt.field — Field to plot: ‘film’ or ‘drop’

  • opt.solveMode — Solve mode: ‘REAL’ or ‘NULL’

  • opt.wall — Wall index(es) to plot

  • opt.zIdx — Axial indices (must include at least 2)

  • opt.tIdx — Time indices (must include at least 2)

  • opt.reverseTime — Logical flag to reverse time axis

  • opt.shading — Surface shading style: ‘faceted’, ‘flat’, or ‘interp’

  • opt.view — View angle for 3D plot [azimuth elevation]

Notes:

  • If ‘ALL’ is passed to opt.display, all supported parameters are plotted.

  • The function validates that at least two axial and time indices are provided to enable meaningful 2D plotting.

  • Each wall is plotted in a separate figure window with appropriate titles and subplot grouping.

class Solvers.FourField.Film(inputSet, fluid)

Bases: Solvers.AbstractFilm

FILM Class for modeling liquid film in four-field solver

Represents the physical and numerical properties of the liquid film in a four-field thermal-hydraulic solver. The film consists of two sub-fields: base film and wave film, which interact with vapor and droplet phases.

This class manages initialization, mass flow distribution, velocity and enthalpy calculations, and provides methods for splitting film flow at the onset of annular flow (OAF) according to model options.

Key responsibilities:

  • Store solver state (time, axial positions, heat flux, evaporation flux)

  • Aggregate flow properties from base and wave components

  • Initialize base and wave fields and apply OAF split logic

  • Support equilibrium calculations for film thickness and flow ratios

Constructor Summary
Film(inputSet, fluid)

FILM Creates a Film object

Initializes the film object with input set and fluid properties required by the four-field solver.

Inputs:

Property Summary
DT = 0

Time step size [s] from Inputs.options.TSTEP

H

(:,:) double {mustBeNumeric} = 1E6 % Enthalpy [J/kg]

HFLUX = 1.

Film heat flux [W/m^2]

ITR

Iteration tracking

MEVAP = -1.

Evaporation mass flux [kg/s/m^2]

NTIME = 0

Number of time steps [-]

NZ = 0

Number of axial steps [-] from Inputs.Model.NNODES

TIDX = 1

Time step index [-]

TIME = 0

Time series [s]

U

(:,:) double {mustBeNumeric} = 1. % Velocity [m/s]

W

(:,:) double {mustBeNumeric} = 1. % Mass flow rate [kg/s]

Z = 1.

Elevation [m]

base = NaN

Solvers.FourField.Base object

wave = NaN

Solvers.FourField.Wave object

Method Summary
EQUIL(film, Wf, zIdx)

EQUIL Equilibrium base/wave split ratio [-]

Calculates the equilibrium split ratio of film mass flow between base and wave using base equilibrium thickness, base velocity, perimeter, and liquid density.

Inputs:

  • film — Solvers.FourField.Film object

  • Wf — Film mass flow rate [kg/s]

  • zIdx — Axial indices to evaluate (optional)

copyFlowProperties(srcObj, targetObj, opts)

COPYFLOWPROPERTIES Copy flow properties between Film objects

Copies flow-related properties (base and wave fields) from a source Film object to one or more target Film objects. Ensures that spatial meshes match before copying.

Inputs:

distributeOAFW(film, Wf, zIdx)

DISTRIBUTEOAFW Distribute film mass flow at OAF [kg/s]

Distributes the film mass flow rate at the onset of annular flow (OAF) between base and wave fields according to the Inputs.Model.OAFFILMSPLIT model option.

Inputs:

  • film — Solvers.FourField.Film object

  • Wf — Film mass flow rate to distribute [kg/s] (per wall or per z-index)

  • zIdx — Axial indices to apply distribution (optional)

initializeBaseAndWave(film, WIN, ITR)

INITIALIZEBASEANDWAVE Initialize base and wave fields

Initializes base and wave arrays given inlet film mass flow rate and iteration settings. Performs OAF film split, computes velocities, enthalpies, and sets initial wave frequency using equilibrium correlations.

Inputs:

  • film — Solvers.FourField.Film object

  • WIN — Inlet film mass flow rate [kg/s]

  • ITR — Iteration control struct

class Solvers.FourField.Wave(film)

Bases: Solvers.AbstractFilm

WAVE Class for modeling liquid waves in four-field solver

This class represents the physical and numerical properties of the wave field in a four-field thermal-hydraulic solver. It handles flow variables, heat transfer, and phase interaction mechanisms between the waves and other fields (base film, drop, vapor).

The class supports initialization from an existing film object and provides methods for computing derived quantities such as mass flow fractions, interfacial fractions, and deposition/evaporation fluxes.

Constructor Summary
Wave(film)

WAVE Constructor for Wave class

Initializes the wave object using properties from an existing film object. Copies solver configuration and geometry data, and allocates flow property arrays for multiple walls.

Inputs:

Property Summary
DT = 0

Time step size [s] from Inputs.options.TSTEP

FREQUENCY = 1

Wave frequency [Hz]

H = 1E6

Enthalpy [J/kg]

HFLUX = 1.

Film heat flux [W/m^2]

ITR

Iteration tracking

NTIME = 0

Number of time steps [-]

NZ = 0

Number of axial steps [-] from Inputs.Model.NNODES

TIDX = 1

Time step index [-]

TIME = 0

Time series [s]

U = 1.

Velocity [m/s]

W = 1.

Mass flow rate [kg/s]

Z = 1.

Elevation [m]

film

Parent Solvers.AbstractFilm object

Method Summary
AMPLITUDE(wave, zIdx)

AMPLITUDE Wave amplitude [m]

Computes wave amplitude based on mass flow, shape factor, fluid density, and frequency. Limited by half the hydraulic diameter.

Inputs:

BETA(wave, zIdx)

BETA Wave film interfacial fraction [-]

Calculates the fraction of the interface occupied by waves.

Inputs:

BETAP(wave, zIdx)

BETAP Wave film heat flux fraction [-]

Determines the fraction of heat flux directed to the waves.

Inputs:

DELTAFREQ(wave, zIdx)

DELTAFREQ Wave frequency source/sink [Hz/m]

Calculates the exchange term for wave frequency relaxation toward equilibrium frequency.

Inputs:

DRAGCOEF(wave, zIdx)

DRAGCOEF Wave drag coefficient [-]

Calculates the drag coefficient for vapor-wave interaction based on Reynolds number and Inputs.Model.WAVEDRAGCOEF model coefficients.

Inputs:

EPSILON(wave, zIdx)

EPSILON Wave mass flow fraction [-]

Computes the fraction of total film mass flow carried by waves.

Inputs:

EQFREQUENCY(wave, zIdx)

EQFREQUENCY Equilibrium wave frequency [Hz]

Computes the equilibrium wave frequency based on Strouhal number and vapor velocity.

Inputs:

EQPERIOD(wave, zIdx)

EQPERIOD Equilibrium wave period [s]

Calculates the inverse of equilibrium wave frequency.

Inputs:

EQSTROUHAL(wave, zIdx)

EQSTROUHAL Equilibrium Strouhal number [-]

Calculates the equilibrium Strouhal number using selected Inputs.Model.EQSTROUHAL model.

Inputs:

ETA(wave, zIdx)

ETA Wave film deposition fraction

Computes the fraction of deposition attributed to waves.

Inputs:

FBASE(wave, drop, zIdx)

FBASE Base interfacial force [N/m^2]

Computes the interfacial force between wave and base film.

Inputs:

FBASEMASS(wave, drop, zIdx)

FBASEMASS Wave mass exchange force [N/m^2]

Calculates the force due to mass exchange between wave and base film.

Inputs:

FDEP(wave, drop, zIdx)

FDEP Droplet interaction force [N/m^2]

Calculates the force due to deposition of droplets onto the wave field.

Inputs:

FDRAG(wave, zIdx)

FDRAG Vapor drag force on wave [N/m^2]

Calculates the shear stress due to vapor drag acting on the wave field.

Inputs:

FSHEAR(wave, zIdx)

FSHEAR Vapor shear stress [N/m^2]

Computes the shear stress exerted by vapor on the wave field using friction factor.

Inputs:

FTOT(wave, drop, zIdx)

FTOT Total interfacial force [N/m^2]

Computes the total interfacial force acting on the wave field, including vapor shear, buoyancy, gravity, droplet interaction, and base film forces.

Inputs:

FVAPOR(wave, zIdx)

FVAPOR Wave vapor shear stress [N/m^2]

Computes the total shear stress exerted by vapor on the wave field.

Inputs:

FWALL(wave, zIdx)

FWALL Wave wall shear stress [N/m^2]

Placeholder method. Wave does not implement wall shear stress.

Inputs:

MBASE(wave, drop, zIdx)

MBASE Mass flux interaction with base film [kg/m^2/s]

Computes the net mass exchange between wave and base film, including turbulent mixing.

Inputs:

MDEP(wave, drop, zIdx)

MDEP Wave deposition mass flux [kg/m^2/s]

Computes the deposition mass flux from waves to droplets.

Inputs:

MENT(wave, zIdx)

MENT Wave entrainment mass flux [kg/m^2/s]

Computes the mass flux associated with entrainment from the wave field.

Inputs:

MEVAP(wave, zIdx)

MEVAP Wave evaporation mass flux [kg/m^2/s]

Calculates the evaporation mass flux from the wave field based on heat flux fraction.

Inputs:

MTOT(wave, drop, zIdx)

MTOT Total mass flux [kg/m^2/s]

Calculates the total mass flux for the wave field, including evaporation, entrainment, deposition, and interactions with base film.

Inputs:

MTURB(wave, zIdx)

MTURB Turbulent mass exchange [kg/m^2/s]

Returns the turbulent mixing mass flux between wave and base film.

Inputs:

N(wave, zIdx)

N Wave number density [1/m]

Computes the number of waves per unit length (inverse of wave spacing).

Inputs:

PERIOD(wave, zIdx)

PERIOD Wave time period [s]

Calculates the inverse of wave frequency.

Inputs:

REV(wave, zIdx)

REV Vapor Reynolds number (relative to waves) [-]

Computes the Reynolds number for vapor flow interacting with waves.

Inputs:

SHAPEFACTOR(wave, zIdx)

SHAPEFACTOR Wave shape factor [-]

Computes the dimensionless shape factor for waves based on Reynolds number and Inputs.Model.SHAPEFACTORCOEF model.

Inputs:

SPACING(wave, zIdx)

SPACING Wave spacing [m]

Computes the axial spacing between waves based on velocity and frequency.

Inputs:

THICK(wave, zIdx)

THICK Wave thickness [m]

Computes the wave thickness as the product of amplitude and interfacial fraction. This parameter represents the thickness that the wave field would have if fully overlayed on top of the base film.

Inputs:

WIDTH(wave, zIdx)

WIDTH Wave width [m]

Calculates wave width as a function of amplitude and shape factor, limited by wave spacing.

Inputs:

WL(wave, zIdx)

WL Wave mass flow rate per unit perimeter [kg/s/m]

Calculates the mass flow carried by the wave field at specified axial locations.

Inputs:

class Solvers.FourField.Base(film)

Bases: Solvers.AbstractFilm

BASE Class for modeling liquid base film in four-field solver

This class represents the physical and numerical properties of the base film field in a four-field thermal-hydraulic solver. It handles flow variables, heat transfer, and phase interaction mechanisms between the base film and other fields (wave, drop, vapor).

The class supports initialization from an existing film object and provides methods for computing derived quantities such as mass flow fractions, interfacial fractions, and deposition/evaporation fluxes.

Constructor Summary
Base(film)

BASE Constructor for Base class

Initializes the base film object using properties from an existing film object. Copies solver configuration and geometry data, and allocates flow property arrays for multiple walls.

Inputs:

Property Summary
DT = 0

Time step size [s] from Inputs.options.TSTEP

H = 1E6

Enthalpy [J/kg]

HFLUX = 1.

Film heat flux [W/m^2]

ITR

Iteration tracking

NTIME = 0

Number of time steps [-]

NZ = 0

Number of axial steps [-] from Inputs.Model.NNODES

TIDX = 1

Time step index [-]

TIME = 0

Time series [s]

U = 1.

Velocity [m/s]

W = 1.

Mass flow rate [kg/s]

Z = 1.

Elevation [m]

Method Summary
BETA(base, zIdx)

BETA Base film interfacial fraction [-]

Computes the fraction of interfacial area attributed to the base film.

Inputs:

BETAP(base, zIdx)

BETAP Base film heat flux fraction [-]

Determines the fraction of heat flux directed to the base film. Applies suppression in pre-annular region to enforce target base film thickness at onset of annular flow (OAF).

Inputs:

Notes:

  • Applies suppression in pre-annular region to enforce target base film thickness at onset of annular flow (OAF). The entire film evaporation hence occurs in the wave field.

EPSILON(base, zIdx)

EPSILON Base film mass flow fraction [-]

Calculates the fraction of total film mass flow carried by the base film at specified axial locations.

Inputs:

EQTHICK(base, zIdx)

EQTHICK Base equilibrium thickness [m]

Computes the equilibrium thickness of the base film based on Inputs.Model.BASEEQTHICK model.

Inputs:

ETA(base, zIdx)

ETA Base film deposition fraction [-]

Computes the fraction of drop deposition allocated to the base film based on mass flow and interfacial fractions.

Inputs:

FBASEVAPOR(base, zIdx)

FBASEVAPOR Vapor shear stress on base film [N/m^2]

Computes the shear force exerted by the vapor phase on the base film, scaled by interfacial fraction.

Inputs:

FDEP(base, drop, zIdx)

FDEP Droplet deposition force [N/m^2]

Calculates the momentum transfer to the base film from droplet deposition, based on velocity difference and deposition flux.

Inputs:

FDRY(base, drop, zIdx)

FDRY Dry fraction [-]

Computes the fraction of the exposure period during which the base film is expected to be dry.

Inputs:

FTOT(base, drop, zIdx)

FTOT Total interfacial force on base film [N/m^2]

Aggregates all force contributions acting on the base film, including wall shear, wave interaction, vapor shear, buoyancy, gravity, and droplet deposition based on Inputs.Model.MOMENTBASE model.

Inputs:

FWAVE(base, drop, zIdx)

FWAVE Wave interfacial force [N/m^2]

Computes the interfacial shear force exerted by the wave field on the base film based on Inputs.Model.WAVEBASEINT model.

Inputs:

FWAVEMASS(base, drop, zIdx)

FWAVEMASS Base-wave mass exchange force [N/m^2]

Calculates the momentum exchange force due to mass transfer between base and wave films, based on velocity difference and MWAVE flux.

Inputs:

MDEP(base, drop, zIdx)

MDEP Base deposition mass flux [kg/m^2/s]

Calculates deposition flux to the base film from drops, scaled by deposition fraction.

Inputs:

MENT(base, zIdx)

MENT Base entrainment mass flux [kg/m^2/s]

Placeholder for future implementation if needed. Currently returns zero entrainment.

Inputs:

MEVAP(base, zIdx)

MEVAP Base evaporation mass flux [kg/m^2/s]

Computes evaporation flux from the base film using heat flux fraction and parent film evaporation rate.

Inputs:

MTOT(base, drop, zIdx)

MTOT Total mass flux of the base film [kg/m^2/s]

Calculates the total mass flux exchange associated with the base film, including evaporation, entrainment, deposition, and wave interaction.

Inputs:

MTURB(base, zIdx)

MTURB Turbulent mass exchange [kg/m^2/s]

Estimates turbulent mass exchange between base and wave films using wave mixing coefficient Inputs.Model.WAVEMIXCOEF.

Inputs:

MWAVE(base, drop, zIdx)

MWAVE Base film mass flux interaction with wave [kg/m^2/s]

Computes the mass flux exchange between the base film and wave field. Includes net exchange term based on relaxation toward equilibrium thickness and turbulent mixing contribution.

Inputs:

Outputs:

  • Mwave — Mass flux transferred to wave field [kg/m^2/s]

  • Mnet — Net exchange term before turbulent adjustment [kg/m^2/s]

TBASE(base, zIdx)

TBASE Base film exposure time [s]

Calculates the time during which the base film is exposed to vapor between successive waves, based on wave spacing and relative velocity (eq 24 of Le Corre [2022]).

Inputs:

TDRY(base, drop, zIdx)

TDRY Base film dry-out time [s]

Estimates the time of base film dry out under current evaporation and deposition conditions. Returns zero if dry-out does not occur.

Inputs:

THICKMIN(base, drop, zIdx)

THICKMIN Minimum base film thickness [m]

Calculates the minimum thickness of the base film at the end of its exposure period, based on mass flow rate, liquid density, velocity, and channel perimeter (eq. 28 of Le Corre [2022]).

Inputs:

WMIN(base, drop, zIdx)

WMIN Minimum base film mass flow rate [kg/s]

Calculates the minimum mass flow rate of the base film at the end of its exposure period (TBASE), accounting for evaporation and deposition.

Inputs:

WMINL(base, drop, zIdx)

WMINL Minimum base film mass flow rate per perimeter [kg/m-s]

Computes the minimum base film mass flow rate normalized by channel perimeter.

Inputs:

copyFlowProperties(srcObj, targetObj, opts)

COPYFLOWPROPERTIES Copy flow property arrays between Base objects

Copies selected flow properties (W, U, H) from a source object to one or more target Base objects. Supports full copy or partial copy (excluding first axial node) based on options.

Inputs:

  • srcObj — Source object containing flow properties

  • targetObj — Array of Solvers.FourField.Base objects to update

  • opts.all — Logical flag; if true, copies all axial nodes, otherwise skips the first node (default: false)

class Solvers.FourField.Drop

Bases: Solvers.ThreeField.Drop

DROP Class for modeling drop field in four-field solver

This class encapsulates the physical and numerical properties of the drop field, including flow variables and phase interactions. It supports multiple solver models and provides methods for computing derived quantities.

All methods are inherited from Solvers.ThreeField.ThreeFieldSolver