banner-transparent.png

Tutorial #2: Mixture solver

Introduction

In this tutorial, you will learn how to use the OpenSTREAM mixture solver, with examples demonstrating various capabilities.
The selected geometry consist 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 Homogeneous Equilibrium Model (HEM) Non-homogeneous models Slip void model Drift flux void model Thermal non-equilibrium models Empirical model Homogeneous Relaxation Model (HRM) Pressure drop models Axial pressure drop components Wall friction Singular pressure losses Transient
% 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.*
We then add 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.
We can examine the physical model inputs by opening the model.inp file or using the command below, then scroll to locate the ID entries labeled TUTORIAL2*.
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
These indicate that the simulations will use 100 nodes and that water is the working fluid. Models not specificed are set to their default values by the code. Default models are listed here and are set within the Inputs.Model class definition file.
In this tutorial , the following model sets will be used:
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. Default options are listed here and are set within the Inputs.Options class definition file.
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 specified in SI units and 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. In this case, the area and perimeter corresponding to a circular test section with a diameter 8.8 mm were specified.
Finally, we can examine the boundary conditions by opening the boundary conditions 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
Defining a single time step indicates that the simulation is steady-state. The axial power distribution is specified using two arrays: WMESH for node size and WPOWER for corresponding wall heat flux. In this case, the configuration represents a 3.5 m uniformly heated section followed by a 2 m adiabatic section.
All inputs are specified in SI units and required except for WMESH and WPOWER. If they are omitted, the test section is assumed to be uniformly heated by default.

Homogeneous Equilibrium Model (HEM)

In this section, we will walk through how to set up and run the mixture solver using the simplest set of physical models available.
By default, the physical models in the mixture solver are consistent with the Homogeneous Equilibrium Model (HEM). That is, hydrodynamic and thermal equilibrium is assumed along with no slip between the phases.
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
inputSetA = InputSet( ...
modelFilePath = [osp '/inputs/models.inp'] , modelID = 'TUTORIAL2A', ...
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');
Warning: /san/reactorsF/job1/2018p4151/mfval/Git/openstream/tutorials/outputs/TUTORIAL1-TUTORIAL2A-DEFAULT was removed.
Warning:
Default values were used for the following variables in INPUTS.MODEL:
PROPERTIES : SATURATED
ANGLE : 0.000000
KLOC : [0.000000 0.000000 ]
WBOILINGXSUB : -0.200000
WBOILINGN : 2.000000
WBOILINGXSAT : 0.000000
OAF : WALLIS
OAFTRANSITION : [0.100000 0.000000 ]
CBT : NONE
CBTMULT : @(z)1
CBTKEFFECT : [0.000000 0.000000 ]
CBTELEVATION : 1.000000
MFBT : NONE
DTMFB : 100.000000
SPMTM : BLASIUS
FRICTION : [0.200000 -0.200000 0.000000 ]
TPFM : HOMOGENEOUS
KLOSS : [0.000000 0.000000 ]
TPKM : HOMOGENEOUS
BTMTM : TPFM
SPHTM : DITTUSBOELTER
DITTUSBOELTERCOEF : [0.023000 0.800000 0.400000 ]
TPHTM : THOM
BTHTM : VAPOR
VOID : HOMOGENEOUS
SLIP : 1.000000
THERMALNONEQ : EQUILIBRIUM
THERMALRELAX : QUALITY
RELAXX : [-0.500000 -0.250000 -0.100000 0.000000 1.000000 ]
RELAXTCOND : [1.000000 0.500000 0.300000 0.100000 0.100000 ]
RELAXTEVAP : [0.300000 0.300000 0.300000 0.300000 0.300000 ]
RELAXCONDCOEF : [0.000100 0.333333 0.050000 0.000000 ]
RELAXEVAPCOEF : [0.000100 0.333333 0.000010 0.000000 ]
KTRELAX : NaN
NEARWALLRATIO : 0.500000
NEARWALLEQOAF :
NEARWALLH : RATIO
NEARWALLHRATIO : 1.000000
NEARWALLRELAX : QUALITY
NEARWALLRELAXX : [-0.500000 -0.250000 -0.100000 0.000000 1.000000 ]
NEARWALLRELAXT : [1.000000 0.500000 0.300000 0.100000 0.100000 ]
NEARWALLRELAXCOEF : [0.000100 0.333333 0.050000 ]
INTLENGTH : CONSTANT
INTAREA : DISPGAS2DISPLIQ
INTLENGTHVCST : 0.002000
INTLENGTHLCST : 0.002000
LOCRELVEL : AREAMEAN
RELVELCST : 1.000000
MOMENTLIQUID : MIXTURE
MOMENTGAS : MIXTURE
DROPDRAG : VISCOUS
DROPDRAGCOEF : 0.450000
BUBBLEDRAG : STOKES
BUBBLEDRAGCOEF : 0.450000
INTNU : RELAXATION
INTNUVCST : 2.000000
INTNULCST : 2.000000
RANZMARSHALLVCST : [2.000000 0.600000 0.500000 0.333333 ]
RANZMARSHALLLCST : [2.000000 0.600000 0.500000 0.333333 ]
INTTRANSH : BULK
POSFILM : 
OAFENTRAINED : EQUILIBRIUM
OAFDROPRATIO : 0.700000
DEPOSITION : OKAWA
DEPENHANCEMENT : NONE
KBLOCKRATIO : [0.000000 0.000000 ]
KTUNING : [0.000000 0.000000 ]
ENTRAINMENT : OKAWA2003
OKAWACOEFS : [320.000000 0.111000 0.000479 1.000000 ]
MOMENTFILM : ALGEBRAIC
VAPORFRIC : SOLVER_DEPENDENT
VAPORFRICCST : 0.005000
THINFILMFRIC : TURBULENT
THINFILMTHICK : 0.000100
MOMENTDROP : SLIP
DROPSLIP : 1.000000
DROPDIAM : 0.001000
OAFFILMSPLIT : EQUILIBRIUM
OAFBASERATIO : 0.500000
BASEEQTHICK : RISO
BASEEQTHICKCOEF : [0.000054 -0.640000 1.210000 ]
BASEYPLUS : 15.000000
RELAXTB : 0.200000
WAVEMIXCOEF : 2.000000
WAVEFREQUENCY : RELAXATION
EQSTROUHAL : RISO
EQSTROUHALCOEF : [0.000112 0.500000 0.000000 ]
RELAXTW : 0.200000
CSTWAVEFREQ : 100.000000
MOMENTBASE : FULLNOP
MOMENTWAVE : FULL
WAVEBASEINT : VAPORSHEAR
SHAPEFACTORCOEF : [132500.000000 2.000000 ]
WAVEDRAGCOEF : [0.020000 135000.000000 0.437000 ]
THINWAVETHICK : 0.000010
Warning:
Default values were used for the following variables in INPUTS.OPTIONS:
AXIALINTERP : next
TIMEINTERP : linear
TSTEP : 0.100000
MAXITER : 100.000000
SSTSTEP : 1.000000
SSMAXITER : 30.000000
ERRORW : 0.001000
ERRORP : 0.100000
ERRORH : 0.100000
SSCONVW : 0.001000
SSCONVP : 1.000000
SSCONVH : 1.000000
RELAXWM : 1.000000
RELAXPM : 1.000000
RELAXHM : 1.000000
RELAXWV : 0.800000
RELAXHV : 0.800000
ERRORU : 0.000100
SSCONVU : 0.001000
RELAXWL : 1.000000
RELAXUL : 0.500000
RELAXUV : 0.500000
RELAXHL : 1.000000
ERRORWF : 0.000100
ERRORUF : 0.010000
ERRORUD : 0.010000
SSCONVWF : 0.000100
SSCONVUF : 0.010000
SSCONVUD : 0.010000
RELAXWF : 0.500000
RELAXUF : 0.200000
RELAXUD : 0.200000
ERRORFW : 0.100000
SSCONVFW : 0.100000
RELAXWB : 0.500000
RELAXUB : 0.200000
RELAXWW : 0.500000
RELAXUW : 0.200000
RELAXFW : 0.500000
% Turn off warnings during input setup
%warning on
Note:
A mixture solver object, mixSolverA, is created from inputSetA (which contains all information related to geometry, boundary conditions, physical models and options).
% Create the mixture solver object
mixSolverA = MixtureSolver(inputSetA);
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
mixSolverA.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: 2.80 sec --------------------------------------------- Mixture solver run completed --------------------------------------------- Output directory: /san/reactorsF/job1/2018p4151/mfval/Git/openstream/tutorials/outputs/TUTORIAL1-TUTORIAL2A-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 mixture solver class includes built-in plotting methods to help visualize the results. Not indicating any argument (mixSolver.plotz()) will generate figures for most relevant parameters. Various name-value arguments are available, such as display to select the parameter(s) to be plotted. Further information how to use the plotz method can be found here.
We start by plotting the axial distribution of wall heat flux along the test section.
% Plot the axial distribution of heat flux
mixSolverA.plotz('display','HFLUX');
The wall heat flux is imposed by the WMESH and WPOWER inputs in the boundary conditions input file. Here we verify that uniform power is applied along the first 3.5 m, followed by a 2 m long adiabatic region.
We can then plot the axial distributions of mixture and phase mass flow rates.
% Plot the axial distribution of mass flow rates
mixSolverA.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).
Finally, we plot the axial distributions of vapor volumetric (void) and mass (quality) ratios.
% Plot the axial distributions of vapor ratio
mixSolverA.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 the solver object
mixSolverA.save('name','mixSolverA','showpath',false)
The solver object can be retrieved later by using the load command.

