LinearFVTemperatureWallFunctionBC

Implements a wall function boundary condition for the effective thermal conductivity at the wall for the linear finite volume discretization.

This boundary condition computes an effective wall thermal conductivity used to represent turbulent heat transfer near solid boundaries without fully resolving the thermal boundary layer. As in LinearFVTurbulentViscosityWallFunctionBC, the near-wall behavior is parameterized in terms of the non-dimensional wall distance $var element = document.getElementById("moose-equation-d347e4cd-bb7d-4f95-a31a-cb778557ec39");katex.render("y^+", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ and the friction velocity $var element = document.getElementById("moose-equation-0a8ba1b3-f323-4713-9021-c1d429711ba7");katex.render("u_{\\tau}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ , using the selected "wall_treatment".

The wall heat flux is modeled through an effective conductivity, $var element = document.getElementById("moose-equation-cbd68f7d-e242-4a54-b1fa-211674bb25f8");katex.render("k_{eff,w}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ , such that

$var element = document.getElementById("moose-equation-0cc4fd84-17fd-472b-b74f-32a2d9342b18");katex.render(" q_w = \\frac{k_{eff,w}}{y_p} \\left( T_p - T_w \\right)\\,,", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$

where:

$var element = document.getElementById("moose-equation-fb6fcf05-ab8d-4d87-9668-9b262d468278");katex.render("q_w", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the wall heat flux,
$var element = document.getElementById("moose-equation-54276d85-10cc-4257-a55c-0450add9733d");katex.render("k_{t,w}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the effective wall thermal conductivity returned by this boundary condition,
$var element = document.getElementById("moose-equation-78ff7943-ae58-46b8-9201-15b1e7122aed");katex.render("T_p", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the temperature at the centroid of the fluid near-wall cell,
$var element = document.getElementById("moose-equation-731b4e46-79ad-4797-b926-73f9c285e09e");katex.render("T_w", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the wall temperature,
$var element = document.getElementById("moose-equation-9c6316d6-b71a-4a96-9a3d-5fdecb5b77f6");katex.render("y_p", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the wall-normal distance from the wall face to the centroid of the near-wall cell.

The wall function computes this effective conductivity as

$var element = document.getElementById("moose-equation-74e6f3ff-4853-4820-a73c-78660f01c5b9");katex.render(" k_{eff,w} = \\frac{\\rho c_p u_{\\tau} y_p}{T^+}\\,,", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$

where:

$var element = document.getElementById("moose-equation-d305af21-8c5c-493e-98c5-bb9d1494e277");katex.render("\\rho", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the density
$var element = document.getElementById("moose-equation-f45f908b-a8b2-452c-8124-3d1890241b89");katex.render("c_p", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the specific heat capacity
$var element = document.getElementById("moose-equation-f40c01f6-dd03-43bb-8149-ddfbd1eeaa0e");katex.render("u_{\\tau}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the friction velocity
$var element = document.getElementById("moose-equation-97e44395-1272-4e77-9a0a-1cd87f87280b");katex.render("T^+", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the dimensionless temperature

The molecular Prandtl number is defined as

$var element = document.getElementById("moose-equation-32b6dca3-866e-4267-907e-ca769d063ec1");katex.render(" \\mathrm{Pr} = \\frac{c_p \\mu}{k}\\,,", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$

where:

$var element = document.getElementById("moose-equation-93e629f5-edb6-4001-b398-49419910e755");katex.render("\\mu", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the molecular dynamic viscosity
$var element = document.getElementById("moose-equation-2aa47dac-e315-4fa2-898f-0e52ea77ed54");katex.render("k", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the molecular thermal conductivity

The turbulent Prandtl number $var element = document.getElementById("moose-equation-d12a6363-edbe-4cba-962b-1f223b948c31");katex.render("\\mathrm{Pr}_t", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is provided through the material property "Pr_t". As in the turbulent viscosity wall function, four different formulations are supported to compute $var element = document.getElementById("moose-equation-aa55b108-272d-4da0-a49b-e13c39a42658");katex.render("u_{\\tau}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ and $var element = document.getElementById("moose-equation-1be746bf-9095-476d-80dc-6e4e8ac7ca78");katex.render("y^+", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ , as defined by the "wall_treatment".

The same recommendations about mesh design and wall-normal spacing given in LinearFVTurbulentViscosityWallFunctionBC apply here as well.

Thermal wall-function formulation

Once $var element = document.getElementById("moose-equation-7c479afc-e9e1-44b8-9646-bc22c58e0b32");katex.render("u_{\\tau}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ and $var element = document.getElementById("moose-equation-9d7b85fc-10e6-4b6d-86a9-184a57eae6bb");katex.render("y^+", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ are computed, the thermal wall function evaluates the dimensionless temperature $var element = document.getElementById("moose-equation-428c90e0-3403-4ade-86f3-a6343fb5c2aa");katex.render("T^+", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ using the thermal sublayer thickness $var element = document.getElementById("moose-equation-325aa573-533d-4dc9-9d86-cf40be8c112f");katex.render("y^+_{t}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ , obtained from a fixed point solve through NS::findYplusThermal, which is dependent on $var element = document.getElementById("moose-equation-719148d2-994c-44ce-af9b-ae683c43de84");katex.render("\\mathrm{Pr}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ and $var element = document.getElementById("moose-equation-da6b4550-83c5-4f82-be19-0c16bcd2da81");katex.render("\\mathrm{Pr}_t", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ .

The thermal sublayer thickness $var element = document.getElementById("moose-equation-23ea0b0f-bece-4ae5-8a14-a6d2a30bc160");katex.render("y^+_{t}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is then used as the parameter to discern whether to employ laminar sublayer thermal wall functions or turbulent logarithmic layer thermal wall functions:

$var element = document.getElementById("moose-equation-a198400a-7841-4ab9-b19a-da48d6b1135c");katex.render(" T^+ = \\begin{cases} y^+ \\mathrm{Pr} & \\text{if } y^+ \\le y^+_{t} \\\\ \\mathrm{Pr}_t \\left[ \\frac{1}{\\kappa}\\ln(E y^+) + P \\right] & \\text{if } y^+ > y^+_{t}, \\end{cases}", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$

with

$var element = document.getElementById("moose-equation-b68160b4-0839-4787-bce0-37f16af010a0");katex.render("P = 9.24 \\left[ \\left( \\frac{\\mathrm{Pr}}{\\mathrm{Pr}_t} \\right)^{0.75} - 1 \\right] \\left( 1 + 0.28 \\exp\\left[-0.007 \\frac{\\mathrm{Pr}}{\\mathrm{Pr}_t}\\right] \\right),", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$

where:

$var element = document.getElementById("moose-equation-72d5d5a0-2686-42c3-9735-cb7aa22dbf46");katex.render("\\kappa = 0.4187", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the von Karman constant,
$var element = document.getElementById("moose-equation-f9e23eff-5925-4335-81eb-1426dde1c7cc");katex.render("E = 9.793", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the wall-function log law offset constant,
$var element = document.getElementById("moose-equation-46fbad89-9c5e-4821-8118-758082481c2a");katex.render("P", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the Jayatilleke temperature offset

The returned boundary value is then

$var element = document.getElementById("moose-equation-79af2593-2b2a-45c0-b440-5add73fb5167");katex.render(" k_{eff,w} = \\frac{\\rho c_p u_{\\tau} y_p}{T^+}.", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$

note

This wall function is formulated for molecular Prandtl numbers of order 1 or greater, as is typical for gases and most liquids. For low-Prandtl-number fluids, such as liquid metals, $var element = document.getElementById("moose-equation-93fbfcc8-29e1-4570-bc48-2e7bb5826011");katex.render("y^+_{t}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ can grow well past the validity range of the underlying equilibrium velocity wall function (i.e. into or beyond the logarithmic layer), making the sub-laminar/logarithmic switch and the resulting $var element = document.getElementById("moose-equation-d0bafcae-7656-4e44-91c0-1e6fc415be5c");katex.render("k_{eff,w}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ unreliable. A dedicated low-Prandtl-number thermal wall treatment is recommended for such fluids.

Input Parameters

boundaryThe list of boundary IDs from the mesh where this object applies
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The list of boundary IDs from the mesh where this object applies
cpThe specific heat at constant pressure. 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 specific heat at constant pressure. 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.
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 boundary condition applies to
C++ Type:LinearVariableName
Unit:(no unit assumed)
Controllable:No
Description:The name of the variable that this boundary condition applies to

Required Parameters

C_mu0.09Coupled turbulent kinetic energy closure.
Default:0.09
C++ Type:Real
Unit:(no unit assumed)
Controllable:No
Description:Coupled turbulent kinetic energy closure.
Pr_t0.7The turbulent Prandtl number. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
Default:0.7
C++ Type:MooseFunctorName
Unit:(no unit assumed)
Controllable:No
Description:The turbulent Prandtl number. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
kThe thermal conductivity. 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 thermal conductivity. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
tkeThe 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:The turbulent kinetic energy. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
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_treatmentneqThe method used for computing the wall functions
Default:neq
C++ Type:MooseEnum
Options:eq_newton, eq_incremental, eq_linearized, neq
Controllable:No
Description:The method used for computing the wall functions

Optional Parameters

absolute_value_vector_tagsThe tags for the vectors this residual object should fill with the absolute value of the residual contribution
C++ Type:std::vector<TagName>
Controllable:No
Description:The tags for the vectors this residual object should fill with the absolute value of the residual contribution
extra_matrix_tagsThe extra tags for the matrices this Kernel should fill
C++ Type:std::vector<TagName>
Controllable:No
Description:The extra tags for the matrices this Kernel should fill
extra_vector_tagsThe extra tags for the vectors this Kernel should fill
C++ Type:std::vector<TagName>
Controllable:No
Description:The extra tags for the vectors this Kernel should fill
matrix_onlyFalseWhether this object is only doing assembly to matrices (no vectors)
Default:False
C++ Type:bool
Controllable:No
Description:Whether this object is only doing assembly to matrices (no vectors)
matrix_tagssystemThe tag for the matrices this Kernel should fill
Default:system
C++ Type:MultiMooseEnum
Options:nontime, system
Controllable:No
Description:The tag for the matrices this Kernel should fill
vector_tagsrhsThe tag for the vectors this Kernel should fill
Default:rhs
C++ Type:MultiMooseEnum
Options:rhs, time
Controllable:No
Description:The tag for the vectors this Kernel should fill

Contribution To Tagged Field Data 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.
implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
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).

Advanced Parameters

Input Files

(modules/navier_stokes/test/tests/finite_volume/ins/turbulence/heated-channel/channel_ERCOFTAC_heated.i)

(modules/navier_stokes/test/tests/finite_volume/ins/turbulence/heated-channel/channel_ERCOFTAC_heated.i)

# Turbulent flow in a heated channel where fluid enters cold and is heated by hot walls.
# Momentum and heat transport are modeled using a k-epsilon turbulence model with
# wall functions for both velocity and temperature.

H = 1 #halfwidth of the channel
L = 100

Re = 13700

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 = 'top bottom'
wall_treatment = 'eq_newton'  # Options: eq_newton, eq_incremental, eq_linearized, neq

### Energy variables for Pr = 1
temperature_init = 300.
cp = 4000.
k = 0.58
Pr_t = 0.7


[Mesh]
  [block_1]
    type = CartesianMeshGenerator
    dim = 2
    dx = '${L}'
    dy = '0.75 0.25'
    ix = '5'
    iy = '3 1'
  []
  [block_2_base]
    type = CartesianMeshGenerator
    dim = 2
    dx = '${L}'
    dy = '0.25 0.75'
    ix = '5'
    iy = '1 3'
  []
  [block_2]
    type = TransformGenerator
    input = block_2_base
    transform = TRANSLATE
    vector_value = '0 -1 0'
  []
  [smg]
    type = StitchedMeshGenerator
    inputs = 'block_1 block_2'
    clear_stitched_boundary_ids = true
    stitch_boundaries_pairs = 'bottom top'
    merge_boundaries_with_same_name = true
  []
  # Prevent sampling diffing on distributed parallel element numbering
  allow_renumbering = false
[]

[Problem]
  linear_sys_names = 'u_system v_system pressure_system TKE_system TKED_system energy_system'
  previous_nl_solution_required = true
[]

[GlobalParams]
  rhie_chow_user_object = 'rc'
  advected_interp_method = ${advected_interp_method}
[]

[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 = ${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 = ${fparse 3*k_init}
  []
  [TKED]
    type = MooseLinearVariableFVReal
    solver_sys = TKED_system
    initial_condition = ${fparse 3*eps_init}
  []
  [temperature]
    type = MooseLinearVariableFVReal
    solver_sys = energy_system
    initial_condition = ${temperature_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}
    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 = 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 = ${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 = LinearFVTKESourceSink
    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 = 10
  []

  [TKED_advection]
    type = LinearFVTurbulentAdvection
    variable = TKED
    walls = ${walls}
    advected_interp_method_name = ${advected_interp_method}
  []
  [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 = LinearFVTKEDSourceSink
    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 = 10
  []

  [temp_advection]
    type = LinearFVEnergyAdvection
    variable = temperature
    advected_quantity = temperature
    cp = ${cp}
  []
  [temp_diffusion]
    type = LinearFVDiffusion
    variable = temperature
    diffusion_coeff = 'k_eff'
  []
[]

[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'
  []
  [inlet-temp]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = temperature
    functor = '300.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
  []
  [walls_temp]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'top bottom'
    variable = temperature
    functor = 400.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
  []
  [outlet_temp]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'right'
    variable = temperature
    use_two_term_expansion = false
  []

  [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}
  []
  [walls_k_t]
    type = LinearFVTemperatureWallFunctionBC
    boundary = 'bottom top'
    variable = 'k_eff'
    u = vel_x
    v = vel_y
    rho = ${rho}
    mu = ${mu}
    cp = ${cp}
    k = ${k}
    Pr_t = ${Pr_t}
    tke = TKE
    wall_treatment = ${wall_treatment}
  []
[]

[AuxVariables]
  [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}'
  []
  [k_eff]
    type = MooseLinearVariableFVReal
  []
  [k_t]
    type = MooseLinearVariableFVReal
  []
[]

[AuxKernels]
  [compute_mu_t]
    type = kEpsilonViscosityAux
    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}
    execute_on = 'NONLINEAR'
    mu_t_ratio_max = 1e20
  []
  [compute_k_eff]
    type = EffectiveTurbulentThermalConductivityAux
    variable = k_eff
    Pr_t = ${Pr_t}
    cp = ${cp}
    k = ${k}
    mu_t = 'mu_t'
    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'
  []
  [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'
  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.35 0.35'
  turbulence_field_relaxation = '0.35 0.35'
  energy_equation_relaxation = 0.95
  num_iterations = 4000
  pressure_absolute_tolerance = 1e-12
  momentum_absolute_tolerance = 1e-12
  turbulence_absolute_tolerance = '1e-12 1e-12'
  energy_absolute_tolerance = 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'
  energy_petsc_options_iname = '-pc_type -pc_hypre_type'
  energy_petsc_options_value = 'hypre boomeramg'

  print_fields = false
  continue_on_max_its = true
[]

[Outputs]
  exodus = true
  csv = true
[]

# sample along walls and in the bulk
[VectorPostprocessors]
  [side_bottom]
    type = SideValueSampler
    boundary = 'bottom'
    variable = 'vel_x vel_y pressure temperature TKE TKED mu_t yplus mu_eff k_eff k_t'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
  [side_top]
    type = SideValueSampler
    boundary = 'top'
    variable = 'vel_x vel_y pressure temperature TKE TKED mu_t yplus mu_eff k_eff k_t'
    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/ix}
    variable = 'vel_x vel_y pressure temperature TKE TKED mu_t yplus mu_eff k_eff k_t'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
  [line_quarter_radius_channel]
    type = LineValueSampler
    start_point = '${fparse 0.125 * L} ${fparse 0.4999 * H} 0'
    end_point = '${fparse 0.875 * L} ${fparse 0.4999 * H} 0'
    num_points =  ${Mesh/block_1/ix}
    variable = 'vel_x vel_y pressure temperature TKE TKED mu_t yplus mu_eff k_eff k_t'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
[]

Thermal wall-function formulation
Input Parameters
Input Files