banner-transparent.png

Tutorial #1: Basic capabilities and functionalities

Introduction

Welcome to OpenSTREAM!
In this tutorial, you will learn how to use the basic features of the code, with provided examples to run a simple case using all four solvers.
The selected geometry consists of a 3.5 m uniformly heated vertical tube with an inner diameter of 8.8 mm, followed by a 2 m adiabatic section (of same diameter). The flow enters slightly subcooled, transitions to annular two-phase flow around 0.7 m, and reaches the critical boiling transition (film dryout) at the end of the heated section. In the adiabatic region, droplets deposit on the wall, regenerating the liquid film at the wall and leading to a fully developed annular flow.
Tutorial1.png
Table of Contents
Introduction Input files Mixture solver Two-fluid solver Three-field solver Four-field solver
% Clear existing variables from workspace
clearvars
Before starting this tutorial, make sure OpenSTREAM is included in your MATLAB search path. If you are running the tutorial directly from the tutorials folder, you can add OpenSTREAM like this:
% Add path to OpenSTREAM
osp = './..'; % When running from the tutorials folder
addpath(osp);

Input files

All the input files you will need for this tutorial are already provided in the openstream/inputs/ folder.
We first start by adding the inputs and solver packages to the import list.
% Import the inputs and solvers packages
import Inputs.*
import Solvers.*
The Inputs package defines a set of classes used to configure and import all the necessary parameters for running OpenSTREAM simulations. The Solvers package defines the OpenSTREAM core solver modules and their associated superclasses.
The path to the physical models, numerical options, test section geometry and boundary condition files are provided and read by InputSet as shown below. Since each file can include various set of inputs (several geometries, etc.), the IDs corresponding to the relevant sets must be indicated. New user-specific files can be created and the paths can be modified by the user, based on the specific working directory and file locations.
The models and options parameters not specified in the input files are set to their default values listed here.
% Turn off warnings during input setup
warning off
% Create the input set using paths to model, options, geometry, and boundary condition files
inputSet = InputSet( ...
modelFilePath = [osp '/inputs/models.inp'] , modelID = 'TUTORIAL1', ...
optionsFilePath = [osp '/inputs/options.inp'], optionsID = 'DEFAULT' , ...
geometryFilePath = [osp '/inputs/geom.inp'] , geometryID = 'TUTORIAL1', ...
bcFilePath = [osp '/inputs/tutorial1.inp'], ...
sessionParentDir = fullfile(pwd,'outputs') , ...
overwriteSessionFiles = true , ...
LOGMODE = 'BOTH');
% Turn off warnings during input setup
warning on
Note:
We can examine the physical model inputs by opening the models.inp file or using the command below, then scroll to locate the ID entries labeled TUTORIAL1.
type('../inputs/models.inp')
ID ! Models ID > TUTORIAL1 NNODES ! Number of axial nodes > 100 FLUID ! Fluid ID > WATER INTNU ! Interfacial heat transfer model > RELAXATION RELAXX ! Interfacial phase change relaxation time Xeq [-] > -0.5 1.0 RELAXTCOND ! Interfacial condensation relaxation time [s] > 0.5 0.5 MOMENTLIQUID ! Liquid momentum conservation model > FULL MOMENTGAS ! Gas momentum conservation model > FULL MOMENTFILM ! Film momentum conservation model > FULL MOMENTDROP ! Drop momentum conservation model > FULL OAFFILMSPLIT ! Film mass flow split model at onset of annular flow > RATIO WAVEBASEINT ! Wave / base film interfacial momentum transfer model > VAPORSHEARDROPMASS OAFBASERATIO ! Base/Film mass ratio at onset of annular flow [-] > 0.2 BASEEQTHICK ! Equlibrium base film thickness model > YPLUS END ID ! Model ID > TUTORIAL2A NNODES ! Number of axial nodes > 100 FLUID ! Fluid ID > WATER END ID ! Model ID > TUTORIAL2B NNODES ! Number of axial nodes > 100 FLUID ! Fluid ID > WATER VOID ! Void fraction model > SLIP SLIP ! Phase velocity ratio [-] > 2 END ID ! Model ID > TUTORIAL2C NNODES ! Number of axial nodes > 100 FLUID ! Fluid ID > WATER VOID ! Void fraction model > BESTION END ID ! Model ID > TUTORIAL2D NNODES ! Number of axial nodes > 100 FLUID ! Fluid ID > WATER VOID ! Void fraction model > BESTION THERMALNONEQ ! Thermal non-equilibrium model > SAHAZUBER END ID ! Model ID > TUTORIAL2E NNODES ! Number of axial nodes > 100 FLUID ! Fluid ID > WATER VOID ! Void fraction model > BESTION THERMALNONEQ ! Thermal non-equilibrium model > HRM END ID ! Model ID > TUTORIAL2F NNODES ! Number of axial nodes > 100 FLUID ! Fluid ID > WATER VOID ! Void fraction model > BESTION THERMALNONEQ ! Thermal non-equilibrium model > SAHAZUBER SPMTM ! Single-phase wall momentum transfer model > BLASIUS FRICTION ! Wall friction coefficients [-] > 0.2 -0.2 0 TPFM ! Two-phase friction multiplier model > EPRI END ID ! Model ID > TUTORIAL2G NNODES ! Number of axial nodes > 100 FLUID ! Fluid ID > WATER VOID ! Void fraction model > BESTION THERMALNONEQ ! Thermal non-equilibrium model > SAHAZUBER KLOC ! Elevation of local perturbations [m] > 1 2 3 4 5 KLOSS ! Pressure loss coefficients of local perturbations [-] > 0.5 0.5 0.5 0.5 0.5 TPKM ! Two-phase local loss multiplier model > HOMOGENEOUS END
This indicates that the simulations will use 100 nodes and that water is the working fluid.
Model inputs are mostly solver-specific and can be reviewed (along with the default models) at the following link. Note that each model is organized by solver within the Inputs.Model class definition file.
A simple and self-explanatory input file format has been defined for use with OpenSTREAM, as shown above (note that any input between ! and > is treated as a comment). Alternatively, a JSON file format can also be used if preferred.
type('../inputs/models.json')
[ { "ID": "TUTORIAL1", "NNODES": 100, "FLUID": "WATER", "INTNU": "RELAXATION", "RELAXX": [ -0.5, 1.0 ], "RELAXTCOND": [ 0.5, 0.5, 0.01 ], "MOMENTLIQUID": "FULL", "MOMENTGAS": "FULL", "MOMENTFILM": "FULL", "MOMENTDROP": "FULL", "WAVEBASEINT": "VAPORSHEARDROPMASS", "OAFFILMSPLIT": "RATIO", "OAFBASERATIO": 0.2, "BASEEQTHICK": "YPLUS" } ]
We can examine the numerical option inputs by opening the options.inp file or using the command below, then scroll to locate the ID entry labeled DEFAULT.
type('../inputs/options.inp')
ID ! Options ID > DEFAULT END
Note that no specific options have been defined, i.e., defaults setting are used for the simulations.
Options inputs are mostly solver-specific and can be reviewed (including the default options) at the following link.
We can examine the geometry inputs by opening the geom.inp file or using the command below, then scroll to locate the ID entry labeled TUTORIAL1.
type('../inputs/geom.inp')
ID ! Channel ID > TUTORIAL1 LENGTH ! Axial length [m] > 5.5 AREA ! Coolant area [m^2] > 6.0821e-05 PERIM ! Perimeters [m] > 0.0276 ANGLE ! Angle [rad] > 0 END
All inputs are required except for ANGLE, which defaults to 0, indicating a vertical (upward) test section. The separate specification of AREA and PERIM allows to define test sections of arbitrary shape. PERIM can be an array to define a multi-wall geometry (e.g., an annulus). In this tutorial, the area and perimeter corresponding to a circular test section with a diameter 8.8 mm were specified.
Geometry inputs are solver-independent and can be reviewed at the following link.
Finally, we can examine the boundary conditions by opening the tutorial1.inp file or using the command below.
type('../inputs/tutorial1.inp')
TIME ! Time [s] > 0.0 PRESSURE ! System pressure [Pa] > 6000000 HIN ! Inlet enthalpy [J/kg] > 1.1394E6 MFLOW ! Mass flow rate [kg/s] > 0.07 POWER ! Total power [W] > 87500 WMESH ! Relative power node size distribution [m] > 3.5 2.0 WPOWER ! Relative power distribution(s) [-] > 1.0 0.0 END
The file defines a single time step, indicating that the simulation is steady-state. The axial power distribution is specified using two arrays: WMESH for the node size array and WPOWER for the corresponding wall heat flux. In this tutorial, the configuration represents a 3.5 m uniformly heated section followed by a 2 m adiabatic section.
All inputs are required except for WMESH and WPOWER. If they are omitted, the test section is assumed to be uniformly heated by default.
Boundary conditions are solver-independent and can be reviewed at the following link.

