Steady-State Griffin-Pronghorn-SAM Coupling
Model link: Steady-State Coupling Model
Griffin, Pronghorn, and SAM are coupled with a domain overlapping approach. The steady-state Griffin-Pronghorn coupled model is described at Griffin-Pronghorn model. The coupling scheme between the Griffin-Pronghorn and SAM models is shown in Figure 1.
Figure 1: Coupling scheme between the different physics solvers
The core power is computed by the Griffin-Pronghorn model. This power is then applied to the core in SAM. Then, the computed cold primary heat exchanger temperature in SAM is applied as the cold source temperature in the Pronghorn model. The main blocks of the coupling are presented here below.
MultiApp coupling Griffin-Pronghorn to SAM
Pronghorn is the main application, whereas SAM is the sub-application. The MultiApp transferring data between Pronghorn and SAM is shown here below. This is a TransientMultiApp with matched time steps between Pronghorn and SAM.
[MultiApps]
[sam_balance_of_plant]
type = TransientMultiApp
app_type = 'SamApp'
input_files = msfr_system_1d.i
max_procs_per_app = 1
execute_on = 'timestep_end'
keep_solution_during_restore = False
# the balance of plant needs to be run separately first for 100 steps
# here we are hardcoding a restart file path, to a file we checked in, but you should
# re-create the checkpoint file and adapt the restart_file_base path for your own runs
cli_args = 'Problem/restart_file_base=restart/msfr_system_1d_checkpoint_cp/0100;Problem/force_restart=true'
[]
[]The two key transfers are shown in the following block. They consist of MultiAppPostprocessorVectorTransfer in which values between post-processors are communicated between both codes.
[Transfers]
# Primary and secondary loops
[send_flow_BCs]
type = MultiAppPostprocessorVectorTransfer
to_multi_app = sam_balance_of_plant
from_postprocessors = 'total_power'
to_postprocessors = 'core_power'
[]
[receive_flow_BCs]
type = MultiAppPostprocessorVectorTransfer
from_multi_app = sam_balance_of_plant
from_postprocessors = 'HX_cold_temp'
to_postprocessors = 'heat_exchanger_T_ambient'
[]
[]Post-processors transferred
Two sets of postprocessors are implemented to transfer data between Pronghorn and SAM. To transfer power, an ElementIntegralVariablePostprocessor computes the total power in Pronghorn (transferred previously from Griffin) and sends it to a Receiver postprocessor in SAM.
[Postprocessors]
[total_power]
type = ElementIntegralVariablePostprocessor
variable = power_density
block = 'fuel pump hx'
[]
[][Postprocessors]
[core_power]
type = Receiver
default = ${core_power_val}
execute_on = 'INITIAL TIMESTEP_BEGIN TIMESTEP_END TRANSFER'
[]
[]To transfer the heat exchanger cold temperature, an ElementAverageValue computes the heat exchanger cold temperature in SAM and sends it to a Receiver postprocessor in Pronghorn.
[Postprocessors]
[HX_cold_temp]
type = ElementAverageValue
variable = temperature
block = 'pipe3'
[]
[][Postprocessors]
[heat_exchanger_T_ambient]
type = Receiver
default = ${T_HX}
execute_on = 'INITIAL TIMESTEP_END TRANSFER'
[]
[](msr/msfr/plant/steady/run_ns.i)
################################################################################
## Molten Salt Fast Reactor - Euratom EVOL + Rosatom MARS Design ##
## Pronghorn Sub-Application input file ##
## Relaxation towards Steady state 3D thermal hydraulics model ##
################################################################################
# If using or referring to this model, please cite as explained in
# https://mooseframework.inl.gov/virtual_test_bed/citing.html
# The flow in this simulation should be initialized with a previous flow
# solution (isothermal, heated or multiphysics) OR using a viscosity rampdown
# An isothermal solution can be generated using 'run_ns_initial.i', and is
# saved in 'restart/run_ns_initial_restart.e'
# A heated solution, with a flat power distribution, can be generated with
# this script 'run_ns.i', and is saved in 'restart/run_ns_restart.e'
# A coupled neutronics-coarse TH solution can be generated with
# 'run_neutronics.i', saved in 'restart/run_neutronics_ns_restart.e'
# Material properties
rho = 4284 # density [kg / m^3] (@1000K)
cp = 1594 # specific heat capacity [J / kg / K]
drho_dT = 0.882 # derivative of density w.r.t. temperature [kg / m^3 / K]
mu = 0.0166 # viscosity [Pa s]
k = 1.7 # thermal conductivity [W / m / K]
# https://www.researchgate.net/publication/337161399_Development_of_a_control-\
# oriented_power_plant_simulator_for_the_molten_salt_fast_reactor/fulltext/5dc9\
# 5c9da6fdcc57503eec39/Development-of-a-control-oriented-power-plant-simulator-\
# for-the-molten-salt-fast-reactor.pdf
von_karman_const = 0.41
# Turbulent properties
Pr_t = 0.9 # turbulent Prandtl number
Sc_t = 1 # turbulent Schmidt number
# Derived material properties
alpha = '${fparse drho_dT / rho}' # thermal expansion coefficient
# Operating parameters
T_HX = 873.15 # heat exchanger temperature [K]
# Mass flow rate tuning, for heat exchanger pressure and temperature drop
friction = 4e3 # [kg / m^4]
pump_force = -20000. # [N / m^3]
# Delayed neutron precursor parameters. Lambda values are decay constants in
# [1 / s]. Beta values are production fractions.
lambda1_m = -0.0133104
lambda2_m = -0.0305427
lambda3_m = -0.115179
lambda4_m = -0.301152
lambda5_m = -0.879376
lambda6_m = -2.91303
beta1 = 8.42817e-05
beta2 = 0.000684616
beta3 = 0.000479796
beta4 = 0.00103883
beta5 = 0.000549185
beta6 = 0.000184087
[GlobalParams]
rhie_chow_user_object = 'ins_rhie_chow_interpolator'
[]
################################################################################
# GEOMETRY
################################################################################
[Mesh]
coord_type = 'RZ'
rz_coord_axis = Y
[restart]
type = FileMeshGenerator
use_for_exodus_restart = true
# Depending on the file chosen, the initialization of variables should be
# adjusted. The following variables can be initalized:
# - vel_x, vel_y, p from isothermal simulation
# file = '../../steady/restart/run_ns_initial_restart.e'
# Below are initialization points created from this input file
# The variable IC should be set from_file_var for temperature and precursors
# - vel_x, vel_y, p, T_fluid, c_i from cosine heated simulation
# file = '../../steady/restart/run_ns_restart.e'
# - adding SAM-coupling
file = 'restart/run_ns_restart.e'
[]
[hx_top]
type = ParsedGenerateSideset
combinatorial_geometry = 'y > 0'
included_subdomains = '3'
included_neighbors = '1'
fixed_normal = true
normal = '0 1 0'
new_sideset_name = 'hx_top'
input = 'restart'
[]
[hx_bot]
type = ParsedGenerateSideset
combinatorial_geometry = 'y <-0.6'
included_subdomains = '3'
included_neighbors = '1'
fixed_normal = true
normal = '0 -1 0'
new_sideset_name = 'hx_bot'
input = 'hx_top'
[]
[]
[Problem]
# Using restart for velocities & pressure initialization
allow_initial_conditions_with_restart = true
[]
################################################################################
# EQUATIONS: VARIABLES, KERNELS, BOUNDARY CONDITIONS
################################################################################
[Modules]
[NavierStokesFV]
# General parameters
compressibility = 'incompressible'
add_energy_equation = true
add_scalar_equation = true
boussinesq_approximation = true
block = 'fuel pump hx'
# Variables, defined below for the Exodus restart
velocity_variable = 'vel_x vel_y'
pressure_variable = 'pressure'
fluid_temperature_variable = 'T_fluid'
# Material properties
density = ${rho}
dynamic_viscosity = ${mu}
thermal_conductivity = ${k}
specific_heat = 'cp'
thermal_expansion = ${alpha}
# Boussinesq parameters
gravity = '0 -9.81 0'
ref_temperature = ${T_HX}
# Boundary conditions
wall_boundaries = 'shield_wall reflector_wall fluid_symmetry'
momentum_wall_types = 'wallfunction wallfunction symmetry'
energy_wall_types = 'heatflux heatflux heatflux'
energy_wall_function = '0 0 0'
# Pressure pin for incompressible flow
pin_pressure = true
pinned_pressure_type = average
pinned_pressure_value = 1e5
# Turbulence parameters
turbulence_handling = 'mixing-length'
turbulent_prandtl = ${Pr_t}
von_karman_const = ${von_karman_const}
mixing_length_delta = 0.1
mixing_length_walls = 'shield_wall reflector_wall'
mixing_length_aux_execute_on = 'initial'
# Numerical scheme
momentum_advection_interpolation = 'upwind'
mass_advection_interpolation = 'upwind'
energy_advection_interpolation = 'upwind'
passive_scalar_advection_interpolation = 'upwind'
# Heat source
external_heat_source = power_density
# Precursor advection, diffusion and source term
passive_scalar_names = 'c1 c2 c3 c4 c5 c6'
passive_scalar_schmidt_number = '${Sc_t} ${Sc_t} ${Sc_t} ${Sc_t} ${Sc_t} ${Sc_t}'
passive_scalar_coupled_source = 'fission_source c1; fission_source c2; fission_source c3;
fission_source c4; fission_source c5; fission_source c6'
passive_scalar_coupled_source_coeff = '${beta1} ${lambda1_m}; ${beta2} ${lambda2_m};
${beta3} ${lambda3_m}; ${beta4} ${lambda4_m};
${beta5} ${lambda5_m}; ${beta6} ${lambda6_m}'
# Heat exchanger
standard_friction_formulation = false
friction_blocks = 'hx'
friction_types = 'FORCHHEIMER'
friction_coeffs = 'friction_coeff_vector'
ambient_convection_blocks = 'hx'
ambient_convection_alpha = '${fparse 600 * 20e3}' # HX specifications
ambient_temperature = 'hx_cold_temp'
[]
[]
[Variables]
[vel_x]
type = INSFVVelocityVariable
block = 'fuel pump hx'
initial_from_file_var = vel_x
[]
[vel_y]
type = INSFVVelocityVariable
block = 'fuel pump hx'
initial_from_file_var = vel_y
[]
[pressure]
type = INSFVPressureVariable
block = 'fuel pump hx'
initial_from_file_var = pressure
[]
[T_fluid]
type = INSFVEnergyVariable
block = 'fuel pump hx'
initial_condition = ${T_HX}
# initial_from_file_var = T_fluid
[]
[c1]
type = MooseVariableFVReal
block = 'fuel pump hx'
# initial_from_file_var = c1
[]
[c2]
type = MooseVariableFVReal
block = 'fuel pump hx'
# initial_from_file_var = c2
[]
[c3]
type = MooseVariableFVReal
block = 'fuel pump hx'
# initial_from_file_var = c3
[]
[c4]
type = MooseVariableFVReal
block = 'fuel pump hx'
# initial_from_file_var = c4
[]
[c5]
type = MooseVariableFVReal
block = 'fuel pump hx'
# initial_from_file_var = c5
[]
[c6]
type = MooseVariableFVReal
block = 'fuel pump hx'
# initial_from_file_var = c6
[]
[]
[FVICs]
[c1]
type = FVFunctionIC
variable = c1
function = 'cosine_guess'
scaling_factor = 0.02
[]
[c2]
type = FVFunctionIC
variable = c2
function = 'cosine_guess'
scaling_factor = 0.1
[]
[c3]
type = FVFunctionIC
variable = c3
function = 'cosine_guess'
scaling_factor = 0.03
[]
[c4]
type = FVFunctionIC
variable = c4
function = 'cosine_guess'
scaling_factor = 0.04
[]
[c5]
type = FVFunctionIC
variable = c5
function = 'cosine_guess'
scaling_factor = 0.01
[]
[c6]
type = FVFunctionIC
variable = c6
function = 'cosine_guess'
scaling_factor = 0.001
[]
# Power density is re-initalized by a transfer from neutronics
[power]
type = FVFunctionIC
variable = power_density
function = 'cosine_guess'
scaling_factor = '${fparse 3e9/2.81543}'
[]
# Fission source is re-initalized by a transfer from neutronics
[fission_source]
type = FVFunctionIC
variable = fission_source
function = 'cosine_guess'
scaling_factor = '${fparse 6.303329e+01/2.81543}'
[]
[]
[AuxVariables]
[power_density]
type = MooseVariableFVReal
block = 'fuel pump hx'
[]
[fission_source]
type = MooseVariableFVReal
block = 'fuel pump hx'
[]
[]
[Functions]
# Guess to have a 3D power distribution
[cosine_guess]
type = ParsedFunction
expression = 'max(0, cos(x*pi/2/1.2))*max(0, cos(y*pi/2/1.1))'
[]
[hx_cold_temp]
type = ParsedFunction
expression = 'Tcold'
symbol_names = 'Tcold'
symbol_values = '${T_HX}'
[]
[]
[FVKernels]
[pump]
type = INSFVBodyForce
variable = vel_y
functor = ${pump_force}
block = 'pump'
momentum_component = 'y'
[]
[]
################################################################################
# MATERIALS
################################################################################
[FunctorMaterials]
# Most of these constants could be specified directly to the action
[mu]
type = ADGenericFunctorMaterial
prop_names = 'mu'
prop_values = '${mu}'
block = 'fuel pump hx'
[]
# [not_used]
# type = ADGenericFunctorMaterial
# prop_names = 'not_used'
# prop_values = 0
# block = 'shield reflector'
# []
[cp]
type = ADGenericFunctorMaterial
prop_names = 'cp'
prop_values = '${cp}'
block = 'fuel pump hx'
[]
[friction_coeff_mat]
type = ADGenericVectorFunctorMaterial
prop_names = 'friction_coeff_vector'
prop_values = '${friction} ${friction} ${friction}'
[]
[]
################################################################################
# EXECUTION / SOLVE
################################################################################
[Functions]
[dts]
type = PiecewiseConstant
x = '0 100'
y = '0.75 2.5'
[]
[]
[Executioner]
type = Transient
# Time stepping parameters
start_time = 0.0
end_time = 20
# end_time will depend on the restart file chosen
# though steady state detection can also be used
# from _initial/no heating : 150 - 200s enough
# from _ns/_ns_coupled/heated: 10s enough
[TimeStepper]
# This time stepper makes the time step grow exponentially
# It can only be used with proper initialization
type = IterationAdaptiveDT
dt = 1 # chosen to obtain convergence with first coupled iteration
growth_factor = 2
[]
# [TimeStepper]
# type = FunctionDT
# function = dts
# []
steady_state_detection = true
steady_state_tolerance = 1e-8
steady_state_start_time = 10
# Time integration scheme
scheme = 'implicit-euler'
# Solver parameters
solve_type = 'NEWTON'
petsc_options_iname = '-pc_type -pc_factor_shift_type -ksp_gmres_restart'
petsc_options_value = 'lu NONZERO 50'
line_search = 'none'
nl_rel_tol = 1e-9
nl_abs_tol = 2e-8
nl_max_its = 15
l_max_its = 50
automatic_scaling = true
# Fixed point iterations are not necessary if the only goal is to obtain
# a steady state solution, no matter the numerical path
# fixed_point_max_its = 10
# fixed_point_abs_tol = 1e-5
# accept_on_max_fixed_point_iteration = true
[]
################################################################################
# MULTIAPPS FOR POWER TRANSFER
################################################################################
[MultiApps]
[sam_balance_of_plant]
type = TransientMultiApp
app_type = 'SamApp'
input_files = msfr_system_1d.i
max_procs_per_app = 1
execute_on = 'timestep_end'
keep_solution_during_restore = False
# the balance of plant needs to be run separately first for 100 steps
# here we are hardcoding a restart file path, to a file we checked in, but you should
# re-create the checkpoint file and adapt the restart_file_base path for your own runs
cli_args = 'Problem/restart_file_base=restart/msfr_system_1d_checkpoint_cp/0100;Problem/force_restart=true'
[]
[]
[Transfers]
# Primary and secondary loops
[send_flow_BCs]
type = MultiAppPostprocessorVectorTransfer
to_multi_app = sam_balance_of_plant
from_postprocessors = 'total_power'
to_postprocessors = 'core_power'
[]
[receive_flow_BCs]
type = MultiAppPostprocessorVectorTransfer
from_multi_app = sam_balance_of_plant
from_postprocessors = 'HX_cold_temp'
to_postprocessors = 'heat_exchanger_T_ambient'
[]
[]
################################################################################
# SIMULATION OUTPUTS
################################################################################
[Outputs]
csv = true
hide = 'flow_hx_bot flow_hx_top min_flow_T max_flow_T'
[restart]
type = Exodus
overwrite = true
[]
# Reduce base output
print_linear_converged_reason = false
print_linear_residuals = false
print_nonlinear_converged_reason = false
[]
[Postprocessors]
[max_v]
type = ElementExtremeValue
variable = vel_x
value_type = max
block = 'fuel pump hx'
[]
# TODO: weakly compressible, switch to mass flow rate
[flow_hx_bot]
type = VolumetricFlowRate
boundary = 'hx_bot'
vel_x = vel_x
vel_y = vel_y
advected_quantity = 1
[]
[flow_hx_top]
type = VolumetricFlowRate
boundary = 'hx_top'
vel_x = vel_x
vel_y = vel_y
advected_quantity = 1
[]
[max_flow_T]
type = VolumetricFlowRate
boundary = 'hx_top'
vel_x = vel_x
vel_y = vel_y
advected_quantity = 'T_fluid'
[]
[min_flow_T]
type = VolumetricFlowRate
boundary = 'hx_bot'
vel_x = vel_x
vel_y = vel_y
advected_quantity = 'T_fluid'
[]
[dT]
type = ParsedPostprocessor
expression = '-max_flow_T / flow_hx_bot + min_flow_T / flow_hx_top'
pp_names = 'max_flow_T min_flow_T flow_hx_bot flow_hx_top'
[]
[total_power]
type = ElementIntegralVariablePostprocessor
variable = power_density
block = 'fuel pump hx'
[]
[total_fission_source]
type = ElementIntegralVariablePostprocessor
variable = fission_source
block = 'fuel pump hx'
[]
[heat_exchanger_T_ambient]
type = Receiver
default = ${T_HX}
execute_on = 'INITIAL TIMESTEP_END TRANSFER'
[]
[]
(msr/msfr/plant/steady/run_ns.i)
################################################################################
## Molten Salt Fast Reactor - Euratom EVOL + Rosatom MARS Design ##
## Pronghorn Sub-Application input file ##
## Relaxation towards Steady state 3D thermal hydraulics model ##
################################################################################
# If using or referring to this model, please cite as explained in
# https://mooseframework.inl.gov/virtual_test_bed/citing.html
# The flow in this simulation should be initialized with a previous flow
# solution (isothermal, heated or multiphysics) OR using a viscosity rampdown
# An isothermal solution can be generated using 'run_ns_initial.i', and is
# saved in 'restart/run_ns_initial_restart.e'
# A heated solution, with a flat power distribution, can be generated with
# this script 'run_ns.i', and is saved in 'restart/run_ns_restart.e'
# A coupled neutronics-coarse TH solution can be generated with
# 'run_neutronics.i', saved in 'restart/run_neutronics_ns_restart.e'
# Material properties
rho = 4284 # density [kg / m^3] (@1000K)
cp = 1594 # specific heat capacity [J / kg / K]
drho_dT = 0.882 # derivative of density w.r.t. temperature [kg / m^3 / K]
mu = 0.0166 # viscosity [Pa s]
k = 1.7 # thermal conductivity [W / m / K]
# https://www.researchgate.net/publication/337161399_Development_of_a_control-\
# oriented_power_plant_simulator_for_the_molten_salt_fast_reactor/fulltext/5dc9\
# 5c9da6fdcc57503eec39/Development-of-a-control-oriented-power-plant-simulator-\
# for-the-molten-salt-fast-reactor.pdf
von_karman_const = 0.41
# Turbulent properties
Pr_t = 0.9 # turbulent Prandtl number
Sc_t = 1 # turbulent Schmidt number
# Derived material properties
alpha = '${fparse drho_dT / rho}' # thermal expansion coefficient
# Operating parameters
T_HX = 873.15 # heat exchanger temperature [K]
# Mass flow rate tuning, for heat exchanger pressure and temperature drop
friction = 4e3 # [kg / m^4]
pump_force = -20000. # [N / m^3]
# Delayed neutron precursor parameters. Lambda values are decay constants in
# [1 / s]. Beta values are production fractions.
lambda1_m = -0.0133104
lambda2_m = -0.0305427
lambda3_m = -0.115179
lambda4_m = -0.301152
lambda5_m = -0.879376
lambda6_m = -2.91303
beta1 = 8.42817e-05
beta2 = 0.000684616
beta3 = 0.000479796
beta4 = 0.00103883
beta5 = 0.000549185
beta6 = 0.000184087
[GlobalParams]
rhie_chow_user_object = 'ins_rhie_chow_interpolator'
[]
################################################################################
# GEOMETRY
################################################################################
[Mesh]
coord_type = 'RZ'
rz_coord_axis = Y
[restart]
type = FileMeshGenerator
use_for_exodus_restart = true
# Depending on the file chosen, the initialization of variables should be
# adjusted. The following variables can be initalized:
# - vel_x, vel_y, p from isothermal simulation
# file = '../../steady/restart/run_ns_initial_restart.e'
# Below are initialization points created from this input file
# The variable IC should be set from_file_var for temperature and precursors
# - vel_x, vel_y, p, T_fluid, c_i from cosine heated simulation
# file = '../../steady/restart/run_ns_restart.e'
# - adding SAM-coupling
file = 'restart/run_ns_restart.e'
[]
[hx_top]
type = ParsedGenerateSideset
combinatorial_geometry = 'y > 0'
included_subdomains = '3'
included_neighbors = '1'
fixed_normal = true
normal = '0 1 0'
new_sideset_name = 'hx_top'
input = 'restart'
[]
[hx_bot]
type = ParsedGenerateSideset
combinatorial_geometry = 'y <-0.6'
included_subdomains = '3'
included_neighbors = '1'
fixed_normal = true
normal = '0 -1 0'
new_sideset_name = 'hx_bot'
input = 'hx_top'
[]
[]
[Problem]
# Using restart for velocities & pressure initialization
allow_initial_conditions_with_restart = true
[]
################################################################################
# EQUATIONS: VARIABLES, KERNELS, BOUNDARY CONDITIONS
################################################################################
[Modules]
[NavierStokesFV]
# General parameters
compressibility = 'incompressible'
add_energy_equation = true
add_scalar_equation = true
boussinesq_approximation = true
block = 'fuel pump hx'
# Variables, defined below for the Exodus restart
velocity_variable = 'vel_x vel_y'
pressure_variable = 'pressure'
fluid_temperature_variable = 'T_fluid'
# Material properties
density = ${rho}
dynamic_viscosity = ${mu}
thermal_conductivity = ${k}
specific_heat = 'cp'
thermal_expansion = ${alpha}
# Boussinesq parameters
gravity = '0 -9.81 0'
ref_temperature = ${T_HX}
# Boundary conditions
wall_boundaries = 'shield_wall reflector_wall fluid_symmetry'
momentum_wall_types = 'wallfunction wallfunction symmetry'
energy_wall_types = 'heatflux heatflux heatflux'
energy_wall_function = '0 0 0'
# Pressure pin for incompressible flow
pin_pressure = true
pinned_pressure_type = average
pinned_pressure_value = 1e5
# Turbulence parameters
turbulence_handling = 'mixing-length'
turbulent_prandtl = ${Pr_t}
von_karman_const = ${von_karman_const}
mixing_length_delta = 0.1
mixing_length_walls = 'shield_wall reflector_wall'
mixing_length_aux_execute_on = 'initial'
# Numerical scheme
momentum_advection_interpolation = 'upwind'
mass_advection_interpolation = 'upwind'
energy_advection_interpolation = 'upwind'
passive_scalar_advection_interpolation = 'upwind'
# Heat source
external_heat_source = power_density
# Precursor advection, diffusion and source term
passive_scalar_names = 'c1 c2 c3 c4 c5 c6'
passive_scalar_schmidt_number = '${Sc_t} ${Sc_t} ${Sc_t} ${Sc_t} ${Sc_t} ${Sc_t}'
passive_scalar_coupled_source = 'fission_source c1; fission_source c2; fission_source c3;
fission_source c4; fission_source c5; fission_source c6'
passive_scalar_coupled_source_coeff = '${beta1} ${lambda1_m}; ${beta2} ${lambda2_m};
${beta3} ${lambda3_m}; ${beta4} ${lambda4_m};
${beta5} ${lambda5_m}; ${beta6} ${lambda6_m}'
# Heat exchanger
standard_friction_formulation = false
friction_blocks = 'hx'
friction_types = 'FORCHHEIMER'
friction_coeffs = 'friction_coeff_vector'
ambient_convection_blocks = 'hx'
ambient_convection_alpha = '${fparse 600 * 20e3}' # HX specifications
ambient_temperature = 'hx_cold_temp'
[]
[]
[Variables]
[vel_x]
type = INSFVVelocityVariable
block = 'fuel pump hx'
initial_from_file_var = vel_x
[]
[vel_y]
type = INSFVVelocityVariable
block = 'fuel pump hx'
initial_from_file_var = vel_y
[]
[pressure]
type = INSFVPressureVariable
block = 'fuel pump hx'
initial_from_file_var = pressure
[]
[T_fluid]
type = INSFVEnergyVariable
block = 'fuel pump hx'
initial_condition = ${T_HX}
# initial_from_file_var = T_fluid
[]
[c1]
type = MooseVariableFVReal
block = 'fuel pump hx'
# initial_from_file_var = c1
[]
[c2]
type = MooseVariableFVReal
block = 'fuel pump hx'
# initial_from_file_var = c2
[]
[c3]
type = MooseVariableFVReal
block = 'fuel pump hx'
# initial_from_file_var = c3
[]
[c4]
type = MooseVariableFVReal
block = 'fuel pump hx'
# initial_from_file_var = c4
[]
[c5]
type = MooseVariableFVReal
block = 'fuel pump hx'
# initial_from_file_var = c5
[]
[c6]
type = MooseVariableFVReal
block = 'fuel pump hx'
# initial_from_file_var = c6
[]
[]
[FVICs]
[c1]
type = FVFunctionIC
variable = c1
function = 'cosine_guess'
scaling_factor = 0.02
[]
[c2]
type = FVFunctionIC
variable = c2
function = 'cosine_guess'
scaling_factor = 0.1
[]
[c3]
type = FVFunctionIC
variable = c3
function = 'cosine_guess'
scaling_factor = 0.03
[]
[c4]
type = FVFunctionIC
variable = c4
function = 'cosine_guess'
scaling_factor = 0.04
[]
[c5]
type = FVFunctionIC
variable = c5
function = 'cosine_guess'
scaling_factor = 0.01
[]
[c6]
type = FVFunctionIC
variable = c6
function = 'cosine_guess'
scaling_factor = 0.001
[]
# Power density is re-initalized by a transfer from neutronics
[power]
type = FVFunctionIC
variable = power_density
function = 'cosine_guess'
scaling_factor = '${fparse 3e9/2.81543}'
[]
# Fission source is re-initalized by a transfer from neutronics
[fission_source]
type = FVFunctionIC
variable = fission_source
function = 'cosine_guess'
scaling_factor = '${fparse 6.303329e+01/2.81543}'
[]
[]
[AuxVariables]
[power_density]
type = MooseVariableFVReal
block = 'fuel pump hx'
[]
[fission_source]
type = MooseVariableFVReal
block = 'fuel pump hx'
[]
[]
[Functions]
# Guess to have a 3D power distribution
[cosine_guess]
type = ParsedFunction
expression = 'max(0, cos(x*pi/2/1.2))*max(0, cos(y*pi/2/1.1))'
[]
[hx_cold_temp]
type = ParsedFunction
expression = 'Tcold'
symbol_names = 'Tcold'
symbol_values = '${T_HX}'
[]
[]
[FVKernels]
[pump]
type = INSFVBodyForce
variable = vel_y
functor = ${pump_force}
block = 'pump'
momentum_component = 'y'
[]
[]
################################################################################
# MATERIALS
################################################################################
[FunctorMaterials]
# Most of these constants could be specified directly to the action
[mu]
type = ADGenericFunctorMaterial
prop_names = 'mu'
prop_values = '${mu}'
block = 'fuel pump hx'
[]
# [not_used]
# type = ADGenericFunctorMaterial
# prop_names = 'not_used'
# prop_values = 0
# block = 'shield reflector'
# []
[cp]
type = ADGenericFunctorMaterial
prop_names = 'cp'
prop_values = '${cp}'
block = 'fuel pump hx'
[]
[friction_coeff_mat]
type = ADGenericVectorFunctorMaterial
prop_names = 'friction_coeff_vector'
prop_values = '${friction} ${friction} ${friction}'
[]
[]
################################################################################
# EXECUTION / SOLVE
################################################################################
[Functions]
[dts]
type = PiecewiseConstant
x = '0 100'
y = '0.75 2.5'
[]
[]
[Executioner]
type = Transient
# Time stepping parameters
start_time = 0.0
end_time = 20
# end_time will depend on the restart file chosen
# though steady state detection can also be used
# from _initial/no heating : 150 - 200s enough
# from _ns/_ns_coupled/heated: 10s enough
[TimeStepper]
# This time stepper makes the time step grow exponentially
# It can only be used with proper initialization
type = IterationAdaptiveDT
dt = 1 # chosen to obtain convergence with first coupled iteration
growth_factor = 2
[]
# [TimeStepper]
# type = FunctionDT
# function = dts
# []
steady_state_detection = true
steady_state_tolerance = 1e-8
steady_state_start_time = 10
# Time integration scheme
scheme = 'implicit-euler'
# Solver parameters
solve_type = 'NEWTON'
petsc_options_iname = '-pc_type -pc_factor_shift_type -ksp_gmres_restart'
petsc_options_value = 'lu NONZERO 50'
line_search = 'none'
nl_rel_tol = 1e-9
nl_abs_tol = 2e-8
nl_max_its = 15
l_max_its = 50
automatic_scaling = true
# Fixed point iterations are not necessary if the only goal is to obtain
# a steady state solution, no matter the numerical path
# fixed_point_max_its = 10
# fixed_point_abs_tol = 1e-5
# accept_on_max_fixed_point_iteration = true
[]
################################################################################
# MULTIAPPS FOR POWER TRANSFER
################################################################################
[MultiApps]
[sam_balance_of_plant]
type = TransientMultiApp
app_type = 'SamApp'
input_files = msfr_system_1d.i
max_procs_per_app = 1
execute_on = 'timestep_end'
keep_solution_during_restore = False
# the balance of plant needs to be run separately first for 100 steps
# here we are hardcoding a restart file path, to a file we checked in, but you should
# re-create the checkpoint file and adapt the restart_file_base path for your own runs
cli_args = 'Problem/restart_file_base=restart/msfr_system_1d_checkpoint_cp/0100;Problem/force_restart=true'
[]
[]
[Transfers]
# Primary and secondary loops
[send_flow_BCs]
type = MultiAppPostprocessorVectorTransfer
to_multi_app = sam_balance_of_plant
from_postprocessors = 'total_power'
to_postprocessors = 'core_power'
[]
[receive_flow_BCs]
type = MultiAppPostprocessorVectorTransfer
from_multi_app = sam_balance_of_plant
from_postprocessors = 'HX_cold_temp'
to_postprocessors = 'heat_exchanger_T_ambient'
[]
[]
################################################################################
# SIMULATION OUTPUTS
################################################################################
[Outputs]
csv = true
hide = 'flow_hx_bot flow_hx_top min_flow_T max_flow_T'
[restart]
type = Exodus
overwrite = true
[]
# Reduce base output
print_linear_converged_reason = false
print_linear_residuals = false
print_nonlinear_converged_reason = false
[]
[Postprocessors]
[max_v]
type = ElementExtremeValue
variable = vel_x
value_type = max
block = 'fuel pump hx'
[]
# TODO: weakly compressible, switch to mass flow rate
[flow_hx_bot]
type = VolumetricFlowRate
boundary = 'hx_bot'
vel_x = vel_x
vel_y = vel_y
advected_quantity = 1
[]
[flow_hx_top]
type = VolumetricFlowRate
boundary = 'hx_top'
vel_x = vel_x
vel_y = vel_y
advected_quantity = 1
[]
[max_flow_T]
type = VolumetricFlowRate
boundary = 'hx_top'
vel_x = vel_x
vel_y = vel_y
advected_quantity = 'T_fluid'
[]
[min_flow_T]
type = VolumetricFlowRate
boundary = 'hx_bot'
vel_x = vel_x
vel_y = vel_y
advected_quantity = 'T_fluid'
[]
[dT]
type = ParsedPostprocessor
expression = '-max_flow_T / flow_hx_bot + min_flow_T / flow_hx_top'
pp_names = 'max_flow_T min_flow_T flow_hx_bot flow_hx_top'
[]
[total_power]
type = ElementIntegralVariablePostprocessor
variable = power_density
block = 'fuel pump hx'
[]
[total_fission_source]
type = ElementIntegralVariablePostprocessor
variable = fission_source
block = 'fuel pump hx'
[]
[heat_exchanger_T_ambient]
type = Receiver
default = ${T_HX}
execute_on = 'INITIAL TIMESTEP_END TRANSFER'
[]
[]
(msr/msfr/plant/steady/run_ns.i)
################################################################################
## Molten Salt Fast Reactor - Euratom EVOL + Rosatom MARS Design ##
## Pronghorn Sub-Application input file ##
## Relaxation towards Steady state 3D thermal hydraulics model ##
################################################################################
# If using or referring to this model, please cite as explained in
# https://mooseframework.inl.gov/virtual_test_bed/citing.html
# The flow in this simulation should be initialized with a previous flow
# solution (isothermal, heated or multiphysics) OR using a viscosity rampdown
# An isothermal solution can be generated using 'run_ns_initial.i', and is
# saved in 'restart/run_ns_initial_restart.e'
# A heated solution, with a flat power distribution, can be generated with
# this script 'run_ns.i', and is saved in 'restart/run_ns_restart.e'
# A coupled neutronics-coarse TH solution can be generated with
# 'run_neutronics.i', saved in 'restart/run_neutronics_ns_restart.e'
# Material properties
rho = 4284 # density [kg / m^3] (@1000K)
cp = 1594 # specific heat capacity [J / kg / K]
drho_dT = 0.882 # derivative of density w.r.t. temperature [kg / m^3 / K]
mu = 0.0166 # viscosity [Pa s]
k = 1.7 # thermal conductivity [W / m / K]
# https://www.researchgate.net/publication/337161399_Development_of_a_control-\
# oriented_power_plant_simulator_for_the_molten_salt_fast_reactor/fulltext/5dc9\
# 5c9da6fdcc57503eec39/Development-of-a-control-oriented-power-plant-simulator-\
# for-the-molten-salt-fast-reactor.pdf
von_karman_const = 0.41
# Turbulent properties
Pr_t = 0.9 # turbulent Prandtl number
Sc_t = 1 # turbulent Schmidt number
# Derived material properties
alpha = '${fparse drho_dT / rho}' # thermal expansion coefficient
# Operating parameters
T_HX = 873.15 # heat exchanger temperature [K]
# Mass flow rate tuning, for heat exchanger pressure and temperature drop
friction = 4e3 # [kg / m^4]
pump_force = -20000. # [N / m^3]
# Delayed neutron precursor parameters. Lambda values are decay constants in
# [1 / s]. Beta values are production fractions.
lambda1_m = -0.0133104
lambda2_m = -0.0305427
lambda3_m = -0.115179
lambda4_m = -0.301152
lambda5_m = -0.879376
lambda6_m = -2.91303
beta1 = 8.42817e-05
beta2 = 0.000684616
beta3 = 0.000479796
beta4 = 0.00103883
beta5 = 0.000549185
beta6 = 0.000184087
[GlobalParams]
rhie_chow_user_object = 'ins_rhie_chow_interpolator'
[]
################################################################################
# GEOMETRY
################################################################################
[Mesh]
coord_type = 'RZ'
rz_coord_axis = Y
[restart]
type = FileMeshGenerator
use_for_exodus_restart = true
# Depending on the file chosen, the initialization of variables should be
# adjusted. The following variables can be initalized:
# - vel_x, vel_y, p from isothermal simulation
# file = '../../steady/restart/run_ns_initial_restart.e'
# Below are initialization points created from this input file
# The variable IC should be set from_file_var for temperature and precursors
# - vel_x, vel_y, p, T_fluid, c_i from cosine heated simulation
# file = '../../steady/restart/run_ns_restart.e'
# - adding SAM-coupling
file = 'restart/run_ns_restart.e'
[]
[hx_top]
type = ParsedGenerateSideset
combinatorial_geometry = 'y > 0'
included_subdomains = '3'
included_neighbors = '1'
fixed_normal = true
normal = '0 1 0'
new_sideset_name = 'hx_top'
input = 'restart'
[]
[hx_bot]
type = ParsedGenerateSideset
combinatorial_geometry = 'y <-0.6'
included_subdomains = '3'
included_neighbors = '1'
fixed_normal = true
normal = '0 -1 0'
new_sideset_name = 'hx_bot'
input = 'hx_top'
[]
[]
[Problem]
# Using restart for velocities & pressure initialization
allow_initial_conditions_with_restart = true
[]
################################################################################
# EQUATIONS: VARIABLES, KERNELS, BOUNDARY CONDITIONS
################################################################################
[Modules]
[NavierStokesFV]
# General parameters
compressibility = 'incompressible'
add_energy_equation = true
add_scalar_equation = true
boussinesq_approximation = true
block = 'fuel pump hx'
# Variables, defined below for the Exodus restart
velocity_variable = 'vel_x vel_y'
pressure_variable = 'pressure'
fluid_temperature_variable = 'T_fluid'
# Material properties
density = ${rho}
dynamic_viscosity = ${mu}
thermal_conductivity = ${k}
specific_heat = 'cp'
thermal_expansion = ${alpha}
# Boussinesq parameters
gravity = '0 -9.81 0'
ref_temperature = ${T_HX}
# Boundary conditions
wall_boundaries = 'shield_wall reflector_wall fluid_symmetry'
momentum_wall_types = 'wallfunction wallfunction symmetry'
energy_wall_types = 'heatflux heatflux heatflux'
energy_wall_function = '0 0 0'
# Pressure pin for incompressible flow
pin_pressure = true
pinned_pressure_type = average
pinned_pressure_value = 1e5
# Turbulence parameters
turbulence_handling = 'mixing-length'
turbulent_prandtl = ${Pr_t}
von_karman_const = ${von_karman_const}
mixing_length_delta = 0.1
mixing_length_walls = 'shield_wall reflector_wall'
mixing_length_aux_execute_on = 'initial'
# Numerical scheme
momentum_advection_interpolation = 'upwind'
mass_advection_interpolation = 'upwind'
energy_advection_interpolation = 'upwind'
passive_scalar_advection_interpolation = 'upwind'
# Heat source
external_heat_source = power_density
# Precursor advection, diffusion and source term
passive_scalar_names = 'c1 c2 c3 c4 c5 c6'
passive_scalar_schmidt_number = '${Sc_t} ${Sc_t} ${Sc_t} ${Sc_t} ${Sc_t} ${Sc_t}'
passive_scalar_coupled_source = 'fission_source c1; fission_source c2; fission_source c3;
fission_source c4; fission_source c5; fission_source c6'
passive_scalar_coupled_source_coeff = '${beta1} ${lambda1_m}; ${beta2} ${lambda2_m};
${beta3} ${lambda3_m}; ${beta4} ${lambda4_m};
${beta5} ${lambda5_m}; ${beta6} ${lambda6_m}'
# Heat exchanger
standard_friction_formulation = false
friction_blocks = 'hx'
friction_types = 'FORCHHEIMER'
friction_coeffs = 'friction_coeff_vector'
ambient_convection_blocks = 'hx'
ambient_convection_alpha = '${fparse 600 * 20e3}' # HX specifications
ambient_temperature = 'hx_cold_temp'
[]
[]
[Variables]
[vel_x]
type = INSFVVelocityVariable
block = 'fuel pump hx'
initial_from_file_var = vel_x
[]
[vel_y]
type = INSFVVelocityVariable
block = 'fuel pump hx'
initial_from_file_var = vel_y
[]
[pressure]
type = INSFVPressureVariable
block = 'fuel pump hx'
initial_from_file_var = pressure
[]
[T_fluid]
type = INSFVEnergyVariable
block = 'fuel pump hx'
initial_condition = ${T_HX}
# initial_from_file_var = T_fluid
[]
[c1]
type = MooseVariableFVReal
block = 'fuel pump hx'
# initial_from_file_var = c1
[]
[c2]
type = MooseVariableFVReal
block = 'fuel pump hx'
# initial_from_file_var = c2
[]
[c3]
type = MooseVariableFVReal
block = 'fuel pump hx'
# initial_from_file_var = c3
[]
[c4]
type = MooseVariableFVReal
block = 'fuel pump hx'
# initial_from_file_var = c4
[]
[c5]
type = MooseVariableFVReal
block = 'fuel pump hx'
# initial_from_file_var = c5
[]
[c6]
type = MooseVariableFVReal
block = 'fuel pump hx'
# initial_from_file_var = c6
[]
[]
[FVICs]
[c1]
type = FVFunctionIC
variable = c1
function = 'cosine_guess'
scaling_factor = 0.02
[]
[c2]
type = FVFunctionIC
variable = c2
function = 'cosine_guess'
scaling_factor = 0.1
[]
[c3]
type = FVFunctionIC
variable = c3
function = 'cosine_guess'
scaling_factor = 0.03
[]
[c4]
type = FVFunctionIC
variable = c4
function = 'cosine_guess'
scaling_factor = 0.04
[]
[c5]
type = FVFunctionIC
variable = c5
function = 'cosine_guess'
scaling_factor = 0.01
[]
[c6]
type = FVFunctionIC
variable = c6
function = 'cosine_guess'
scaling_factor = 0.001
[]
# Power density is re-initalized by a transfer from neutronics
[power]
type = FVFunctionIC
variable = power_density
function = 'cosine_guess'
scaling_factor = '${fparse 3e9/2.81543}'
[]
# Fission source is re-initalized by a transfer from neutronics
[fission_source]
type = FVFunctionIC
variable = fission_source
function = 'cosine_guess'
scaling_factor = '${fparse 6.303329e+01/2.81543}'
[]
[]
[AuxVariables]
[power_density]
type = MooseVariableFVReal
block = 'fuel pump hx'
[]
[fission_source]
type = MooseVariableFVReal
block = 'fuel pump hx'
[]
[]
[Functions]
# Guess to have a 3D power distribution
[cosine_guess]
type = ParsedFunction
expression = 'max(0, cos(x*pi/2/1.2))*max(0, cos(y*pi/2/1.1))'
[]
[hx_cold_temp]
type = ParsedFunction
expression = 'Tcold'
symbol_names = 'Tcold'
symbol_values = '${T_HX}'
[]
[]
[FVKernels]
[pump]
type = INSFVBodyForce
variable = vel_y
functor = ${pump_force}
block = 'pump'
momentum_component = 'y'
[]
[]
################################################################################
# MATERIALS
################################################################################
[FunctorMaterials]
# Most of these constants could be specified directly to the action
[mu]
type = ADGenericFunctorMaterial
prop_names = 'mu'
prop_values = '${mu}'
block = 'fuel pump hx'
[]
# [not_used]
# type = ADGenericFunctorMaterial
# prop_names = 'not_used'
# prop_values = 0
# block = 'shield reflector'
# []
[cp]
type = ADGenericFunctorMaterial
prop_names = 'cp'
prop_values = '${cp}'
block = 'fuel pump hx'
[]
[friction_coeff_mat]
type = ADGenericVectorFunctorMaterial
prop_names = 'friction_coeff_vector'
prop_values = '${friction} ${friction} ${friction}'
[]
[]
################################################################################
# EXECUTION / SOLVE
################################################################################
[Functions]
[dts]
type = PiecewiseConstant
x = '0 100'
y = '0.75 2.5'
[]
[]
[Executioner]
type = Transient
# Time stepping parameters
start_time = 0.0
end_time = 20
# end_time will depend on the restart file chosen
# though steady state detection can also be used
# from _initial/no heating : 150 - 200s enough
# from _ns/_ns_coupled/heated: 10s enough
[TimeStepper]
# This time stepper makes the time step grow exponentially
# It can only be used with proper initialization
type = IterationAdaptiveDT
dt = 1 # chosen to obtain convergence with first coupled iteration
growth_factor = 2
[]
# [TimeStepper]
# type = FunctionDT
# function = dts
# []
steady_state_detection = true
steady_state_tolerance = 1e-8
steady_state_start_time = 10
# Time integration scheme
scheme = 'implicit-euler'
# Solver parameters
solve_type = 'NEWTON'
petsc_options_iname = '-pc_type -pc_factor_shift_type -ksp_gmres_restart'
petsc_options_value = 'lu NONZERO 50'
line_search = 'none'
nl_rel_tol = 1e-9
nl_abs_tol = 2e-8
nl_max_its = 15
l_max_its = 50
automatic_scaling = true
# Fixed point iterations are not necessary if the only goal is to obtain
# a steady state solution, no matter the numerical path
# fixed_point_max_its = 10
# fixed_point_abs_tol = 1e-5
# accept_on_max_fixed_point_iteration = true
[]
################################################################################
# MULTIAPPS FOR POWER TRANSFER
################################################################################
[MultiApps]
[sam_balance_of_plant]
type = TransientMultiApp
app_type = 'SamApp'
input_files = msfr_system_1d.i
max_procs_per_app = 1
execute_on = 'timestep_end'
keep_solution_during_restore = False
# the balance of plant needs to be run separately first for 100 steps
# here we are hardcoding a restart file path, to a file we checked in, but you should
# re-create the checkpoint file and adapt the restart_file_base path for your own runs
cli_args = 'Problem/restart_file_base=restart/msfr_system_1d_checkpoint_cp/0100;Problem/force_restart=true'
[]
[]
[Transfers]
# Primary and secondary loops
[send_flow_BCs]
type = MultiAppPostprocessorVectorTransfer
to_multi_app = sam_balance_of_plant
from_postprocessors = 'total_power'
to_postprocessors = 'core_power'
[]
[receive_flow_BCs]
type = MultiAppPostprocessorVectorTransfer
from_multi_app = sam_balance_of_plant
from_postprocessors = 'HX_cold_temp'
to_postprocessors = 'heat_exchanger_T_ambient'
[]
[]
################################################################################
# SIMULATION OUTPUTS
################################################################################
[Outputs]
csv = true
hide = 'flow_hx_bot flow_hx_top min_flow_T max_flow_T'
[restart]
type = Exodus
overwrite = true
[]
# Reduce base output
print_linear_converged_reason = false
print_linear_residuals = false
print_nonlinear_converged_reason = false
[]
[Postprocessors]
[max_v]
type = ElementExtremeValue
variable = vel_x
value_type = max
block = 'fuel pump hx'
[]
# TODO: weakly compressible, switch to mass flow rate
[flow_hx_bot]
type = VolumetricFlowRate
boundary = 'hx_bot'
vel_x = vel_x
vel_y = vel_y
advected_quantity = 1
[]
[flow_hx_top]
type = VolumetricFlowRate
boundary = 'hx_top'
vel_x = vel_x
vel_y = vel_y
advected_quantity = 1
[]
[max_flow_T]
type = VolumetricFlowRate
boundary = 'hx_top'
vel_x = vel_x
vel_y = vel_y
advected_quantity = 'T_fluid'
[]
[min_flow_T]
type = VolumetricFlowRate
boundary = 'hx_bot'
vel_x = vel_x
vel_y = vel_y
advected_quantity = 'T_fluid'
[]
[dT]
type = ParsedPostprocessor
expression = '-max_flow_T / flow_hx_bot + min_flow_T / flow_hx_top'
pp_names = 'max_flow_T min_flow_T flow_hx_bot flow_hx_top'
[]
[total_power]
type = ElementIntegralVariablePostprocessor
variable = power_density
block = 'fuel pump hx'
[]
[total_fission_source]
type = ElementIntegralVariablePostprocessor
variable = fission_source
block = 'fuel pump hx'
[]
[heat_exchanger_T_ambient]
type = Receiver
default = ${T_HX}
execute_on = 'INITIAL TIMESTEP_END TRANSFER'
[]
[]
(msr/msfr/plant/steady/msfr_system_1d.i)
################################################################################
## Molten Salt Fast Reactor - Euratom EVOL + Rosatom MARS Design ##
## SAM Sub-Sub-Application input file ##
## Relaxation towards steady state balance of plant model ##
################################################################################
# If using or referring to this model, please cite as explained in
# https://mooseframework.inl.gov/virtual_test_bed/citing.html
core_power_val = 3e9
[Application]
type = SamApp
[]
[GlobalParams]
global_init_P = 1e5
global_init_V = 0.01
global_init_T = 898.15
Tsolid_sf = 1e-3
scaling_factor_var = '1 1e-3 1e-6'
[]
[Functions]
[fuel_salt_rho_func] # Linear fitting used by Rouch et al. 2014
type = PiecewiseLinear
x = '800 1200'
y = '4277.96 3925.16'
[]
[fuel_salt_cp_func]
type = PiecewiseLinear
x = '800 1200'
y = '1113.00 2225.00'
[]
[fuel_salt_k_func]
type = PiecewiseLinear
x = '800 1200'
y = '0.995200 1.028800'
[]
[fuel_salt_mu_func] # Nonlinear fitting used by Rouch et al. 2014
type = PiecewiseLinear
x = ' 800 820 840 860 880 900
920 940 960 980 1000 1020
1040 1060 1080 1100 1120 1140
1160 1180 1200'
y = '2.3887E-02 2.1258E-02 1.9020E-02 1.7102E-02 1.5449E-02 1.4015E-02
1.2767E-02 1.1673E-02 1.0711E-02 9.8609E-03 9.1066E-03 8.4347E-03
7.8340E-03 7.2951E-03 6.8099E-03 6.3719E-03 5.9751E-03 5.6147E-03
5.2865E-03 4.9868E-03 4.7124E-03'
[]
[fuel_salt_enthalpy_func] # approximated by cp*T
type = PiecewiseLinear
x = ' 800 820 840 860 880 900
920 940 960 980 1000 1020
1040 1060 1080 1100 1120 1140
1160 1180 1200'
y = '890400 913216 937144 962184 988336 1015600
1043976 1073464 1104064 1135776 1168600 1202536
1237584 1273744 1311016 1349400 1388896 1429504
1471224 1514056 1558000'
[]
[TimeStepperFunc]
type = PiecewiseLinear
#x = '-200 -180 -160 -140 -100 -60 -50 -20 -10. 0.0 100'
#y = ' 0.2 0.5 1.0 2.0 5.0 10.0 50. 10. 1.0 0.1 0.5'
x = '-200 -180 -160 -140 -100 -60 -50 -20 -10. 0.0 100'
y = ' 1.0 2.0 5.0 10.0 20.0 20.0 50. 10. 1.0 0.1 0.5'
[]
[core_specific_power]
type = ParsedFunction
expression = '${fparse core_power_val / 9.27}'
[]
[pump_head_fun]
type = PiecewiseConstant
xy_data = '0.0 156010.45
5.0 78005.225'
direction = 'left'
[]
[]
[EOS]
[fuel_salt_eos]
type = PTFunctionsEOS
rho = fuel_salt_rho_func
# beta = 2.1051E-04 # thermal expansion coefficient - approximated by -(drho/dT)/rho
cp = fuel_salt_cp_func
mu = fuel_salt_mu_func
k = fuel_salt_k_func
enthalpy = fuel_salt_enthalpy_func
[]
[hx_salt_eos]
type = SaltEquationOfState
salt_type = Flibe
[]
[Helium]
type = HeEquationOfState
[]
[]
[MaterialProperties]
[alloy-mat] # Based on Hastelloy N alloy
type = SolidMaterialProps
k = 23.6
Cp = 578
rho = 8.86e3
[]
[]
[Components]
[MSFR_core]
type = PBOneDFluidComponent
eos = fuel_salt_eos
orientation = '0 0 1'
position = '0 0 0'
A = 3.4636
Dh = 2.1
length = 2.65
n_elems = 25
heat_source = core_specific_power
[]
[Core_Out_Branch]
type = PBBranch
inputs = 'MSFR_core(out)'
outputs = 'pipe1(in) loop_pbc(in)'
eos = fuel_salt_eos
K = '4.25 4.25 500.'
Area = 2.24
[]
[pipe1] # Horizontal hot channel
type = PBOneDFluidComponent
eos = fuel_salt_eos
position = '0 0 2.65'
orientation = '0 1 0'
A = 2.24
Dh = 1.6888
length = 2.0
n_elems = 12
[]
[loop_pbc]
type = PBOneDFluidComponent
eos = fuel_salt_eos
position = '0 0 2.65'
orientation = '0 0 1'
A = 2.24
Dh = 1.6888
length = 0.025
n_elems = 1
[]
[TDV1]
type = PBTDV
input = 'loop_pbc(out)'
eos = fuel_salt_eos
p_bc = 1.0e5
[]
#
# ====== pump ======
#
[pump]
type = PBPump
Area = 2.24
K = '0.15 0.1'
eos = fuel_salt_eos
inputs = 'pipe1(out)'
outputs = 'pipe2(in)'
initial_P = 1.0e5
Head = 156010.45
[]
[pipe2] # Vertical hot channel from pump to HX
type = PBOneDFluidComponent
eos = fuel_salt_eos
position = '0 2.0 2.65'
orientation = '0 0 -1'
A = 2.24
Dh = 1.6888
length = 0.25
n_elems = 3
[]
#
# ====== Heat Exchanger ======
#
[IHX1]
type = PBHeatExchanger
eos = fuel_salt_eos
eos_secondary = hx_salt_eos
position = '0 2.0 2.4'
orientation = '0 0 -1'
A = 2.24
A_secondary = 3.60
Dh = 0.02
Dh_secondary = 0.0125
length = 2.4
n_elems = 24
SC_HTC = 4
HTC_geometry_type = Pipe
HTC_geometry_type_secondary = Pipe
HT_surface_area_density = 800.0
HT_surface_area_density_secondary = 497.78
initial_V_secondary = -2.67
Twall_init = 898.15
wall_thickness = 0.001
dim_wall = 1
material_wall = alloy-mat
n_wall_elems = 2
[]
[J_P2_IHX1]
type = PBBranch
inputs = 'pipe2(out)'
outputs = 'IHX1(primary_in) '
eos = fuel_salt_eos
K = '1. 1.'
Area = 2.24
[]
#
# ====== Heat Exchanger Secondary Side ======
#
[IHX1_P3]
type = PBBranch
inputs = 'IHX1(primary_out)'
outputs = 'pipe3(in) '
eos = fuel_salt_eos
K = '1. 1.'
Area = 2.24
[]
[pipe3] # Horizontal cold channel
type = PBOneDFluidComponent
eos = fuel_salt_eos
position = '0 2.0 0.'
orientation = '0 -1 0'
A = 2.24
Dh = 1.6888
length = 2.0
n_elems = 12
[]
[Core_In_Branch]
type = PBBranch
inputs = 'pipe3(out)'
outputs = 'MSFR_core(in)'
eos = fuel_salt_eos
K = '4.25 4.25'
Area = 2.24
[]
#
# ====== Intermediate circuit connected to HX1 ======
#
[IHX1_P4]
type = PBBranch
inputs = 'IHX1(secondary_out)'
outputs = 'pipe4(in), loop2_pbc(in)'
eos = hx_salt_eos
K = '4.25 4.25 500.'
Area = 3.6
[]
[loop2_pbc]
type = PBOneDFluidComponent
eos = hx_salt_eos
position = '0 2.0 2.4'
orientation = '1 0 0'
A = 3.6
Dh = 2.14
length = 0.025
n_elems = 1
[]
[loop2_TDV1]
type = PBTDV
input = 'loop2_pbc(out)'
eos = hx_salt_eos
p_bc = 1.0e5
[]
[pipe4] # Horizontal intermediate hot leg
type = PBOneDFluidComponent
eos = hx_salt_eos
position = '0 2.0 2.4'
orientation = '0 1 0'
A = 3.6
Dh = 2.14
length = 2.0
n_elems = 10
[]
[J_P4_P5]
type = PBBranch
inputs = 'pipe4(out)'
outputs = 'pipe5(in)'
eos = hx_salt_eos
K = '1. 1.'
Area = 3.6
[]
[pipe5] # Vertical intermediate hot leg
type = PBOneDFluidComponent
eos = hx_salt_eos
position = '0 4.0 2.4'
orientation = '0 0 -1'
A = 3.6
Dh = 2.14
length = 0.2
n_elems = 2
[]
[IHX2]
type = PBHeatExchanger
eos = hx_salt_eos
eos_secondary = Helium
position = '0 4.0 2.2'
orientation = '0 0 -1'
A = 2.4
A_secondary = 7.2
Dh = 0.04
Dh_secondary = 0.01
length = 3.2
n_elems = 24
SC_HTC = 4
HTC_geometry_type = Pipe
HTC_geometry_type_secondary = Pipe
HT_surface_area_density = 510.0
HT_surface_area_density_secondary = 170.0
initial_V_secondary = -50
Twall_init = 900
wall_thickness = 0.002
dim_wall = 1
material_wall = alloy-mat
n_wall_elems = 2
[]
[J_P5_IHX2]
type = PBBranch
inputs = 'pipe5(out)'
outputs = 'IHX2(primary_in)'
eos = hx_salt_eos
K = '1. 1.'
Area = 3.6
[]
[IHX2_S_In]
type = PBTDJ
input = 'IHX2(secondary_in)'
eos = Helium
v_bc = -66.65
T_bc = 673.15
[]
[IHX2_S_Out]
type = PBTDV
input = 'IHX2(secondary_out)'
eos = Helium
p_bc = 7.5e6
[]
[J_IHX2_P6]
type = PBBranch
inputs = 'IHX2(primary_out)'
outputs = 'pipe6(in)'
eos = hx_salt_eos
K = '1. 1.'
Area = 3.6
[]
[pipe6] # Intermediate cold leg
type = PBOneDFluidComponent
eos = hx_salt_eos
position = '0 4.0 -1.0'
orientation = '0 -1 0'
A = 3.6
Dh = 2.14
length = 1.0
n_elems = 5
[]
[J_P6_P7]
type = PBBranch
inputs = 'pipe6(out)'
outputs = 'pipe7(in)'
eos = hx_salt_eos
K = '1. 1.'
Area = 3.6
[]
[pipe7] # Intermediate cold leg
type = PBOneDFluidComponent
eos = hx_salt_eos
position = '0 3.0 -1.0'
orientation = '0 0 1'
A = 3.6
Dh = 2.14
length = 1.0
n_elems = 5
[]
[pump2]
type = PBPump
Area = 3.6
K = '0.15 0.1'
eos = hx_salt_eos
inputs = 'pipe7(out)'
outputs = 'pipe8(in)'
initial_P = 1.0e5
Head = 1.76e5
[]
[pipe8] # Intermediate cold leg
type = PBOneDFluidComponent
eos = hx_salt_eos
position = '0 3.0 0.0'
orientation = '0 -1 0'
A = 3.6
Dh = 2.14
length = 1.0
n_elems = 5
[]
[P8_IHX1]
type = PBBranch
inputs = 'pipe8(out)'
outputs = 'IHX1(secondary_in)'
eos = hx_salt_eos
K = '1. 1.'
Area = 3.6
[]
[]
[Postprocessors]
[Core_P_out]
type = ComponentBoundaryVariableValue
variable = pressure
input = MSFR_core(out)
[]
[Core_T_out]
type = ComponentBoundaryVariableValue
variable = temperature
input = MSFR_core(out)
[]
[Fuel_mass_flow] # Output mass flow rate at inlet of CH1
type = ComponentBoundaryFlow
input = MSFR_core(in)
[]
[Core_T_in]
type = ComponentBoundaryVariableValue
variable = temperature
input = MSFR_core(in)
[]
[Core_P_in]
type = ComponentBoundaryVariableValue
variable = pressure
input = MSFR_core(in)
[]
[HX_Tin_p]
type = ComponentBoundaryVariableValue
variable = temperature
input = IHX1(primary_in)
[]
[HX_Tout_p]
type = ComponentBoundaryVariableValue
variable = temperature
input = IHX1(primary_out)
[]
[HX_Tout_s]
type = ComponentBoundaryVariableValue
variable = temperature
input = IHX1(secondary_out)
[]
[HX_Uout_s]
type = ComponentBoundaryVariableValue
variable = velocity
input = IHX1(secondary_out)
[]
[FliBe_mass_flow]
type = ComponentBoundaryFlow
input = IHX1(secondary_in)
[]
[HX_Pin_s]
type = ComponentBoundaryVariableValue
variable = pressure
input = IHX1(secondary_in)
[]
[FliBe_mass_flow2]
type = ComponentBoundaryFlow
input = IHX2(primary_in)
[]
[Gas_mass_flow]
type = ComponentBoundaryFlow
input = IHX2(secondary_in)
[]
[HX2_Tin_p]
type = ComponentBoundaryVariableValue
variable = temperature
input = IHX2(primary_in)
[]
[HX2_Tout_p]
type = ComponentBoundaryVariableValue
variable = temperature
input = IHX2(primary_out)
[]
[HX2_Tout_s]
type = ComponentBoundaryVariableValue
variable = temperature
input = IHX2(secondary_out)
[]
[HX_cold_temp]
type = ElementAverageValue
variable = temperature
block = 'pipe3'
[]
[core_power]
type = Receiver
default = ${core_power_val}
execute_on = 'INITIAL TIMESTEP_BEGIN TIMESTEP_END TRANSFER'
[]
[]
[Preconditioning]
[SMP_PJFNK]
type = SMP
full = true
solve_type = 'PJFNK'
petsc_options_iname = '-pc_type -ksp_gmres_restart'
petsc_options_value = 'lu 101'
[]
[]
[Executioner]
type = Transient
[TimeStepper]
type = FunctionDT
function = TimeStepperFunc
min_dt = 1e-3
[]
#[TimeStepper]
# type = AB2PredictorCorrector
# dt = .01
# e_max = 10
# e_tol = 1
#[]
nl_rel_tol = 1e-8
nl_abs_tol = 1e-6
nl_max_its = 20
l_tol = 1e-6
l_max_its = 200
# the start_time can be set earlier to initialize SAM on its own
# sub_cycling must be set to true in the parent application to do that
start_time = -200.0
end_time = 200.0
num_steps = 20000
[Quadrature]
type = TRAP
order = FIRST
[]
[]
[Outputs]
print_linear_residuals = false
perf_graph = true
[out_displaced]
type = Exodus
use_displaced = true
execute_on = 'initial timestep_end'
sequence = false
[]
[csv]
type = CSV
[]
[checkpoint]
type = Checkpoint
num_files = 2
[]
[console]
type = Console
execute_scalars_on = 'none'
[]
[]
(msr/msfr/plant/steady/msfr_system_1d.i)
################################################################################
## Molten Salt Fast Reactor - Euratom EVOL + Rosatom MARS Design ##
## SAM Sub-Sub-Application input file ##
## Relaxation towards steady state balance of plant model ##
################################################################################
# If using or referring to this model, please cite as explained in
# https://mooseframework.inl.gov/virtual_test_bed/citing.html
core_power_val = 3e9
[Application]
type = SamApp
[]
[GlobalParams]
global_init_P = 1e5
global_init_V = 0.01
global_init_T = 898.15
Tsolid_sf = 1e-3
scaling_factor_var = '1 1e-3 1e-6'
[]
[Functions]
[fuel_salt_rho_func] # Linear fitting used by Rouch et al. 2014
type = PiecewiseLinear
x = '800 1200'
y = '4277.96 3925.16'
[]
[fuel_salt_cp_func]
type = PiecewiseLinear
x = '800 1200'
y = '1113.00 2225.00'
[]
[fuel_salt_k_func]
type = PiecewiseLinear
x = '800 1200'
y = '0.995200 1.028800'
[]
[fuel_salt_mu_func] # Nonlinear fitting used by Rouch et al. 2014
type = PiecewiseLinear
x = ' 800 820 840 860 880 900
920 940 960 980 1000 1020
1040 1060 1080 1100 1120 1140
1160 1180 1200'
y = '2.3887E-02 2.1258E-02 1.9020E-02 1.7102E-02 1.5449E-02 1.4015E-02
1.2767E-02 1.1673E-02 1.0711E-02 9.8609E-03 9.1066E-03 8.4347E-03
7.8340E-03 7.2951E-03 6.8099E-03 6.3719E-03 5.9751E-03 5.6147E-03
5.2865E-03 4.9868E-03 4.7124E-03'
[]
[fuel_salt_enthalpy_func] # approximated by cp*T
type = PiecewiseLinear
x = ' 800 820 840 860 880 900
920 940 960 980 1000 1020
1040 1060 1080 1100 1120 1140
1160 1180 1200'
y = '890400 913216 937144 962184 988336 1015600
1043976 1073464 1104064 1135776 1168600 1202536
1237584 1273744 1311016 1349400 1388896 1429504
1471224 1514056 1558000'
[]
[TimeStepperFunc]
type = PiecewiseLinear
#x = '-200 -180 -160 -140 -100 -60 -50 -20 -10. 0.0 100'
#y = ' 0.2 0.5 1.0 2.0 5.0 10.0 50. 10. 1.0 0.1 0.5'
x = '-200 -180 -160 -140 -100 -60 -50 -20 -10. 0.0 100'
y = ' 1.0 2.0 5.0 10.0 20.0 20.0 50. 10. 1.0 0.1 0.5'
[]
[core_specific_power]
type = ParsedFunction
expression = '${fparse core_power_val / 9.27}'
[]
[pump_head_fun]
type = PiecewiseConstant
xy_data = '0.0 156010.45
5.0 78005.225'
direction = 'left'
[]
[]
[EOS]
[fuel_salt_eos]
type = PTFunctionsEOS
rho = fuel_salt_rho_func
# beta = 2.1051E-04 # thermal expansion coefficient - approximated by -(drho/dT)/rho
cp = fuel_salt_cp_func
mu = fuel_salt_mu_func
k = fuel_salt_k_func
enthalpy = fuel_salt_enthalpy_func
[]
[hx_salt_eos]
type = SaltEquationOfState
salt_type = Flibe
[]
[Helium]
type = HeEquationOfState
[]
[]
[MaterialProperties]
[alloy-mat] # Based on Hastelloy N alloy
type = SolidMaterialProps
k = 23.6
Cp = 578
rho = 8.86e3
[]
[]
[Components]
[MSFR_core]
type = PBOneDFluidComponent
eos = fuel_salt_eos
orientation = '0 0 1'
position = '0 0 0'
A = 3.4636
Dh = 2.1
length = 2.65
n_elems = 25
heat_source = core_specific_power
[]
[Core_Out_Branch]
type = PBBranch
inputs = 'MSFR_core(out)'
outputs = 'pipe1(in) loop_pbc(in)'
eos = fuel_salt_eos
K = '4.25 4.25 500.'
Area = 2.24
[]
[pipe1] # Horizontal hot channel
type = PBOneDFluidComponent
eos = fuel_salt_eos
position = '0 0 2.65'
orientation = '0 1 0'
A = 2.24
Dh = 1.6888
length = 2.0
n_elems = 12
[]
[loop_pbc]
type = PBOneDFluidComponent
eos = fuel_salt_eos
position = '0 0 2.65'
orientation = '0 0 1'
A = 2.24
Dh = 1.6888
length = 0.025
n_elems = 1
[]
[TDV1]
type = PBTDV
input = 'loop_pbc(out)'
eos = fuel_salt_eos
p_bc = 1.0e5
[]
#
# ====== pump ======
#
[pump]
type = PBPump
Area = 2.24
K = '0.15 0.1'
eos = fuel_salt_eos
inputs = 'pipe1(out)'
outputs = 'pipe2(in)'
initial_P = 1.0e5
Head = 156010.45
[]
[pipe2] # Vertical hot channel from pump to HX
type = PBOneDFluidComponent
eos = fuel_salt_eos
position = '0 2.0 2.65'
orientation = '0 0 -1'
A = 2.24
Dh = 1.6888
length = 0.25
n_elems = 3
[]
#
# ====== Heat Exchanger ======
#
[IHX1]
type = PBHeatExchanger
eos = fuel_salt_eos
eos_secondary = hx_salt_eos
position = '0 2.0 2.4'
orientation = '0 0 -1'
A = 2.24
A_secondary = 3.60
Dh = 0.02
Dh_secondary = 0.0125
length = 2.4
n_elems = 24
SC_HTC = 4
HTC_geometry_type = Pipe
HTC_geometry_type_secondary = Pipe
HT_surface_area_density = 800.0
HT_surface_area_density_secondary = 497.78
initial_V_secondary = -2.67
Twall_init = 898.15
wall_thickness = 0.001
dim_wall = 1
material_wall = alloy-mat
n_wall_elems = 2
[]
[J_P2_IHX1]
type = PBBranch
inputs = 'pipe2(out)'
outputs = 'IHX1(primary_in) '
eos = fuel_salt_eos
K = '1. 1.'
Area = 2.24
[]
#
# ====== Heat Exchanger Secondary Side ======
#
[IHX1_P3]
type = PBBranch
inputs = 'IHX1(primary_out)'
outputs = 'pipe3(in) '
eos = fuel_salt_eos
K = '1. 1.'
Area = 2.24
[]
[pipe3] # Horizontal cold channel
type = PBOneDFluidComponent
eos = fuel_salt_eos
position = '0 2.0 0.'
orientation = '0 -1 0'
A = 2.24
Dh = 1.6888
length = 2.0
n_elems = 12
[]
[Core_In_Branch]
type = PBBranch
inputs = 'pipe3(out)'
outputs = 'MSFR_core(in)'
eos = fuel_salt_eos
K = '4.25 4.25'
Area = 2.24
[]
#
# ====== Intermediate circuit connected to HX1 ======
#
[IHX1_P4]
type = PBBranch
inputs = 'IHX1(secondary_out)'
outputs = 'pipe4(in), loop2_pbc(in)'
eos = hx_salt_eos
K = '4.25 4.25 500.'
Area = 3.6
[]
[loop2_pbc]
type = PBOneDFluidComponent
eos = hx_salt_eos
position = '0 2.0 2.4'
orientation = '1 0 0'
A = 3.6
Dh = 2.14
length = 0.025
n_elems = 1
[]
[loop2_TDV1]
type = PBTDV
input = 'loop2_pbc(out)'
eos = hx_salt_eos
p_bc = 1.0e5
[]
[pipe4] # Horizontal intermediate hot leg
type = PBOneDFluidComponent
eos = hx_salt_eos
position = '0 2.0 2.4'
orientation = '0 1 0'
A = 3.6
Dh = 2.14
length = 2.0
n_elems = 10
[]
[J_P4_P5]
type = PBBranch
inputs = 'pipe4(out)'
outputs = 'pipe5(in)'
eos = hx_salt_eos
K = '1. 1.'
Area = 3.6
[]
[pipe5] # Vertical intermediate hot leg
type = PBOneDFluidComponent
eos = hx_salt_eos
position = '0 4.0 2.4'
orientation = '0 0 -1'
A = 3.6
Dh = 2.14
length = 0.2
n_elems = 2
[]
[IHX2]
type = PBHeatExchanger
eos = hx_salt_eos
eos_secondary = Helium
position = '0 4.0 2.2'
orientation = '0 0 -1'
A = 2.4
A_secondary = 7.2
Dh = 0.04
Dh_secondary = 0.01
length = 3.2
n_elems = 24
SC_HTC = 4
HTC_geometry_type = Pipe
HTC_geometry_type_secondary = Pipe
HT_surface_area_density = 510.0
HT_surface_area_density_secondary = 170.0
initial_V_secondary = -50
Twall_init = 900
wall_thickness = 0.002
dim_wall = 1
material_wall = alloy-mat
n_wall_elems = 2
[]
[J_P5_IHX2]
type = PBBranch
inputs = 'pipe5(out)'
outputs = 'IHX2(primary_in)'
eos = hx_salt_eos
K = '1. 1.'
Area = 3.6
[]
[IHX2_S_In]
type = PBTDJ
input = 'IHX2(secondary_in)'
eos = Helium
v_bc = -66.65
T_bc = 673.15
[]
[IHX2_S_Out]
type = PBTDV
input = 'IHX2(secondary_out)'
eos = Helium
p_bc = 7.5e6
[]
[J_IHX2_P6]
type = PBBranch
inputs = 'IHX2(primary_out)'
outputs = 'pipe6(in)'
eos = hx_salt_eos
K = '1. 1.'
Area = 3.6
[]
[pipe6] # Intermediate cold leg
type = PBOneDFluidComponent
eos = hx_salt_eos
position = '0 4.0 -1.0'
orientation = '0 -1 0'
A = 3.6
Dh = 2.14
length = 1.0
n_elems = 5
[]
[J_P6_P7]
type = PBBranch
inputs = 'pipe6(out)'
outputs = 'pipe7(in)'
eos = hx_salt_eos
K = '1. 1.'
Area = 3.6
[]
[pipe7] # Intermediate cold leg
type = PBOneDFluidComponent
eos = hx_salt_eos
position = '0 3.0 -1.0'
orientation = '0 0 1'
A = 3.6
Dh = 2.14
length = 1.0
n_elems = 5
[]
[pump2]
type = PBPump
Area = 3.6
K = '0.15 0.1'
eos = hx_salt_eos
inputs = 'pipe7(out)'
outputs = 'pipe8(in)'
initial_P = 1.0e5
Head = 1.76e5
[]
[pipe8] # Intermediate cold leg
type = PBOneDFluidComponent
eos = hx_salt_eos
position = '0 3.0 0.0'
orientation = '0 -1 0'
A = 3.6
Dh = 2.14
length = 1.0
n_elems = 5
[]
[P8_IHX1]
type = PBBranch
inputs = 'pipe8(out)'
outputs = 'IHX1(secondary_in)'
eos = hx_salt_eos
K = '1. 1.'
Area = 3.6
[]
[]
[Postprocessors]
[Core_P_out]
type = ComponentBoundaryVariableValue
variable = pressure
input = MSFR_core(out)
[]
[Core_T_out]
type = ComponentBoundaryVariableValue
variable = temperature
input = MSFR_core(out)
[]
[Fuel_mass_flow] # Output mass flow rate at inlet of CH1
type = ComponentBoundaryFlow
input = MSFR_core(in)
[]
[Core_T_in]
type = ComponentBoundaryVariableValue
variable = temperature
input = MSFR_core(in)
[]
[Core_P_in]
type = ComponentBoundaryVariableValue
variable = pressure
input = MSFR_core(in)
[]
[HX_Tin_p]
type = ComponentBoundaryVariableValue
variable = temperature
input = IHX1(primary_in)
[]
[HX_Tout_p]
type = ComponentBoundaryVariableValue
variable = temperature
input = IHX1(primary_out)
[]
[HX_Tout_s]
type = ComponentBoundaryVariableValue
variable = temperature
input = IHX1(secondary_out)
[]
[HX_Uout_s]
type = ComponentBoundaryVariableValue
variable = velocity
input = IHX1(secondary_out)
[]
[FliBe_mass_flow]
type = ComponentBoundaryFlow
input = IHX1(secondary_in)
[]
[HX_Pin_s]
type = ComponentBoundaryVariableValue
variable = pressure
input = IHX1(secondary_in)
[]
[FliBe_mass_flow2]
type = ComponentBoundaryFlow
input = IHX2(primary_in)
[]
[Gas_mass_flow]
type = ComponentBoundaryFlow
input = IHX2(secondary_in)
[]
[HX2_Tin_p]
type = ComponentBoundaryVariableValue
variable = temperature
input = IHX2(primary_in)
[]
[HX2_Tout_p]
type = ComponentBoundaryVariableValue
variable = temperature
input = IHX2(primary_out)
[]
[HX2_Tout_s]
type = ComponentBoundaryVariableValue
variable = temperature
input = IHX2(secondary_out)
[]
[HX_cold_temp]
type = ElementAverageValue
variable = temperature
block = 'pipe3'
[]
[core_power]
type = Receiver
default = ${core_power_val}
execute_on = 'INITIAL TIMESTEP_BEGIN TIMESTEP_END TRANSFER'
[]
[]
[Preconditioning]
[SMP_PJFNK]
type = SMP
full = true
solve_type = 'PJFNK'
petsc_options_iname = '-pc_type -ksp_gmres_restart'
petsc_options_value = 'lu 101'
[]
[]
[Executioner]
type = Transient
[TimeStepper]
type = FunctionDT
function = TimeStepperFunc
min_dt = 1e-3
[]
#[TimeStepper]
# type = AB2PredictorCorrector
# dt = .01
# e_max = 10
# e_tol = 1
#[]
nl_rel_tol = 1e-8
nl_abs_tol = 1e-6
nl_max_its = 20
l_tol = 1e-6
l_max_its = 200
# the start_time can be set earlier to initialize SAM on its own
# sub_cycling must be set to true in the parent application to do that
start_time = -200.0
end_time = 200.0
num_steps = 20000
[Quadrature]
type = TRAP
order = FIRST
[]
[]
[Outputs]
print_linear_residuals = false
perf_graph = true
[out_displaced]
type = Exodus
use_displaced = true
execute_on = 'initial timestep_end'
sequence = false
[]
[csv]
type = CSV
[]
[checkpoint]
type = Checkpoint
num_files = 2
[]
[console]
type = Console
execute_scalars_on = 'none'
[]
[]
(msr/msfr/plant/steady/run_ns.i)
################################################################################
## Molten Salt Fast Reactor - Euratom EVOL + Rosatom MARS Design ##
## Pronghorn Sub-Application input file ##
## Relaxation towards Steady state 3D thermal hydraulics model ##
################################################################################
# If using or referring to this model, please cite as explained in
# https://mooseframework.inl.gov/virtual_test_bed/citing.html
# The flow in this simulation should be initialized with a previous flow
# solution (isothermal, heated or multiphysics) OR using a viscosity rampdown
# An isothermal solution can be generated using 'run_ns_initial.i', and is
# saved in 'restart/run_ns_initial_restart.e'
# A heated solution, with a flat power distribution, can be generated with
# this script 'run_ns.i', and is saved in 'restart/run_ns_restart.e'
# A coupled neutronics-coarse TH solution can be generated with
# 'run_neutronics.i', saved in 'restart/run_neutronics_ns_restart.e'
# Material properties
rho = 4284 # density [kg / m^3] (@1000K)
cp = 1594 # specific heat capacity [J / kg / K]
drho_dT = 0.882 # derivative of density w.r.t. temperature [kg / m^3 / K]
mu = 0.0166 # viscosity [Pa s]
k = 1.7 # thermal conductivity [W / m / K]
# https://www.researchgate.net/publication/337161399_Development_of_a_control-\
# oriented_power_plant_simulator_for_the_molten_salt_fast_reactor/fulltext/5dc9\
# 5c9da6fdcc57503eec39/Development-of-a-control-oriented-power-plant-simulator-\
# for-the-molten-salt-fast-reactor.pdf
von_karman_const = 0.41
# Turbulent properties
Pr_t = 0.9 # turbulent Prandtl number
Sc_t = 1 # turbulent Schmidt number
# Derived material properties
alpha = '${fparse drho_dT / rho}' # thermal expansion coefficient
# Operating parameters
T_HX = 873.15 # heat exchanger temperature [K]
# Mass flow rate tuning, for heat exchanger pressure and temperature drop
friction = 4e3 # [kg / m^4]
pump_force = -20000. # [N / m^3]
# Delayed neutron precursor parameters. Lambda values are decay constants in
# [1 / s]. Beta values are production fractions.
lambda1_m = -0.0133104
lambda2_m = -0.0305427
lambda3_m = -0.115179
lambda4_m = -0.301152
lambda5_m = -0.879376
lambda6_m = -2.91303
beta1 = 8.42817e-05
beta2 = 0.000684616
beta3 = 0.000479796
beta4 = 0.00103883
beta5 = 0.000549185
beta6 = 0.000184087
[GlobalParams]
rhie_chow_user_object = 'ins_rhie_chow_interpolator'
[]
################################################################################
# GEOMETRY
################################################################################
[Mesh]
coord_type = 'RZ'
rz_coord_axis = Y
[restart]
type = FileMeshGenerator
use_for_exodus_restart = true
# Depending on the file chosen, the initialization of variables should be
# adjusted. The following variables can be initalized:
# - vel_x, vel_y, p from isothermal simulation
# file = '../../steady/restart/run_ns_initial_restart.e'
# Below are initialization points created from this input file
# The variable IC should be set from_file_var for temperature and precursors
# - vel_x, vel_y, p, T_fluid, c_i from cosine heated simulation
# file = '../../steady/restart/run_ns_restart.e'
# - adding SAM-coupling
file = 'restart/run_ns_restart.e'
[]
[hx_top]
type = ParsedGenerateSideset
combinatorial_geometry = 'y > 0'
included_subdomains = '3'
included_neighbors = '1'
fixed_normal = true
normal = '0 1 0'
new_sideset_name = 'hx_top'
input = 'restart'
[]
[hx_bot]
type = ParsedGenerateSideset
combinatorial_geometry = 'y <-0.6'
included_subdomains = '3'
included_neighbors = '1'
fixed_normal = true
normal = '0 -1 0'
new_sideset_name = 'hx_bot'
input = 'hx_top'
[]
[]
[Problem]
# Using restart for velocities & pressure initialization
allow_initial_conditions_with_restart = true
[]
################################################################################
# EQUATIONS: VARIABLES, KERNELS, BOUNDARY CONDITIONS
################################################################################
[Modules]
[NavierStokesFV]
# General parameters
compressibility = 'incompressible'
add_energy_equation = true
add_scalar_equation = true
boussinesq_approximation = true
block = 'fuel pump hx'
# Variables, defined below for the Exodus restart
velocity_variable = 'vel_x vel_y'
pressure_variable = 'pressure'
fluid_temperature_variable = 'T_fluid'
# Material properties
density = ${rho}
dynamic_viscosity = ${mu}
thermal_conductivity = ${k}
specific_heat = 'cp'
thermal_expansion = ${alpha}
# Boussinesq parameters
gravity = '0 -9.81 0'
ref_temperature = ${T_HX}
# Boundary conditions
wall_boundaries = 'shield_wall reflector_wall fluid_symmetry'
momentum_wall_types = 'wallfunction wallfunction symmetry'
energy_wall_types = 'heatflux heatflux heatflux'
energy_wall_function = '0 0 0'
# Pressure pin for incompressible flow
pin_pressure = true
pinned_pressure_type = average
pinned_pressure_value = 1e5
# Turbulence parameters
turbulence_handling = 'mixing-length'
turbulent_prandtl = ${Pr_t}
von_karman_const = ${von_karman_const}
mixing_length_delta = 0.1
mixing_length_walls = 'shield_wall reflector_wall'
mixing_length_aux_execute_on = 'initial'
# Numerical scheme
momentum_advection_interpolation = 'upwind'
mass_advection_interpolation = 'upwind'
energy_advection_interpolation = 'upwind'
passive_scalar_advection_interpolation = 'upwind'
# Heat source
external_heat_source = power_density
# Precursor advection, diffusion and source term
passive_scalar_names = 'c1 c2 c3 c4 c5 c6'
passive_scalar_schmidt_number = '${Sc_t} ${Sc_t} ${Sc_t} ${Sc_t} ${Sc_t} ${Sc_t}'
passive_scalar_coupled_source = 'fission_source c1; fission_source c2; fission_source c3;
fission_source c4; fission_source c5; fission_source c6'
passive_scalar_coupled_source_coeff = '${beta1} ${lambda1_m}; ${beta2} ${lambda2_m};
${beta3} ${lambda3_m}; ${beta4} ${lambda4_m};
${beta5} ${lambda5_m}; ${beta6} ${lambda6_m}'
# Heat exchanger
standard_friction_formulation = false
friction_blocks = 'hx'
friction_types = 'FORCHHEIMER'
friction_coeffs = 'friction_coeff_vector'
ambient_convection_blocks = 'hx'
ambient_convection_alpha = '${fparse 600 * 20e3}' # HX specifications
ambient_temperature = 'hx_cold_temp'
[]
[]
[Variables]
[vel_x]
type = INSFVVelocityVariable
block = 'fuel pump hx'
initial_from_file_var = vel_x
[]
[vel_y]
type = INSFVVelocityVariable
block = 'fuel pump hx'
initial_from_file_var = vel_y
[]
[pressure]
type = INSFVPressureVariable
block = 'fuel pump hx'
initial_from_file_var = pressure
[]
[T_fluid]
type = INSFVEnergyVariable
block = 'fuel pump hx'
initial_condition = ${T_HX}
# initial_from_file_var = T_fluid
[]
[c1]
type = MooseVariableFVReal
block = 'fuel pump hx'
# initial_from_file_var = c1
[]
[c2]
type = MooseVariableFVReal
block = 'fuel pump hx'
# initial_from_file_var = c2
[]
[c3]
type = MooseVariableFVReal
block = 'fuel pump hx'
# initial_from_file_var = c3
[]
[c4]
type = MooseVariableFVReal
block = 'fuel pump hx'
# initial_from_file_var = c4
[]
[c5]
type = MooseVariableFVReal
block = 'fuel pump hx'
# initial_from_file_var = c5
[]
[c6]
type = MooseVariableFVReal
block = 'fuel pump hx'
# initial_from_file_var = c6
[]
[]
[FVICs]
[c1]
type = FVFunctionIC
variable = c1
function = 'cosine_guess'
scaling_factor = 0.02
[]
[c2]
type = FVFunctionIC
variable = c2
function = 'cosine_guess'
scaling_factor = 0.1
[]
[c3]
type = FVFunctionIC
variable = c3
function = 'cosine_guess'
scaling_factor = 0.03
[]
[c4]
type = FVFunctionIC
variable = c4
function = 'cosine_guess'
scaling_factor = 0.04
[]
[c5]
type = FVFunctionIC
variable = c5
function = 'cosine_guess'
scaling_factor = 0.01
[]
[c6]
type = FVFunctionIC
variable = c6
function = 'cosine_guess'
scaling_factor = 0.001
[]
# Power density is re-initalized by a transfer from neutronics
[power]
type = FVFunctionIC
variable = power_density
function = 'cosine_guess'
scaling_factor = '${fparse 3e9/2.81543}'
[]
# Fission source is re-initalized by a transfer from neutronics
[fission_source]
type = FVFunctionIC
variable = fission_source
function = 'cosine_guess'
scaling_factor = '${fparse 6.303329e+01/2.81543}'
[]
[]
[AuxVariables]
[power_density]
type = MooseVariableFVReal
block = 'fuel pump hx'
[]
[fission_source]
type = MooseVariableFVReal
block = 'fuel pump hx'
[]
[]
[Functions]
# Guess to have a 3D power distribution
[cosine_guess]
type = ParsedFunction
expression = 'max(0, cos(x*pi/2/1.2))*max(0, cos(y*pi/2/1.1))'
[]
[hx_cold_temp]
type = ParsedFunction
expression = 'Tcold'
symbol_names = 'Tcold'
symbol_values = '${T_HX}'
[]
[]
[FVKernels]
[pump]
type = INSFVBodyForce
variable = vel_y
functor = ${pump_force}
block = 'pump'
momentum_component = 'y'
[]
[]
################################################################################
# MATERIALS
################################################################################
[FunctorMaterials]
# Most of these constants could be specified directly to the action
[mu]
type = ADGenericFunctorMaterial
prop_names = 'mu'
prop_values = '${mu}'
block = 'fuel pump hx'
[]
# [not_used]
# type = ADGenericFunctorMaterial
# prop_names = 'not_used'
# prop_values = 0
# block = 'shield reflector'
# []
[cp]
type = ADGenericFunctorMaterial
prop_names = 'cp'
prop_values = '${cp}'
block = 'fuel pump hx'
[]
[friction_coeff_mat]
type = ADGenericVectorFunctorMaterial
prop_names = 'friction_coeff_vector'
prop_values = '${friction} ${friction} ${friction}'
[]
[]
################################################################################
# EXECUTION / SOLVE
################################################################################
[Functions]
[dts]
type = PiecewiseConstant
x = '0 100'
y = '0.75 2.5'
[]
[]
[Executioner]
type = Transient
# Time stepping parameters
start_time = 0.0
end_time = 20
# end_time will depend on the restart file chosen
# though steady state detection can also be used
# from _initial/no heating : 150 - 200s enough
# from _ns/_ns_coupled/heated: 10s enough
[TimeStepper]
# This time stepper makes the time step grow exponentially
# It can only be used with proper initialization
type = IterationAdaptiveDT
dt = 1 # chosen to obtain convergence with first coupled iteration
growth_factor = 2
[]
# [TimeStepper]
# type = FunctionDT
# function = dts
# []
steady_state_detection = true
steady_state_tolerance = 1e-8
steady_state_start_time = 10
# Time integration scheme
scheme = 'implicit-euler'
# Solver parameters
solve_type = 'NEWTON'
petsc_options_iname = '-pc_type -pc_factor_shift_type -ksp_gmres_restart'
petsc_options_value = 'lu NONZERO 50'
line_search = 'none'
nl_rel_tol = 1e-9
nl_abs_tol = 2e-8
nl_max_its = 15
l_max_its = 50
automatic_scaling = true
# Fixed point iterations are not necessary if the only goal is to obtain
# a steady state solution, no matter the numerical path
# fixed_point_max_its = 10
# fixed_point_abs_tol = 1e-5
# accept_on_max_fixed_point_iteration = true
[]
################################################################################
# MULTIAPPS FOR POWER TRANSFER
################################################################################
[MultiApps]
[sam_balance_of_plant]
type = TransientMultiApp
app_type = 'SamApp'
input_files = msfr_system_1d.i
max_procs_per_app = 1
execute_on = 'timestep_end'
keep_solution_during_restore = False
# the balance of plant needs to be run separately first for 100 steps
# here we are hardcoding a restart file path, to a file we checked in, but you should
# re-create the checkpoint file and adapt the restart_file_base path for your own runs
cli_args = 'Problem/restart_file_base=restart/msfr_system_1d_checkpoint_cp/0100;Problem/force_restart=true'
[]
[]
[Transfers]
# Primary and secondary loops
[send_flow_BCs]
type = MultiAppPostprocessorVectorTransfer
to_multi_app = sam_balance_of_plant
from_postprocessors = 'total_power'
to_postprocessors = 'core_power'
[]
[receive_flow_BCs]
type = MultiAppPostprocessorVectorTransfer
from_multi_app = sam_balance_of_plant
from_postprocessors = 'HX_cold_temp'
to_postprocessors = 'heat_exchanger_T_ambient'
[]
[]
################################################################################
# SIMULATION OUTPUTS
################################################################################
[Outputs]
csv = true
hide = 'flow_hx_bot flow_hx_top min_flow_T max_flow_T'
[restart]
type = Exodus
overwrite = true
[]
# Reduce base output
print_linear_converged_reason = false
print_linear_residuals = false
print_nonlinear_converged_reason = false
[]
[Postprocessors]
[max_v]
type = ElementExtremeValue
variable = vel_x
value_type = max
block = 'fuel pump hx'
[]
# TODO: weakly compressible, switch to mass flow rate
[flow_hx_bot]
type = VolumetricFlowRate
boundary = 'hx_bot'
vel_x = vel_x
vel_y = vel_y
advected_quantity = 1
[]
[flow_hx_top]
type = VolumetricFlowRate
boundary = 'hx_top'
vel_x = vel_x
vel_y = vel_y
advected_quantity = 1
[]
[max_flow_T]
type = VolumetricFlowRate
boundary = 'hx_top'
vel_x = vel_x
vel_y = vel_y
advected_quantity = 'T_fluid'
[]
[min_flow_T]
type = VolumetricFlowRate
boundary = 'hx_bot'
vel_x = vel_x
vel_y = vel_y
advected_quantity = 'T_fluid'
[]
[dT]
type = ParsedPostprocessor
expression = '-max_flow_T / flow_hx_bot + min_flow_T / flow_hx_top'
pp_names = 'max_flow_T min_flow_T flow_hx_bot flow_hx_top'
[]
[total_power]
type = ElementIntegralVariablePostprocessor
variable = power_density
block = 'fuel pump hx'
[]
[total_fission_source]
type = ElementIntegralVariablePostprocessor
variable = fission_source
block = 'fuel pump hx'
[]
[heat_exchanger_T_ambient]
type = Receiver
default = ${T_HX}
execute_on = 'INITIAL TIMESTEP_END TRANSFER'
[]
[]