Three-field solver

The Solvers.ThreeField module in OpenSTREAM provides classes and tools for simulating annular two-phase flows using a three-field approach. In this method, the two liquid fields (film and droplets) are treated separately, 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 explicit modeling of liquid–liquid interactions, slip velocities, and hydrodynamic non-equilibrium effects, making it well-suited for applications involving annular boiling two-phase flow regimes where field separation and interfacial dynamics are significant.

This module includes:

Solvers.ThreeField.ThreeFieldSolver.solve(tfSolver)

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

Runs the full solution process for the Solvers.ThreeField.ThreeFieldSolver 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.ThreeField.ThreeFieldSolver.ThreeFieldSolver(inputSet, mixSolver)

Bases: Solvers.AbstractSolver

THREEFIELDSOLVER Solver for initializing, solving, and visualizing three-field flow.

The ThreeFieldSolver class handles the setup, execution, and visualization of thermal-hydraulic simulations using a two-phase three-field approach.

Responsibilities:

Key Components:

  • Liquid field construction (film/drop)

  • 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
ThreeFieldSolver(inputSet, mixSolver)

THREEFIELDSOLVER Constructor for the ThreeFieldSolver class.

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

Input:

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

Notes:

  • This constructor assumes Inputs.InputSet is fully validated.

  • Solver initialization includes liquid and vapor objects setup.

  • Time/space discretization and boundary condition interpolation initialized by Solvers.Mixture.MixtureSolver are used

Property Summary
DT = 0

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

DZ = 0

Axial step size [m]

NTIME = 0

Number of time steps [-]

NZ = 0

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

TIME = 0

Time series [s]

Z = 1.

Elevation [m]

boundaryConditions

Boundary conditions object Inputs.BoundaryConditions

drop

Drop object

dropInit

Null transient drop object

film

Film object

filmInit

Null transient film object

fluid

Fluid object Inputs.FluidProperties

inputSet

Input set object Inputs.InputSet

mixSolver

Mixture object Solvers.Mixture.Mixture

Method Summary
EQUIL(tfSolver, 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:

  • Enhanced drop deposition is disabled

  • Uses interpolation for improved convergence

  • Falls back to ad-hoc update if interpolation fails

  • Method update film and drop mass flow rates. Use flm.copy() and drp.copy() as input arguments if used in post-process

initializeSolver(tfSolver)

INITIALIZESOLVER Initializes solver parameters and constructs film/drop objects.

Sets up the solver’s internal state using the provided Inputs.InputSet and Solvers.ThreeField.ThreeFieldSolver, This includes initialization of data structures for transient and steady-state simulations based on the mixture solver solutions.

Operations performed:

  • Initializes film and drop arrays for each time step

  • Initializes fluid property objects

  • Sets up mixture object

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

Notes:

  • The function assumes uniform axial discretization.

  • Wall heat flux is computed from interpolated boundary conditions.

plott(tfSolver, zIdx, opt)

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

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

Inputs:

  • tfSolver — Solvers.ThreeField.ThreeFieldSolver 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(tfSolver, tIdx, opts)

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

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

Inputs:

  • tfSolver — Solvers.ThreeField.ThreeFieldSolver 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(tfSolver, opt)

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

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

Inputs:

  • tfSolver — Solvers.ThreeField.ThreeFieldSolver 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.ThreeField.Film

Bases: Solvers.AbstractFilm

FILM Class for modeling liquid film in three-field solver

This class encapsulates the physical and numerical properties of the film field, including flow variables and phase interactions.

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

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 = 1.

Velocity [m/s]

W = 1.

Mass flow rate [kg/s]

Z = 1.

Elevation [m]

class Solvers.ThreeField.Drop(inputSet, fluid)

Bases: Solvers.AbstractField

DROP Class for modeling drop field in three-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.

Constructor Summary
Drop(inputSet, fluid)

DROP Constructor for Drop class

Initializes the drop field object with input configuration and fluid properties.

Inputs:

Property Summary
DT = 0

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

H = 1E6

Enthalpy [J/kg]

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]

mix = NaN

Solvers.Mixture.Mixture object

Method Summary
AI(drop, zIdx)

AI Drop volumetric interfacial area [m^-1]

Calculates the interfacial area per unit volume for drops.