Non-homogeneous models

In this section, two non-homogeneous models are demonstrated where the phases have different velocities. This is achieved by calculating the void fraction based either on the slip void model or a drift flux void model, instead of the homogeneous void model.

Slip void model

In this section, we will walk through how to set up and run the mixture solver using a simple slip void model and a phase velocity ratio set to 2.
The only difference with the previous case is that input VOID is set to SLIP and input SLIP is set to 2 in the model set.
The path to the physical models, numerical options, test section geometry and boundary condition files are provided and read by InputSet. The models ID corresponding to the void slip model (TUTORIAL2B) is indicated.
% Turn off warnings during input setup
warning off
% Create the input set using paths to model, options, geometry, and boundary condition files
inputSetB = InputSet( ...
modelFilePath = [osp '/inputs/models.inp'] , modelID = 'TUTORIAL2B', ...
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
A mixture solver object, mixSolverB, is created from inputSetB.
% Create the mixture solver object
mixSolverB = MixtureSolver(inputSetB);
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
mixSolverB.solve();
--------------------------------------------- Mixture solver run initiated --------------------------------------------- Solve steady-state ... Time 1.00 [s] max point iter = 5 in node 6, max errors: W = 0.0000000 [kg/s], P = 0.07355 [Pa], H = 0.08337 [J/kg] Time 2.00 [s] max point iter = 4 in node 5, max errors: W = 0.0000023 [kg/s], P = 0.09317 [Pa], H = 0.04083 [J/kg] Time 3.00 [s] max point iter = 4 in node 57, max errors: W = 0.0000001 [kg/s], P = 0.09325 [Pa], H = 0.03948 [J/kg] Time 4.00 [s] max point iter = 3 in node 5, max errors: W = 0.0000000 [kg/s], P = 0.01212 [Pa], H = 0.00181 [J/kg] Time 5.00 [s] max point iter = 3 in node 8, max errors: W = 0.0000000 [kg/s], P = 0.08790 [Pa], H = 0.02859 [J/kg] Time 6.00 [s] max point iter = 3 in node 15, max errors: W = 0.0000000 [kg/s], P = 0.08555 [Pa], H = 0.08811 [J/kg] Time 7.00 [s] max point iter = 3 in node 30, max errors: W = 0.0000000 [kg/s], P = 0.08933 [Pa], H = 0.07321 [J/kg] Time 8.00 [s] max point iter = 3 in node 60, max errors: W = 0.0000000 [kg/s], P = 0.09924 [Pa], H = 0.09113 [J/kg] Time 9.00 [s] max point iter = 2 in node 42, max errors: W = 0.0000000 [kg/s], P = 0.02878 [Pa], H = 0.09626 [J/kg] STEADY-STATE CONVERGED max errors: W = 0.0000001 [kg/s], P = 0.07850 [Pa], H = 0.83469 [J/kg] Elapsed time: 2.61 sec --------------------------------------------- Mixture solver run completed --------------------------------------------- Output directory: /san/reactorsF/job1/2018p4151/mfval/Git/openstream/tutorials/outputs/TUTORIAL1-TUTORIAL2B-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.
We plot the axial distributions of vapor volumetric (void) and mass (quality) ratios.
% Plot the axial distributions of vapor ratio
mixSolverB.plotz('display','VR');
While the quality is the same as in the previous simulation (thermodynamic equilibrium assumption), the void is slightly lower, characteristic of the imposed velocity slip between the vapor and liquid phases.
This is verified by plotting the axial distributions of phase velocity.
% Plot the axial distributions of phase velocity
mixSolverB.plotz('display','U');
The calculated vapor velocity is observed to be twice higher than the liquid velocity, as specified by the SLIP void model.

Drift flux void model

In this section, we will walk through how to set up and run the mixture solver using a drift flux void model.
The only difference with the previous case is that input VOID is set to BESTION in the model set. This model is one of several drift flux void model available in the literature, which has the merit to be simple and accurate for many applications.
The path to the physical models, numerical options, test section geometry and boundary condition files are provided and read by InputSet. The models ID corresponding to the Bestion drift flux model (TUTORIAL2C) is indicated.
% Turn off warnings during input setup
warning off
% Create the input set using paths to model, options, geometry, and boundary condition files
inputSetC = InputSet( ...
modelFilePath = [osp '/inputs/models.inp'] , modelID = 'TUTORIAL2C', ...
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
The mixture solver object, mixSolverC, is created from inputSetC.
% Create the mixture solver object
mixSolverC = MixtureSolver(inputSetC);
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
mixSolverC.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.08279 [Pa], H = 0.09085 [J/kg] Time 2.00 [s] max point iter = 4 in node 5, max errors: W = 0.0000001 [kg/s], P = 0.07251 [Pa], H = 0.01908 [J/kg] Time 3.00 [s] max point iter = 4 in node 46, max errors: W = 0.0000002 [kg/s], P = 0.09745 [Pa], H = 0.05819 [J/kg] Time 4.00 [s] max point iter = 3 in node 5, max errors: W = 0.0000000 [kg/s], P = 0.00945 [Pa], H = 0.00283 [J/kg] Time 5.00 [s] max point iter = 3 in node 7, max errors: W = 0.0000000 [kg/s], P = 0.07799 [Pa], H = 0.02859 [J/kg] Time 6.00 [s] max point iter = 3 in node 15, max errors: W = 0.0000000 [kg/s], P = 0.08429 [Pa], H = 0.09584 [J/kg] Time 7.00 [s] max point iter = 3 in node 40, max errors: W = 0.0000000 [kg/s], P = 0.09510 [Pa], H = 0.09358 [J/kg] Time 8.00 [s] max point iter = 2 in node 25, max errors: W = 0.0000000 [kg/s], P = 0.06086 [Pa], H = 0.09582 [J/kg] Time 9.00 [s] max point iter = 2 in node 64, max errors: W = 0.0000000 [kg/s], P = 0.02491 [Pa], H = 0.09815 [J/kg] STEADY-STATE CONVERGED max errors: W = 0.0000000 [kg/s], P = 0.03292 [Pa], H = 0.25962 [J/kg] Elapsed time: 2.67 sec --------------------------------------------- Mixture solver run completed --------------------------------------------- Output directory: /san/reactorsF/job1/2018p4151/mfval/Git/openstream/tutorials/outputs/TUTORIAL1-TUTORIAL2C-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.
We plot the axial distributions of vapor volumetric (void) and mass (quality) ratios.
% Plot the axial distributions of vapor ratio
mixSolverC.plotz('display','VR');
While the quality is the same as in the previous simulation (thermodynamic equilibrium assumption), the void is slightly different due to different and varying velocity differences between the phase along the channel imposed by the drift flux model.
This is verified by plotting the axial distributions of phase velocity.
% Plot the axial distributions of phase velocity
mixSolverC.plotz('display','U');
The calculated vapor velocity is observed to be higher than the liquid velocity (with a non-constant slip ratio), as specified by the BESTION drift flux void model.

Thermal non-equilibrium models

In this section, two thermal non-equilibrium models are demonstrated where the phases have different temperatures. This is achieved by calculating the vapor mass quality based either on an empirical model or a physical model representative of interfacial mass and heat transfer, instead of the thermal equilibrium assumption.

Empirical model

In this section, we will walk through how to set up and run the mixture solver using an empirical thermal non-equilibrium model to simulate subcooled boiling.
The only difference with the previous case is that input THERMALNONEQ is set to SAHAZUBER in the models file.
The path to the physical models, numerical options, test section geometry and boundary condition files are provided and read by InputSet. The models ID corresponding to the Saha-Zuber subcooled boiling model (TUTORIAL2D) is indicated.
% Turn off warnings during input setup
warning off
% Create the input set using paths to model, options, geometry, and boundary condition files
inputSetD = InputSet( ...
modelFilePath = [osp '/inputs/models.inp'] , modelID = 'TUTORIAL2D', ...
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
The mixture solver object, mixSolverD, is created from inputSetD.
% Create the mixture solver object
mixSolverD = MixtureSolver(inputSetD);
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
mixSolverD.solve();
--------------------------------------------- Mixture solver run initiated --------------------------------------------- Solve steady-state ... Time 1.00 [s] max point iter = 5 in node 4, max errors: W = 0.0000023 [kg/s], P = 0.07521 [Pa], H = 0.07250 [J/kg] Time 2.00 [s] max point iter = 4 in node 4, max errors: W = 0.0000007 [kg/s], P = 0.08997 [Pa], H = 0.07848 [J/kg] Time 3.00 [s] max point iter = 4 in node 46, max errors: W = 0.0000000 [kg/s], P = 0.09609 [Pa], H = 0.00691 [J/kg] Time 4.00 [s] max point iter = 3 in node 4, max errors: W = 0.0000000 [kg/s], P = 0.03266 [Pa], H = 0.00180 [J/kg] Time 5.00 [s] max point iter = 3 in node 8, max errors: W = 0.0000000 [kg/s], P = 0.06573 [Pa], H = 0.03066 [J/kg] Time 6.00 [s] max point iter = 3 in node 19, max errors: W = 0.0000000 [kg/s], P = 0.09802 [Pa], H = 0.05823 [J/kg] Time 7.00 [s] max point iter = 3 in node 48, max errors: W = 0.0000000 [kg/s], P = 0.09564 [Pa], H = 0.08215 [J/kg] Time 8.00 [s] max point iter = 2 in node 31, max errors: W = 0.0000000 [kg/s], P = 0.04116 [Pa], H = 0.09545 [J/kg] Time 9.00 [s] max point iter = 2 in node 84, max errors: W = 0.0000000 [kg/s], P = 0.03762 [Pa], H = 0.09856 [J/kg] STEADY-STATE CONVERGED max errors: W = 0.0000000 [kg/s], P = 0.04045 [Pa], H = 0.16051 [J/kg] Elapsed time: 3.26 sec --------------------------------------------- Mixture solver run completed --------------------------------------------- Output directory: /san/reactorsF/job1/2018p4151/mfval/Git/openstream/tutorials/outputs/TUTORIAL1-TUTORIAL2D-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.
We plot the axial distributions of vapor volumetric (void) and mass (quality) ratios.
% Plot the axial distribution of vapor ratios
mixSolverD.plotz('display','VR');
The vapor mass quality begins to increase near the inlet, even though the equilibrium thermodynamic quality remains negative. Consequently, the void fraction also increases in that region. This behavior is characteristic of subcooled boiling, as captured by the Saha-Zuber model.

Homogeneous Relaxation Model (HRM)

In this section, we will walk through how to set up and run the mixture solver using a physical and generic thermal non-equilibrium model that can simulate both subcooled boiling and post-boiling transition (i.e., superheated vapor).
The only difference with the previous case is that input THERMALNONEQ is set to HRM in the models file. This Homogeneous Relaxation Model (HRM) is based on two additional conservation equations for the vapor phase (mass and energy) allowing to physically describe thermal non-equilibrium exchanges (evaporation and condensation) across the interface. In this tutorial, the time relaxation model is set to its default value.
The path to the physical models, numerical options, test section geometry and boundary condition files are provided and read by InputSet. The models ID corresponding to the homogeneous relaxation model (TUTORIAL2E) is indicated.
% Turn off warnings during input setup
warning off
% Create the input set using paths to model, options, geometry, and boundary condition files
inputSetE = InputSet( ...
modelFilePath = [osp '/inputs/models.inp'] , modelID = 'TUTORIAL2E', ...
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
The mixture solver object, mixSolverE, is created from inputSetE.
% Create the mixture solver object
mixSolverE = MixtureSolver(inputSetE);
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
mixSolverE.solve();
--------------------------------------------- Mixture solver run initiated --------------------------------------------- Solve steady-state ... Time 1.00 [s] max point iter = 12 in node 72, max errors: W = 0.0000000 [kg/s], Wv = 0.0009702 [kg/s], P = 0.09918 [Pa], H = 0.01855 [J/kg], Hv = 0.00000 [J/kg] Time 2.00 [s] max point iter = 9 in node 19, max errors: W = 0.0000000 [kg/s], Wv = 0.0003605 [kg/s], P = 0.09315 [Pa], H = 0.00549 [J/kg], Hv = 0.00000 [J/kg] Time 3.00 [s] max point iter = 8 in node 27, max errors: W = 0.0000000 [kg/s], Wv = 0.0000326 [kg/s], P = 0.09835 [Pa], H = 0.00095 [J/kg], Hv = 0.00000 [J/kg] Time 4.00 [s] max point iter = 6 in node 16, max errors: W = 0.0000000 [kg/s], Wv = 0.0000020 [kg/s], P = 0.05159 [Pa], H = 0.00065 [J/kg], Hv = 0.00000 [J/kg] Time 5.00 [s] max point iter = 5 in node 37, max errors: W = 0.0000000 [kg/s], Wv = 0.0000001 [kg/s], P = 0.09956 [Pa], H = 0.02586 [J/kg], Hv = 0.00000 [J/kg] Time 6.00 [s] max point iter = 4 in node 55, max errors: W = 0.0000000 [kg/s], Wv = 0.0000000 [kg/s], P = 0.09937 [Pa], H = 0.06437 [J/kg], Hv = 0.00000 [J/kg] Time 7.00 [s] max point iter = 3 in node 72, max errors: W = 0.0000000 [kg/s], Wv = 0.0000000 [kg/s], P = 0.09799 [Pa], H = 0.04341 [J/kg], Hv = 0.00000 [J/kg] Time 8.00 [s] max point iter = 2 in node 37, max errors: W = 0.0000000 [kg/s], Wv = 0.0000000 [kg/s], P = 0.09799 [Pa], H = 0.05806 [J/kg], Hv = 0.00000 [J/kg] Time 9.00 [s] max point iter = 2 in node 76, max errors: W = 0.0000000 [kg/s], Wv = 0.0000000 [kg/s], P = 0.09723 [Pa], H = 0.05152 [J/kg], Hv = 0.00000 [J/kg] STEADY-STATE CONVERGED max errors: W = 0.0000000 [kg/s], Wv = 0.0000000 [kg/s], P = 0.17283 [Pa], H = 0.09033 [J/kg], Hv = 0.00000 [J/kg] Elapsed time: 4.42 sec --------------------------------------------- Mixture solver run completed --------------------------------------------- Output directory: /san/reactorsF/job1/2018p4151/mfval/Git/openstream/tutorials/outputs/TUTORIAL1-TUTORIAL2E-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.
We plot the axial distributions of vapor volumetric (void) and mass (quality) ratios.
% Plot the axial distribution of vapor ratios
mixSolverE.plotz('display','VR');
The vapor mass quality begins to increase near the inlet, even though the equilibrium thermodynamic quality remains negative. Consequently, the void fraction also increases in that region. This behavior is characteristic of subcooled boiling, as captured by the HRM model.
A future tutorial will describe in details the features and capabilities of the Homogeneous Relaxation Model, including the simulation of thermal non-equilibrium in post-boiling transition.

Pressure drop models

In this section, the axial pressure drop models and simulation results are discussed. The axial pressure drop components (gravitational, wall shear and accelerations) are resolved from the mixture momentum conservation equation (which solves for the total pressure drop). The code saves the pressure drop components individually for easy retrieval.

Axial pressure drop components

Using an already solved mixture solver object (mixSolverD, based on Saha-Zuber subcooled boiling model), the cumulative pressure drop components can be plotted as follows.
% Plot the axial distribution of pressure drop contributors
mixSolverD.plotz('display','DP');
The references pressure is based on the inlet and, while the cumulative pressure drop is negative along the axial lengths, it is represented as a positive value on the figure.
For this case, both the local and temporal acceleration pressure drop components are 0. This is expected for a steady-state case in a channel with no obstruction. In addition, the spacial acceleration component remains constant along the adiabatic length, as expected.
As can be observed from the figure for two-phase convective boiling at high mass flux, the pressure drop is dominated by the wall friction, followed by the spatial acceleration (from boiling and related density decrease).

Wall friction

The wall fraction pressure drop is due to the shear of the flowing mixture on the channel boundaries. For one-dimensional analysis, correlations are typically used in term of a wall friction factor, from which the wall shear and wall friction pressure drop are derived. For two-phase flow applications, a wall friction two-phase multiplier is used to account for the increase of wall shear. Both wall friction factor and wall friction two-phase multiplier are generally derived from experimental data. A few typical models are implemented in OpenSTREAM. In this section, a Blasius type wall friction factor is set, along with the EPRI wall friction two-phase multiplier.
The path to the physical models, numerical options, test section geometry and boundary condition files are provided and read by InputSet. The models ID (TUTORIAL2F) corresponding to a Blasius type wall friction model () along with the EPRI two-phase flow multiplier is indicated.
% Turn off warnings during input setup
warning off
% Create the input set using paths to model, options, geometry, and boundary condition files
inputSetF = InputSet( ...
modelFilePath = [osp '/inputs/models.inp'] , modelID = 'TUTORIAL2F', ...
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
The mixture solver object, mixSolverF, is created from inputSetF.
% Create the mixture solver object
mixSolverF = MixtureSolver(inputSetF);
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
mixSolverF.solve();
--------------------------------------------- Mixture solver run initiated --------------------------------------------- Solve steady-state ... Time 1.00 [s] max point iter = 5 in node 4, max errors: W = 0.0000023 [kg/s], P = 0.07628 [Pa], H = 0.07250 [J/kg] Time 2.00 [s] max point iter = 4 in node 4, max errors: W = 0.0000007 [kg/s], P = 0.08903 [Pa], H = 0.07848 [J/kg] Time 3.00 [s] max point iter = 4 in node 46, max errors: W = 0.0000000 [kg/s], P = 0.09638 [Pa], H = 0.00691 [J/kg] Time 4.00 [s] max point iter = 3 in node 4, max errors: W = 0.0000000 [kg/s], P = 0.03593 [Pa], H = 0.00180 [J/kg] Time 5.00 [s] max point iter = 3 in node 8, max errors: W = 0.0000000 [kg/s], P = 0.06882 [Pa], H = 0.03066 [J/kg] Time 6.00 [s] max point iter = 3 in node 18, max errors: W = 0.0000000 [kg/s], P = 0.08498 [Pa], H = 0.05823 [J/kg] Time 7.00 [s] max point iter = 3 in node 48, max errors: W = 0.0000000 [kg/s], P = 0.09760 [Pa], H = 0.08215 [J/kg] Time 8.00 [s] max point iter = 2 in node 31, max errors: W = 0.0000000 [kg/s], P = 0.04243 [Pa], H = 0.09545 [J/kg] Time 9.00 [s] max point iter = 2 in node 84, max errors: W = 0.0000000 [kg/s], P = 0.03756 [Pa], H = 0.09856 [J/kg] STEADY-STATE CONVERGED max errors: W = 0.0000000 [kg/s], P = 0.03918 [Pa], H = 0.16051 [J/kg] Elapsed time: 2.44 sec --------------------------------------------- Mixture solver run completed --------------------------------------------- Output directory: /san/reactorsF/job1/2018p4151/mfval/Git/openstream/tutorials/outputs/TUTORIAL1-TUTORIAL2F-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.
We plot the axial cumulative distributions of pressure drop components.
% Plot the axial distribution of pressure drop contributors
mixSolverF.plotz('display','DP');
The change in wall friction model (from default) leads to an increase in wall friction, and hence total pressure drop.

Singular pressure losses

In an obstructed channel (e.g., due to a structural component, a flow mixer, a valve, etc.), the pressure drop can increase significantly across the obstruction. This is represented in the code by a singular pressure drop scaled with a pressure drop coefficient. For two-phase flow applications, a local two-phase multiplier is used to account for the increase of pressure drop. Both pressure drop coefficient and local pressure drop two-phase multiplier are generally derived from experimental data. A few typical models are implemented in OpenSTREAM. In this section, constant local pressure drop coefficients are set, along with the homogeneous local pressure drop two-phase multiplier.
The path to the physical models, numerical options, test section geometry and boundary condition files are provided and read by InputSet. The models ID (TUTORIAL2G) corresponding to channel with 5 obstructions (at 1, 2, 3, 4 and 5 m) of loss coefficient 0.5 is indicated. The homogeneous model for the two-phase local pressure drop multiplier is selected.
% Turn off warnings during input setup
warning off
% Create the input set using paths to model, options, geometry, and boundary condition files
inputSetG = InputSet( ...
modelFilePath = [osp '/inputs/models.inp'] , modelID = 'TUTORIAL2G', ...
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
The mixture solver object, mixSolverG, is created from inputSetG.
% Create the mixture solver object
mixSolverG = MixtureSolver(inputSetG);
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
mixSolverG.solve();
--------------------------------------------- Mixture solver run initiated --------------------------------------------- Solve steady-state ... Time 1.00 [s] max point iter = 5 in node 4, max errors: W = 0.0000023 [kg/s], P = 0.09284 [Pa], H = 0.07250 [J/kg] Time 2.00 [s] max point iter = 4 in node 4, max errors: W = 0.0000007 [kg/s], P = 0.08997 [Pa], H = 0.07848 [J/kg] Time 3.00 [s] max point iter = 4 in node 46, max errors: W = 0.0000000 [kg/s], P = 0.09609 [Pa], H = 0.00691 [J/kg] Time 4.00 [s] max point iter = 3 in node 4, max errors: W = 0.0000000 [kg/s], P = 0.03266 [Pa], H = 0.00180 [J/kg] Time 5.00 [s] max point iter = 3 in node 8, max errors: W = 0.0000000 [kg/s], P = 0.06573 [Pa], H = 0.03066 [J/kg] Time 6.00 [s] max point iter = 3 in node 19, max errors: W = 0.0000000 [kg/s], P = 0.09802 [Pa], H = 0.05823 [J/kg] Time 7.00 [s] max point iter = 3 in node 48, max errors: W = 0.0000000 [kg/s], P = 0.09564 [Pa], H = 0.08215 [J/kg] Time 8.00 [s] max point iter = 2 in node 31, max errors: W = 0.0000000 [kg/s], P = 0.04371 [Pa], H = 0.09545 [J/kg] Time 9.00 [s] max point iter = 2 in node 84, max errors: W = 0.0000000 [kg/s], P = 0.03895 [Pa], H = 0.09856 [J/kg] STEADY-STATE CONVERGED max errors: W = 0.0000000 [kg/s], P = 0.04241 [Pa], H = 0.16051 [J/kg] Elapsed time: 4.46 sec --------------------------------------------- Mixture solver run completed --------------------------------------------- Output directory: /san/reactorsF/job1/2018p4151/mfval/Git/openstream/tutorials/outputs/TUTORIAL1-TUTORIAL2G-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.
We plot the axial cumulative distributions of pressure drop components.
% Plot the axial distribution of pressure drop contributors
mixSolverG.plotz('display','DP');
The addition of local perturbations (with associated pressure loss coefficients) leads to clear local increases in pressure drop, and hence total pressure drop.

Transient

In this section a transient case is simulated, representative of a power decrease (by half of the nominal power over 1 sec). This example demonstrates the dynamic capabilities of OpenSTREAM for a simple case.
The path to the physical models, numerical options, test section geometry and boundary condition files are provided and read by InputSet.
The transient boundary conditions are specified in the tutorial2_transient.inp file. We can examine the file by using the command below.
type('../inputs/tutorial2_transient.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 TIME ! Time [s] > 1.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] > 43750 WMESH ! Relative power node size distribution [m] > 3.5 2.0 WPOWER ! Relative power distribution(s) [-] > 1.0 0.0 END TIME ! Time [s] > 5.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] > 43750 WMESH ! Relative power node size distribution [m] > 3.5 2.0 WPOWER ! Relative power distribution(s) [-] > 1.0 0.0 END
The format of the file is similar to setting a steady-state case, except that several set of inputs are provided at various (arbitrary) time, defining the transient. Any input parameter in the boundary conditons file can vary with time, that is:
In this tutorial, only the power is varied.
Regarding the model file, the models ID previously used to set the Bestion drift flux void model and the Saha-Zuber subcooled boiling model (TUTORIAL2D) is indicated.
% Turn off warnings during input setup
warning off
% Create the input set using paths to model, options, geometry, and boundary condition files
inputSetH = InputSet( ...
modelFilePath = [osp '/inputs/models.inp'] , modelID = 'TUTORIAL2D', ...
optionsFilePath = [osp '/inputs/options.inp'], optionsID = 'DEFAULT' , ...
geometryFilePath = [osp '/inputs/geom.inp'] , geometryID = 'TUTORIAL1', ...
bcFilePath = [osp '/inputs/tutorial2_transient.inp'], ...
sessionParentDir = fullfile(pwd,'outputs') , ...
overwriteSessionFiles = true , ...
LOGMODE = 'BOTH');
% Turn off warnings during input setup
warning on
The mixture solver object, mixSolverH, is created from inputSetH.
% Create the mixture solver object
mixSolverH = MixtureSolver(inputSetH);
To solve the flow fields, simply call the solve method (no input arguments are needed). Initially, the solver uses a pseudo-transient approach with large time steps to drive the solution toward steady-state. After steady-state convergence, the solver proceeds to simulate the transient scenario specified in the boundary conditions file.
% Solve the flow fields
mixSolverH.solve();
--------------------------------------------- Mixture solver run initiated --------------------------------------------- Solve steady-state ... Time 1.00 [s] max point iter = 5 in node 4, max errors: W = 0.0000023 [kg/s], P = 0.07521 [Pa], H = 0.07250 [J/kg] Time 2.00 [s] max point iter = 4 in node 4, max errors: W = 0.0000007 [kg/s], P = 0.08997 [Pa], H = 0.07848 [J/kg] Time 3.00 [s] max point iter = 4 in node 46, max errors: W = 0.0000000 [kg/s], P = 0.09609 [Pa], H = 0.00691 [J/kg] Time 4.00 [s] max point iter = 3 in node 4, max errors: W = 0.0000000 [kg/s], P = 0.03266 [Pa], H = 0.00180 [J/kg] Time 5.00 [s] max point iter = 3 in node 8, max errors: W = 0.0000000 [kg/s], P = 0.06573 [Pa], H = 0.03066 [J/kg] Time 6.00 [s] max point iter = 3 in node 19, max errors: W = 0.0000000 [kg/s], P = 0.09802 [Pa], H = 0.05823 [J/kg] Time 7.00 [s] max point iter = 3 in node 48, max errors: W = 0.0000000 [kg/s], P = 0.09564 [Pa], H = 0.08215 [J/kg] Time 8.00 [s] max point iter = 2 in node 31, max errors: W = 0.0000000 [kg/s], P = 0.04116 [Pa], H = 0.09545 [J/kg] Time 9.00 [s] max point iter = 2 in node 84, max errors: W = 0.0000000 [kg/s], P = 0.03762 [Pa], H = 0.09856 [J/kg] STEADY-STATE CONVERGED max errors: W = 0.0000000 [kg/s], P = 0.04045 [Pa], H = 0.16051 [J/kg] Elapsed time: 2.66 sec Solve transient ... Time 0.10 [s] max point iter = 6 in node 2, max errors: W = 0.0000004 [kg/s], P = 0.03952 [Pa], H = 0.09406 [J/kg] Time 0.20 [s] max point iter = 7 in node 4, max errors: W = 0.0000004 [kg/s], P = 0.03380 [Pa], H = 0.09482 [J/kg] Time 0.30 [s] max point iter = 7 in node 4, max errors: W = 0.0000004 [kg/s], P = 0.03390 [Pa], H = 0.09844 [J/kg] Time 0.40 [s] max point iter = 7 in node 4, max errors: W = 0.0000005 [kg/s], P = 0.03223 [Pa], H = 0.09766 [J/kg] Time 0.50 [s] max point iter = 7 in node 4, max errors: W = 0.0000005 [kg/s], P = 0.03169 [Pa], H = 0.09776 [J/kg] Time 0.60 [s] max point iter = 7 in node 4, max errors: W = 0.0000005 [kg/s], P = 0.03033 [Pa], H = 0.09139 [J/kg] Time 0.70 [s] max point iter = 7 in node 4, max errors: W = 0.0000005 [kg/s], P = 0.03194 [Pa], H = 0.09989 [J/kg] Time 0.80 [s] max point iter = 7 in node 5, max errors: W = 0.0000004 [kg/s], P = 0.03189 [Pa], H = 0.09636 [J/kg] Time 0.90 [s] max point iter = 7 in node 11, max errors: W = 0.0000006 [kg/s], P = 0.03645 [Pa], H = 0.09646 [J/kg] Time 1.00 [s] max point iter = 7 in node 9, max errors: W = 0.0000006 [kg/s], P = 0.03812 [Pa], H = 0.09655 [J/kg] Time 1.10 [s] max point iter = 6 in node 4, max errors: W = 0.0000009 [kg/s], P = 0.02457 [Pa], H = 0.09936 [J/kg] Time 1.20 [s] max point iter = 5 in node 7, max errors: W = 0.0000005 [kg/s], P = 0.07811 [Pa], H = 0.09973 [J/kg] Time 1.30 [s] max point iter = 5 in node 7, max errors: W = 0.0000005 [kg/s], P = 0.06471 [Pa], H = 0.09975 [J/kg] Time 1.40 [s] max point iter = 5 in node 10, max errors: W = 0.0000007 [kg/s], P = 0.08810 [Pa], H = 0.09787 [J/kg] Time 1.50 [s] max point iter = 4 in node 7, max errors: W = 0.0000005 [kg/s], P = 0.09983 [Pa], H = 0.09716 [J/kg] Time 1.60 [s] max point iter = 4 in node 9, max errors: W = 0.0000006 [kg/s], P = 0.07780 [Pa], H = 0.09848 [J/kg] Time 1.70 [s] max point iter = 4 in node 13, max errors: W = 0.0000007 [kg/s], P = 0.06492 [Pa], H = 0.09668 [J/kg] Time 1.80 [s] max point iter = 3 in node 6, max errors: W = 0.0000004 [kg/s], P = 0.04988 [Pa], H = 0.05803 [J/kg] Time 1.90 [s] max point iter = 3 in node 7, max errors: W = 0.0000002 [kg/s], P = 0.04125 [Pa], H = 0.06142 [J/kg] Time 2.00 [s] max point iter = 3 in node 9, max errors: W = 0.0000001 [kg/s], P = 0.09642 [Pa], H = 0.09237 [J/kg] Time 2.10 [s] max point iter = 3 in node 10, max errors: W = 0.0000001 [kg/s], P = 0.07636 [Pa], H = 0.02840 [J/kg] Time 2.20 [s] max point iter = 3 in node 12, max errors: W = 0.0000000 [kg/s], P = 0.09926 [Pa], H = 0.03546 [J/kg] Time 2.30 [s] max point iter = 3 in node 13, max errors: W = 0.0000000 [kg/s], P = 0.06584 [Pa], H = 0.03756 [J/kg] Time 2.40 [s] max point iter = 3 in node 16, max errors: W = 0.0000000 [kg/s], P = 0.09386 [Pa], H = 0.08582 [J/kg] Time 2.50 [s] max point iter = 3 in node 18, max errors: W = 0.0000000 [kg/s], P = 0.07649 [Pa], H = 0.06638 [J/kg] Time 2.60 [s] max point iter = 3 in node 22, max errors: W = 0.0000000 [kg/s], P = 0.09669 [Pa], H = 0.09094 [J/kg] Time 2.70 [s] max point iter = 3 in node 26, max errors: W = 0.0000000 [kg/s], P = 0.09887 [Pa], H = 0.05911 [J/kg] Time 2.80 [s] max point iter = 3 in node 30, max errors: W = 0.0000000 [kg/s], P = 0.08802 [Pa], H = 0.09394 [J/kg] Time 2.90 [s] max point iter = 3 in node 36, max errors: W = 0.0000000 [kg/s], P = 0.09382 [Pa], H = 0.07830 [J/kg] Time 3.00 [s] max point iter = 3 in node 43, max errors: W = 0.0000000 [kg/s], P = 0.09620 [Pa], H = 0.08169 [J/kg] Time 3.10 [s] max point iter = 3 in node 51, max errors: W = 0.0000000 [kg/s], P = 0.09432 [Pa], H = 0.09420 [J/kg] Time 3.20 [s] max point iter = 3 in node 61, max errors: W = 0.0000000 [kg/s], P = 0.09523 [Pa], H = 0.09018 [J/kg] Time 3.30 [s] max point iter = 3 in node 75, max errors: W = 0.0000000 [kg/s], P = 0.09908 [Pa], H = 0.09021 [J/kg] Time 3.40 [s] max point iter = 3 in node 89, max errors: W = 0.0000000 [kg/s], P = 0.09629 [Pa], H = 0.09036 [J/kg] Time 3.50 [s] max point iter = 2 in node 48, max errors: W = 0.0000000 [kg/s], P = 0.08836 [Pa], H = 0.09932 [J/kg] Time 3.60 [s] max point iter = 2 in node 56, max errors: W = 0.0000000 [kg/s], P = 0.04038 [Pa], H = 0.09216 [J/kg] Time 3.70 [s] max point iter = 2 in node 67, max errors: W = 0.0000000 [kg/s], P = 0.01805 [Pa], H = 0.09434 [J/kg] Time 3.80 [s] max point iter = 2 in node 81, max errors: W = 0.0000000 [kg/s], P = 0.01345 [Pa], H = 0.09721 [J/kg] Time 3.90 [s] max point iter = 2 in node 95, max errors: W = 0.0000000 [kg/s], P = 0.01697 [Pa], H = 0.09941 [J/kg] Time 4.00 [s] max point iter = 1 in node 2, max errors: W = 0.0000000 [kg/s], P = 0.01868 [Pa], H = 0.06475 [J/kg] Time 4.10 [s] max point iter = 1 in node 2, max errors: W = 0.0000000 [kg/s], P = 0.01230 [Pa], H = 0.02658 [J/kg] Time 4.20 [s] max point iter = 1 in node 2, max errors: W = 0.0000000 [kg/s], P = 0.00496 [Pa], H = 0.01068 [J/kg] Time 4.30 [s] max point iter = 1 in node 2, max errors: W = 0.0000000 [kg/s], P = 0.00196 [Pa], H = 0.00421 [J/kg] Time 4.40 [s] max point iter = 1 in node 2, max errors: W = 0.0000000 [kg/s], P = 0.00076 [Pa], H = 0.00162 [J/kg] Time 4.50 [s] max point iter = 1 in node 2, max errors: W = 0.0000000 [kg/s], P = 0.00029 [Pa], H = 0.00061 [J/kg] Time 4.60 [s] max point iter = 1 in node 2, max errors: W = 0.0000000 [kg/s], P = 0.00011 [Pa], H = 0.00023 [J/kg] Time 4.70 [s] max point iter = 1 in node 2, max errors: W = 0.0000000 [kg/s], P = 0.00004 [Pa], H = 0.00008 [J/kg] Time 4.80 [s] max point iter = 1 in node 2, max errors: W = 0.0000000 [kg/s], P = 0.00001 [Pa], H = 0.00003 [J/kg] Time 4.90 [s] max point iter = 1 in node 2, max errors: W = 0.0000000 [kg/s], P = 0.00001 [Pa], H = 0.00001 [J/kg] Time 5.00 [s] max point iter = 1 in node 2, max errors: W = 0.0000000 [kg/s], P = 0.00000 [Pa], H = 0.00000 [J/kg] Elapsed time: 14.42 sec --------------------------------------------- Mixture solver run completed --------------------------------------------- Output directory: /san/reactorsF/job1/2018p4151/mfval/Git/openstream/tutorials/outputs/TUTORIAL1-TUTORIAL2D-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). During the pseudo-transient simulation, 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.
We plot the axial distributions of heat flux (a boundary conditions) and of vapor volumetric (void) and mass (quality) ratios.
% Plot the axial distribution of heat fluc and vapor ratios
mixSolverH.plotz('display',{'HFLUX','VR'},'arrangement','horizontal','resize',1);
The initial solution is the same as the steady-state solution of case 2D, as expected. The plotz method of the mixtureSolver has the capability of showing the axial distribution of the selected parameters for each time step and play the entire transient (once or in a loop) by using the command button on the bottom left of the figure window. A video of the transient can be generated by using the animation drop dow menu
.