- epsilonCoupled turbulent kinetic energy dissipation rate. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
C++ Type:MooseFunctorName
Unit:(no unit assumed)
Controllable:No
Description:Coupled turbulent kinetic energy dissipation rate. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
- muDynamic viscosity. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
C++ Type:MooseFunctorName
Unit:(no unit assumed)
Controllable:No
Description:Dynamic viscosity. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
- rhoDensity. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
C++ Type:MooseFunctorName
Unit:(no unit assumed)
Controllable:No
Description:Density. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
- tkeCoupled turbulent kinetic energy. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
C++ Type:MooseFunctorName
Unit:(no unit assumed)
Controllable:No
Description:Coupled turbulent kinetic energy. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
- uThe velocity in the x direction. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
C++ Type:MooseFunctorName
Unit:(no unit assumed)
Controllable:No
Description:The velocity in the x direction. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
- variableThe name of the variable that this object applies to
C++ Type:AuxVariableName
Unit:(no unit assumed)
Controllable:No
Description:The name of the variable that this object applies to
kEpsilonViscosity
kEpsilonViscosity is an auxiliary kernel that computes the turbulent dynamic viscosity used in the k– family of turbulence models. It forms the closure term for the Reynolds stresses in the momentum equations:
where the turbulent time scale and the effective coefficient depend on the selected k– variant and the local flow state.
This object implements:
Standard, Low-Re, Two-Layer, and Realizable viscosity formulations,
Several wall treatments (Newton, incremental, linearized, non-equilibrium),
Two-layer near-wall blending methods (Wolfstein (1969), NORRIS III (1975), Xu et al. (1998)),
Low-Re damping functions for the StandardLowRe model,
Realizable variable ,
Optional scale limiting for the time scale,
Bulk and near-wall formulations that match Menter-type corrections Menter (1994).
Turbulent viscosity in the bulk region
Away from walls (or when [bulk_wall_treatment=false]), the turbulent viscosity is computed as:
where:
is the turbulent kinetic energy,
is the fluid density,
depends on the k– variant,
is the turbulent time scale (possibly limited).
Turbulent time scale
The base time scale is:
If [scale_limiter=standard], the time scale is limited using:
where:
is the molecular kinematic viscosity,
is a model constant,
The second term becomes important near walls or when becomes too large.
If [scale_limiter=none], then .
Effective coefficient
The meaning of varies based on the selected turbulence variant:
Standard k–
StandardLowRe
Low-Re models use a damping function :
where is computed as:
with:
and is the wall distance (supplied via a functor).
Realizable k–
The realizable model computes dynamically as a function of strain and rotation:
where:
,
.
This value replaces the standard constant everywhere in the model.
Two-layer models (StandardTwoLayer and RealizableTwoLayer)
Two-layer models blend between:
an outer k– viscosity,
an inner near-wall viscosity based on two-layer relations.
Let:
,
= near-wall turbulent viscosity from the chosen flavor (Wolfstein, Norris–Reynolds, or Xu), see the theory section for more details.
Using the wall-distance Reynolds number , the two-layer blending function is:
with:
,
chosen such that at .
The final viscosity is:
Near-wall bulk treatment
When [bulk_wall_treatment=true], the kernel applies wall functions to compute in wall-bounded cells.
The procedure:
Identify the minimum wall distance in the cell.
Compute the tangential velocity and friction velocity .
Compute the nondimensional wall distance:
Classify into:
Viscous sublayer: → ,
Log‑layer: → (see theory for more details),
Buffer layer: linear blending between the two regimes.
Four wall-function types are supported:
eq_newtoneq_incrementaleq_linearizedneq(non-equilibrium using local )
Each provides a different method for computing and the corresponding .
This wall-function implementation is fully consistent with the finite-volume INSFVTurbulentViscosity models used in the Navier–Stokes FV physics.
Wall distance requirements
The following variants require a wall-distance functor:
StandardTwoLayerRealizableTwoLayerStandardLowRe
If a required functor is missing, the kernel throws an error during construction.
Interaction with k and kernels
kEpsilonViscosity is tightly coupled with the k and kernels:
It uses
kand to compute viscosity.The k and kernels require to compute production and destruction terms.
Realizable and nonlinear models require strain/rotation invariants also used by k and kernels.
Two-layer models require consistency between: - -equation two-layer region, - viscosity two-layer region.
This object must therefore be included in every k– turbulence simulation.
Input Parameters
- C_mu0.09Base turbulent kinetic energy closure constant C_mu.
Default:0.09
C++ Type:Real
Unit:(no unit assumed)
Controllable:No
Description:Base turbulent kinetic energy closure constant C_mu.
- Ca00.667Realizable k-epsilon coefficient Ca0 used in the C_mu expression.
Default:0.667
C++ Type:Real
Unit:(no unit assumed)
Controllable:No
Description:Realizable k-epsilon coefficient Ca0 used in the C_mu expression.
- Ca11.25Realizable k-epsilon coefficient Ca1 used in the C_mu expression.
Default:1.25
C++ Type:Real
Unit:(no unit assumed)
Controllable:No
Description:Realizable k-epsilon coefficient Ca1 used in the C_mu expression.
- Ca21Realizable k-epsilon coefficient Ca2 used in the C_mu expression.
Default:1
C++ Type:Real
Unit:(no unit assumed)
Controllable:No
Description:Realizable k-epsilon coefficient Ca2 used in the C_mu expression.
- Ca30.9Realizable k-epsilon coefficient Ca3 used in the C_mu expression.
Default:0.9
C++ Type:Real
Unit:(no unit assumed)
Controllable:No
Description:Realizable k-epsilon coefficient Ca3 used in the C_mu expression.
- Cd00.091Low-Re coefficient Cd0 used in the f_mu damping function for the Standard k-epsilon Low-Re model.
Default:0.091
C++ Type:Real
Unit:(no unit assumed)
Controllable:No
Description:Low-Re coefficient Cd0 used in the f_mu damping function for the Standard k-epsilon Low-Re model.
- Cd10.0042Low-Re coefficient Cd1 used in the f_mu damping function for the Standard k-epsilon Low-Re model.
Default:0.0042
C++ Type:Real
Unit:(no unit assumed)
Controllable:No
Description:Low-Re coefficient Cd1 used in the f_mu damping function for the Standard k-epsilon Low-Re model.
- Cd20.00011Low-Re coefficient Cd2 used in the f_mu damping function for the Standard k-epsilon Low-Re model.
Default:0.00011
C++ Type:Real
Unit:(no unit assumed)
Controllable:No
Description:Low-Re coefficient Cd2 used in the f_mu damping function for the Standard k-epsilon Low-Re model.
- Ct1Time-scale constant Ct used in the k-epsilon time scale T = max(Te, Ct*sqrt(nu/epsilon)) when scale_limiter='standard'.
Default:1
C++ Type:Real
Unit:(no unit assumed)
Controllable:No
Description:Time-scale constant Ct used in the k-epsilon time scale T = max(Te, Ct*sqrt(nu/epsilon)) when scale_limiter='standard'.
- blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
- boundaryThe list of boundaries (ids or names) from the mesh where this object applies
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The list of boundaries (ids or names) from the mesh where this object applies
- bulk_wall_treatmentFalseIf true, use classical wall functions for near-wall cells.
Default:False
C++ Type:bool
Controllable:No
Description:If true, use classical wall functions for near-wall cells.
- check_boundary_restrictedTrueWhether to check for multiple element sides on the boundary in the case of a boundary restricted, element aux variable. Setting this to false will allow contribution to a single element's elemental value(s) from multiple boundary sides on the same element (example: when the restricted boundary exists on two or more sides of an element, such as at a corner of a mesh
Default:True
C++ Type:bool
Controllable:No
Description:Whether to check for multiple element sides on the boundary in the case of a boundary restricted, element aux variable. Setting this to false will allow contribution to a single element's elemental value(s) from multiple boundary sides on the same element (example: when the restricted boundary exists on two or more sides of an element, such as at a corner of a mesh
- execute_onLINEAR TIMESTEP_ENDThe list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html.
Default:LINEAR TIMESTEP_END
C++ Type:ExecFlagEnum
Options:NONE, INITIAL, LINEAR, LINEAR_CONVERGENCE, NONLINEAR, NONLINEAR_CONVERGENCE, POSTCHECK, TIMESTEP_END, TIMESTEP_BEGIN, MULTIAPP_FIXED_POINT_END, MULTIAPP_FIXED_POINT_BEGIN, MULTIAPP_FIXED_POINT_CONVERGENCE, MULTISYSTEM_FIXED_POINT_ITERATION_END, FINAL, CUSTOM, PRE_DISPLACE
Controllable:No
Description:The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html.
- k_epsilon_variantStandardk-epsilon model variant used for viscosity.
Default:Standard
C++ Type:MooseEnum
Options:Standard, StandardLowRe, StandardTwoLayer, Realizable, RealizableTwoLayer
Controllable:No
Description:k-epsilon model variant used for viscosity.
- mu_t_ratio_max100000Maximum allowable mu_t_ratio : mu/mu_t.
Default:100000
C++ Type:Real
Unit:(no unit assumed)
Controllable:No
Description:Maximum allowable mu_t_ratio : mu/mu_t.
- scale_limiterstandardThe method used to limit the k-epsilon time scale: 'none' or 'standard' (max(T_e, sqrt(nu/epsilon))).
Default:standard
C++ Type:MooseEnum
Options:none, standard
Controllable:No
Description:The method used to limit the k-epsilon time scale: 'none' or 'standard' (max(T_e, sqrt(nu/epsilon))).
- two_layer_flavorWolfsteinTwo-layer formulation to use for 2L variants (Wolfstein, NorrisReynolds, Xu).
Default:Wolfstein
C++ Type:MooseEnum
Options:Wolfstein, NorrisReynolds, Xu
Controllable:No
Description:Two-layer formulation to use for 2L variants (Wolfstein, NorrisReynolds, Xu).
- vThe velocity in the y direction. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
C++ Type:MooseFunctorName
Unit:(no unit assumed)
Controllable:No
Description:The velocity in the y direction. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
- wThe velocity in the z direction. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
C++ Type:MooseFunctorName
Unit:(no unit assumed)
Controllable:No
Description:The velocity in the z direction. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
- wall_distanceDistance to the closest wall; required for two-layer and Low-Re k-epsilon variants. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
C++ Type:MooseFunctorName
Unit:(no unit assumed)
Controllable:No
Description:Distance to the closest wall; required for two-layer and Low-Re k-epsilon variants. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
- wall_treatmenteq_newtonThe method used for computing the wall functions.
Default:eq_newton
C++ Type:MooseEnum
Options:eq_newton, eq_incremental, eq_linearized, neq
Controllable:No
Description:The method used for computing the wall functions.
- wallsBoundaries that correspond to solid walls.
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:Boundaries that correspond to solid walls.
Optional Parameters
- control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
- enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:Yes
Description:Set the enabled status of the MooseObject.
- search_methodnearest_node_connected_sidesChoice of search algorithm. All options begin by finding the nearest node in the primary boundary to a query point in the secondary boundary. In the default nearest_node_connected_sides algorithm, primary boundary elements are searched iff that nearest node is one of their nodes. This is fast to determine via a pregenerated node-to-elem map and is robust on conforming meshes. In the optional all_proximate_sides algorithm, primary boundary elements are searched iff they touch that nearest node, even if they are not topologically connected to it. This is more CPU-intensive but is necessary for robustness on any boundary surfaces which has disconnections (such as Flex IGA meshes) or non-conformity (such as hanging nodes in adaptively h-refined meshes).
Default:nearest_node_connected_sides
C++ Type:MooseEnum
Options:nearest_node_connected_sides, all_proximate_sides
Controllable:No
Description:Choice of search algorithm. All options begin by finding the nearest node in the primary boundary to a query point in the secondary boundary. In the default nearest_node_connected_sides algorithm, primary boundary elements are searched iff that nearest node is one of their nodes. This is fast to determine via a pregenerated node-to-elem map and is robust on conforming meshes. In the optional all_proximate_sides algorithm, primary boundary elements are searched iff they touch that nearest node, even if they are not topologically connected to it. This is more CPU-intensive but is necessary for robustness on any boundary surfaces which has disconnections (such as Flex IGA meshes) or non-conformity (such as hanging nodes in adaptively h-refined meshes).
- seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Controllable:No
Description:The seed for the master random number generator
- use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:False
C++ Type:bool
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Advanced Parameters
- prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
- use_interpolated_stateFalseFor the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Default:False
C++ Type:bool
Controllable:No
Description:For the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Material Property Retrieval Parameters
Input Files
- (test/tests/kEpsilon/special-cases/channel_ERCOFTAC_buoyant.i)
- (test/tests/kEpsilon/lid-driven-cavity/lid-driven.i)
- (test/tests/kEpsilon/special-cases/lid-driven-curvature-corrected.i)
- (validation/free_flow/isothermal/single-stream-jet/RTL10000/jet_RTL10000.i)
- (test/tests/kEpsilon/special-cases/channel_ERCOFTAC_low_Re.i)
- (test/tests/kEpsilon/special-cases/channel_ERCOFTAC.i)
- (test/tests/kEpsilon/special-cases/channel_ERCOFTAC_bulk_treatment.i)
- (validation/free_flow/isothermal/ercoftac_030_bfs/bfs_input.i)
- (test/tests/kEpsilon/channel/channel_ERCOFTAC.i)
References
- Florian R Menter.
Two-equation eddy-viscosity turbulence models for engineering applications.
AIAA journal, 32(8):1598–1605, 1994.[Export]
BibTeX
@article{menter1994two, author = "Menter, Florian R", title = "Two-equation eddy-viscosity turbulence models for engineering applications", journal = "AIAA journal", volume = "32", number = "8", pages = "1598--1605", year = "1994" }RIS
TY - JOUR AU - Menter, Florian R TI - Two-equation eddy-viscosity turbulence models for engineering applications JO - AIAA journal PY - 1994 VL - 32 IS - 8 SP - 1598 EP - 1605 ER -Plain Text
Florian R Menter. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA journal, 32(8):1598–1605, 1994. - H LEE NORRIS III.
Turbulent channel flow with a moving wavy boundary.
Stanford University, 1975.[Export]
BibTeX
@book{norris1975oneequation, author = "NORRIS III, H LEE", title = "Turbulent channel flow with a moving wavy boundary.", year = "1975", publisher = "Stanford University" }RIS
TY - BOOK AU - NORRIS III, H LEE TI - Turbulent channel flow with a moving wavy boundary. PY - 1975 PB - Stanford University ER -Plain Text
H LEE NORRIS III. Turbulent channel flow with a moving wavy boundary. Stanford University, 1975. - M. Wolfstein.
The velocity and temperature distribution in one-dimensional flow with turbulence augmentation and pressure gradient.
International Journal of Heat and Mass Transfer, 12(3):301–318, 1969.[Export]
BibTeX
@article{wolfstein1969velocity, author = "Wolfstein, M.", title = "The velocity and temperature distribution in one-dimensional flow with turbulence augmentation and pressure gradient", journal = "International Journal of Heat and Mass Transfer", volume = "12", number = "3", pages = "301--318", year = "1969" }RIS
TY - JOUR AU - Wolfstein, M. TI - The velocity and temperature distribution in one-dimensional flow with turbulence augmentation and pressure gradient JO - International Journal of Heat and Mass Transfer PY - 1969 VL - 12 IS - 3 SP - 301 EP - 318 ER -Plain Text
M. Wolfstein. The velocity and temperature distribution in one-dimensional flow with turbulence augmentation and pressure gradient. International Journal of Heat and Mass Transfer, 12(3):301–318, 1969. - Wen Xu, Q Chen, and FTM Nieuwstadt.
A new turbulence model for near-wall natural convection.
International Journal of Heat and Mass Transfer, 41(21):3161–3176, 1998.[Export]
BibTeX
@article{xu1998new, author = "Xu, Wen and Chen, Q and Nieuwstadt, FTM", title = "A new turbulence model for near-wall natural convection", journal = "International Journal of Heat and Mass Transfer", volume = "41", number = "21", pages = "3161--3176", year = "1998" }RIS
TY - JOUR AU - Xu, Wen AU - Chen, Q AU - Nieuwstadt, FTM TI - A new turbulence model for near-wall natural convection JO - International Journal of Heat and Mass Transfer PY - 1998 VL - 41 IS - 21 SP - 3161 EP - 3176 ER -Plain Text
Wen Xu, Q Chen, and FTM Nieuwstadt. A new turbulence model for near-wall natural convection. International Journal of Heat and Mass Transfer, 41(21):3161–3176, 1998.
(test/tests/kEpsilon/special-cases/channel_ERCOFTAC_buoyant.i)
##########################################################
# ERCOFTAC test case for turbulent channel flow
# Case Number: 032
# Author: Dr. Mauricio Tano & Hailey Tran Kieu
# Last Update: November, 2023 & 2025
# Turbulent model using:
# k-epsilon
# Equilibrium + Newton wall treatment
# SIMPLE solve
##########################################################
### Problem Parameters ###
H = 1 # half-width of the channel
L = 120
Re = 14000
rho_0 = 1
bulk_u = 1
mu = '${fparse rho_0 * bulk_u * 2 * H / Re}'
cp = 4186.0
k = 1.0
alpha_b = 1e-3
T_0 = 300.0
c = 1480.0
advected_interp_method = 'upwind'
### k-epsilon Closure Parameters ###
sigma_k = 1.0
sigma_eps = 1.3
C1_eps = 1.44
C2_eps = 1.92
C_mu = 0.09
### Initial and Boundary Conditions ###
intensity = 0.01
k_init = '${fparse 1.5*(intensity * bulk_u)^2}'
eps_init = '${fparse C_mu^0.75 * k_init^1.5 / (2*H)}'
### Modeling parameters ###
bulk_wall_treatment = false
walls = 'bottom top'
wall_treatment = 'eq_newton' # Options: eq_newton, eq_incremental, eq_linearized, neq
# Turbulence-model knobs (optional)
k_epsilon_variant = 'RealizableTwoLayer' # Standard | StandardLowRe | StandardTwoLayer | Realizable | RealizableTwoLayer
two_layer_flavor = 'Wolfstein' # Wolfstein | NorrisReynolds | Xu (only used for *TwoLayer variants)
use_buoyancy = true
use_compressibility = true
nonlinear_model = 'none'
curvature_model = 'none'
use_yap = false
use_low_re_Gprime = false
[Mesh]
[block_1]
type = GeneratedMeshGenerator
dim = 2
xmin = 0
xmax = ${L}
ymin = 0
ymax = ${H}
nx = 10
ny = 5
bias_y = 0.7
[]
[block_2]
type = GeneratedMeshGenerator
dim = 2
xmin = 0
xmax = ${L}
ymin = ${fparse -H}
ymax = 0
nx = 10
ny = 5
bias_y = ${fparse 1/0.7}
[]
[smg]
type = StitchMeshGenerator
inputs = 'block_1 block_2'
clear_stitched_boundary_ids = true
stitch_boundaries_pairs = 'bottom top'
merge_boundaries_with_same_name = true
[]
# Prevent test diffing on distributed parallel element numbering
allow_renumbering = false
[]
[Problem]
linear_sys_names = 'u_system v_system pressure_system energy_system TKE_system TKED_system'
previous_nl_solution_required = true
[]
[GlobalParams]
rhie_chow_user_object = 'rc'
[]
[UserObjects]
[rc]
type = RhieChowMassFlux
u = vel_x
v = vel_y
pressure = pressure
rho = 'rho_functor'
p_diffusion_kernel = p_diffusion
pressure_projection_method = 'consistent'
[]
[]
[Variables]
[vel_x]
type = MooseLinearVariableFVReal
initial_condition = ${bulk_u}
solver_sys = u_system
[]
[vel_y]
type = MooseLinearVariableFVReal
initial_condition = 0
solver_sys = v_system
[]
[pressure]
type = MooseLinearVariableFVReal
initial_condition = 1e-8
solver_sys = pressure_system
[]
[TKE]
type = MooseLinearVariableFVReal
solver_sys = TKE_system
initial_condition = ${k_init}
[]
[TKED]
type = MooseLinearVariableFVReal
solver_sys = TKED_system
initial_condition = ${eps_init}
[]
[T]
type = MooseLinearVariableFVReal
solver_sys = energy_system
initial_condition = 300.0
[]
[]
[FVInterpolationMethods]
[upwind]
type = FVAdvectedUpwind
[]
[]
[LinearFVKernels]
[u_advection_stress]
type = LinearWCNSFVMomentumFlux
variable = vel_x
advected_interp_method_name = ${advected_interp_method}
mu = 'mu_t'
u = vel_x
v = vel_y
momentum_component = 'x'
rhie_chow_user_object = 'rc'
use_nonorthogonal_correction = false
use_deviatoric_terms = yes
[]
[u_diffusion]
type = LinearFVDiffusion
variable = vel_x
diffusion_coeff = '${mu}'
[]
[u_pressure]
type = LinearFVMomentumPressure
variable = vel_x
pressure = pressure
momentum_component = 'x'
[]
[u_buoyancy]
type = LinearFVMomentumBuoyancy
variable = vel_x
rho = 'rho_functor'
reference_rho = ${rho_0}
gravity = '-9.81 0 0'
momentum_component = 'x'
[]
[v_advection_stress]
type = LinearWCNSFVMomentumFlux
variable = vel_y
advected_interp_method_name = ${advected_interp_method}
mu = 'mu_t'
u = vel_x
v = vel_y
momentum_component = 'y'
rhie_chow_user_object = 'rc'
use_nonorthogonal_correction = false
use_deviatoric_terms = yes
[]
[v_diffusion]
type = LinearFVDiffusion
variable = vel_y
diffusion_coeff = '${mu}'
[]
[v_pressure]
type = LinearFVMomentumPressure
variable = vel_y
pressure = pressure
momentum_component = 'y'
[]
[v_buoyancy]
type = LinearFVMomentumBuoyancy
variable = vel_y
rho = 'rho_functor'
reference_rho = ${rho_0}
gravity = '-9.81 0 0'
momentum_component = 'y'
[]
[p_diffusion]
type = LinearFVPressureCorrectionDiffusion
variable = pressure
diffusion_tensor = Ainv
use_nonorthogonal_correction = false
[]
[HbyA_divergence]
type = LinearFVDivergence
variable = pressure
face_flux = HbyA
force_boundary_execution = true
[]
[TKE_advection]
type = LinearFVTurbulentAdvection
variable = TKE
advected_interp_method_name = ${advected_interp_method}
[]
[TKE_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKE
diffusion_coeff = ${mu}
use_nonorthogonal_correction = false
[]
[TKE_turb_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKE
diffusion_coeff = 'mu_t'
scaling_coeff = ${sigma_k}
use_nonorthogonal_correction = false
[]
[TKE_source_sink]
type = kEpsilonTKESourceSink
variable = TKE
u = vel_x
v = vel_y
epsilon = TKED
rho = 'rho_functor'
mu = ${mu}
mu_t = 'mu_t'
walls = ${walls}
wall_treatment = ${wall_treatment}
C_pl = 1e10
# NEW (optional)
k_epsilon_variant = ${k_epsilon_variant} # 'Standard', 'Realizable', etc.
use_buoyancy = ${use_buoyancy}
use_compressibility = ${use_compressibility}
nonlinear_model = ${nonlinear_model}
curvature_model = ${curvature_model}
Pr_t = 0.9
C_M = 1.0
gravity = '-9.81 0 0'
# if/when you have these fields:
temperature = T
beta = ${alpha_b}
speed_of_sound = ${c}
# nonlinear_production = Gnl
# curvature_factor = fc
[]
[TKED_advection]
type = LinearFVTurbulentAdvection
variable = TKED
advected_interp_method_name = ${advected_interp_method}
walls = ${walls}
[]
[TKED_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKED
diffusion_coeff = ${mu}
use_nonorthogonal_correction = false
walls = ${walls}
[]
[TKED_turb_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKED
diffusion_coeff = 'mu_t'
scaling_coeff = ${sigma_eps}
use_nonorthogonal_correction = false
walls = ${walls}
[]
[TKED_source_sink]
type = kEpsilonTKEDSourceSink
variable = TKED
u = vel_x
v = vel_y
tke = TKE
rho = 'rho_functor'
mu = ${mu}
mu_t = 'mu_t'
C1_eps = ${C1_eps}
C2_eps = ${C2_eps}
walls = ${walls}
wall_treatment = ${wall_treatment}
C_pl = 1e10
# NEW (optional)
k_epsilon_variant = ${k_epsilon_variant}
use_buoyancy = ${use_buoyancy}
use_compressibility = ${use_compressibility}
nonlinear_model = ${nonlinear_model}
curvature_model = ${curvature_model}
use_yap = ${use_yap}
use_low_re_Gprime = ${use_low_re_Gprime}
Pr_t = 0.9
C_M = 1.0
gravity = '-9.81 0 0'
# same functors as for TKE if you use them:
temperature = T
beta = ${alpha_b}
speed_of_sound = ${c}
# nonlinear_production = Gnl
# curvature_factor = fc
wall_distance = wall_distance # for low-Re / two-layer Yap / G' terms
[]
[heat_advection]
type = LinearFVEnergyAdvection
variable = T
advected_quantity = temperature
cp = ${cp}
[]
[conduction]
type = LinearFVDiffusion
variable = T
diffusion_coeff = ${k}
[]
[source]
type = LinearFVSource
variable = T
source_density = 1e3
[]
[]
[LinearFVBCs]
[inlet-u]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'left'
variable = vel_x
functor = '${bulk_u}'
[]
[inlet-v]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'left'
variable = vel_y
functor = '0.0'
[]
[walls-u]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'top bottom'
variable = vel_x
functor = 0.0
[]
[walls-v]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'top bottom'
variable = vel_y
functor = 0.0
[]
[outlet_u]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = 'right'
variable = vel_x
use_two_term_expansion = false
[]
[outlet_v]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = 'right'
variable = vel_y
use_two_term_expansion = false
[]
[outlet_p]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'right'
variable = pressure
functor = 0.0
[]
[inlet_TKE]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'left'
variable = TKE
functor = '${k_init}'
[]
[outlet_TKE]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = 'right'
variable = TKE
use_two_term_expansion = false
[]
[inlet_TKED]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'left'
variable = TKED
functor = '${eps_init}'
[]
[outlet_TKED]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = 'right'
variable = TKED
use_two_term_expansion = false
[]
[walls_mu_t]
type = LinearFVTurbulentViscosityWallFunctionBC
boundary = 'bottom top'
variable = 'mu_t'
u = vel_x
v = vel_y
rho = 'rho_functor'
mu = ${mu}
tke = TKE
wall_treatment = ${wall_treatment}
[]
[inlet_and_wall_T]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
variable = T
functor = 300.0
boundary = 'left top bottom'
[]
[outlet_T]
type = LinearFVAdvectionDiffusionOutflowBC
variable = T
use_two_term_expansion = false
boundary = 'right'
[]
[]
[AuxVariables]
[wall_distance]
type = MooseVariableFVReal
initial_condition = 1.0
[]
[mu_t]
type = MooseLinearVariableFVReal
initial_condition = '${fparse rho_0 * C_mu * ${k_init}^2 / eps_init}'
[]
[yplus]
type = MooseVariableFVReal
two_term_boundary_expansion = false
[]
[]
[AuxKernels]
[compute_wall_distance]
type = WallDistanceAux
variable = wall_distance
walls = ${walls}
execute_on = 'INITIAL NONLINEAR'
[]
[compute_mu_t]
type = kEpsilonViscosity
variable = mu_t
C_mu = ${C_mu}
tke = TKE
epsilon = TKED
mu = ${mu}
rho = 'rho_functor'
u = vel_x
v = vel_y
bulk_wall_treatment = ${bulk_wall_treatment}
walls = ${walls}
wall_treatment = ${wall_treatment}
mu_t_ratio_max = 1e20
execute_on = 'NONLINEAR'
# NEW (optional) – choose model and options
k_epsilon_variant = ${k_epsilon_variant} # e.g. 'Standard' or 'Realizable'
two_layer_flavor = ${two_layer_flavor} # ignored unless *TwoLayer variants
wall_distance = wall_distance # only needed for LowRe/TwoLayer (see below)
[]
[compute_y_plus]
type = RANSYPlusAux
variable = yplus
tke = TKE
mu = ${mu}
rho = 'rho_functor'
u = vel_x
v = vel_y
walls = ${walls}
wall_treatment = ${wall_treatment}
execute_on = 'NONLINEAR'
[]
[]
[FunctorMaterials]
[rho_function]
type = ParsedFunctorMaterial
property_name = 'rho_functor'
functor_names = 'T'
expression = '${rho_0}*(1-${alpha_b}*(T-${T_0})) '
[]
[]
[Executioner]
type = SIMPLE
rhie_chow_user_object = 'rc'
momentum_systems = 'u_system v_system'
pressure_system = 'pressure_system'
turbulence_systems = 'TKE_system TKED_system'
energy_system = 'energy_system'
momentum_l_abs_tol = 1e-14
pressure_l_abs_tol = 1e-14
turbulence_l_abs_tol = 1e-14
energy_l_abs_tol = 1e-14
momentum_l_tol = 1e-14
pressure_l_tol = 1e-14
turbulence_l_tol = 1e-14
energy_l_tol = 1e-14
momentum_equation_relaxation = 0.7
pressure_variable_relaxation = 0.3
turbulence_equation_relaxation = '0.2 0.2'
turbulence_field_relaxation = '0.2 0.2'
energy_equation_relaxation = 0.9
num_iterations = 1000
pressure_absolute_tolerance = 1e-7
momentum_absolute_tolerance = 1e-7
turbulence_absolute_tolerance = '1e-7 1e-7'
energy_absolute_tolerance = 1e-7
momentum_petsc_options_iname = '-pc_type -pc_hypre_type'
momentum_petsc_options_value = 'hypre boomeramg'
pressure_petsc_options_iname = '-pc_type -pc_hypre_type'
pressure_petsc_options_value = 'hypre boomeramg'
turbulence_petsc_options_iname = '-pc_type -pc_hypre_type'
turbulence_petsc_options_value = 'hypre boomeramg'
momentum_l_max_its = 300
pressure_l_max_its = 300
turbulence_l_max_its = 30
energy_l_max_its = 300
print_fields = false
continue_on_max_its = true
[]
[Outputs]
exodus = false
[csv]
type = CSV
execute_on = FINAL
[]
[]
variables_to_sample = 'vel_x vel_y pressure T TKE TKED'
[VectorPostprocessors]
[side_bottom]
type = SideValueSampler
boundary = 'bottom'
variable = ${variables_to_sample}
sort_by = 'x'
execute_on = 'timestep_end'
[]
[side_top]
type = SideValueSampler
boundary = 'top'
variable = ${variables_to_sample}
sort_by = 'x'
execute_on = 'timestep_end'
[]
[line_center_channel]
type = LineValueSampler
start_point = '${fparse 0.125 * L} ${fparse 0.0001} 0'
end_point = '${fparse 0.875 * L} ${fparse 0.0001} 0'
num_points = ${Mesh/block_1/nx}
variable = ${variables_to_sample}
sort_by = 'x'
execute_on = 'timestep_end'
[]
[line_quarter_radius_channel]
type = LineValueSampler
start_point = '${fparse 0.125 * L} ${fparse 0.5 * H} 0'
end_point = '${fparse 0.875 * L} ${fparse 0.5 * H} 0'
num_points = ${Mesh/block_1/nx}
variable = ${variables_to_sample}
sort_by = 'x'
execute_on = 'timestep_end'
[]
[]
(test/tests/kEpsilon/lid-driven-cavity/lid-driven.i)
### Thermophysical Properties ###
mu = 1e-3
rho = 1.0
### Operation Conditions ###
lid_velocity = 1.0
side_length = 0.1
### Initial Conditions ###
intensity = 0.01
k_init = '${fparse 1.5*(intensity * lid_velocity)^2}'
eps_init = '${fparse C_mu^0.75 * k_init^1.5 / side_length}'
### k-epsilon Closure Parameters ###
sigma_k = 1.0
sigma_eps = 1.3
C1_eps = 1.44
C2_eps = 1.92
C_mu = 0.09
### Modeling parameters ###
bulk_wall_treatment = false
walls = 'left top right bottom'
wall_treatment = 'neq' # Options: eq_newton, eq_incremental, eq_linearized, neq
[GlobalParams]
rhie_chow_user_object = 'rc'
[]
[Mesh]
[gen]
type = GeneratedMeshGenerator
dim = 2
xmin = 0
xmax = ${side_length}
ymin = 0
ymax = ${side_length}
nx = 12
ny = 12
[]
# Prevent test diffing on distributed parallel element numbering
allow_renumbering = false
[]
[Problem]
linear_sys_names = 'u_system v_system pressure_system TKE_system TKED_system'
previous_nl_solution_required = true
[]
[UserObjects]
[rc]
type = RhieChowMassFlux
u = vel_x
v = vel_y
pressure = pressure
rho = ${rho}
p_diffusion_kernel = p_diffusion
[]
[]
[Variables]
[vel_x]
type = MooseLinearVariableFVReal
initial_condition = ${lid_velocity}
solver_sys = u_system
[]
[vel_y]
type = MooseLinearVariableFVReal
initial_condition = 0
solver_sys = v_system
[]
[pressure]
type = MooseLinearVariableFVReal
initial_condition = 1e-8
solver_sys = pressure_system
[]
[TKE]
type = MooseLinearVariableFVReal
solver_sys = TKE_system
initial_condition = ${k_init}
[]
[TKED]
type = MooseLinearVariableFVReal
solver_sys = TKED_system
initial_condition = ${eps_init}
[]
[]
[FVInterpolationMethods]
[upwind]
type = FVAdvectedUpwind
[]
[]
[LinearFVKernels]
[u_advection_stress]
type = LinearWCNSFVMomentumFlux
variable = vel_x
advected_interp_method_name = upwind
mu = 'mu_t'
u = vel_x
v = vel_y
momentum_component = 'x'
rhie_chow_user_object = 'rc'
use_nonorthogonal_correction = false
use_deviatoric_terms = yes
[]
[u_diffusion]
type = LinearFVDiffusion
variable = vel_x
diffusion_coeff = ${mu}
use_nonorthogonal_correction = false
[]
[u_pressure]
type = LinearFVMomentumPressure
variable = vel_x
pressure = pressure
momentum_component = 'x'
[]
[v_advection_stress]
type = LinearWCNSFVMomentumFlux
variable = vel_y
advected_interp_method_name = upwind
mu = 'mu_t'
u = vel_x
v = vel_y
momentum_component = 'y'
rhie_chow_user_object = 'rc'
use_nonorthogonal_correction = false
use_deviatoric_terms = yes
[]
[v_diffusion]
type = LinearFVDiffusion
variable = vel_y
diffusion_coeff = ${mu}
use_nonorthogonal_correction = false
[]
[v_pressure]
type = LinearFVMomentumPressure
variable = vel_y
pressure = pressure
momentum_component = 'y'
[]
[p_diffusion]
type = LinearFVPressureCorrectionDiffusion
variable = pressure
diffusion_tensor = Ainv
use_nonorthogonal_correction = false
[]
[HbyA_divergence]
type = LinearFVDivergence
variable = pressure
face_flux = HbyA
force_boundary_execution = true
[]
[TKE_advection]
type = LinearFVTurbulentAdvection
variable = TKE
advected_interp_method_name = upwind
[]
[TKE_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKE
diffusion_coeff = ${mu}
use_nonorthogonal_correction = false
[]
[TKE_turb_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKE
diffusion_coeff = 'mu_t'
scaling_coeff = ${sigma_k}
use_nonorthogonal_correction = false
[]
[TKE_source_sink]
type = kEpsilonTKESourceSink
variable = TKE
u = vel_x
v = vel_y
epsilon = TKED
rho = ${rho}
mu = ${mu}
mu_t = 'mu_t'
walls = ${walls}
wall_treatment = ${wall_treatment}
C_pl = 1e10
# if/when you have these fields:
# temperature = T
# beta = beta
# speed_of_sound = c
# nonlinear_production = Gnl
# curvature_factor = fc
[]
[TKED_advection]
type = LinearFVTurbulentAdvection
variable = TKED
advected_interp_method_name = upwind
walls = ${walls}
[]
[TKED_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKED
diffusion_coeff = ${mu}
use_nonorthogonal_correction = false
walls = ${walls}
[]
[TKED_turb_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKED
diffusion_coeff = 'mu_t'
scaling_coeff = ${sigma_eps}
use_nonorthogonal_correction = false
walls = ${walls}
[]
[TKED_source_sink]
type = kEpsilonTKEDSourceSink
variable = TKED
u = vel_x
v = vel_y
tke = TKE
rho = ${rho}
mu = ${mu}
mu_t = 'mu_t'
C1_eps = ${C1_eps}
C2_eps = ${C2_eps}
walls = ${walls}
wall_treatment = ${wall_treatment}
# C_pl = 1e10
# same functors as for TKE if you use them:
# temperature = T
# beta = beta
# speed_of_sound = c
# nonlinear_production = Gnl
# curvature_factor = fc
wall_distance = wall_distance # for low-Re / two-layer Yap / G' terms
[]
[]
[LinearFVBCs]
[top_x]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
variable = vel_x
boundary = 'top'
functor = 1
[]
[no_slip_x]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
variable = vel_x
boundary = 'left right bottom'
functor = 0
[]
[no_slip_y]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
variable = vel_y
boundary = 'left right top bottom'
functor = 0
[]
[pressure-extrapolation]
type = LinearFVExtrapolatedPressureBC
boundary = 'left right top bottom'
variable = pressure
use_two_term_expansion = true
[]
[walls_mu_t]
type = LinearFVTurbulentViscosityWallFunctionBC
boundary = 'bottom top'
variable = 'mu_t'
u = vel_x
v = vel_y
rho = ${rho}
mu = ${mu}
tke = TKE
wall_treatment = ${wall_treatment}
[]
[]
[AuxVariables]
[wall_distance]
type = MooseVariableFVReal
initial_condition = 1.0
[]
[mu_t]
type = MooseLinearVariableFVReal
initial_condition = '${fparse rho * C_mu * ${k_init}^2 / eps_init}'
[]
[yplus]
type = MooseLinearVariableFVReal
[]
[mu_eff]
type = MooseLinearVariableFVReal
initial_condition = '${fparse rho * C_mu * ${k_init}^2 / eps_init}'
[]
[]
[AuxKernels]
[compute_wall_distance]
type = WallDistanceAux
variable = wall_distance
walls = ${walls}
execute_on = 'INITIAL NONLINEAR'
[]
[compute_mu_t]
type = kEpsilonViscosity
variable = mu_t
C_mu = ${C_mu}
tke = TKE
epsilon = TKED
mu = ${mu}
rho = ${rho}
u = vel_x
v = vel_y
bulk_wall_treatment = ${bulk_wall_treatment}
walls = ${walls}
wall_treatment = ${wall_treatment}
mu_t_ratio_max = 1e20
execute_on = 'NONLINEAR'
wall_distance = wall_distance # only needed for LowRe/TwoLayer (see below)
[]
[compute_y_plus]
type = RANSYPlusAux
variable = yplus
tke = TKE
mu = ${mu}
rho = ${rho}
u = vel_x
v = vel_y
walls = ${walls}
wall_treatment = ${wall_treatment}
execute_on = 'NONLINEAR'
[]
[compute_mu_eff]
type = ParsedAux
variable = 'mu_eff'
coupled_variables = 'mu_t'
expression = 'mu_t + ${mu}'
execute_on = 'NONLINEAR'
[]
[]
[Executioner]
type = SIMPLE
rhie_chow_user_object = 'rc'
momentum_systems = 'u_system v_system'
pressure_system = 'pressure_system'
turbulence_systems = 'TKE_system TKED_system'
momentum_l_abs_tol = 1e-14
pressure_l_abs_tol = 1e-14
turbulence_l_abs_tol = 1e-14
momentum_l_tol = 1e-14
pressure_l_tol = 1e-14
turbulence_l_tol = 1e-14
momentum_equation_relaxation = 0.7
pressure_variable_relaxation = 0.3
turbulence_equation_relaxation = '0.5 0.5'
num_iterations = 500
pressure_absolute_tolerance = 1e-10
momentum_absolute_tolerance = 1e-10
turbulence_absolute_tolerance = '1e-10 1e-10'
momentum_petsc_options_iname = '-pc_type -pc_hypre_type'
momentum_petsc_options_value = 'hypre boomeramg'
pressure_petsc_options_iname = '-pc_type -pc_hypre_type'
pressure_petsc_options_value = 'hypre boomeramg'
turbulence_petsc_options_iname = '-pc_type -pc_hypre_type'
turbulence_petsc_options_value = 'hypre boomeramg'
print_fields = false
continue_on_max_its = true
pin_pressure = true
pressure_pin_value = 0.0
pressure_pin_point = '0.01 0.099 0.0'
[]
[Outputs]
exodus = false
[csv]
type = CSV
execute_on = FINAL
[]
[]
[VectorPostprocessors]
[side_bottom]
type = SideValueSampler
boundary = 'bottom'
variable = 'vel_x vel_y pressure TKE TKED'
sort_by = 'x'
execute_on = 'timestep_end'
[]
[side_top]
type = SideValueSampler
boundary = 'top'
variable = 'vel_x vel_y pressure TKE TKED'
sort_by = 'x'
execute_on = 'timestep_end'
[]
[side_left]
type = SideValueSampler
boundary = 'left'
variable = 'vel_x vel_y pressure TKE TKED'
sort_by = 'y'
execute_on = 'timestep_end'
[]
[side_right]
type = SideValueSampler
boundary = 'right'
variable = 'vel_x vel_y pressure TKE TKED'
sort_by = 'y'
execute_on = 'timestep_end'
[]
[horizontal_center]
type = LineValueSampler
start_point = '${fparse 0.01 * side_length} ${fparse 0.499 * side_length} 0'
end_point = '${fparse 0.99 * side_length} ${fparse 0.499 * side_length} 0'
num_points = ${Mesh/gen/nx}
variable = 'vel_x vel_y pressure TKE TKED'
sort_by = 'x'
execute_on = 'timestep_end'
[]
[vertical_center]
type = LineValueSampler
start_point = '${fparse 0.499 * side_length} ${fparse 0.01 * side_length} 0'
end_point = '${fparse 0.499 * side_length} ${fparse 0.99 * side_length} 0'
num_points = ${Mesh/gen/ny}
variable = 'vel_x vel_y pressure TKE TKED'
sort_by = 'y'
execute_on = 'timestep_end'
[]
[]
(test/tests/kEpsilon/special-cases/lid-driven-curvature-corrected.i)
### Thermophysical Properties ###
mu = 1e-3
rho = 1.0
### Operation Conditions ###
lid_velocity = 1.0
side_length = 0.1
### Initial Conditions ###
intensity = 0.01
k_init = '${fparse 1.5*(intensity * lid_velocity)^2}'
eps_init = '${fparse C_mu^0.75 * k_init^1.5 / side_length}'
### k-epsilon Closure Parameters ###
sigma_k = 1.0
sigma_eps = 1.3
C1_eps = 1.44
C2_eps = 1.92
C_mu = 0.09
### Modeling parameters ###
bulk_wall_treatment = false
walls = 'left top right bottom'
wall_treatment = 'neq' # Options: eq_newton, eq_incremental, eq_linearized, neq
# Turbulence-model knobs (optional)
k_epsilon_variant = 'Realizable' # Standard | StandardLowRe | StandardTwoLayer | Realizable | RealizableTwoLayer
use_buoyancy = false
use_compressibility = false
nonlinear_model = 'none'
curvature_model = 'standard'
use_yap = false
use_low_re_Gprime = false
[GlobalParams]
rhie_chow_user_object = 'rc'
[]
[Mesh]
[gen]
type = GeneratedMeshGenerator
dim = 2
xmin = 0
xmax = ${side_length}
ymin = 0
ymax = ${side_length}
nx = 8
ny = 8
[]
# Prevent test diffing on distributed parallel element numbering
allow_renumbering = false
[]
[Problem]
linear_sys_names = 'u_system v_system pressure_system TKE_system TKED_system'
previous_nl_solution_required = true
[]
[UserObjects]
[rc]
type = RhieChowMassFlux
u = vel_x
v = vel_y
pressure = pressure
rho = ${rho}
p_diffusion_kernel = p_diffusion
[]
[]
[Variables]
[vel_x]
type = MooseLinearVariableFVReal
initial_condition = ${lid_velocity}
solver_sys = u_system
[]
[vel_y]
type = MooseLinearVariableFVReal
initial_condition = 0
solver_sys = v_system
[]
[pressure]
type = MooseLinearVariableFVReal
initial_condition = 1e-8
solver_sys = pressure_system
[]
[TKE]
type = MooseLinearVariableFVReal
solver_sys = TKE_system
initial_condition = ${k_init}
[]
[TKED]
type = MooseLinearVariableFVReal
solver_sys = TKED_system
initial_condition = ${eps_init}
[]
[]
[FVInterpolationMethods]
[upwind]
type = FVAdvectedUpwind
[]
[]
[LinearFVKernels]
[u_advection_stress]
type = LinearWCNSFVMomentumFlux
variable = vel_x
advected_interp_method_name = upwind
mu = 'mu_t'
u = vel_x
v = vel_y
momentum_component = 'x'
rhie_chow_user_object = 'rc'
use_nonorthogonal_correction = false
use_deviatoric_terms = yes
[]
[u_diffusion]
type = LinearFVDiffusion
variable = vel_x
diffusion_coeff = ${mu}
use_nonorthogonal_correction = false
[]
[u_pressure]
type = LinearFVMomentumPressure
variable = vel_x
pressure = pressure
momentum_component = 'x'
[]
[v_advection_stress]
type = LinearWCNSFVMomentumFlux
variable = vel_y
advected_interp_method_name = upwind
mu = 'mu_t'
u = vel_x
v = vel_y
momentum_component = 'y'
rhie_chow_user_object = 'rc'
use_nonorthogonal_correction = false
use_deviatoric_terms = yes
[]
[v_diffusion]
type = LinearFVDiffusion
variable = vel_y
diffusion_coeff = ${mu}
use_nonorthogonal_correction = false
[]
[v_pressure]
type = LinearFVMomentumPressure
variable = vel_y
pressure = pressure
momentum_component = 'y'
[]
[p_diffusion]
type = LinearFVPressureCorrectionDiffusion
variable = pressure
diffusion_tensor = Ainv
use_nonorthogonal_correction = false
[]
[HbyA_divergence]
type = LinearFVDivergence
variable = pressure
face_flux = HbyA
force_boundary_execution = true
[]
[TKE_advection]
type = LinearFVTurbulentAdvection
variable = TKE
advected_interp_method_name = upwind
[]
[TKE_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKE
diffusion_coeff = ${mu}
use_nonorthogonal_correction = false
[]
[TKE_turb_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKE
diffusion_coeff = 'mu_t'
scaling_coeff = ${sigma_k}
use_nonorthogonal_correction = false
[]
[TKE_source_sink]
type = kEpsilonTKESourceSink
variable = TKE
u = vel_x
v = vel_y
epsilon = TKED
rho = ${rho}
mu = ${mu}
mu_t = 'mu_t'
walls = ${walls}
wall_treatment = ${wall_treatment}
C_pl = 1e10
# NEW (optional)
k_epsilon_variant = ${k_epsilon_variant} # 'Standard', 'Realizable', etc.
use_buoyancy = ${use_buoyancy}
use_compressibility = ${use_compressibility}
nonlinear_model = ${nonlinear_model}
use_curvature_correction = true
curvature_model = ${curvature_model}
[]
[TKED_advection]
type = LinearFVTurbulentAdvection
variable = TKED
advected_interp_method_name = upwind
walls = ${walls}
[]
[TKED_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKED
diffusion_coeff = ${mu}
use_nonorthogonal_correction = false
walls = ${walls}
[]
[TKED_turb_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKED
diffusion_coeff = 'mu_t'
scaling_coeff = ${sigma_eps}
use_nonorthogonal_correction = false
walls = ${walls}
[]
[TKED_source_sink]
type = kEpsilonTKEDSourceSink
variable = TKED
u = vel_x
v = vel_y
tke = TKE
rho = ${rho}
mu = ${mu}
mu_t = 'mu_t'
C1_eps = ${C1_eps}
C2_eps = ${C2_eps}
walls = ${walls}
wall_treatment = ${wall_treatment}
# C_pl = 1e10
# NEW (optional)
k_epsilon_variant = ${k_epsilon_variant}
use_buoyancy = ${use_buoyancy}
use_compressibility = ${use_compressibility}
nonlinear_model = ${nonlinear_model}
use_curvature_correction = true
curvature_model = ${curvature_model}
use_yap = ${use_yap}
use_low_re_Gprime = ${use_low_re_Gprime}
wall_distance = wall_distance # for low-Re / two-layer Yap / G' terms
[]
[]
[LinearFVBCs]
[top_x]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
variable = vel_x
boundary = 'top'
functor = 1
[]
[no_slip_x]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
variable = vel_x
boundary = 'left right bottom'
functor = 0
[]
[no_slip_y]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
variable = vel_y
boundary = 'left right top bottom'
functor = 0
[]
[pressure-extrapolation]
type = LinearFVExtrapolatedPressureBC
boundary = 'left right top bottom'
variable = pressure
use_two_term_expansion = true
[]
[walls_mu_t]
type = LinearFVTurbulentViscosityWallFunctionBC
boundary = 'bottom top'
variable = 'mu_t'
u = vel_x
v = vel_y
rho = ${rho}
mu = ${mu}
tke = TKE
wall_treatment = ${wall_treatment}
[]
[]
[AuxVariables]
[wall_distance]
type = MooseVariableFVReal
initial_condition = 1.0
[]
[mu_t]
type = MooseLinearVariableFVReal
initial_condition = '${fparse rho * C_mu * ${k_init}^2 / eps_init}'
[]
[yplus]
type = MooseLinearVariableFVReal
[]
[mu_eff]
type = MooseLinearVariableFVReal
initial_condition = '${fparse rho * C_mu * ${k_init}^2 / eps_init}'
[]
[]
[AuxKernels]
[compute_wall_distance]
type = WallDistanceAux
variable = wall_distance
walls = ${walls}
execute_on = 'INITIAL NONLINEAR'
[]
[compute_mu_t]
type = kEpsilonViscosity
variable = mu_t
C_mu = ${C_mu}
tke = TKE
epsilon = TKED
mu = ${mu}
rho = ${rho}
u = vel_x
v = vel_y
bulk_wall_treatment = ${bulk_wall_treatment}
walls = ${walls}
wall_treatment = ${wall_treatment}
mu_t_ratio_max = 1e20
execute_on = 'NONLINEAR'
# NEW (optional) – choose model and options
k_epsilon_variant = ${k_epsilon_variant} # e.g. 'Standard' or 'Realizable'
wall_distance = wall_distance # only needed for LowRe/TwoLayer (see below)
[]
[compute_y_plus]
type = RANSYPlusAux
variable = yplus
tke = TKE
mu = ${mu}
rho = ${rho}
u = vel_x
v = vel_y
walls = ${walls}
wall_treatment = ${wall_treatment}
execute_on = 'NONLINEAR'
[]
[compute_mu_eff]
type = ParsedAux
variable = 'mu_eff'
coupled_variables = 'mu_t'
expression = 'mu_t + ${mu}'
execute_on = 'NONLINEAR'
[]
[]
[Executioner]
type = SIMPLE
rhie_chow_user_object = 'rc'
momentum_systems = 'u_system v_system'
pressure_system = 'pressure_system'
turbulence_systems = 'TKE_system TKED_system'
momentum_l_abs_tol = 1e-10
pressure_l_abs_tol = 1e-10
turbulence_l_abs_tol = 1e-10
momentum_l_tol = 1e-10
pressure_l_tol = 1e-10
turbulence_l_tol = 1e-10
momentum_equation_relaxation = 0.7
pressure_variable_relaxation = 0.3
turbulence_equation_relaxation = '0.5 0.5'
num_iterations = 500
pressure_absolute_tolerance = 1e-10
momentum_absolute_tolerance = 1e-10
turbulence_absolute_tolerance = '1e-10 1e-10'
momentum_petsc_options_iname = '-pc_type -pc_hypre_type'
momentum_petsc_options_value = 'hypre boomeramg'
pressure_petsc_options_iname = '-pc_type -pc_hypre_type'
pressure_petsc_options_value = 'hypre boomeramg'
turbulence_petsc_options_iname = '-pc_type -pc_hypre_type'
turbulence_petsc_options_value = 'hypre boomeramg'
print_fields = false
continue_on_max_its = true
pin_pressure = true
pressure_pin_value = 0.0
pressure_pin_point = '0.01 0.099 0.0'
[]
[Outputs]
exodus = false
[csv]
type = CSV
execute_on = FINAL
[]
[]
[VectorPostprocessors]
[side_bottom]
type = SideValueSampler
boundary = 'bottom'
variable = 'vel_x vel_y pressure TKE TKED'
sort_by = 'x'
execute_on = 'timestep_end'
[]
[side_top]
type = SideValueSampler
boundary = 'top'
variable = 'vel_x vel_y pressure TKE TKED'
sort_by = 'x'
execute_on = 'timestep_end'
[]
[side_left]
type = SideValueSampler
boundary = 'left'
variable = 'vel_x vel_y pressure TKE TKED'
sort_by = 'y'
execute_on = 'timestep_end'
[]
[side_right]
type = SideValueSampler
boundary = 'right'
variable = 'vel_x vel_y pressure TKE TKED'
sort_by = 'y'
execute_on = 'timestep_end'
[]
[horizontal_center]
type = LineValueSampler
start_point = '${fparse 0.01 * side_length} ${fparse 0.499 * side_length} 0'
end_point = '${fparse 0.99 * side_length} ${fparse 0.499 * side_length} 0'
num_points = ${Mesh/gen/nx}
variable = 'vel_x vel_y pressure TKE TKED'
sort_by = 'x'
execute_on = 'timestep_end'
[]
[vertical_center]
type = LineValueSampler
start_point = '${fparse 0.499 * side_length} ${fparse 0.01 * side_length} 0'
end_point = '${fparse 0.499 * side_length} ${fparse 0.99 * side_length} 0'
num_points = ${Mesh/gen/ny}
variable = 'vel_x vel_y pressure TKE TKED'
sort_by = 'y'
execute_on = 'timestep_end'
[]
[]
(validation/free_flow/isothermal/single-stream-jet/RTL10000/jet_RTL10000.i)
# ================================================================================================================
# Model description
# ================================================================================================================
# Open-Pronghorn simulation of a steady state incompressible turbulent flow of a jet into a big space/tank.
# Reference: https://pubs.aip.org/aip/pof/article/36/11/115133/3319195/Calibration-of-the-Reynolds-stress-model-for
# Author: Dr. Kyriakopoulos Vasileios, Dr. Mauricio Tano
# Idaho Falls, INL, April 10, 2026
# ================================================================================================================
# MODEL PARAMETERS
# ================================================================================================================
# Geometry
# ================================================================================================================
H = 0.2 # Height/width of the tank (m).
L = 0.3 # Total tank length (m)
Dh = 0.01 # Inlet Pipe diameter (m)
nz = 150 # number of cells in the axial direction
n1 = 8 # number of radial cells per ring in the cross-section
n2 = 8 # number of radial cells per ring in the cross-section
n3 = 8 # number of radial cells per ring in the cross-section
n4 = 8 # number of radial cells per ring in the cross-section
n5 = 8 # number of radial cells per ring in the cross-section
ny = ${fparse n1 + n2 + n3 + n4 + n5}
# A = '${fparse pi*Dh*Dh/4}' # Cross-section area of the inlet (m2)
# ================================================================================================================
# Air properties
# ================================================================================================================
rho = 1.184 # Density (kg/m3)
# ================================================================================================================
# Flow and Boundary conditions
# ================================================================================================================
Re = 10000.0 # Reynolds number
mu = 1.844e-5 # Pa*s
Vin = '${fparse Re * mu / (rho * Dh)}' # Inlet velocity (m/s) # Uniform velocity profile Faghani et al. 2011
intensity = 0.05 # Turutoglu et al 2024
k_init = '${fparse 1.5*(intensity * Vin)^2}'
eps_init = '${fparse C_mu^0.75 * k_init^1.5 / (2*0.07*Dh)}'
advected_interp_method = 'upwind'
# ================================================================================================================
# k-epsilon Closure Parameters
# ================================================================================================================
sigma_k = 1.0
sigma_eps = 1.3
C1_eps = 1.44
C2_eps = 1.92
C_mu = 0.09
# The front boundary is the solid inlet-plane area outside the jet opening;
# the remaining outer tank boundaries are pressure outlets.
walls = 'front'
wall_treatment = 'neq' # Options: eq_newton, eq_incremental, eq_linearized, neq
eps = 1e-4
# ================================================================================================================
# Modeling parameters: Turbulence-model knobs (optional)
# ================================================================================================================
k_epsilon_variant = 'RealizableTwoLayer' # Standard | StandardLowRe | StandardTwoLayer | Realizable | RealizableTwoLayer
two_layer_flavor = 'Wolfstein' # Wolfstein | NorrisReynolds | Xu (only used for *TwoLayer variants)
use_buoyancy = false
use_compressibility = false
nonlinear_model = 'none'
curvature_model = 'none'
use_yap = false
use_low_re_Gprime = false
bulk_wall_treatment = false
outlet = 'back bottom top left right'
[Mesh]
[front_face]
type = ConcentricCircleMeshGenerator
num_sectors = 8
radii = '${fparse Dh/2.0} ${Dh} ${fparse Dh*2.0} ${fparse Dh*4.0}'
rings = '${n1} ${n2} ${n3} ${n4} ${n5}'
has_outer_square = on
pitch = ${H}
preserve_volumes = on
smoothing_max_it = 3
[]
[extrude]
type = AdvancedExtruderGenerator
direction = '0 0 1'
input = front_face
heights = '${L}'
num_layers = ${nz}
[]
[inlet]
type = ParsedGenerateSideset
input = extrude
combinatorial_geometry = 'abs(z) < 1e-6'
included_subdomains = '1'
normal = '0 0 -1'
new_sideset_name = inlet
[]
[front]
type = ParsedGenerateSideset
input = inlet
combinatorial_geometry = 'abs(z) < 1e-6'
included_subdomains = '2 3 4 5'
normal = '0 0 -1'
new_sideset_name = front
[]
[back]
type = ParsedGenerateSideset
input = front
combinatorial_geometry = 'abs(z) > ${fparse L - 1e-6}'
included_subdomains = '1 2 3 4 5'
normal = '0 0 1'
new_sideset_name = back
[]
[delete_surfaces]
type = BoundaryDeletionGenerator
input = back
operation = remove
boundary_names = '5 6'
[]
[rename]
type = RenameBlockGenerator
input = delete_surfaces
old_block = '1 2 3 4 5'
new_block = 'air air air air air'
[]
[]
[Problem]
linear_sys_names = 'u_system v_system w_system pressure_system TKE_system TKED_system'
previous_nl_solution_required = true
[]
[GlobalParams]
rhie_chow_user_object = 'rc'
[]
[UserObjects]
[rc]
type = RhieChowMassFlux
u = vel_x
v = vel_y
w = vel_z
pressure = pressure
rho = ${rho}
p_diffusion_kernel = p_diffusion
pressure_projection_method = 'consistent'
[]
[]
[Variables]
[vel_x]
type = MooseLinearVariableFVReal
initial_condition = 0
solver_sys = u_system
[]
[vel_y]
type = MooseLinearVariableFVReal
initial_condition = 0
solver_sys = v_system
[]
[vel_z]
type = MooseLinearVariableFVReal
initial_condition = 0
solver_sys = w_system
[]
[pressure]
type = MooseLinearVariableFVReal
initial_condition = 1e-8
solver_sys = pressure_system
[]
[TKE]
type = MooseLinearVariableFVReal
solver_sys = TKE_system
initial_condition = ${k_init}
[]
[TKED]
type = MooseLinearVariableFVReal
solver_sys = TKED_system
initial_condition = ${eps_init}
[]
[]
[FVInterpolationMethods]
[upwind]
type = FVAdvectedUpwind
[]
[]
[LinearFVKernels]
[u_advection_stress]
type = LinearWCNSFVMomentumFlux
variable = vel_x
advected_interp_method_name = ${advected_interp_method}
mu = 'mu_eff'
u = vel_x
v = vel_y
w = vel_z
momentum_component = 'x'
rhie_chow_user_object = 'rc'
use_nonorthogonal_correction = true
use_deviatoric_terms = yes
[]
[u_pressure]
type = LinearFVMomentumPressure
variable = vel_x
pressure = pressure
momentum_component = 'x'
[]
[v_advection_stress]
type = LinearWCNSFVMomentumFlux
variable = vel_y
advected_interp_method_name = ${advected_interp_method}
mu = 'mu_eff'
u = vel_x
v = vel_y
w = vel_z
momentum_component = 'y'
rhie_chow_user_object = 'rc'
use_nonorthogonal_correction = true
use_deviatoric_terms = yes
[]
[v_pressure]
type = LinearFVMomentumPressure
variable = vel_y
pressure = pressure
momentum_component = 'y'
[]
[w_advection_stress]
type = LinearWCNSFVMomentumFlux
variable = vel_z
advected_interp_method_name = ${advected_interp_method}
mu = 'mu_eff'
u = vel_x
v = vel_y
w = vel_z
momentum_component = 'z'
rhie_chow_user_object = 'rc'
use_nonorthogonal_correction = true
use_deviatoric_terms = yes
[]
[w_pressure]
type = LinearFVMomentumPressure
variable = vel_z
pressure = pressure
momentum_component = 'z'
[]
[p_diffusion]
type = LinearFVPressureCorrectionDiffusion
variable = pressure
diffusion_tensor = Ainv
use_nonorthogonal_correction = true
[]
[HbyA_divergence]
type = LinearFVDivergence
variable = pressure
face_flux = HbyA
force_boundary_execution = true
[]
[TKE_advection]
type = LinearFVTurbulentAdvection
variable = TKE
advected_interp_method_name = ${advected_interp_method}
[]
[TKE_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKE
diffusion_coeff = 'mu_eff_k'
use_nonorthogonal_correction = true
[]
[TKE_source_sink]
type = kEpsilonTKESourceSink
variable = TKE
u = vel_x
v = vel_y
w = vel_z
epsilon = TKED
rho = ${rho}
mu = ${mu}
mu_t = 'mu_t'
walls = ${walls}
wall_treatment = ${wall_treatment}
C_pl = 1e10
# optional
k_epsilon_variant = ${k_epsilon_variant} # 'Standard', 'Realizable', etc.
use_buoyancy = ${use_buoyancy}
use_compressibility = ${use_compressibility}
nonlinear_model = ${nonlinear_model}
curvature_model = ${curvature_model}
[]
[TKED_advection]
type = LinearFVTurbulentAdvection
variable = TKED
advected_interp_method_name = ${advected_interp_method}
walls = ${walls}
[]
[TKED_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKED
diffusion_coeff = 'mu_eff_eps'
use_nonorthogonal_correction = true
walls = ${walls}
[]
[TKED_source_sink]
type = kEpsilonTKEDSourceSink
variable = TKED
u = vel_x
v = vel_y
w = vel_z
tke = TKE
rho = ${rho}
mu = ${mu}
mu_t = 'mu_t'
C1_eps = ${C1_eps}
C2_eps = ${C2_eps}
walls = ${walls}
wall_treatment = ${wall_treatment}
# optional
k_epsilon_variant = ${k_epsilon_variant}
use_buoyancy = ${use_buoyancy}
use_compressibility = ${use_compressibility}
nonlinear_model = ${nonlinear_model}
curvature_model = ${curvature_model}
use_yap = ${use_yap}
use_low_re_Gprime = ${use_low_re_Gprime}
wall_distance = wall_distance # for low-Re / two-layer Yap / G' terms
[]
[]
[LinearFVBCs]
[inlet-u]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'inlet'
variable = vel_x
functor = '0.0'
[]
[inlet-v]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'inlet'
variable = vel_y
functor = '0.0'
[]
[inlet-w]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'inlet'
variable = vel_z
functor = ${Vin}
[]
[walls-u]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = ${walls}
variable = vel_x
functor = 0.0
[]
[walls-v]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = ${walls}
variable = vel_y
functor = 0.0
[]
[walls-w]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = ${walls}
variable = vel_z
functor = 0.0
[]
[outlet_u]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = ${outlet}
variable = vel_x
use_two_term_expansion = false
[]
[outlet_v]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = ${outlet}
variable = vel_y
use_two_term_expansion = false
[]
[outlet_w]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = ${outlet}
variable = vel_z
use_two_term_expansion = false
[]
[outlet_p]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = ${outlet}
variable = pressure
functor = 0.0
[]
[inlet_TKE]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'inlet'
variable = TKE
functor = '${k_init}'
[]
[outlet_TKE]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = ${outlet}
variable = TKE
use_two_term_expansion = false
[]
[inlet_TKED]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'inlet'
variable = TKED
functor = '${eps_init}'
[]
[outlet_TKED]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = ${outlet}
variable = TKED
use_two_term_expansion = false
[]
[walls_mu_t]
type = LinearFVTurbulentViscosityWallFunctionBC
boundary = ${walls}
variable = 'mu_t'
u = vel_x
v = vel_y
w = vel_z
rho = ${rho}
mu = ${mu}
tke = TKE
wall_treatment = ${wall_treatment}
[]
[]
[AuxVariables]
[wall_distance]
type = MooseVariableFVReal
initial_condition = 1.0
[]
[mu_t]
type = MooseLinearVariableFVReal
initial_condition = '${fparse rho * C_mu * ${k_init}^2 / eps_init}'
[]
[yplus]
type = MooseVariableFVReal
two_term_boundary_expansion = false
[]
[]
[AuxKernels]
[compute_wall_distance]
type = WallDistanceAux
variable = wall_distance
walls = ${walls}
execute_on = 'INITIAL'
[]
[compute_mu_t]
type = kEpsilonViscosity
variable = mu_t
C_mu = ${C_mu}
tke = TKE
epsilon = TKED
mu = ${mu}
rho = ${rho}
u = vel_x
v = vel_y
w = vel_z
bulk_wall_treatment = ${bulk_wall_treatment}
walls = ${walls}
wall_treatment = ${wall_treatment}
mu_t_ratio_max = 1e20
execute_on = 'NONLINEAR'
# NEW (optional) – choose model and options
k_epsilon_variant = ${k_epsilon_variant} # e.g. 'Standard' or 'Realizable'
two_layer_flavor = ${two_layer_flavor} # ignored unless *TwoLayer variants
# Cd0 = 0.091 # defaults, can omit if you keep Standard
# Cd1 = 0.0042
# Cd2 = 0.00011
Ca0 = 0.667 # Realizable C_mu coefficients
Ca1 = 1.25
Ca2 = 1.0
Ca3 = 0.9
wall_distance = wall_distance # only needed for LowRe/TwoLayer (see below)
[]
[compute_y_plus]
type = RANSYPlusAux
variable = yplus
tke = TKE
mu = ${mu}
rho = ${rho}
u = vel_x
v = vel_y
w = vel_z
walls = ${walls}
wall_treatment = ${wall_treatment}
execute_on = 'NONLINEAR'
[]
[]
[FunctorMaterials]
[compute_mu_eff]
type = FunctorEffectiveDynamicViscosity
property_name = mu_eff
mu = ${mu}
mu_t = mu_t
mu_t_inverse_factor = 1
execute_on = 'ALWAYS'
[]
[compute_mu_eff_k]
type = FunctorEffectiveDynamicViscosity
property_name = mu_eff_k
mu = ${mu}
mu_t = mu_t
mu_t_inverse_factor = ${sigma_k}
execute_on = 'ALWAYS'
[]
[compute_mu_eff_eps]
type = FunctorEffectiveDynamicViscosity
property_name = mu_eff_eps
mu = ${mu}
mu_t = mu_t
mu_t_inverse_factor = ${sigma_eps}
execute_on = 'ALWAYS'
[]
[]
[VectorPostprocessors]
[centerline]
type = LineValueSampler
start_point = '${eps} ${eps} ${fparse 0.5 * L / nz}'
end_point = '${eps} ${eps} ${fparse L - 0.5 * L / nz}'
num_points = ${nz}
variable = vel_z
sort_by = 'z'
execute_on = 'FINAL'
[]
!include radial.i
[]
[Executioner]
type = SIMPLE
rhie_chow_user_object = 'rc'
momentum_systems = 'u_system v_system w_system'
pressure_system = 'pressure_system'
turbulence_systems = 'TKE_system TKED_system'
momentum_l_abs_tol = 1e-6
pressure_l_abs_tol = 1e-6
turbulence_l_abs_tol = 1e-6
momentum_l_tol = 1e-6
pressure_l_tol = 1e-6
turbulence_l_tol = 1e-6
momentum_equation_relaxation = 0.6
pressure_variable_relaxation = 0.25
turbulence_equation_relaxation = '0.4 0.4'
turbulence_field_relaxation = '0.4 0.4'
num_iterations = 1500
pressure_absolute_tolerance = 1e-7
momentum_absolute_tolerance = 1e-7
turbulence_absolute_tolerance = '1e-7 1e-7'
momentum_petsc_options_iname = '-u_system_pc_type -u_system_pc_hypre_type -v_system_pc_type -v_system_pc_hypre_type'
momentum_petsc_options_value = 'hypre boomeramg hypre boomeramg'
pressure_petsc_options_iname = '-pc_type -pc_hypre_type'
pressure_petsc_options_value = 'hypre boomeramg'
turbulence_petsc_options_iname = '-TKE_system_pc_type -TKE_system_pc_hypre_type -TKED_system_pc_type -TKED_system_pc_hypre_type'
turbulence_petsc_options_value = 'hypre boomeramg hypre boomeramg'
momentum_l_max_its = 300
pressure_l_max_its = 300
turbulence_l_max_its = 300
print_fields = false
continue_on_max_its = true
[]
[Outputs]
exodus = true
[out]
type = CSV
execute_on = FINAL
[]
[]
(test/tests/kEpsilon/special-cases/channel_ERCOFTAC_low_Re.i)
##########################################################
# ERCOFTAC test case for turbulent channel flow
# Case Number: 032
# Author: Dr. Mauricio Tano & Hailey Tran Kieu
# Last Update: November, 2023 & 2025
# Turbulent model using:
# k-epsilon
# Equilibrium + Newton wall treatment
# SIMPLE solve
##########################################################
### Problem Parameters ###
H = 1 # half-width of the channel
L = 120
Re = 2000
rho = 1
bulk_u = 1
mu = '${fparse rho * bulk_u * 2 * H / Re}'
advected_interp_method = 'upwind'
### k-epsilon Closure Parameters ###
sigma_k = 1.0
sigma_eps = 1.3
C1_eps = 1.44
C2_eps = 1.92
C_mu = 0.09
### Initial and Boundary Conditions ###
intensity = 0.01
k_init = '${fparse 1.5*(intensity * bulk_u)^2}'
eps_init = '${fparse C_mu^0.75 * k_init^1.5 / (2*H)}'
### Modeling parameters ###
bulk_wall_treatment = false
walls = 'bottom top'
wall_treatment = 'eq_newton' # Options: eq_newton, eq_incremental, eq_linearized, neq
# Turbulence-model knobs (optional)
k_epsilon_variant = 'RealizableTwoLayer' # Standard | StandardLowRe | StandardTwoLayer | Realizable | RealizableTwoLayer
two_layer_flavor = 'Wolfstein' # Wolfstein | NorrisReynolds | Xu (only used for *TwoLayer variants)
use_buoyancy = false
use_compressibility = false
nonlinear_model = 'none'
curvature_model = 'none'
use_yap = false
[Mesh]
[block_1]
type = GeneratedMeshGenerator
dim = 2
xmin = 0
xmax = ${L}
ymin = 0
ymax = ${H}
nx = 10
ny = 5
bias_y = 0.9
[]
[block_2]
type = GeneratedMeshGenerator
dim = 2
xmin = 0
xmax = ${L}
ymin = ${fparse -H}
ymax = 0
nx = 10
ny = 5
bias_y = ${fparse 1/0.9}
[]
[smg]
type = StitchMeshGenerator
inputs = 'block_1 block_2'
clear_stitched_boundary_ids = true
stitch_boundaries_pairs = 'bottom top'
merge_boundaries_with_same_name = true
[]
# Prevent test diffing on distributed parallel element numbering
allow_renumbering = false
[]
[Problem]
linear_sys_names = 'u_system v_system pressure_system TKE_system TKED_system'
previous_nl_solution_required = true
[]
[GlobalParams]
rhie_chow_user_object = 'rc'
[]
[UserObjects]
[rc]
type = RhieChowMassFlux
u = vel_x
v = vel_y
pressure = pressure
rho = ${rho}
p_diffusion_kernel = p_diffusion
pressure_projection_method = 'consistent'
[]
[]
[Variables]
[vel_x]
type = MooseLinearVariableFVReal
initial_condition = ${bulk_u}
solver_sys = u_system
[]
[vel_y]
type = MooseLinearVariableFVReal
initial_condition = 0
solver_sys = v_system
[]
[pressure]
type = MooseLinearVariableFVReal
initial_condition = 1e-8
solver_sys = pressure_system
[]
[TKE]
type = MooseLinearVariableFVReal
solver_sys = TKE_system
initial_condition = ${k_init}
[]
[TKED]
type = MooseLinearVariableFVReal
solver_sys = TKED_system
initial_condition = ${eps_init}
[]
[]
[FVInterpolationMethods]
[upwind]
type = FVAdvectedUpwind
[]
[]
[LinearFVKernels]
[u_advection_stress]
type = LinearWCNSFVMomentumFlux
variable = vel_x
advected_interp_method_name = ${advected_interp_method}
mu = 'mu_t'
u = vel_x
v = vel_y
momentum_component = 'x'
rhie_chow_user_object = 'rc'
use_nonorthogonal_correction = false
use_deviatoric_terms = yes
[]
[u_diffusion]
type = LinearFVDiffusion
variable = vel_x
diffusion_coeff = '${mu}'
[]
[u_pressure]
type = LinearFVMomentumPressure
variable = vel_x
pressure = pressure
momentum_component = 'x'
[]
[v_advection_stress]
type = LinearWCNSFVMomentumFlux
variable = vel_y
advected_interp_method_name = ${advected_interp_method}
mu = 'mu_t'
u = vel_x
v = vel_y
momentum_component = 'y'
rhie_chow_user_object = 'rc'
use_nonorthogonal_correction = false
use_deviatoric_terms = yes
[]
[v_diffusion]
type = LinearFVDiffusion
variable = vel_y
diffusion_coeff = '${mu}'
[]
[v_pressure]
type = LinearFVMomentumPressure
variable = vel_y
pressure = pressure
momentum_component = 'y'
[]
[p_diffusion]
type = LinearFVPressureCorrectionDiffusion
variable = pressure
diffusion_tensor = Ainv
use_nonorthogonal_correction = false
[]
[HbyA_divergence]
type = LinearFVDivergence
variable = pressure
face_flux = HbyA
force_boundary_execution = true
[]
[TKE_advection]
type = LinearFVTurbulentAdvection
variable = TKE
advected_interp_method_name = ${advected_interp_method}
[]
[TKE_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKE
diffusion_coeff = ${mu}
use_nonorthogonal_correction = false
[]
[TKE_turb_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKE
diffusion_coeff = 'mu_t'
scaling_coeff = ${sigma_k}
use_nonorthogonal_correction = false
[]
[TKE_source_sink]
type = kEpsilonTKESourceSink
variable = TKE
u = vel_x
v = vel_y
epsilon = TKED
rho = ${rho}
mu = ${mu}
mu_t = 'mu_t'
walls = ${walls}
wall_treatment = ${wall_treatment}
C_pl = 1e10
# NEW (optional)
k_epsilon_variant = ${k_epsilon_variant} # 'Standard', 'Realizable', etc.
use_buoyancy = ${use_buoyancy}
use_compressibility = ${use_compressibility}
nonlinear_model = ${nonlinear_model}
curvature_model = ${curvature_model}
[]
[TKED_advection]
type = LinearFVTurbulentAdvection
variable = TKED
advected_interp_method_name = ${advected_interp_method}
walls = ${walls}
[]
[TKED_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKED
diffusion_coeff = ${mu}
use_nonorthogonal_correction = false
walls = ${walls}
[]
[TKED_turb_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKED
diffusion_coeff = 'mu_t'
scaling_coeff = ${sigma_eps}
use_nonorthogonal_correction = false
walls = ${walls}
[]
[TKED_source_sink]
type = kEpsilonTKEDSourceSink
variable = TKED
u = vel_x
v = vel_y
tke = TKE
rho = ${rho}
mu = ${mu}
mu_t = 'mu_t'
C1_eps = ${C1_eps}
C2_eps = ${C2_eps}
walls = ${walls}
wall_treatment = ${wall_treatment}
C_pl = 1e10
# NEW (optional)
k_epsilon_variant = ${k_epsilon_variant}
use_buoyancy = ${use_buoyancy}
use_compressibility = ${use_compressibility}
nonlinear_model = ${nonlinear_model}
curvature_model = ${curvature_model}
use_yap = ${use_yap}
wall_distance = wall_distance # for low-Re / two-layer Yap / G' terms
[]
[]
[LinearFVBCs]
[inlet-u]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'left'
variable = vel_x
functor = '${bulk_u}'
[]
[inlet-v]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'left'
variable = vel_y
functor = '0.0'
[]
[walls-u]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'top bottom'
variable = vel_x
functor = 0.0
[]
[walls-v]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'top bottom'
variable = vel_y
functor = 0.0
[]
[outlet_u]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = 'right'
variable = vel_x
use_two_term_expansion = false
[]
[outlet_v]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = 'right'
variable = vel_y
use_two_term_expansion = false
[]
[outlet_p]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'right'
variable = pressure
functor = 0.0
[]
[inlet_TKE]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'left'
variable = TKE
functor = '${k_init}'
[]
[outlet_TKE]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = 'right'
variable = TKE
use_two_term_expansion = false
[]
[inlet_TKED]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'left'
variable = TKED
functor = '${eps_init}'
[]
[outlet_TKED]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = 'right'
variable = TKED
use_two_term_expansion = false
[]
[walls_mu_t]
type = LinearFVTurbulentViscosityWallFunctionBC
boundary = 'bottom top'
variable = 'mu_t'
u = vel_x
v = vel_y
rho = ${rho}
mu = ${mu}
tke = TKE
wall_treatment = ${wall_treatment}
[]
[]
[AuxVariables]
[wall_distance]
type = MooseVariableFVReal
initial_condition = 1.0
[]
[mu_t]
type = MooseLinearVariableFVReal
initial_condition = '${fparse rho * C_mu * ${k_init}^2 / eps_init}'
[]
[yplus]
type = MooseVariableFVReal
two_term_boundary_expansion = false
[]
[]
[AuxKernels]
[compute_wall_distance]
type = WallDistanceAux
variable = wall_distance
walls = ${walls}
execute_on = 'INITIAL NONLINEAR'
[]
[compute_mu_t]
type = kEpsilonViscosity
variable = mu_t
C_mu = ${C_mu}
tke = TKE
epsilon = TKED
mu = ${mu}
rho = ${rho}
u = vel_x
v = vel_y
bulk_wall_treatment = ${bulk_wall_treatment}
walls = ${walls}
wall_treatment = ${wall_treatment}
mu_t_ratio_max = 1e20
execute_on = 'NONLINEAR'
# NEW (optional) – choose model and options
k_epsilon_variant = ${k_epsilon_variant} # e.g. 'Standard' or 'Realizable'
two_layer_flavor = ${two_layer_flavor} # ignored unless *TwoLayer variants
wall_distance = wall_distance # only needed for LowRe/TwoLayer (see below)
[]
[compute_y_plus]
type = RANSYPlusAux
variable = yplus
tke = TKE
mu = ${mu}
rho = ${rho}
u = vel_x
v = vel_y
walls = ${walls}
wall_treatment = ${wall_treatment}
execute_on = 'NONLINEAR'
[]
[]
[Executioner]
type = SIMPLE
rhie_chow_user_object = 'rc'
momentum_systems = 'u_system v_system'
pressure_system = 'pressure_system'
turbulence_systems = 'TKE_system TKED_system'
momentum_l_abs_tol = 1e-14
pressure_l_abs_tol = 1e-14
turbulence_l_abs_tol = 1e-14
momentum_l_tol = 1e-14
pressure_l_tol = 1e-14
turbulence_l_tol = 1e-14
momentum_equation_relaxation = 0.7
pressure_variable_relaxation = 0.3
turbulence_equation_relaxation = '0.2 0.2'
turbulence_field_relaxation = '0.2 0.2'
num_iterations = 1000
pressure_absolute_tolerance = 1e-7
momentum_absolute_tolerance = 1e-7
turbulence_absolute_tolerance = '1e-7 1e-7'
momentum_petsc_options_iname = '-pc_type -pc_hypre_type'
momentum_petsc_options_value = 'hypre boomeramg'
pressure_petsc_options_iname = '-pc_type -pc_hypre_type'
pressure_petsc_options_value = 'hypre boomeramg'
turbulence_petsc_options_iname = '-pc_type -pc_hypre_type'
turbulence_petsc_options_value = 'hypre boomeramg'
momentum_l_max_its = 300
pressure_l_max_its = 300
turbulence_l_max_its = 30
print_fields = false
continue_on_max_its = true
[]
[Outputs]
exodus = false
[csv]
type = CSV
execute_on = FINAL
[]
[]
variables_to_sample = 'vel_x vel_y pressure TKE TKED'
[VectorPostprocessors]
[side_bottom]
type = SideValueSampler
boundary = 'bottom'
variable = ${variables_to_sample}
sort_by = 'x'
execute_on = 'timestep_end'
[]
[side_top]
type = SideValueSampler
boundary = 'top'
variable = ${variables_to_sample}
sort_by = 'x'
execute_on = 'timestep_end'
[]
[line_center_channel]
type = LineValueSampler
start_point = '${fparse 0.125 * L} ${fparse 0.0001} 0'
end_point = '${fparse 0.875 * L} ${fparse 0.0001} 0'
num_points = ${Mesh/block_1/nx}
variable = ${variables_to_sample}
sort_by = 'x'
execute_on = 'timestep_end'
[]
[line_quarter_radius_channel]
type = LineValueSampler
start_point = '${fparse 0.125 * L} ${fparse 0.5 * H} 0'
end_point = '${fparse 0.875 * L} ${fparse 0.5 * H} 0'
num_points = ${Mesh/block_1/nx}
variable = ${variables_to_sample}
sort_by = 'x'
execute_on = 'timestep_end'
[]
[]
(test/tests/kEpsilon/special-cases/channel_ERCOFTAC.i)
##########################################################
# ERCOFTAC test case for turbulent channel flow
# Case Number: 032
# Author: Dr. Mauricio Tano & Hailey Tran Kieu
# Last Update: November, 2023 & 2025
# Turbulent model using:
# k-epsilon
# Equilibrium + Newton wall treatment
# SIMPLE solve
##########################################################
### Problem Parameters ###
H = 1 # half-width of the channel
L = 120
Re = 14000
rho = 1
bulk_u = 1
mu = '${fparse rho * bulk_u * 2 * H / Re}'
advected_interp_method = 'upwind'
### k-epsilon Closure Parameters ###
sigma_k = 1.0
sigma_eps = 1.3
C1_eps = 1.44
C2_eps = 1.92
C_mu = 0.09
### Initial and Boundary Conditions ###
intensity = 0.01
k_init = '${fparse 1.5*(intensity * bulk_u)^2}'
eps_init = '${fparse C_mu^0.75 * k_init^1.5 / (2*H)}'
### Modeling parameters ###
bulk_wall_treatment = false
walls = 'bottom top'
wall_treatment = 'eq_newton' # Options: eq_newton, eq_incremental, eq_linearized, neq
# NOTE:
# The k-epsilon kernels/aux-kernels now include strict parameter applicability
# checks. To keep this input file usable across all test permutations, we do
# *not* hard-code model-variant-specific options here. Instead, the test harness
# sets per-object options via CLI overrides (see create_tests.py).
[Mesh]
[block_1]
type = GeneratedMeshGenerator
dim = 2
xmin = 0
xmax = ${L}
ymin = 0
ymax = ${H}
nx = 10
ny = 5
bias_y = 0.7
[]
[block_2]
type = GeneratedMeshGenerator
dim = 2
xmin = 0
xmax = ${L}
ymin = ${fparse -H}
ymax = 0
nx = 10
ny = 5
bias_y = ${fparse 1/0.7}
[]
[smg]
type = StitchMeshGenerator
inputs = 'block_1 block_2'
clear_stitched_boundary_ids = true
stitch_boundaries_pairs = 'bottom top'
merge_boundaries_with_same_name = true
[]
# Prevent test diffing on distributed parallel element numbering
allow_renumbering = false
[]
[Problem]
linear_sys_names = 'u_system v_system pressure_system TKE_system TKED_system'
previous_nl_solution_required = true
[]
[GlobalParams]
rhie_chow_user_object = 'rc'
[]
[UserObjects]
[rc]
type = RhieChowMassFlux
u = vel_x
v = vel_y
pressure = pressure
rho = ${rho}
p_diffusion_kernel = p_diffusion
pressure_projection_method = 'consistent'
[]
[]
[Variables]
[vel_x]
type = MooseLinearVariableFVReal
initial_condition = ${bulk_u}
solver_sys = u_system
[]
[vel_y]
type = MooseLinearVariableFVReal
initial_condition = 0
solver_sys = v_system
[]
[pressure]
type = MooseLinearVariableFVReal
initial_condition = 1e-8
solver_sys = pressure_system
[]
[TKE]
type = MooseLinearVariableFVReal
solver_sys = TKE_system
initial_condition = ${k_init}
[]
[TKED]
type = MooseLinearVariableFVReal
solver_sys = TKED_system
initial_condition = ${eps_init}
[]
[]
[FVInterpolationMethods]
[upwind]
type = FVAdvectedUpwind
[]
[]
[LinearFVKernels]
[u_advection_stress]
type = LinearWCNSFVMomentumFlux
variable = vel_x
advected_interp_method_name = ${advected_interp_method}
mu = 'mu_t'
u = vel_x
v = vel_y
momentum_component = 'x'
rhie_chow_user_object = 'rc'
use_nonorthogonal_correction = false
use_deviatoric_terms = yes
[]
[u_diffusion]
type = LinearFVDiffusion
variable = vel_x
diffusion_coeff = '${mu}'
[]
[u_pressure]
type = LinearFVMomentumPressure
variable = vel_x
pressure = pressure
momentum_component = 'x'
[]
[v_advection_stress]
type = LinearWCNSFVMomentumFlux
variable = vel_y
advected_interp_method_name = ${advected_interp_method}
mu = 'mu_t'
u = vel_x
v = vel_y
momentum_component = 'y'
rhie_chow_user_object = 'rc'
use_nonorthogonal_correction = false
use_deviatoric_terms = yes
[]
[v_diffusion]
type = LinearFVDiffusion
variable = vel_y
diffusion_coeff = '${mu}'
[]
[v_pressure]
type = LinearFVMomentumPressure
variable = vel_y
pressure = pressure
momentum_component = 'y'
[]
[p_diffusion]
type = LinearFVPressureCorrectionDiffusion
variable = pressure
diffusion_tensor = Ainv
use_nonorthogonal_correction = false
[]
[HbyA_divergence]
type = LinearFVDivergence
variable = pressure
face_flux = HbyA
force_boundary_execution = true
[]
[TKE_advection]
type = LinearFVTurbulentAdvection
variable = TKE
advected_interp_method_name = ${advected_interp_method}
[]
[TKE_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKE
diffusion_coeff = ${mu}
use_nonorthogonal_correction = false
[]
[TKE_turb_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKE
diffusion_coeff = 'mu_t'
scaling_coeff = ${sigma_k}
use_nonorthogonal_correction = false
[]
[TKE_source_sink]
type = kEpsilonTKESourceSink
variable = TKE
u = vel_x
v = vel_y
epsilon = TKED
rho = ${rho}
mu = ${mu}
mu_t = 'mu_t'
walls = ${walls}
wall_treatment = ${wall_treatment}
C_pl = 1e10
[]
[TKED_advection]
type = LinearFVTurbulentAdvection
variable = TKED
advected_interp_method_name = ${advected_interp_method}
walls = ${walls}
[]
[TKED_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKED
diffusion_coeff = ${mu}
use_nonorthogonal_correction = false
walls = ${walls}
[]
[TKED_turb_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKED
diffusion_coeff = 'mu_t'
scaling_coeff = ${sigma_eps}
use_nonorthogonal_correction = false
walls = ${walls}
[]
[TKED_source_sink]
type = kEpsilonTKEDSourceSink
variable = TKED
u = vel_x
v = vel_y
tke = TKE
rho = ${rho}
mu = ${mu}
mu_t = 'mu_t'
C1_eps = ${C1_eps}
C2_eps = ${C2_eps}
walls = ${walls}
wall_treatment = ${wall_treatment}
C_pl = 1e10
[]
[]
[LinearFVBCs]
[inlet-u]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'left'
variable = vel_x
functor = '${bulk_u}'
[]
[inlet-v]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'left'
variable = vel_y
functor = '0.0'
[]
[walls-u]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'top bottom'
variable = vel_x
functor = 0.0
[]
[walls-v]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'top bottom'
variable = vel_y
functor = 0.0
[]
[outlet_u]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = 'right'
variable = vel_x
use_two_term_expansion = false
[]
[outlet_v]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = 'right'
variable = vel_y
use_two_term_expansion = false
[]
[outlet_p]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'right'
variable = pressure
functor = 0.0
[]
[inlet_TKE]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'left'
variable = TKE
functor = '${k_init}'
[]
[outlet_TKE]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = 'right'
variable = TKE
use_two_term_expansion = false
[]
[inlet_TKED]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'left'
variable = TKED
functor = '${eps_init}'
[]
[outlet_TKED]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = 'right'
variable = TKED
use_two_term_expansion = false
[]
[walls_mu_t]
type = LinearFVTurbulentViscosityWallFunctionBC
boundary = 'bottom top'
variable = 'mu_t'
u = vel_x
v = vel_y
rho = ${rho}
mu = ${mu}
tke = TKE
wall_treatment = ${wall_treatment}
[]
[]
[AuxVariables]
[mu_t]
type = MooseLinearVariableFVReal
initial_condition = '${fparse rho * C_mu * ${k_init}^2 / eps_init}'
[]
[yplus]
type = MooseVariableFVReal
two_term_boundary_expansion = false
[]
[]
[AuxKernels]
[compute_mu_t]
type = kEpsilonViscosity
variable = mu_t
C_mu = ${C_mu}
tke = TKE
epsilon = TKED
mu = ${mu}
rho = ${rho}
u = vel_x
v = vel_y
bulk_wall_treatment = ${bulk_wall_treatment}
walls = ${walls}
wall_treatment = ${wall_treatment}
mu_t_ratio_max = 1e20
execute_on = 'NONLINEAR'
[]
[compute_y_plus]
type = RANSYPlusAux
variable = yplus
tke = TKE
mu = ${mu}
rho = ${rho}
u = vel_x
v = vel_y
walls = ${walls}
wall_treatment = ${wall_treatment}
execute_on = 'NONLINEAR'
[]
[]
[Executioner]
type = SIMPLE
rhie_chow_user_object = 'rc'
momentum_systems = 'u_system v_system'
pressure_system = 'pressure_system'
turbulence_systems = 'TKE_system TKED_system'
momentum_l_abs_tol = 1e-14
pressure_l_abs_tol = 1e-14
turbulence_l_abs_tol = 1e-14
momentum_l_tol = 1e-14
pressure_l_tol = 1e-14
turbulence_l_tol = 1e-14
momentum_equation_relaxation = 0.7
pressure_variable_relaxation = 0.3
turbulence_equation_relaxation = '0.2 0.2'
turbulence_field_relaxation = '0.2 0.2'
num_iterations = 1000
pressure_absolute_tolerance = 1e-7
momentum_absolute_tolerance = 1e-7
turbulence_absolute_tolerance = '1e-7 1e-7'
momentum_petsc_options_iname = '-pc_type -pc_hypre_type'
momentum_petsc_options_value = 'hypre boomeramg'
pressure_petsc_options_iname = '-pc_type -pc_hypre_type'
pressure_petsc_options_value = 'hypre boomeramg'
turbulence_petsc_options_iname = '-pc_type -pc_hypre_type'
turbulence_petsc_options_value = 'hypre boomeramg'
momentum_l_max_its = 300
pressure_l_max_its = 300
turbulence_l_max_its = 30
print_fields = false
continue_on_max_its = true
[]
[Outputs]
exodus = false
[csv]
type = CSV
execute_on = FINAL
[]
[]
variables_to_sample = 'vel_x vel_y pressure TKE TKED'
[VectorPostprocessors]
[side_bottom]
type = SideValueSampler
boundary = 'bottom'
variable = ${variables_to_sample}
sort_by = 'x'
execute_on = 'timestep_end'
[]
[side_top]
type = SideValueSampler
boundary = 'top'
variable = ${variables_to_sample}
sort_by = 'x'
execute_on = 'timestep_end'
[]
[line_center_channel]
type = LineValueSampler
start_point = '${fparse 0.125 * L} ${fparse 0.0001} 0'
end_point = '${fparse 0.875 * L} ${fparse 0.0001} 0'
num_points = ${Mesh/block_1/nx}
variable = ${variables_to_sample}
sort_by = 'x'
execute_on = 'timestep_end'
[]
[line_quarter_radius_channel]
type = LineValueSampler
start_point = '${fparse 0.125 * L} ${fparse 0.5 * H} 0'
end_point = '${fparse 0.875 * L} ${fparse 0.5 * H} 0'
num_points = ${Mesh/block_1/nx}
variable = ${variables_to_sample}
sort_by = 'x'
execute_on = 'timestep_end'
[]
[]
(test/tests/kEpsilon/special-cases/channel_ERCOFTAC_bulk_treatment.i)
##########################################################
# ERCOFTAC test case for turbulent channel flow
# Case Number: 032
# Author: Dr. Mauricio Tano & Hailey Tran Kieu
# Last Update: November, 2023 & 2025
# Turbulent model using:
# k-epsilon
# Equilibrium + Newton wall treatment
# SIMPLE solve
##########################################################
### Problem Parameters ###
H = 1 # half-width of the channel
L = 120
Re = 14000
rho = 1
bulk_u = 1
mu = '${fparse rho * bulk_u * 2 * H / Re}'
advected_interp_method = 'upwind'
### k-epsilon Closure Parameters ###
sigma_k = 1.0
sigma_eps = 1.3
C1_eps = 1.44
C2_eps = 1.92
C_mu = 0.09
### Initial and Boundary Conditions ###
intensity = 0.01
k_init = '${fparse 1.5*(intensity * bulk_u)^2}'
eps_init = '${fparse C_mu^0.75 * k_init^1.5 / (2*H)}'
### Modeling parameters ###
bulk_wall_treatment = true
walls = 'bottom top'
wall_treatment = 'eq_newton' # Options: eq_newton, eq_incremental, eq_linearized, neq
# Turbulence-model knobs (optional)
k_epsilon_variant = 'Standard'
use_buoyancy = false
use_compressibility = false
nonlinear_model = 'none'
curvature_model = 'none'
use_yap = false
use_low_re_Gprime = false
[Mesh]
[block_1]
type = GeneratedMeshGenerator
dim = 2
xmin = 0
xmax = ${L}
ymin = 0
ymax = ${H}
nx = 10
ny = 3
bias_y = 1.0
[]
[block_2]
type = GeneratedMeshGenerator
dim = 2
xmin = 0
xmax = ${L}
ymin = ${fparse -H}
ymax = 0
nx = 10
ny = 3
bias_y = 1.0
[]
[smg]
type = StitchMeshGenerator
inputs = 'block_1 block_2'
clear_stitched_boundary_ids = true
stitch_boundaries_pairs = 'bottom top'
merge_boundaries_with_same_name = true
[]
# Prevent test diffing on distributed parallel element numbering
allow_renumbering = false
[]
[Problem]
linear_sys_names = 'u_system v_system pressure_system TKE_system TKED_system'
previous_nl_solution_required = true
[]
[GlobalParams]
rhie_chow_user_object = 'rc'
[]
[UserObjects]
[rc]
type = RhieChowMassFlux
u = vel_x
v = vel_y
pressure = pressure
rho = ${rho}
p_diffusion_kernel = p_diffusion
pressure_projection_method = 'consistent'
[]
[]
[Variables]
[vel_x]
type = MooseLinearVariableFVReal
initial_condition = ${bulk_u}
solver_sys = u_system
[]
[vel_y]
type = MooseLinearVariableFVReal
initial_condition = 0
solver_sys = v_system
[]
[pressure]
type = MooseLinearVariableFVReal
initial_condition = 1e-8
solver_sys = pressure_system
[]
[TKE]
type = MooseLinearVariableFVReal
solver_sys = TKE_system
initial_condition = ${k_init}
[]
[TKED]
type = MooseLinearVariableFVReal
solver_sys = TKED_system
initial_condition = ${eps_init}
[]
[]
[FVInterpolationMethods]
[upwind]
type = FVAdvectedUpwind
[]
[]
[LinearFVKernels]
[u_advection_stress]
type = LinearWCNSFVMomentumFlux
variable = vel_x
advected_interp_method_name = ${advected_interp_method}
mu = 'mu_t'
u = vel_x
v = vel_y
momentum_component = 'x'
rhie_chow_user_object = 'rc'
use_nonorthogonal_correction = false
use_deviatoric_terms = yes
[]
[u_diffusion]
type = LinearFVDiffusion
variable = vel_x
diffusion_coeff = '${mu}'
[]
[u_pressure]
type = LinearFVMomentumPressure
variable = vel_x
pressure = pressure
momentum_component = 'x'
[]
[v_advection_stress]
type = LinearWCNSFVMomentumFlux
variable = vel_y
advected_interp_method_name = ${advected_interp_method}
mu = 'mu_t'
u = vel_x
v = vel_y
momentum_component = 'y'
rhie_chow_user_object = 'rc'
use_nonorthogonal_correction = false
use_deviatoric_terms = yes
[]
[v_diffusion]
type = LinearFVDiffusion
variable = vel_y
diffusion_coeff = '${mu}'
[]
[v_pressure]
type = LinearFVMomentumPressure
variable = vel_y
pressure = pressure
momentum_component = 'y'
[]
[p_diffusion]
type = LinearFVPressureCorrectionDiffusion
variable = pressure
diffusion_tensor = Ainv
use_nonorthogonal_correction = false
[]
[HbyA_divergence]
type = LinearFVDivergence
variable = pressure
face_flux = HbyA
force_boundary_execution = true
[]
[TKE_advection]
type = LinearFVTurbulentAdvection
variable = TKE
advected_interp_method_name = ${advected_interp_method}
[]
[TKE_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKE
diffusion_coeff = ${mu}
use_nonorthogonal_correction = false
[]
[TKE_turb_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKE
diffusion_coeff = 'mu_t'
scaling_coeff = ${sigma_k}
use_nonorthogonal_correction = false
[]
[TKE_source_sink]
type = kEpsilonTKESourceSink
variable = TKE
u = vel_x
v = vel_y
epsilon = TKED
rho = ${rho}
mu = ${mu}
mu_t = 'mu_t'
walls = ${walls}
wall_treatment = ${wall_treatment}
C_pl = 1e10
# NEW (optional)
k_epsilon_variant = ${k_epsilon_variant} # 'Standard', 'Realizable', etc.
use_buoyancy = ${use_buoyancy}
use_compressibility = ${use_compressibility}
nonlinear_model = ${nonlinear_model}
curvature_model = ${curvature_model}
[]
[TKED_advection]
type = LinearFVTurbulentAdvection
variable = TKED
advected_interp_method_name = ${advected_interp_method}
walls = ${walls}
[]
[TKED_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKED
diffusion_coeff = ${mu}
use_nonorthogonal_correction = false
walls = ${walls}
[]
[TKED_turb_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKED
diffusion_coeff = 'mu_t'
scaling_coeff = ${sigma_eps}
use_nonorthogonal_correction = false
walls = ${walls}
[]
[TKED_source_sink]
type = kEpsilonTKEDSourceSink
variable = TKED
u = vel_x
v = vel_y
tke = TKE
rho = ${rho}
mu = ${mu}
mu_t = 'mu_t'
C1_eps = ${C1_eps}
C2_eps = ${C2_eps}
walls = ${walls}
wall_treatment = ${wall_treatment}
C_pl = 1e10
# NEW (optional)
k_epsilon_variant = ${k_epsilon_variant}
use_buoyancy = ${use_buoyancy}
use_compressibility = ${use_compressibility}
nonlinear_model = ${nonlinear_model}
curvature_model = ${curvature_model}
use_yap = ${use_yap}
use_low_re_Gprime = ${use_low_re_Gprime}
wall_distance = wall_distance # for low-Re / two-layer Yap / G' terms
[]
[]
[LinearFVBCs]
[inlet-u]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'left'
variable = vel_x
functor = '${bulk_u}'
[]
[inlet-v]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'left'
variable = vel_y
functor = '0.0'
[]
[walls-u]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'top bottom'
variable = vel_x
functor = 0.0
[]
[walls-v]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'top bottom'
variable = vel_y
functor = 0.0
[]
[outlet_u]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = 'right'
variable = vel_x
use_two_term_expansion = false
[]
[outlet_v]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = 'right'
variable = vel_y
use_two_term_expansion = false
[]
[outlet_p]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'right'
variable = pressure
functor = 0.0
[]
[inlet_TKE]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'left'
variable = TKE
functor = '${k_init}'
[]
[outlet_TKE]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = 'right'
variable = TKE
use_two_term_expansion = false
[]
[inlet_TKED]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'left'
variable = TKED
functor = '${eps_init}'
[]
[outlet_TKED]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = 'right'
variable = TKED
use_two_term_expansion = false
[]
[]
[AuxVariables]
[wall_distance]
type = MooseVariableFVReal
initial_condition = 1.0
[]
[mu_t]
type = MooseLinearVariableFVReal
initial_condition = '${fparse rho * C_mu * ${k_init}^2 / eps_init}'
[]
[yplus]
type = MooseVariableFVReal
two_term_boundary_expansion = false
[]
[]
[AuxKernels]
[compute_wall_distance]
type = WallDistanceAux
variable = wall_distance
walls = ${walls}
execute_on = 'INITIAL NONLINEAR'
[]
[compute_mu_t]
type = kEpsilonViscosity
variable = mu_t
C_mu = ${C_mu}
tke = TKE
epsilon = TKED
mu = ${mu}
rho = ${rho}
u = vel_x
v = vel_y
bulk_wall_treatment = ${bulk_wall_treatment}
walls = ${walls}
wall_treatment = ${wall_treatment}
mu_t_ratio_max = 1e20
execute_on = 'NONLINEAR'
# NEW (optional) – choose model and options
k_epsilon_variant = ${k_epsilon_variant} # e.g. 'Standard' or 'Realizable'
wall_distance = wall_distance # only needed for LowRe/TwoLayer (see below)
[]
[compute_y_plus]
type = RANSYPlusAux
variable = yplus
tke = TKE
mu = ${mu}
rho = ${rho}
u = vel_x
v = vel_y
walls = ${walls}
wall_treatment = ${wall_treatment}
execute_on = 'NONLINEAR'
[]
[]
[Executioner]
type = SIMPLE
rhie_chow_user_object = 'rc'
momentum_systems = 'u_system v_system'
pressure_system = 'pressure_system'
turbulence_systems = 'TKE_system TKED_system'
momentum_l_abs_tol = 1e-14
pressure_l_abs_tol = 1e-14
turbulence_l_abs_tol = 1e-14
momentum_l_tol = 1e-14
pressure_l_tol = 1e-14
turbulence_l_tol = 1e-14
momentum_equation_relaxation = 0.7
pressure_variable_relaxation = 0.3
turbulence_equation_relaxation = '0.2 0.2'
turbulence_field_relaxation = '0.2 0.2'
num_iterations = 1000
pressure_absolute_tolerance = 1e-7
momentum_absolute_tolerance = 1e-7
turbulence_absolute_tolerance = '1e-7 1e-7'
momentum_petsc_options_iname = '-pc_type -pc_hypre_type'
momentum_petsc_options_value = 'hypre boomeramg'
pressure_petsc_options_iname = '-pc_type -pc_hypre_type'
pressure_petsc_options_value = 'hypre boomeramg'
turbulence_petsc_options_iname = '-pc_type -pc_hypre_type'
turbulence_petsc_options_value = 'hypre boomeramg'
momentum_l_max_its = 300
pressure_l_max_its = 300
turbulence_l_max_its = 30
print_fields = false
continue_on_max_its = true
[]
[Outputs]
exodus = false
[csv]
type = CSV
execute_on = FINAL
[]
[]
variables_to_sample = 'vel_x vel_y pressure TKE TKED'
[VectorPostprocessors]
[side_bottom]
type = SideValueSampler
boundary = 'bottom'
variable = ${variables_to_sample}
sort_by = 'x'
execute_on = 'timestep_end'
[]
[side_top]
type = SideValueSampler
boundary = 'top'
variable = ${variables_to_sample}
sort_by = 'x'
execute_on = 'timestep_end'
[]
[line_center_channel]
type = LineValueSampler
start_point = '${fparse 0.125 * L} ${fparse 0.0001} 0'
end_point = '${fparse 0.875 * L} ${fparse 0.0001} 0'
num_points = ${Mesh/block_1/nx}
variable = ${variables_to_sample}
sort_by = 'x'
execute_on = 'timestep_end'
[]
[line_quarter_radius_channel]
type = LineValueSampler
start_point = '${fparse 0.125 * L} ${fparse 0.5 * H} 0'
end_point = '${fparse 0.875 * L} ${fparse 0.5 * H} 0'
num_points = ${Mesh/block_1/nx}
variable = ${variables_to_sample}
sort_by = 'x'
execute_on = 'timestep_end'
[]
[]
(validation/free_flow/isothermal/ercoftac_030_bfs/bfs_input.i)
########################################################################
# BFS simulation
########################################################################
# Modeling Parameters
rho = 1.18415
bulk_u = 44.2
D = 0.1016
mu = 1.8551e-5
advected_interp_method = 'upwind'
### k-epsilon Closure Parameters ###
sigma_k = 1.0
sigma_eps = 1.3
C1_eps = 1.44
C2_eps = 1.92
C_mu = 0.09
walls = 'wall'
wall_treatment = 'neq' # Options: eq_newton, eq_incremental, eq_linearized, neq
### Initial and Boundary Conditions ###
intensity = 0.01
k_init = '${fparse 1.5*(intensity * bulk_u)^2}'
eps_init = '${fparse C_mu^0.75 * k_init^1.5 / D}'
### Postprocessing
y_first_cell_1 = 5e-4
y_first_cell_2 = ${fparse (0.0127 - 0.0119) / 2}
# Turbulence-model knobs (optional)
k_epsilon_variant = 'Standard' # Standard | StandardLowRe | StandardTwoLayer | Realizable | RealizableTwoLayer
# two_layer_flavor = 'Wolfstein' # Wolfstein | NorrisReynolds | Xu (only used for *TwoLayer variants)
use_buoyancy = false
use_compressibility = false
nonlinear_model = 'none'
curvature_model = 'none'
use_yap = false
use_low_re_Gprime = false
bulk_wall_treatment = false
[Mesh]
[load_mesh]
type = FileMeshGenerator
file = 'BFS_mesh_fine1.e'
[]
[]
[Problem]
linear_sys_names = 'u_system v_system pressure_system TKE_system TKED_system'
previous_nl_solution_required = true
[]
[GlobalParams]
rhie_chow_user_object = 'rc'
[]
[UserObjects]
[rc]
type = RhieChowMassFlux
u = vel_x
v = vel_y
pressure = pressure
rho = ${rho}
p_diffusion_kernel = p_diffusion
pressure_projection_method = 'consistent'
[]
[]
[Variables]
[vel_x]
type = MooseLinearVariableFVReal
solver_sys = u_system
[]
[vel_y]
type = MooseLinearVariableFVReal
solver_sys = v_system
[]
[pressure]
type = MooseLinearVariableFVReal
solver_sys = pressure_system
[]
[TKE]
type = MooseLinearVariableFVReal
solver_sys = TKE_system
[]
[TKED]
type = MooseLinearVariableFVReal
solver_sys = TKED_system
[]
[]
[FVICs]
[vel_x_ic]
type = FVFunctionIC
variable = 'vel_x'
function = fully_developed_velocity
[]
[vel_y_ic]
type = FVFunctionIC
variable = 'vel_y'
function = 0
[]
[pressure_ic]
type = FVFunctionIC
variable = 'pressure'
function = 0
[]
[TKE_ic]
type = FVFunctionIC
variable = 'TKE'
function = fully_developed_tke
[]
[TKED_ic]
type = FVFunctionIC
variable = 'TKED'
function = fully_developed_tked
[]
[]
[FVInterpolationMethods]
[upwind]
type = FVAdvectedUpwind
[]
[]
[LinearFVKernels]
[u_advection_stress]
type = LinearWCNSFVMomentumFlux
variable = vel_x
advected_interp_method_name = ${advected_interp_method}
mu = 'mu_t'
u = vel_x
v = vel_y
momentum_component = 'x'
rhie_chow_user_object = 'rc'
use_nonorthogonal_correction = false
use_deviatoric_terms = true
[]
[u_diffusion]
type = LinearFVDiffusion
variable = vel_x
diffusion_coeff = ${mu}
use_nonorthogonal_correction = false
[]
[u_pressure]
type = LinearFVMomentumPressure
variable = vel_x
pressure = pressure
momentum_component = 'x'
[]
[v_advection_stress]
type = LinearWCNSFVMomentumFlux
variable = vel_y
advected_interp_method_name = ${advected_interp_method}
mu = 'mu_t'
u = vel_x
v = vel_y
momentum_component = 'y'
rhie_chow_user_object = 'rc'
use_nonorthogonal_correction = false
use_deviatoric_terms = true
[]
[v_diffusion]
type = LinearFVDiffusion
variable = vel_y
diffusion_coeff = ${mu}
use_nonorthogonal_correction = false
[]
[v_pressure]
type = LinearFVMomentumPressure
variable = vel_y
pressure = pressure
momentum_component = 'y'
[]
[p_diffusion]
type = LinearFVPressureCorrectionDiffusion
variable = pressure
diffusion_tensor = Ainv
use_nonorthogonal_correction = false
[]
[HbyA_divergence]
type = LinearFVDivergence
variable = pressure
face_flux = HbyA
force_boundary_execution = true
[]
[TKE_advection]
type = LinearFVTurbulentAdvection
variable = TKE
advected_interp_method_name = ${advected_interp_method}
[]
[TKE_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKE
diffusion_coeff = ${mu}
use_nonorthogonal_correction = false
[]
[TKE_turb_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKE
diffusion_coeff = 'mu_t'
scaling_coeff = ${sigma_k}
use_nonorthogonal_correction = false
[]
[TKE_source_sink]
type = kEpsilonTKESourceSink
variable = TKE
u = vel_x
v = vel_y
epsilon = TKED
rho = ${rho}
mu = ${mu}
mu_t = 'mu_t'
walls = ${walls}
wall_treatment = ${wall_treatment}
C_pl = 2
# NEW (optional)
k_epsilon_variant = ${k_epsilon_variant} # 'Standard', 'Realizable', etc.
use_buoyancy = ${use_buoyancy}
use_compressibility = ${use_compressibility}
nonlinear_model = ${nonlinear_model}
curvature_model = ${curvature_model}
[]
[TKED_advection]
type = LinearFVTurbulentAdvection
variable = TKED
advected_interp_method_name = ${advected_interp_method}
walls = ${walls}
[]
[TKED_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKED
diffusion_coeff = ${mu}
use_nonorthogonal_correction = false
walls = ${walls}
[]
[TKED_turb_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKED
diffusion_coeff = 'mu_t'
scaling_coeff = ${sigma_eps}
use_nonorthogonal_correction = false
walls = ${walls}
[]
[TKED_source_sink]
type = kEpsilonTKEDSourceSink
variable = TKED
u = vel_x
v = vel_y
tke = TKE
rho = ${rho}
mu = ${mu}
mu_t = 'mu_t'
C1_eps = ${C1_eps}
C2_eps = ${C2_eps}
walls = ${walls}
wall_treatment = ${wall_treatment}
C_pl = 2
# NEW (optional)
k_epsilon_variant = ${k_epsilon_variant}
use_buoyancy = ${use_buoyancy}
use_compressibility = ${use_compressibility}
nonlinear_model = ${nonlinear_model}
curvature_model = ${curvature_model}
use_yap = ${use_yap}
use_low_re_Gprime = ${use_low_re_Gprime}
wall_distance = distance
[]
[]
[LinearFVBCs]
[inlet_u]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'inlet'
variable = vel_x
functor = 'fully_developed_velocity'
[]
[inlet-v]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'inlet'
variable = vel_y
functor = 0
[]
[inlet_TKE]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'inlet'
variable = TKE
functor = 'fully_developed_tke'
[]
# Alternative constant inlet BC
# [inlet_TKE]
# type = INSFVInletIntensityTKEBC
# boundary = 'inlet'
# variable = TKE
# u = vel_x
# v = vel_y
# w = vel_z
# intensity = ${intensity}
# []
[inlet_TKED]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'inlet'
variable = TKED
functor = 'fully_developed_tked'
[]
# Alternative constant inlet BC
# [inlet_TKED]
# type = INSFVMixingLengthTKEDBC
# boundary = 'inlet'
# variable = TKED
# k = TKE
# characteristic_length = '${D}'
# []
[outlet_u]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = 'outlet'
variable = vel_x
use_two_term_expansion = false
[]
[outlet_v]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = 'outlet'
variable = vel_y
use_two_term_expansion = false
[]
[outlet_tke]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = 'outlet'
variable = TKE
use_two_term_expansion = false
[]
[outlet_tked]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = 'outlet'
variable = TKED
use_two_term_expansion = false
[]
[outlet_p]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'outlet'
variable = pressure
functor = 0.0
[]
[walls-u]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = ${walls}
variable = vel_x
functor = 0
[]
[walls-v]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = ${walls}
variable = vel_y
functor = 0
[]
[walls_mu_t]
type = LinearFVTurbulentViscosityWallFunctionBC
boundary = ${walls}
variable = 'mu_t'
u = vel_x
v = vel_y
rho = ${rho}
mu = ${mu}
tke = TKE
wall_treatment = ${wall_treatment}
[]
[]
[AuxVariables]
[mu_t]
type = MooseLinearVariableFVReal
initial_condition = '${fparse rho * C_mu * ${k_init}^2 / eps_init}'
[]
# For computing c_f
[mu_t_wall]
type = MooseLinearVariableFVReal
[]
[yplus]
type = MooseLinearVariableFVReal
[]
[mu_eff]
type = MooseLinearVariableFVReal
initial_condition = '${fparse mu + rho * C_mu * ${k_init}^2 / eps_init}'
[]
[distance]
type = MooseVariableFVReal
[]
[]
[AuxKernels]
[compute_mu_t]
type = kEpsilonViscosity
variable = mu_t
C_mu = ${C_mu}
tke = TKE
epsilon = TKED
mu = ${mu}
rho = ${rho}
u = vel_x
v = vel_y
bulk_wall_treatment = ${bulk_wall_treatment}
walls = ${walls}
wall_treatment = ${wall_treatment}
mu_t_ratio_max = 1e20
execute_on = 'NONLINEAR'
# NEW (optional) – choose model and options
k_epsilon_variant = ${k_epsilon_variant} # e.g. 'Standard' or 'Realizable'
# two_layer_flavor = ${two_layer_flavor} # ignored unless *TwoLayer variants
# Cd0 = 0.091 # defaults, can omit if you keep Standard
# Cd1 = 0.0042
# Cd2 = 0.00011
# Ca0 = 0.667 # Realizable C_mu coefficients
# Ca1 = 1.25
# Ca2 = 1.0
# Ca3 = 0.9
wall_distance = distance # only needed for LowRe/TwoLayer (see below)
[]
[compute_mu_t_wall_neq]
type = ParsedAux
variable = 'mu_t_wall'
expression = '${mu} * (0.4187 * yplus / log(max(9.793 * yplus, 1.0 + 1e-4)) - 1)'
coupled_variables = 'yplus'
execute_on = 'TIMESTEP_END'
[]
# Equivalent way of computing mu_t_wall
# [compute_mu_t_wall]
# type = kEpsilonViscosityAux
# variable = mu_t_wall
# C_mu = ${C_mu}
# tke = TKE
# epsilon = TKED
# mu = ${mu}
# rho = ${rho}
# u = vel_x
# v = vel_y
# bulk_wall_treatment = true
# walls = ${walls}
# wall_treatment = ${wall_treatment}
# execute_on = 'TIMESTEP_END'
# mu_t_ratio_max = 1e20
# []
[compute_y_plus]
type = RANSYPlusAux
variable = yplus
tke = TKE
mu = ${mu}
rho = ${rho}
u = vel_x
v = vel_y
walls = ${walls}
wall_treatment = ${wall_treatment}
execute_on = 'NONLINEAR'
[]
[compute_mu_eff]
type = ParsedAux
variable = 'mu_eff'
coupled_variables = 'mu_t'
expression = 'mu_t + ${mu}'
execute_on = 'NONLINEAR'
[]
[distance_aux]
type = WallDistanceMixingLengthAux
walls = 'wall'
variable = 'distance'
execute_on = 'initial'
von_karman_const = 1.0
[]
[]
[UserObjects]
[read_recycling]
type = PropertyReadFile
prop_file_name = 'BFS_FDprofile.csv'
read_type = 'voronoi'
nprop = 7
execute_on = TIMESTEP_BEGIN
nvoronoi = 60
[]
[]
[Functions]
[fully_developed_velocity]
type = PiecewiseConstantFromCSV
read_prop_user_object = 'read_recycling'
read_type = 'voronoi'
column_number = '6'
[]
[fully_developed_tke]
type = PiecewiseConstantFromCSV
read_prop_user_object = 'read_recycling'
read_type = 'voronoi'
column_number = '3'
[]
[fully_developed_tked]
type = PiecewiseConstantFromCSV
read_prop_user_object = 'read_recycling'
read_type = 'voronoi'
column_number = '4'
[]
[]
[Executioner]
type = SIMPLE
rhie_chow_user_object = 'rc'
momentum_systems = 'u_system v_system'
pressure_system = 'pressure_system'
turbulence_systems = 'TKE_system TKED_system'
momentum_l_abs_tol = 1e-8
pressure_l_abs_tol = 1e-8
turbulence_l_abs_tol = 1e-8
momentum_l_tol = 1e-10
pressure_l_tol = 1e-10
turbulence_l_tol = 1e-10
momentum_equation_relaxation = 0.7
pressure_variable_relaxation = 0.3
turbulence_equation_relaxation = '0.4 0.4'
pressure_absolute_tolerance = 1e-8
momentum_absolute_tolerance = 1e-8
turbulence_absolute_tolerance = '1e-8 1e-8'
momentum_petsc_options_iname = '-pc_type -pc_hypre_type'
momentum_petsc_options_value = 'hypre boomeramg'
pressure_petsc_options_iname = '-pc_type -pc_hypre_type'
pressure_petsc_options_value = 'hypre boomeramg'
turbulence_petsc_options_iname = '-pc_type -pc_hypre_type'
turbulence_petsc_options_value = 'hypre boomeramg'
momentum_l_max_its = 300
pressure_l_max_its = 300
turbulence_l_max_its = 30
print_fields = false
num_iterations = 8000
continue_on_max_its = true
[]
[VectorPostprocessors]
[inlet_channel_wall_sampler]
type = LineValueSampler
variable = 'distance vel_x pressure mu_t_wall'
start_point = '${fparse -0.048} ${y_first_cell_1} 0.0'
end_point = '${fparse -1e-3} ${y_first_cell_1} 0.0'
sort_by = 'x'
num_points = 100
[]
[outlet_channel_wall_sampler]
type = LineValueSampler
variable = 'distance vel_x pressure mu_t_wall'
start_point = '6.1e-4 ${fparse -0.0127 + y_first_cell_2} 0.0'
end_point = '0.2492 ${fparse -0.0127 + y_first_cell_2} 0.0'
sort_by = 'x'
num_points = 100
[]
[vel_x_xoH_1_sampler]
type = LineValueSampler
variable = 'vel_x'
start_point = '${fparse 0.0127*1.0} -0.0127 0.0'
end_point = '${fparse 0.0127*1.0} ${fparse 0.0127*8.2-0.0127} 0.0'
sort_by = 'y'
num_points = 100
[]
[vel_x_xoH_4_sampler]
type = LineValueSampler
variable = 'vel_x'
start_point = '${fparse 0.0127*4.0} -0.0127 0.0'
end_point = '${fparse 0.0127*4.0} ${fparse 0.0127*8.2-0.0127} 0.0'
sort_by = 'y'
num_points = 100
[]
[vel_x_xoH_6_sampler]
type = LineValueSampler
variable = 'vel_x'
start_point = '${fparse 0.0127*6.0} -0.0127 0.0'
end_point = '${fparse 0.0127*6.0} ${fparse 0.0127*8.2-0.0127} 0.0'
sort_by = 'y'
num_points = 100
[]
[vel_x_xoH_10_sampler]
type = LineValueSampler
variable = 'vel_x'
start_point = '${fparse 0.0127*10.0} -0.0127 0.0'
end_point = '${fparse 0.0127*10.0} ${fparse 0.0127*8.2-0.0127} 0.0'
sort_by = 'y'
num_points = 100
[]
[]
[Postprocessors]
[mdot]
type = VolumetricFlowRate
vel_x = vel_x
vel_y = vel_y
advected_interp_method = ${advected_interp_method}
advected_quantity = ${rho}
boundary = 'inlet'
execute_on = 'TIMESTEP_END'
[]
[bulk_u]
type = SideExtremeValue
boundary = 'inlet'
variable = 'vel_x'
execute_on = 'TIMESTEP_END'
[]
[Reynolds_in]
type = ParsedPostprocessor
expression = 'bulk_u / mu * rho * D'
pp_names = 'bulk_u'
constant_expressions = '${mu} ${rho} ${D}'
constant_names = 'mu rho D'
execute_on = 'TIMESTEP_END'
[]
[]
[Outputs]
[exo]
type = Exodus
[]
[csv]
type = CSV
execute_on = FINAL
[]
[]
(test/tests/kEpsilon/channel/channel_ERCOFTAC.i)
##########################################################
# ERCOFTAC test case for turbulent channel flow
# Case Number: 032
# Author: Dr. Mauricio Tano & Hailey Tran Kieu
# Last Update: November, 2023 & 2025
# Turbulent model using:
# k-epsilon
# Equilibrium + Newton wall treatment
# SIMPLE solve
##########################################################
### Problem Parameters ###
H = 1 # half-width of the channel
L = 120
Re = 14000
rho = 1
bulk_u = 1
mu = '${fparse rho * bulk_u * 2 * H / Re}'
advected_interp_method = 'upwind'
### k-epsilon Closure Parameters ###
sigma_k = 1.0
sigma_eps = 1.3
C1_eps = 1.44
C2_eps = 1.92
C_mu = 0.09
### Initial and Boundary Conditions ###
intensity = 0.01
k_init = '${fparse 1.5*(intensity * bulk_u)^2}'
eps_init = '${fparse C_mu^0.75 * k_init^1.5 / (2*H)}'
### Modeling parameters ###
bulk_wall_treatment = false
walls = 'bottom top'
wall_treatment = 'eq_newton' # Options: eq_newton, eq_incremental, eq_linearized, neq
# NOTE:
# The k-epsilon kernels/aux-kernels now include strict parameter applicability
# checks. To keep this input file usable across all test permutations, we do
# *not* hard-code model-variant-specific options here. Instead, the test harness
# sets per-object options via CLI overrides (see create_tests.py).
[Mesh]
[block_1]
type = GeneratedMeshGenerator
dim = 2
xmin = 0
xmax = ${L}
ymin = 0
ymax = ${H}
nx = 10
ny = 5
bias_y = 0.7
[]
[block_2]
type = GeneratedMeshGenerator
dim = 2
xmin = 0
xmax = ${L}
ymin = ${fparse -H}
ymax = 0
nx = 10
ny = 5
bias_y = ${fparse 1/0.7}
[]
[smg]
type = StitchMeshGenerator
inputs = 'block_1 block_2'
clear_stitched_boundary_ids = true
stitch_boundaries_pairs = 'bottom top'
merge_boundaries_with_same_name = true
[]
# Prevent test diffing on distributed parallel element numbering
allow_renumbering = false
[]
[Problem]
linear_sys_names = 'u_system v_system pressure_system TKE_system TKED_system'
previous_nl_solution_required = true
[]
[GlobalParams]
rhie_chow_user_object = 'rc'
[]
[UserObjects]
[rc]
type = RhieChowMassFlux
u = vel_x
v = vel_y
pressure = pressure
rho = ${rho}
p_diffusion_kernel = p_diffusion
pressure_projection_method = 'consistent'
[]
[]
[Variables]
[vel_x]
type = MooseLinearVariableFVReal
initial_condition = ${bulk_u}
solver_sys = u_system
[]
[vel_y]
type = MooseLinearVariableFVReal
initial_condition = 0
solver_sys = v_system
[]
[pressure]
type = MooseLinearVariableFVReal
initial_condition = 1e-8
solver_sys = pressure_system
[]
[TKE]
type = MooseLinearVariableFVReal
solver_sys = TKE_system
initial_condition = ${k_init}
[]
[TKED]
type = MooseLinearVariableFVReal
solver_sys = TKED_system
initial_condition = ${eps_init}
[]
[]
[FVInterpolationMethods]
[upwind]
type = FVAdvectedUpwind
[]
[]
[LinearFVKernels]
[u_advection_stress]
type = LinearWCNSFVMomentumFlux
variable = vel_x
advected_interp_method_name = ${advected_interp_method}
mu = 'mu_t'
u = vel_x
v = vel_y
momentum_component = 'x'
rhie_chow_user_object = 'rc'
use_nonorthogonal_correction = false
use_deviatoric_terms = yes
[]
[u_diffusion]
type = LinearFVDiffusion
variable = vel_x
diffusion_coeff = '${mu}'
[]
[u_pressure]
type = LinearFVMomentumPressure
variable = vel_x
pressure = pressure
momentum_component = 'x'
[]
[v_advection_stress]
type = LinearWCNSFVMomentumFlux
variable = vel_y
advected_interp_method_name = ${advected_interp_method}
mu = 'mu_t'
u = vel_x
v = vel_y
momentum_component = 'y'
rhie_chow_user_object = 'rc'
use_nonorthogonal_correction = false
use_deviatoric_terms = yes
[]
[v_diffusion]
type = LinearFVDiffusion
variable = vel_y
diffusion_coeff = '${mu}'
[]
[v_pressure]
type = LinearFVMomentumPressure
variable = vel_y
pressure = pressure
momentum_component = 'y'
[]
[p_diffusion]
type = LinearFVPressureCorrectionDiffusion
variable = pressure
diffusion_tensor = Ainv
use_nonorthogonal_correction = false
[]
[HbyA_divergence]
type = LinearFVDivergence
variable = pressure
face_flux = HbyA
force_boundary_execution = true
[]
[TKE_advection]
type = LinearFVTurbulentAdvection
variable = TKE
advected_interp_method_name = ${advected_interp_method}
[]
[TKE_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKE
diffusion_coeff = ${mu}
use_nonorthogonal_correction = false
[]
[TKE_turb_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKE
diffusion_coeff = 'mu_t'
scaling_coeff = ${sigma_k}
use_nonorthogonal_correction = false
[]
[TKE_source_sink]
type = kEpsilonTKESourceSink
variable = TKE
u = vel_x
v = vel_y
epsilon = TKED
rho = ${rho}
mu = ${mu}
mu_t = 'mu_t'
walls = ${walls}
wall_treatment = ${wall_treatment}
C_pl = 1e10
[]
[TKED_advection]
type = LinearFVTurbulentAdvection
variable = TKED
advected_interp_method_name = ${advected_interp_method}
walls = ${walls}
[]
[TKED_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKED
diffusion_coeff = ${mu}
use_nonorthogonal_correction = false
walls = ${walls}
[]
[TKED_turb_diffusion]
type = LinearFVTurbulentDiffusion
variable = TKED
diffusion_coeff = 'mu_t'
scaling_coeff = ${sigma_eps}
use_nonorthogonal_correction = false
walls = ${walls}
[]
[TKED_source_sink]
type = kEpsilonTKEDSourceSink
variable = TKED
u = vel_x
v = vel_y
tke = TKE
rho = ${rho}
mu = ${mu}
mu_t = 'mu_t'
C1_eps = ${C1_eps}
C2_eps = ${C2_eps}
walls = ${walls}
wall_treatment = ${wall_treatment}
C_pl = 1e10
wall_distance = wall_distance # for low-Re / two-layer Yap / G' terms
[]
[]
[LinearFVBCs]
[inlet-u]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'left'
variable = vel_x
functor = '${bulk_u}'
[]
[inlet-v]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'left'
variable = vel_y
functor = '0.0'
[]
[walls-u]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'top bottom'
variable = vel_x
functor = 0.0
[]
[walls-v]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'top bottom'
variable = vel_y
functor = 0.0
[]
[outlet_u]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = 'right'
variable = vel_x
use_two_term_expansion = false
[]
[outlet_v]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = 'right'
variable = vel_y
use_two_term_expansion = false
[]
[outlet_p]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'right'
variable = pressure
functor = 0.0
[]
[inlet_TKE]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'left'
variable = TKE
functor = '${k_init}'
[]
[outlet_TKE]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = 'right'
variable = TKE
use_two_term_expansion = false
[]
[inlet_TKED]
type = LinearFVAdvectionDiffusionFunctorDirichletBC
boundary = 'left'
variable = TKED
functor = '${eps_init}'
[]
[outlet_TKED]
type = LinearFVAdvectionDiffusionOutflowBC
boundary = 'right'
variable = TKED
use_two_term_expansion = false
[]
[walls_mu_t]
type = LinearFVTurbulentViscosityWallFunctionBC
boundary = 'bottom top'
variable = 'mu_t'
u = vel_x
v = vel_y
rho = ${rho}
mu = ${mu}
tke = TKE
wall_treatment = ${wall_treatment}
[]
[]
[AuxVariables]
[wall_distance]
type = MooseVariableFVReal
initial_condition = 1.0
[]
[mu_t]
type = MooseLinearVariableFVReal
initial_condition = '${fparse rho * C_mu * ${k_init}^2 / eps_init}'
[]
[yplus]
type = MooseVariableFVReal
two_term_boundary_expansion = false
[]
[]
[AuxKernels]
[compute_wall_distance]
type = WallDistanceAux
variable = wall_distance
walls = ${walls}
execute_on = 'INITIAL NONLINEAR'
[]
[compute_mu_t]
type = kEpsilonViscosity
variable = mu_t
C_mu = ${C_mu}
tke = TKE
epsilon = TKED
mu = ${mu}
rho = ${rho}
u = vel_x
v = vel_y
bulk_wall_treatment = ${bulk_wall_treatment}
walls = ${walls}
wall_treatment = ${wall_treatment}
mu_t_ratio_max = 1e20
execute_on = 'NONLINEAR'
wall_distance = wall_distance # only needed for LowRe/TwoLayer (see below)
[]
[compute_y_plus]
type = RANSYPlusAux
variable = yplus
tke = TKE
mu = ${mu}
rho = ${rho}
u = vel_x
v = vel_y
walls = ${walls}
wall_treatment = ${wall_treatment}
execute_on = 'NONLINEAR'
[]
[]
[Executioner]
type = SIMPLE
rhie_chow_user_object = 'rc'
momentum_systems = 'u_system v_system'
pressure_system = 'pressure_system'
turbulence_systems = 'TKE_system TKED_system'
momentum_l_abs_tol = 1e-14
pressure_l_abs_tol = 1e-14
turbulence_l_abs_tol = 1e-14
momentum_l_tol = 1e-14
pressure_l_tol = 1e-14
turbulence_l_tol = 1e-14
momentum_equation_relaxation = 0.7
pressure_variable_relaxation = 0.3
turbulence_equation_relaxation = '0.2 0.2'
turbulence_field_relaxation = '0.2 0.2'
num_iterations = 1000
pressure_absolute_tolerance = 1e-7
momentum_absolute_tolerance = 1e-7
turbulence_absolute_tolerance = '1e-7 1e-7'
momentum_petsc_options_iname = '-pc_type -pc_hypre_type'
momentum_petsc_options_value = 'hypre boomeramg'
pressure_petsc_options_iname = '-pc_type -pc_hypre_type'
pressure_petsc_options_value = 'hypre boomeramg'
turbulence_petsc_options_iname = '-pc_type -pc_hypre_type'
turbulence_petsc_options_value = 'hypre boomeramg'
momentum_l_max_its = 300
pressure_l_max_its = 300
turbulence_l_max_its = 30
print_fields = false
continue_on_max_its = true
[]
[Outputs]
exodus = false
[csv]
type = CSV
execute_on = FINAL
[]
[]
variables_to_sample = 'vel_x vel_y pressure TKE TKED'
[VectorPostprocessors]
[side_bottom]
type = SideValueSampler
boundary = 'bottom'
variable = ${variables_to_sample}
sort_by = 'x'
execute_on = 'timestep_end'
[]
[side_top]
type = SideValueSampler
boundary = 'top'
variable = ${variables_to_sample}
sort_by = 'x'
execute_on = 'timestep_end'
[]
[line_center_channel]
type = LineValueSampler
start_point = '${fparse 0.125 * L} ${fparse 0.0001} 0'
end_point = '${fparse 0.875 * L} ${fparse 0.0001} 0'
num_points = ${Mesh/block_1/nx}
variable = ${variables_to_sample}
sort_by = 'x'
execute_on = 'timestep_end'
[]
[line_quarter_radius_channel]
type = LineValueSampler
start_point = '${fparse 0.125 * L} ${fparse 0.5 * H} 0'
end_point = '${fparse 0.875 * L} ${fparse 0.5 * H} 0'
num_points = ${Mesh/block_1/nx}
variable = ${variables_to_sample}
sort_by = 'x'
execute_on = 'timestep_end'
[]
[]
- Turbulent viscosity in the bulk region
- Near - wall bulk treatment
- Wall distance requirements
- Interaction with k and kernels var element = document.getElementById("moose-equation-6e9ddf69-ac2a-4481-a942-3307b4118b2f");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});
- Input Parameters
- Input Files
- References