Mixture solver

In this section, we will walk through how to set up and run the OpenSTREAM mixture solver using a simple example.
The physical models specified in the model.inp file do not affect any of the solver’s default settings for this solver, which hence reduces to a Homogeneous Equilibrium Model (HEM).
We first start by adding the mixture solver module to the import list.
% Import the mixture solver module
import Solvers.Mixture.*
The mixture solver module in OpenSTREAM provides classes and tools for simulating two-phase flows using a mixture approach.
A mixture solver object, mixSolver, is created from the inputSet (which contains all information related to geometry, boundary conditions, physical models and options).
% Create the mixture solver object
mixSolver = MixtureSolver(inputSet);
While the mixture solver object is initialized at construction. Here, it is explicitly initialized for clarity.
% Initialize solver
mixSolver.initializeSolver();
The mixture solver object contains several properties, including the inputSet object, boundaryConditions object, fluid object, and unsolved mixture object.
disp(mixSolver)
MixtureSolver with properties: NZ: 101 NTIME: 1 TIME: 0 DT: 0.1000 Z: [101x1 double] DZ: 0.0550 fluid: [1x1 Inputs.FluidProperties] boundaryConditions: [1x1 struct] mixtureInit: [1x30 Solvers.Mixture.Mixture] mixture: [1x1 Solvers.Mixture.Mixture] inputSet: [1x1 Inputs.InputSet] STATE: UNSOLVED
To solve the flow fields, simply call the solve method (no input arguments are needed). The solver uses a pseudo-transient approach with large time steps to drive the solution toward steady-state.
% Solve the flow fields
mixSolver.solve();
--------------------------------------------- Mixture solver run initiated --------------------------------------------- Solve steady-state ... Time 1.00 [s] max point iter = 5 in node 6, max errors: W = 0.0000001 [kg/s], P = 0.07547 [Pa], H = 0.08701 [J/kg] Time 2.00 [s] max point iter = 4 in node 5, max errors: W = 0.0000014 [kg/s], P = 0.09865 [Pa], H = 0.05141 [J/kg] Time 3.00 [s] max point iter = 4 in node 80, max errors: W = 0.0000002 [kg/s], P = 0.09691 [Pa], H = 0.06371 [J/kg] Time 4.00 [s] max point iter = 3 in node 5, max errors: W = 0.0000000 [kg/s], P = 0.00273 [Pa], H = 0.00318 [J/kg] Time 5.00 [s] max point iter = 3 in node 7, max errors: W = 0.0000000 [kg/s], P = 0.09601 [Pa], H = 0.02859 [J/kg] Time 6.00 [s] max point iter = 3 in node 14, max errors: W = 0.0000001 [kg/s], P = 0.08306 [Pa], H = 0.09930 [J/kg] Time 7.00 [s] max point iter = 3 in node 39, max errors: W = 0.0000000 [kg/s], P = 0.09984 [Pa], H = 0.09096 [J/kg] Time 8.00 [s] max point iter = 2 in node 26, max errors: W = 0.0000000 [kg/s], P = 0.05851 [Pa], H = 0.09187 [J/kg] Time 9.00 [s] max point iter = 2 in node 76, max errors: W = 0.0000000 [kg/s], P = 0.05091 [Pa], H = 0.09916 [J/kg] STEADY-STATE CONVERGED max errors: W = 0.0000000 [kg/s], P = 0.05626 [Pa], H = 0.17789 [J/kg] Elapsed time: 3.16 sec --------------------------------------------- Mixture solver run completed --------------------------------------------- Output directory: /san/reactorsF/job1/2018p4151/mfval/Git/openstream/tutorials/outputs/TUTORIAL1-TUTORIAL1-DEFAULT
At each time step, the solver prints the results of the point-wise iterations (number of iterations and the residuals for mixture mass flow rate, pressure, and enthalpy). It then checks whether the solution has converged for that time step. If not, the solver continues to the next time step. If steady-state convergence is reached, it prints a summary of the time step convergence.
The relevant properties of the mixture object (W, P, H) are now populated with the converged solution. For reference, the mixtureInit object contains the transient iterations toward steady-state convergence of the first time step. That is, in a steady-state simulation, the mixture object contains the last time step of the mixtureInit object.
disp(mixSolver)
MixtureSolver with properties: NZ: 101 NTIME: 1 TIME: 0 DT: 0.1000 Z: [101x1 double] DZ: 0.0550 fluid: [1x1 Inputs.FluidProperties] boundaryConditions: [1x1 struct] mixtureInit: [1x10 Solvers.Mixture.Mixture] mixture: [1x1 Solvers.Mixture.Mixture] inputSet: [1x1 Inputs.InputSet] STATE: INITIALSTEPCONVERGED
We can now look at the solved code outputs. Any properties or method parameters from the mixture object can be easily accessed. For instance, the void fraction array can be printed as follows:
% Display the void fraction
disp(mixSolver.mixture.VF')
0 0 0 0 0.0707 0.2820 0.4178 0.5124 0.5822 0.6357 0.6781 0.7124 0.7409 0.7648 0.7852 0.8028 0.8182 0.8317 0.8437 0.8544 0.8639 0.8726 0.8805 0.8876 0.8942 0.9002 0.9058 0.9109 0.9157 0.9201 0.9242 0.9281 0.9317 0.9351 0.9383 0.9413 0.9441 0.9468 0.9494 0.9518 0.9541 0.9562 0.9583 0.9603 0.9622 0.9640 0.9657 0.9673 0.9689 0.9704 0.9719 0.9733 0.9746 0.9759 0.9771 0.9783 0.9795 0.9806 0.9817 0.9827 0.9837 0.9847 0.9856 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865 0.9865
In addition, the mixture solver class includes built-in plotting methods to help visualize the results. Here are a few examples you can try:
% Plot the axial distribution of mass flow rates
mixSolver.plotz('display','W');
Vapor generation begins near the inlet of the test section as the mixture reaches saturation. From there, vapor mass increases linearly along the heated length due to the uniform heat flux, and then remains constant throughout the adiabatic section. This is a typical behavior expected for a homogeneous equilibrium model (HEM).
% Plot the axial distribution of vapor ratios
mixSolver.plotz('display','VR');
The figure confirms that the mixture becomes saturated () as the steam mass (and void fraction) starts increasing.
Results can be saved if desired.
% Save results
mixSolver.save('name','mixSolver','showpath',false);

Two-fluid solver

In this section, we will walk through how to set up and run the OpenSTREAM two-fluid solver using a simple example.
The physical models specified in the model.inp file (displayed previously) set the phase velocity calculations to a full momentum conservation approach (MOMENTGAS and MOMENTLIQUID set to FULL). In addition, the interfacial transfer model is simplified by a applying a time relaxation approach (INTNU set to RELAXATION).
We first start by adding the two-fluid solver module to the import list.
% Import the two-fluid solver module
import Solvers.TwoFluid.*
The two-fluid solver module in OpenSTREAM provides classes and tools for simulating two-phase flows using a two-fluid approach.
Since the input parameters are already loaded, we can reuse the same input set without re-reading the files.
We create a two-fluid solver object, twfSolver, using the inputSet and the previously solved mixSolver object. This lets the two-fluid solver starts from a physically reasonable state and inherits the pressure gradient already computed by the mixture solver.
% Create the two-fluid solver object
twfSolver = TwoFluidSolver(inputSet, mixSolver);
In case a solved mixSolver object is not specified, it will be generated automatically.
The two-fluid solver object contain several properties, including the inputSet object, boundaryConditions object, fluid object, solved mixSolver object, and unsolved liquid and vapor objects.
disp(twfSolver)
TwoFluidSolver with properties: NZ: 101 NTIME: 1 TIME: 0 DT: 0.1000 Z: [101x1 double] DZ: 0.0550 fluid: [1x1 Inputs.FluidProperties] boundaryConditions: [1x1 struct] liquidInit: [1x30 Solvers.TwoFluid.Liquid] vaporInit: [1x30 Solvers.TwoFluid.Vapor] liquid: [1x1 Solvers.TwoFluid.Liquid] vapor: [1x1 Solvers.TwoFluid.Vapor] mixSolver: [1x1 Solvers.Mixture.MixtureSolver] inputSet: [1x1 Inputs.InputSet] STATE: UNSOLVED
To solve the flow fields, simply call the solve method (no input arguments are needed). The solver uses a pseudo-transient approach with large time steps to drive the solution toward steady-state.
% Solve the flow fields
twfSolver.solve();
-------------------------------------------- Two-fluid solver run initiated -------------------------------------------- Solve steady-state ... Time 1.00 [s] max point iter = 37 in node 14, max errors: Wl = 0.0000000 [kg/s], Wv = 0.0000001 [kg/s], Ul = 0.0000999 [m/s], Uv = 0.0000911 [m/s], Hl = 0.08613 [J/kg], Hv = 0.00000 [J/kg] Time 2.00 [s] max point iter = 28 in node 14, max errors: Wl = 0.0000002 [kg/s], Wv = 0.0000001 [kg/s], Ul = 0.0000994 [m/s], Uv = 0.0000961 [m/s], Hl = 0.04708 [J/kg], Hv = 0.00000 [J/kg] Time 3.00 [s] max point iter = 19 in node 14, max errors: Wl = 0.0000001 [kg/s], Wv = 0.0000001 [kg/s], Ul = 0.0000988 [m/s], Uv = 0.0000963 [m/s], Hl = 0.08815 [J/kg], Hv = 0.00000 [J/kg] Time 4.00 [s] max point iter = 10 in node 65, max errors: Wl = 0.0000000 [kg/s], Wv = 0.0000001 [kg/s], Ul = 0.0000987 [m/s], Uv = 0.0000981 [m/s], Hl = 0.07355 [J/kg], Hv = 0.00000 [J/kg] Time 5.00 [s] max point iter = 8 in node 96, max errors: Wl = 0.0000000 [kg/s], Wv = 0.0000000 [kg/s], Ul = 0.0000992 [m/s], Uv = 0.0000967 [m/s], Hl = 0.05329 [J/kg], Hv = 0.00000 [J/kg] Time 6.00 [s] max point iter = 4 in node 78, max errors: Wl = 0.0000000 [kg/s], Wv = 0.0000001 [kg/s], Ul = 0.0000989 [m/s], Uv = 0.0000996 [m/s], Hl = 0.08611 [J/kg], Hv = 0.00000 [J/kg] Time 7.00 [s] max point iter = 2 in node 33, max errors: Wl = 0.0000002 [kg/s], Wv = 0.0000000 [kg/s], Ul = 0.0000979 [m/s], Uv = 0.0000987 [m/s], Hl = 0.09880 [J/kg], Hv = 0.00000 [J/kg] Time 8.00 [s] max point iter = 2 in node 51, max errors: Wl = 0.0000001 [kg/s], Wv = 0.0000000 [kg/s], Ul = 0.0000960 [m/s], Uv = 0.0000962 [m/s], Hl = 0.09986 [J/kg], Hv = 0.00000 [J/kg] Time 9.00 [s] max point iter = 2 in node 48, max errors: Wl = 0.0000000 [kg/s], Wv = 0.0000000 [kg/s], Ul = 0.0000938 [m/s], Uv = 0.0000936 [m/s], Hl = 0.09890 [J/kg], Hv = 0.00000 [J/kg] Time 10.00 [s] max point iter = 2 in node 51, max errors: Wl = 0.0000000 [kg/s], Wv = 0.0000000 [kg/s], Ul = 0.0000920 [m/s], Uv = 0.0000915 [m/s], Hl = 0.09909 [J/kg], Hv = 0.00000 [J/kg] STEADY-STATE CONVERGED max errors: Wl = 0.0000000 [kg/s], Wv = 0.0000001 [kg/s], Ul = 0.0002758 [m/s], Uv = 0.0002769 [m/s], Hl = 0.69621 [J/kg], Hv = 0.00000 [J/kg] Elapsed time: 6.81 sec -------------------------------------------- Two-fluid solver run completed -------------------------------------------- Output directory: /san/reactorsF/job1/2018p4151/mfval/Git/openstream/tutorials/outputs/TUTORIAL1-TUTORIAL1-DEFAULT
At each time step, the solver prints the results of the point-wise iterations (number of iterations and the residuals for phase mass flow rates, velocities, and enthalpies. It then checks whether the solution has converged for that time step. If not, the solver continues to the next time step. If steady-state convergence is reached, it prints a summary of the time step convergence.
The relevant properties of the liquid and vapor objects (W, U, H) are now populated with the converged solution. For reference, the liquidInit and vaporInit objects contain the transient iterations toward steady-state convergence of the first time step. That is, in a steady-state simulation, the liquid and vapor objects contains the last time step of the mixtureInit and vaporInit objects.
disp(twfSolver)
TwoFluidSolver with properties: NZ: 101 NTIME: 1 TIME: 0 DT: 0.1000 Z: [101x1 double] DZ: 0.0550 fluid: [1x1 Inputs.FluidProperties] boundaryConditions: [1x1 struct] liquidInit: [1x11 Solvers.TwoFluid.Liquid] vaporInit: [1x11 Solvers.TwoFluid.Vapor] liquid: [1x1 Solvers.TwoFluid.Liquid] vapor: [1x1 Solvers.TwoFluid.Vapor] mixSolver: [1x1 Solvers.Mixture.MixtureSolver] inputSet: [1x1 Inputs.InputSet] STATE: INITIALSTEPCONVERGED
The two-fluid solver class includes built-in plotting methods to help visualize the results. Here are a few examples you can try:
% Plot the axial distribution of mass flow rates
twfSolver.plotz('display','W');
The solution is the same to what the mixture solver provides, but with one key difference: vapor generation starts further upstream in the subcooled region due to non-equilibrium wall boiling and interfacial condensation.
% Plot the axial distribution of vapor ratios
twfSolver.plotz('display','VR');
The observation from the previous plot is confirmed. Near the inlet of the test section, the vapor mass quality begins to increase even though the thermodynamic equilibrium quality is still negative, characteristic of subcooled boiling.
Results can be saved if desired.
% Save results
twfSolver.save('name','twfSolver','showpath',false);

Three-field solver

In this section, we will walk through how to set up and run the OpenSTREAM three-field solver using a simple example.
The physical models specified in the model.inp file (displayed previously) set the drop velocity calculations to a full momentum conservation model (MOMENTDROP set to FULL). Note that thermal equilibrium is assumed, hence interfacial heat transfer models are not relevant.
We first start by adding the three-field solver module to the import list.
% Import the three-field solver module
import Solvers.ThreeField.*
The three-field solver module in OpenSTREAM provides classes and tools for simulating annular two-phase flows using a three-field approach.
Since the input parameters are already loaded, we can reuse the same input set without re-reading the files.
We create a three-field solver object, tfSolver, using inputSet and the previously solved mixSolver object. This lets the three-field solver starts from a physically reasonable state and inherits the pressure gradient already computed by the mixture solver.
% Create the three-field solver object
tfSolver = ThreeFieldSolver(inputSet, mixSolver);
In case a solved mixSolver object is not specified, it will be generated it automatically.
The three-field solver object contains several properties, including the inputSet object, boundaryConditions object, fluid object, solved mixSolver object, and unsolved liquid film and drop objects. The vapor solution is the same as for the mixture solver.
disp(tfSolver)
ThreeFieldSolver with properties: NZ: 101 NTIME: 1 TIME: 0 DT: 0.1000 Z: [101x1 double] DZ: 0.0550 fluid: [1x1 Inputs.FluidProperties] boundaryConditions: [1x1 struct] filmInit: [1x30 Solvers.ThreeField.Film] dropInit: [1x30 Solvers.ThreeField.Drop] film: [1x1 Solvers.ThreeField.Film] drop: [1x1 Solvers.ThreeField.Drop] mixSolver: [1x1 Solvers.Mixture.MixtureSolver] inputSet: [1x1 Inputs.InputSet] STATE: UNSOLVED
To solve the flow fields, simply call the solve method (no input arguments are needed). The solver uses a pseudo-transient approach with large time steps to drive the solution toward steady-state.
% Solve the flow fields
tfSolver.solve();
------------------------------------------- Three-field solver run initiated ------------------------------------------- Solve steady-state ... Time 1.00 [s] max point iter = 16 in node 63, max errors: Wf = 0.0000985 [kg/s/m], Uf = 0.00629 [m/s], Ud = 0.00999 [m/s] Time 2.00 [s] max point iter = 8 in node 18, max errors: Wf = 0.0000986 [kg/s/m], Uf = 0.00020 [m/s], Ud = 0.00950 [m/s] Time 3.00 [s] max point iter = 5 in node 24, max errors: Wf = 0.0000981 [kg/s/m], Uf = 0.00010 [m/s], Ud = 0.00769 [m/s] Time 4.00 [s] max point iter = 2 in node 23, max errors: Wf = 0.0000975 [kg/s/m], Uf = 0.00007 [m/s], Ud = 0.00725 [m/s] Time 5.00 [s] max point iter = 1 in node 2, max errors: Wf = 0.0000843 [kg/s/m], Uf = 0.00007 [m/s], Ud = 0.00681 [m/s] STEADY-STATE CONVERGED max errors: Wf = 0.0000843 [kg/s/m], Uf = 0.00007 [m/s], Ud = 0.00681 [m/s] Elapsed time: 2.13 sec ------------------------------------------- Three-field solver run completed ------------------------------------------- Output directory: /san/reactorsF/job1/2018p4151/mfval/Git/openstream/tutorials/outputs/TUTORIAL1-TUTORIAL1-DEFAULT
At each time step, the solver prints the results of the point-wise iterations (number of iterations and the residuals for liquid film mass flow rate, liquid film velocity, and droplet velocity. It then checks whether the solution has converged for that time step. If not, the solver continues to the next time step. If steady-state convergence is reached, it prints a summary of the time step convergence.
The relevant properties of the liquid film and drop objects (W, U) are now populated with the converged solution. For reference, the filmInit and dropInit objects contain the transient iterations toward steady-state convergence of the first time step. That is, in a steady-state simulation, the liquid film and drop objects contains the last time step of the filmInit and dropInit objects.
disp(tfSolver)
ThreeFieldSolver with properties: NZ: 101 NTIME: 1 TIME: 0 DT: 0.1000 Z: [101x1 double] DZ: 0.0550 fluid: [1x1 Inputs.FluidProperties] boundaryConditions: [1x1 struct] filmInit: [1x6 Solvers.ThreeField.Film] dropInit: [1x6 Solvers.ThreeField.Drop] film: [1x1 Solvers.ThreeField.Film] drop: [1x1 Solvers.ThreeField.Drop] mixSolver: [1x1 Solvers.Mixture.MixtureSolver] inputSet: [1x1 Inputs.InputSet] STATE: INITIALSTEPCONVERGED
The three-field solver class includes built-in plotting methods to help visualize the results. Here are a few examples you can try:
% Plot the axial distribution of mass flow rates
tfSolver.plotz('display','W');
The three-field solution is shown only downstream the onset of annular flow conditions (marked by a vertical dashed line).
The solution is the same to what the mixture solver provides, but with one key difference: the liquid mass flow rate is split into film and droplets.
As mentioned in the tutorial description, the liquid film mass flow rate is observed to reduce to 0 at the end of heated section (3.5 m). Further downstream, the film is gradually regenerated as droplet deposit onto the wall along the adiabatic section.
% Plot the axial distribution of velocities
tfSolver.plotz('display','U');
As expected, the film velocity is significantly lower than the droplet velocity. The droplet velocity is slightly lower than the vapor velocity.
Results can be saved if desired.
% Save results
tfSolver.save('name','tfSolver','showpath',false);

Four-field solver

In this section, we will walk through how to set up and run the OpenSTREAM four-field solver using a simple example.
The physical models specified in the model.inp file (displayed previously) modify several default settings related to base film/wave mass exchange, film/wave distribution at the onset of annular flow, and interfacial momentum exchange for the base film. These parameters will be covered in more detail in advanced tutorials dedicated to the four-field solver.
We first start by adding the four-field solver module to the import list.
% Import the four-field solver module
import Solvers.FourField.*
The four-field solver module in OpenSTREAM provides tools and classes for simulating annular two-phase flows using a four-field approach.
The input parameters are unchanged, hence the input set does not need to be read again.
We create a four-field solver object, ffSolver, using inputSet and the previously solved mixSolver object. This lets the four-field solver starts from a physically reasonable state and inherits the pressure gradient already computed by the mixture solver.
% Create the four-field solver object
ffSolver = FourFieldSolver(inputSet, mixSolver);
In case a solved mixSolver object is not specified, it will be generated it automatically.
The four-field solver object contain several properties, including the inputSet object, boundaryConditions object, fluid object, solved mixSolver object, and unsolved liquid film and drop objects. The film object contains the wave and base objects (describing the disturbance waves and base film, respectively). The vapor solution is the same as for the mixture solver.
disp(ffSolver)
FourFieldSolver with properties: NZ: 101 NTIME: 1 TIME: 0 DT: 0.1000 Z: [101x1 double] DZ: 0.0550 fluid: [1x1 Inputs.FluidProperties] boundaryConditions: [1x1 struct] filmInit: [1x30 Solvers.FourField.Film] dropInit: [1x30 Solvers.FourField.Drop] film: [1x1 Solvers.FourField.Film] drop: [1x1 Solvers.FourField.Drop] mixSolver: [1x1 Solvers.Mixture.MixtureSolver] inputSet: [1x1 Inputs.InputSet] STATE: UNSOLVED
To solve the flow fields, simply call the solve method (no input arguments are needed). The solver uses a pseudo-transient approach with large time steps to drive the solution toward steady-state.
% Solve the flow fields
ffSolver.solve();
------------------------------------------- Four-field solver run initiated -------------------------------------------- Solve steady-state ... Time 1.00 [s] max point iter = 100 in node 14, max errors: Wb = 0.0007327 [kg/s/m], Ub = 0.67180 [m/s], Ww = 0.0005019 [kg/s/m], Uw = 0.00981 [m/s], Fw = 0.00649 [Hz], Ud = 0.00998 [m/s] Time 2.00 [s] max point iter = 100 in node 12, max errors: Wb = 0.0014390 [kg/s/m], Ub = 0.76035 [m/s], Ww = 0.0009916 [kg/s/m], Uw = 0.00578 [m/s], Fw = 0.09903 [Hz], Ud = 0.00964 [m/s] Time 3.00 [s] max point iter = 100 in node 14, max errors: Wb = 0.0008920 [kg/s/m], Ub = 0.74878 [m/s], Ww = 0.0006128 [kg/s/m], Uw = 0.00440 [m/s], Fw = 0.09970 [Hz], Ud = 0.00777 [m/s] Time 4.00 [s] max point iter = 100 in node 14, max errors: Wb = 0.0008376 [kg/s/m], Ub = 0.73718 [m/s], Ww = 0.0005745 [kg/s/m], Uw = 0.00372 [m/s], Fw = 0.09320 [Hz], Ud = 0.00728 [m/s] Time 5.00 [s] max point iter = 100 in node 14, max errors: Wb = 0.0008239 [kg/s/m], Ub = 0.73393 [m/s], Ww = 0.0005648 [kg/s/m], Uw = 0.00343 [m/s], Fw = 0.08577 [Hz], Ud = 0.00680 [m/s] Time 6.00 [s] max point iter = 100 in node 14, max errors: Wb = 0.0008225 [kg/s/m], Ub = 0.73360 [m/s], Ww = 0.0005638 [kg/s/m], Uw = 0.00325 [m/s], Fw = 0.07873 [Hz], Ud = 0.00636 [m/s] Time 7.00 [s] max point iter = 100 in node 14, max errors: Wb = 0.0008224 [kg/s/m], Ub = 0.73357 [m/s], Ww = 0.0005637 [kg/s/m], Uw = 0.00300 [m/s], Fw = 0.07216 [Hz], Ud = 0.00593 [m/s] STEADY-STATE CONVERGED max errors: Wb = 0.0000964 [kg/s/m], Ub = 0.00032 [m/s], Ww = 0.0000549 [kg/s/m], Uw = 0.00300 [m/s], Fw = 0.07216 [Hz], Ud = 0.00593 [m/s] Elapsed time: 16.22 sec ------------------------------------------- Four-field solver run completed -------------------------------------------- Output directory: /san/reactorsF/job1/2018p4151/mfval/Git/openstream/tutorials/outputs/TUTORIAL1-TUTORIAL1-DEFAULT
At each time step, the solver prints the results of the point-wise iterations (number of iterations and the residuals for liquid base film mass flow rate and velocity, wave mass flow rate and velocity, and drop velocity. It then checks whether the solution has converged for that time step. If not, the solver continues to the next time step. If steady-state convergence is reached, it prints a summary of the time step convergence.
The relevant properties of the wave, base and drop objects (W, U) are now populated with the converged solution. For reference, the filmInit and dropInit objects contain the transient iterations toward steady-state convergence of the first time step. That is, in a steady-state simulation, the liquid film and drop objects contains the last time step of the filmInit and dropInit objects.
disp(ffSolver)
FourFieldSolver with properties: NZ: 101 NTIME: 1 TIME: 0 DT: 0.1000 Z: [101x1 double] DZ: 0.0550 fluid: [1x1 Inputs.FluidProperties] boundaryConditions: [1x1 struct] filmInit: [1x8 Solvers.FourField.Film] dropInit: [1x8 Solvers.FourField.Drop] film: [1x1 Solvers.FourField.Film] drop: [1x1 Solvers.FourField.Drop] mixSolver: [1x1 Solvers.Mixture.MixtureSolver] inputSet: [1x1 Inputs.InputSet] STATE: INITIALSTEPCONVERGED
The four-field solver class includes built-in plotting methods to help visualize the results. Here are a few examples you can try:
% Plot the axial distribution of mass flow rates
ffSolver.plotz('display','W');
The four-field solution is shown only downstream the onset of annular flow conditions (marked by a vertical dashed line).
The solution is the same to what the three-field solver provides, but with one key difference: the liquid film mass flow rate is split into base film and disturbance waves. Away from the film dryout region, it can be observed that most the the film mass is transported by the waves (consistent with experimental observations).
% Plot the axial distribution of velocities
ffSolver.plotz('display','U');
As expected, the disturbance wave velocity is higher than the base film velocity, except near film dryout (where the velocities are equal). The droplet velocity is slightly lower than the vapor velocity.
% Plot the axial distribution of wave frequency
ffSolver.plotz('display','FREQUENCY');
The figure shows the predicted disturbance wave frequency, alongside the result from the reference equilibrium model. As expected, the wave frequency initially lags behind its equilibrium value and gradually approaches it along the adiabatic section.
Results can be saved if desired.
% Save results
ffSolver.save('name','ffSolver','showpath',false);