Inputs:

AREA(drop, zIdx)

AREA Drop interfacial area [m^2]

Calculates the surface area of a spherical drop.

Inputs:

CONC(drop, zIdx)

CONC Drop concentration [kg/m^3]

Calculates the local drop-field concentration based on mass flow and vapor-phase properties.

Inputs:

DENSITY(drop, zIdx)

DENSITY Drop number density [m^-3]

Calculates the number of drops per unit volume based on flow properties and geometry.

Inputs:

DIAM(drop, zIdx)

DIAM Drop diameter [m]

Returns the drop diameter based on Inputs.Model.DROPDIAM.

Inputs:

DRAG(drop, zIdx)

DRAG Drop drag coefficient [-]

Calculates drag coefficient using selected correlation based on Inputs.Model.DROPDRAG.

Inputs:

FBUOY(drop, zIdx)

FBUOY Drop buoyancy force [N/m^3]

Calculates buoyancy force based on pressure gradient.

Inputs:

FDRAG(drop, zIdx)

FDRAG Vapor drag force on drops [N/m^3]

Calculates shear force exerted by vapor on drops.

Inputs:

FENT(drop, film, zIdx)

FENT Film entrainment shear force [N/m^3]

Calculates shear force due to film entrainment acting on drops.

Inputs:

FGRAV(drop, zIdx)

FGRAV Gravitational force on drops [N/m^3]

Calculates gravity force acting on drops based on channel inclination angle.

Inputs:

FTOT(drop, film, zIdx, opt)

FTOT Total force acting on drops [N/m^3]

Calculates combined forces (drag, buoyancy, gravity, entrainment) based on selected option.

Inputs:

KENH(drop, zIdx)

KENH Drop deposition enhancement factor [-]

Calculates axial variation of deposition enhancement using Windecker correlation or other models based on Inputs.Model.DEPENHANCEMENT.

Inputs:

MDEP(drop, zIdx, enhanced)

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

Implements various deposition models selected based on Inputs.Model.DEPOSITION. Applies enhancement factors, and restricts deposition to annular flow region. The enhancement flag allow to enable/disable deposition enhancement

Inputs:

  • drop — Solvers.ThreeField.Drop object

  • zIdx — Axial indices to evaluate (optional)

  • enhanced — Flag indicating enhanced deposition downstream obstructions (optional, default = true)

RE(drop, zIdx)

RE Drop Reynolds number [-]

Calculates the Reynolds number for the drop field based on mass flow rate, viscosity, and channel perimeter.

Inputs:

REV(drop, zIdx)

REV Drop Reynolds number with respect to vapor properties [-]

Calculates Reynolds number using vapor density, viscosity, and relative phase velocity.

Inputs:

UALGEBR(drop, film, zIdx)

UALGEBR Algebraic drop velocity [m/s]

Calculates drop velocity consistent with mixture model based on void fraction and geometry.

Inputs:

UEQUIL(drop, film, zIdx, opt)

UEQUIL Drop velocity based on equilibrium model (Ftot = 0) [m/s]

Iteratively solves for drop velocity such that the total force acting on the drop phase is zero (force balance). Uses secant method with relaxation for convergence.

Inputs:

  • drop — Solvers.ThreeField.Drop object

  • film — Solvers.ThreeField.Film object

  • zIdx — Axial indices to evaluate (optional)

  • opt — Option flag for force balance:

    0 = full force balance (drag, buoyancy, gravity, entrainment) 1 = simplified (drag + buoyancy only)

USLIP(drop, zIdx)

USLIP Slip velocity model for drops [m/s]

Calculates drop velocity using slip ratio model based on Inputs.Model.DROPSLIP.

Inputs:

VOLUME(drop, zIdx)

VOLUME Drop volume [m^3]

Calculates the volume of a spherical drop.

Inputs:

VR(drop, zIdx)

VR Drop local relative velocity [m/s]

Calculates the velocity difference between vapor and drop fields.

Inputs:

Notes:

  • Only AREAMEAN option is implemented.

X(drop, zIdx)

X Vapor quality in the bulk (drop) region [-]

Vapor quality in the bulk region, considering only the droplet for the liquid phase. The bulk region is complementary of the near-wall region, defined in the mixture solver.

Inputs: