kEpsilonTKESourceSink

kEpsilonTKESourceSink is a finite volume elemental kernel that computes the turbulent source and sink term for the turbulent kinetic energy equation in the k– $var element = document.getElementById("moose-equation-af8d9e7f-de13-48c3-825a-4f07c969e4d5");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ family of models.

note

The explanations in this kernel documentation are a summary. The reader is referred to the theory for more details if needed.

The turbulent kinetic energy (TKE) equation is written in conservative form as

$var element = document.getElementById("moose-equation-8e9377c3-f76e-40ce-b39a-7fbaa778d8a7");katex.render("\\frac{\\partial \\rho k}{\\partial t} + \\nabla \\cdot (\\rho \\mathbf{u} k) = \\rho (P_k + G_b + G_{\\text{nl}} - \\gamma_M - \\epsilon),", element, {displayMode:true,throwOnError:false});$

where:

$var element = document.getElementById("moose-equation-dade8072-9691-4202-9f1c-2b60c7297c9f");katex.render("k", element, {displayMode:false,throwOnError:false});$ is the turbulent kinetic energy,
$var element = document.getElementById("moose-equation-f052a23a-7157-4ea2-b596-12396a4e4224");katex.render("\\mathbf{u}", element, {displayMode:false,throwOnError:false});$ is the mean velocity,
$var element = document.getElementById("moose-equation-ec55567b-de4a-4b17-a849-68ae6bbafd3a");katex.render("P_k", element, {displayMode:false,throwOnError:false});$ is the shear production,
$var element = document.getElementById("moose-equation-a9997436-3137-47d7-9f29-2b01e4faf848");katex.render("G_b", element, {displayMode:false,throwOnError:false});$ is the buoyancy production,
$var element = document.getElementById("moose-equation-022160e3-4df5-4323-b393-f248cea48439");katex.render("G_{\\text{nl}}", element, {displayMode:false,throwOnError:false});$ is the production due to nonlinear Reynolds stresses,
$var element = document.getElementById("moose-equation-c57e05f1-37ee-468b-9574-a5e5bb1a0fc7");katex.render("\\gamma_M", element, {displayMode:false,throwOnError:false});$ is a compressibility correction, and
$var element = document.getElementById("moose-equation-b9a7d4cc-3c72-4f59-96f5-519c39842566");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ is the dissipation rate supplied by the $var element = document.getElementById("moose-equation-71beabf4-9d3b-4d07-a9fb-4505d2d2f505");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ -equation.

kEpsilonTKESourceSink forms the right-hand side source term for the k-equation by computing the combination

$var element = document.getElementById("moose-equation-98ba3cd2-3a12-4010-bde6-e0047b160f3e");katex.render("\\text{source} = P_k + G_b + G_{\\text{nl}} - \\gamma_M - \\epsilon,", element, {displayMode:true,throwOnError:false});$

with expressions that depend on the chosen k– $var element = document.getElementById("moose-equation-6609d7e9-400f-45b1-bcb2-f5ee2938611b");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ variant and options.

The kernel is designed to be used together with:

kEpsilonTKEDSourceSink for the $var element = document.getElementById("moose-equation-e275a592-d9e2-4d21-8a82-a6d8de22a0a9");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ -equation, and
kEpsilonViscosity for the turbulent viscosity $var element = document.getElementById("moose-equation-17231e69-d722-4c9d-ab64-fc498237f446");katex.render("\\mu_t", element, {displayMode:false,throwOnError:false});$ .

Model variants

The kernel supports several members of the k– $var element = document.getElementById("moose-equation-63753a11-f4e3-42fe-b8e0-89b0f00e7724");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ family via the "k_epsilon_variant" parameter:

Standard — classical high-Reynolds-number k– $var element = document.getElementById("moose-equation-e6dd33cf-9a38-4b5e-8dae-a96acfa7521f");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ model.
StandardLowRe — low-Reynolds-number k– $var element = document.getElementById("moose-equation-56d911f0-5918-450a-8c8b-90e2bf6a0289");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ model with damping functions.
StandardTwoLayer — two-layer k– $var element = document.getElementById("moose-equation-d75547bc-b6a0-4c5e-bdcd-e58db27eb85e");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ formulation using a blending between outer k– $var element = document.getElementById("moose-equation-b62d9c6b-8d8c-4990-b9cd-8bfbbd2e2e05");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ and near-wall length scales.
Realizable — realizable k– $var element = document.getElementById("moose-equation-ea9bfd4c-6a9e-4982-b06b-c61ce2cfb552");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ model with variable $var element = document.getElementById("moose-equation-af54fba1-0aba-48d9-b872-de6ef27b6b34");katex.render("C_\\mu", element, {displayMode:false,throwOnError:false});$ depending on strain/rotation invariants.
RealizableTwoLayer — realizable model with a two-layer near-wall treatment.

The choice of variant affects:

how the shear production is formed,
whether a curvature correction factor is applied,
how near-wall behavior is modeled (two-layer vs. high-Re),
whether low-Re damping functions and extra terms are used (via $var element = document.getElementById("moose-equation-8d2738fb-2029-4ff9-8c3a-ac49f411fa4a");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ -equation).

Bulk formulation

In bulk cells (not explicitly treated as near-wall by this kernel), the source term is computed as

$var element = document.getElementById("moose-equation-bdedcb03-b91d-4fed-b289-5997d92b5bf0");katex.render("\\text{source}_{\\text{bulk}} = P_k + G_b + G_{\\text{nl}} - \\gamma_M - \\epsilon.", element, {displayMode:true,throwOnError:false});$

The individual contributions are detailed below.

Strain-rate tensor and invariants

The mean velocity gradient $var element = document.getElementById("moose-equation-684cc433-c11d-4688-855f-e19362514f39");katex.render("\\nabla \\mathbf{u}", element, {displayMode:false,throwOnError:false});$ is used to form the symmetric strain-rate tensor $var element = document.getElementById("moose-equation-3916e941-e065-4d0d-89a2-281f7cdcc94e");katex.render("\\mathbf{S}", element, {displayMode:false,throwOnError:false});$ and the antisymmetric rotation tensor $var element = document.getElementById("moose-equation-be003abb-bb97-474a-8e3f-ab41e4db2094");katex.render("\\mathbf{W}", element, {displayMode:false,throwOnError:false});$ :

$var element = document.getElementById("moose-equation-53fba840-6900-4a93-875f-f1a84f8ac55f");katex.render("\\mathbf{S} = \\frac{1}{2} \\left(\\nabla \\mathbf{u} + \\nabla \\mathbf{u}^T \\right)", element, {displayMode:true,throwOnError:false});$

$var element = document.getElementById("moose-equation-a91b12bc-de97-47cd-a31f-738c12a71de4");katex.render("\\mathbf{W} = \\frac{1}{2} \\left(\\nabla \\mathbf{u} - \\nabla \\mathbf{u}^T \\right)", element, {displayMode:true,throwOnError:false});$

The invariants are defined as:

$var element = document.getElementById("moose-equation-860264ca-54f1-4eff-aa2f-034514655888");katex.render("S^2 = 2 S_{ij} S_{ij}", element, {displayMode:false,throwOnError:false});$ ,
$var element = document.getElementById("moose-equation-0b2b0082-41ae-422f-bff2-9497f54eab39");katex.render("W^2 = 2 W_{ij} W_{ij}", element, {displayMode:false,throwOnError:false});$ ,
$var element = document.getElementById("moose-equation-346eba33-0f98-4500-8ac2-4438db3dce7e");katex.render("\\nabla \\cdot \\mathbf{u}", element, {displayMode:false,throwOnError:false});$ .

These invariants are shared by the turbulence viscosity, nonlinear models, and curvature correction.

Shear production $var element = document.getElementById("moose-equation-26504802-ec99-42fb-bdc0-ae04228526fd");katex.render("G_k", element, {displayMode:false,throwOnError:false});$

The base formulation of the shear production is

$var element = document.getElementById("moose-equation-64356ad5-2c0c-4771-88e1-941659c14ebb");katex.render("G_k = \\mu_t S^2,", element, {displayMode:true,throwOnError:false});$

where

$var element = document.getElementById("moose-equation-34222078-f700-4d15-8fff-9bfc8b75ce3f");katex.render("\\mu_t", element, {displayMode:false,throwOnError:false});$ is the turbulent viscosity provided as a functor, and
$var element = document.getElementById("moose-equation-d4c18932-77b0-4df5-91ed-9f54f2f69038");katex.render("S = \\sqrt{2 \\mathbf{S}:\\mathbf{S}} = \\sqrt{S^2}", element, {displayMode:false,throwOnError:false});$ .

For compressible flows, the implementation allows inclusion of the extra trace terms

$var element = document.getElementById("moose-equation-c5a68ac1-1e84-4ec8-9bc1-5dcb4f016e45");katex.render("G_k = \\mu_t S^2 - \\frac{2}{3}\\left(\\rho k \\nabla \\cdot \\mathbf{u} + \\mu_t (\\nabla\\cdot\\mathbf{u})^2\\right),", element, {displayMode:true,throwOnError:false});$

consistent with the constitutive relation for the Reynolds stress tensor in compressible flow.

The shear production is then modified depending on the k– $var element = document.getElementById("moose-equation-9907e18b-4667-4657-ad29-17d3c46090c5");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ variant and curvature correction model (see below).

Curvature / rotation correction $var element = document.getElementById("moose-equation-8ff96247-5c6a-46ae-ba43-20184f91776b");katex.render("f_c", element, {displayMode:false,throwOnError:false});$

For realizable variants, an optional curvature correction factor $var element = document.getElementById("moose-equation-6fd750b2-1362-49d1-84e1-12e61c3c1cdd");katex.render("f_c", element, {displayMode:false,throwOnError:false});$ may be applied via "curvature_model". The current implementation supports:

curvature_model = none — no curvature correction, $var element = document.getElementById("moose-equation-813e72fd-67c4-47de-9768-53db8848bf9d");katex.render("f_c = 1", element, {displayMode:false,throwOnError:false});$ ;
curvature_model = standard — Spalart–Shur–type rotation/curvature correction based on $var element = document.getElementById("moose-equation-29a98f7b-38ec-407c-8929-10691c268009");katex.render("S^2", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-ea73f0e9-164a-4d1e-a8bd-03d692ad3650");katex.render("W^2", element, {displayMode:false,throwOnError:false});$ .

The curvature correction builds a correction $var element = document.getElementById("moose-equation-7f5583e7-2253-4596-868f-38f0ed1e4723");katex.render("f_c = f_c(S^2, W^2)", element, {displayMode:false,throwOnError:false});$ that enhances production in stabilizing curvature and reduces it in destabilizing regions. The corrected shear production is

Standard, StandardLowRe, StandardTwoLayer: $var element = document.getElementById("moose-equation-caad84ca-3c71-4d4d-9508-47d51586c53a");katex.render("P_k = G_k;", element, {displayMode:true,throwOnError:false});$
Realizable, RealizableTwoLayer: $var element = document.getElementById("moose-equation-0557b0fb-379f-4434-ac21-74298250407e");katex.render("P_k = f_c G_k.", element, {displayMode:true,throwOnError:false});$

Buoyancy production $var element = document.getElementById("moose-equation-4785b38c-edca-44d3-8911-8a1dafebb261");katex.render("G_b", element, {displayMode:false,throwOnError:false});$

If "use_buoyancy" is true, buoyancy production is added via

$var element = document.getElementById("moose-equation-95acfe33-fbe4-474d-a6c8-7bee3bf607a1");katex.render("G_b = \\beta \\frac{\\mu_t}{Pr_t} \\left(\\nabla T \\cdot \\mathbf{g}\\right),", element, {displayMode:true,throwOnError:false});$

where

$var element = document.getElementById("moose-equation-b14f5408-d47e-45c2-a66d-7957fed43e9b");katex.render("\\beta", element, {displayMode:false,throwOnError:false});$ is the thermal expansion coefficient,
$var element = document.getElementById("moose-equation-993445b9-0931-441f-b6ce-73acd1833fbf");katex.render("Pr_t", element, {displayMode:false,throwOnError:false});$ is the turbulent Prandtl number,
$var element = document.getElementById("moose-equation-f5faa66a-bd79-4ca7-8f00-f42a300327d5");katex.render("T", element, {displayMode:false,throwOnError:false});$ is the temperature,
$var element = document.getElementById("moose-equation-51c059ad-9741-4125-b92a-4c5c58e2d104");katex.render("\\mathbf{g}", element, {displayMode:false,throwOnError:false});$ is the gravitational acceleration.

The gradient $var element = document.getElementById("moose-equation-38205f56-1bc3-4449-9be9-168f0dcf45e5");katex.render("\\nabla T", element, {displayMode:false,throwOnError:false});$ and gravity vector are supplied as functors / parameters. Buoyancy can either augment or reduce total production depending on the sign of the dot product $var element = document.getElementById("moose-equation-258e6792-da2d-4842-a838-8b8f31b6f528");katex.render("\\nabla T\\cdot\\mathbf{g}", element, {displayMode:false,throwOnError:false});$ .

nonlinear production $var element = document.getElementById("moose-equation-ea140f4a-a26b-438b-926f-877c55a6f8c1");katex.render("G_{\\text{nl}}", element, {displayMode:false,throwOnError:false});$

kEpsilonTKESourceSink supports optional nonlinear constitutive relations for the Reynolds stress tensor. The choice is controlled by "nonlinear_model":

nonlinear_model = none — no nonlinear production, $var element = document.getElementById("moose-equation-d2d41ecd-5bba-436e-a990-2ba9dea51025");katex.render("G_{\\text{nl}} = 0", element, {displayMode:false,throwOnError:false});$ ;
nonlinear_model = quadratic — quadratic constitutive relation;
nonlinear_model = cubic — quadratic + cubic terms.

The nonlinear Reynolds stress contribution is computed in TurbulenceMethods as

$var element = document.getElementById("moose-equation-021e1b01-d6bd-4d26-ab9b-d43a298f7f09");katex.render("G_{\\text{nl}} = \\mathbf{T}_{\\text{RANS,NL}} : \\nabla \\mathbf{u},", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-b377ba90-8b2f-4ffb-942d-a3d3068a5504");katex.render("\\mathbf{T}_{\\text{RANS,NL}}", element, {displayMode:false,throwOnError:false});$ depends on:

the velocity gradient $var element = document.getElementById("moose-equation-3425200a-2a65-4b2f-b2a9-5b5fb9d39d49");katex.render("\\nabla \\mathbf{u}", element, {displayMode:false,throwOnError:false});$ (through $var element = document.getElementById("moose-equation-38b9f249-b85b-45de-b24b-9f601dee7d84");katex.render("\\mathbf{S}", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-35f37267-2f31-4e87-8b58-4d2dd7f971c8");katex.render("\\mathbf{W}", element, {displayMode:false,throwOnError:false});$ ),
the invariants $var element = document.getElementById("moose-equation-3daf1ac0-fd52-4689-8cf6-3154de21c09a");katex.render("S^2", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-015df84b-2446-4345-969a-86e13eb8601a");katex.render("W^2", element, {displayMode:false,throwOnError:false});$ ,
$var element = document.getElementById("moose-equation-ae36335f-0f07-4e79-bb2b-17a5e234a955");katex.render("\\mu_t", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-bd788970-440c-4e84-adef-69a5000ea1f0");katex.render("k", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-b0e41832-ee2d-481e-9049-44f1ca19198d");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ ,
fixed model coefficients $var element = document.getElementById("moose-equation-34e34fc9-d181-4206-8763-877a571e694c");katex.render("C_{\\text{NL}1\\ldots 7}", element, {displayMode:false,throwOnError:false});$ and a variable $var element = document.getElementById("moose-equation-b5df0a19-42f6-4a09-9cbc-1ec36c9673af");katex.render("C_\\mu", element, {displayMode:false,throwOnError:false});$ .

In all k– $var element = document.getElementById("moose-equation-8b331e42-af79-40a7-a9b8-96675b39c059");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ variants, kEpsilonTKESourceSink adds $var element = document.getElementById("moose-equation-b769607e-a8cd-40b4-9394-eb607818fd31");katex.render("G_{\\text{nl}}", element, {displayMode:false,throwOnError:false});$ to the k-equation’s production:

$var element = document.getElementById("moose-equation-daec48f9-fb7e-4f42-9153-5be790bd925d");katex.render("\\text{source} = (P_k + G_b + G_{\\text{nl}} - \\gamma_M) - \\epsilon.", element, {displayMode:true,throwOnError:false});$

Compressibility correction $var element = document.getElementById("moose-equation-d10bdada-b632-4336-a844-5609bc75a1e4");katex.render("\\gamma_M", element, {displayMode:false,throwOnError:false});$

If "use_compressibility" is true, a compressibility correction is added to the sink term:

$var element = document.getElementById("moose-equation-1847687e-0f93-4dd6-8ff0-87d5e3d5a569");katex.render("\\gamma_M = \\rho C_M \\frac{k \\epsilon}{c^2},", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-4b1de6a4-2e41-4f08-bd78-89a0fb389f8c");katex.render("C_M", element, {displayMode:false,throwOnError:false});$ is a model coefficient and $var element = document.getElementById("moose-equation-d7a95e16-3eb7-457a-b472-ce3d34bb2e1f");katex.render("c", element, {displayMode:false,throwOnError:false});$ is the speed of sound. This term is computed by NS::computeGammaM and is subtracted from total production, representing the transfer of TKE into acoustic modes in compressible flow.

Production limiter

To avoid excessive TKE generation in stagnation regions, a production limiter is applied:

$var element = document.getElementById("moose-equation-e35800f3-1302-4a77-9448-5169cbbd23b1");katex.render("G_k \\le C_{PL}\\,\\rho\\,\\epsilon.", element, {displayMode:true,throwOnError:false});$

The limiter coefficient $var element = document.getElementById("moose-equation-ecad3762-1255-4b62-b585-d932c45c35a5");katex.render("C_{PL}", element, {displayMode:false,throwOnError:false});$ is available as "C_pl". This limiter follows the standard approach from Menter (1994) and is applied in the bulk region.

Summary of bulk source

In non-wall-bounded cells, the kernel assembles the following source term for the residual:

$var element = document.getElementById("moose-equation-bda2e2d7-6a27-40fd-817a-a598083cba5e");katex.render("\\text{source}_{\\text{bulk}} = \\underbrace{P_k}_{\\text{shear (varies by variant)}} + \\underbrace{G_b}_{\\text{buoyancy}} + \\underbrace{G_{\\text{nl}}}_{\\text{nonlinear model}} - \\underbrace{\\gamma_M}_{\\text{compressibility}} - \\underbrace{\\epsilon}_{\\text{dissipation}}.", element, {displayMode:true,throwOnError:false});$

Wall formulation

Cells adjacent to boundaries listed in "walls" are treated as next-to-wall cells with a different expression for production and destruction, similar in spirit to LinearFVTKESourceSink. Wall distances are typically computed using WallDistanceAux, which returns the distance from each cell centroid to the nearest wall.

The near-wall treatment is implemented in the matrix contribution for this kernel, where a distinction is made between the viscous sublayer and the logarithmic region based on the non-dimensional wall distance

$var element = document.getElementById("moose-equation-318ee2a8-d174-4c41-bdbd-68256233385d");katex.render("y^+ = \\frac{\\rho y_p u_\\tau}{\\mu},", element, {displayMode:true,throwOnError:false});$

with

$var element = document.getElementById("moose-equation-857dfc60-d653-49b6-b9fb-3ccc6cf3f186");katex.render("y_p", element, {displayMode:false,throwOnError:false});$ the distance from the cell center to the wall,
$var element = document.getElementById("moose-equation-1477a874-e44c-47fc-8aff-0b1e5e657f4f");katex.render("u_\\tau", element, {displayMode:false,throwOnError:false});$ the friction velocity obtained from a wall-function model,
$var element = document.getElementById("moose-equation-96c6a371-f660-4acc-9158-2f6651d29a05");katex.render("\\mu", element, {displayMode:false,throwOnError:false});$ the molecular viscosity.

The same wall-function options as used for the turbulent viscosity are available via "wall_treatment" (e.g. eq_newton, eq_incremental, eq_linearized, neq).

Viscous sublayer ( $var element = document.getElementById("moose-equation-71fb82ec-709c-4723-9157-c8e1699de07b");katex.render("y^+ < 11.25", element, {displayMode:false,throwOnError:false});$ )

For cells in the viscous sublayer (sub-laminar region), turbulent production is negligible and the model sets

$var element = document.getElementById("moose-equation-b2ec5e6c-681e-4b8f-bf93-0bc0f66c151b");katex.render("G_k = 0.", element, {displayMode:true,throwOnError:false});$

Destruction is modeled using a near-wall expression of the form

$var element = document.getElementById("moose-equation-6bcf3beb-8b74-4d90-acc6-08c8e2a0cc58");katex.render("\\epsilon = \\frac{2\\mu k}{y_p^2},", element, {displayMode:true,throwOnError:false});$ which enforces rapid dissipation very close to the wall.

Logarithmic region ( $var element = document.getElementById("moose-equation-bcf88d1c-09f8-44e3-8c2c-c7e3c85d2437");katex.render("y^+ \\ge 11.25", element, {displayMode:false,throwOnError:false});$ )

For cells in the log layer, production and destruction take the classical wall-function forms:

Production: $var element = document.getElementById("moose-equation-bf98d459-eb85-4e3b-8186-7521476e56cd");katex.render("G_k = \\tau_w \\|\\nabla \\mathbf{u}\\| = \\mu_w \\|\\nabla \\mathbf{u}\\| \\frac{C_\\mu^{1/4}\\sqrt{k}}{\\kappa y_p},", element, {displayMode:true,throwOnError:false});$
where
- $var element = document.getElementById("moose-equation-d60bb052-18a6-40fd-9e13-c4eb5caef869");katex.render("\\tau_w", element, {displayMode:false,throwOnError:false});$ is the wall shear stress,
- $var element = document.getElementById("moose-equation-93feac7d-b667-4ced-a4b8-98c34a1aee69");katex.render("\\mu_w", element, {displayMode:false,throwOnError:false});$ is the effective wall viscosity (including turbulent contributions),
- $var element = document.getElementById("moose-equation-280df3de-3300-4616-bfd0-d413ea9ef6b1");katex.render("\\|\\nabla \\mathbf{u}\\|", element, {displayMode:false,throwOnError:false});$ is the wall-normal velocity gradient norm,
- $var element = document.getElementById("moose-equation-65441dc4-e409-47dc-b135-b730c75d2099");katex.render("C_\\mu", element, {displayMode:false,throwOnError:false});$ is the k– $var element = document.getElementById("moose-equation-0159dc8f-f483-4637-aaa8-d3a2439818fc");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ closure coefficient,
- $var element = document.getElementById("moose-equation-fa5fef2d-efc0-4a3e-bcb0-709a3c2b5d3c");katex.render("\\kappa \\approx 0.41", element, {displayMode:false,throwOnError:false});$ is the von Kármán constant.
Destruction: $var element = document.getElementById("moose-equation-b4138104-7944-4948-8787-519bb5ca9085");katex.render("\\epsilon = C_\\mu^{3/4} \\frac{\\rho k^{3/2}}{\\kappa y_p}.", element, {displayMode:true,throwOnError:false});$

In the logarithmic region, production and destruction are balanced in such a way that the near-wall TKE profile is consistent with the log-law of the wall.

note

When the near-wall treatment is handled by this kernel, the wall boundary condition for $var element = document.getElementById("moose-equation-5a6582c8-4610-4caa-9bd2-b148ade17cda");katex.render("k", element, {displayMode:false,throwOnError:false});$ is effectively provided by the model and explicit Dirichlet boundary conditions for $var element = document.getElementById("moose-equation-d3e26592-8632-4b6f-96e8-4f6d3f937b6d");katex.render("k", element, {displayMode:false,throwOnError:false});$ on those walls are typically unnecessary.

Interaction with $var element = document.getElementById("moose-equation-9abfb639-887f-4d25-84d6-250894a9d73d");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ -equation and viscosity

kEpsilonTKESourceSink only forms the k-equation source terms. For a complete k– $var element = document.getElementById("moose-equation-f6c37a40-98ef-40a6-b38b-ec38ca1320a8");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ model, this kernel must be used together with:

kEpsilonTKEDSourceSink, which provides the $var element = document.getElementById("moose-equation-6ceee49c-0860-4415-b42e-608a61282ff9");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ -equation sources,

including:

realizable production terms,
nonlinear contributions,
ow-Re extra production $var element = document.getElementById("moose-equation-89458e35-a6b0-4b84-aef0-532f54e6b570");katex.render("G'", element, {displayMode:false,throwOnError:false});$ ,
Yap correction,
buoyancy source for $var element = document.getElementById("moose-equation-a9fbcf99-4ffb-4790-a93c-63b1e6468fb5");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ ;

kEpsilonViscosity, which computes the turbulent viscosity $var element = document.getElementById("moose-equation-41b8212d-eb35-4c26-857f-0826ff0a6fa1");katex.render("\\mu_t", element, {displayMode:false,throwOnError:false});$ from the current $var element = document.getElementById("moose-equation-3f9d55e7-c309-493b-af29-0188b5b5ee70");katex.render("k", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-a5212721-777c-4c6f-b14e-cee3917f3cc7");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ , including:

low-Re damping for SKE-LRe,
two-layer viscosity blending,
realizable variable $var element = document.getElementById("moose-equation-7e7b627b-c0a4-49cf-906f-3e721bf28686");katex.render("C_\\mu", element, {displayMode:false,throwOnError:false});$ ,
bulk wall treatment for $var element = document.getElementById("moose-equation-d1c1c1d8-502c-47ff-9e72-047b72cb746f");katex.render("\\mu_t", element, {displayMode:false,throwOnError:false});$ (if requested).

The three objects together implement the full k– $var element = document.getElementById("moose-equation-490d9065-d610-4f5a-b020-87f76d9edf36");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ turbulence model family.

Input Parameters

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.
mu_tTurbulent 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:Turbulent viscosity. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
rhoFluid density. 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:Fluid 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 whose linear system this object contributes to
C++ Type:LinearVariableName
Unit:(no unit assumed)
Controllable:No
Description:The name of the variable whose linear system this object contributes to

Required Parameters

C_M1Coefficient used in the compressibility correction gamma_M.
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Coefficient used in the compressibility correction gamma_M.
C_mu0.09Turbulent kinetic energy closure constant C_mu.
Default:0.09
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Turbulent kinetic energy closure constant C_mu.
C_pl10Production limiter multiplier C_pl.
Default:10
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Production limiter multiplier C_pl.
Pr_t0.9Turbulent Prandtl number used in the buoyancy production term.
Default:0.9
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Turbulent Prandtl number used in the buoyancy production term.
betaThermal expansion coefficient functor used for buoyancy production (beta(T)). 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:Thermal expansion coefficient functor used for buoyancy production (beta(T)). A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
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
curvature_modelnoneCurvature / rotation correction model: 'none' (default) or 'standard' (Spalart–Shur style).
Default:none
C++ Type:MooseEnum
Options:none, standard
Controllable:No
Description:Curvature / rotation correction model: 'none' (default) or 'standard' (Spalart–Shur style).
gravity0 0 0Gravity vector used for buoyancy production.
Default:0 0 0
C++ Type:libMesh::VectorValue<double>
Unit:(no unit assumed)
Controllable:No
Description:Gravity vector used for buoyancy production.
k_epsilon_variantStandardk-epsilon model variant used for the bulk production term.
Default:Standard
C++ Type:MooseEnum
Options:Standard, StandardLowRe, StandardTwoLayer, Realizable, RealizableTwoLayer
Controllable:No
Description:k-epsilon model variant used for the bulk production term.
nonlinear_modelnoneNon-linear constitutive relation for standard k-epsilon.
Default:none
C++ Type:MooseEnum
Options:none, quadratic, cubic
Controllable:No
Description:Non-linear constitutive relation for standard k-epsilon.
speed_of_soundSpeed of sound functor used for compressibility correction gamma_M. 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:Speed of sound functor used for compressibility correction gamma_M. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
temperatureTemperature functor used for buoyancy production. 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:Temperature functor used for buoyancy production. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
use_buoyancyFalseInclude buoyancy production Gb in the bulk k-equation.
Default:False
C++ Type:bool
Controllable:No
Description:Include buoyancy production Gb in the bulk k-equation.
use_compressibilityFalseInclude compressibility correction gamma_M in the bulk k-equation.
Default:False
C++ Type:bool
Controllable:No
Description:Include compressibility correction gamma_M in the bulk k-equation.
use_curvature_correctionFalseApply a curvature correction factor f_c to the shear production.
Default:False
C++ Type:bool
Controllable:No
Description:Apply a curvature correction factor f_c to the shear production.
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_treatmentneqMethod used for wall functions.
Default:neq
C++ Type:MooseEnum
Options:eq_newton, eq_incremental, eq_linearized, neq
Controllable:No
Description:Method used for wall functions.
wallsBoundaries that correspond to solid walls.
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:Boundaries that correspond to solid walls.

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
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

ghost_layers1The number of layers of elements to ghost.
Default:1
C++ Type:unsigned short
Controllable:No
Description:The number of layers of elements to ghost.
use_point_neighborsFalseWhether to use point neighbors, which introduces additional ghosting to that used for simple face neighbors.
Default:False
C++ Type:bool
Controllable:No
Description:Whether to use point neighbors, which introduces additional ghosting to that used for simple face neighbors.

Parallel Ghosting Parameters

Input Files

(test/tests/kEpsilon/special-cases/channel_ERCOFTAC.i)
(validation/free_flow/isothermal/ercoftac_030_bfs/bfs_input.i)
(test/tests/kEpsilon/channel/channel_ERCOFTAC.i)
(test/tests/kEpsilon/special-cases/channel_ERCOFTAC_buoyant.i)
(test/tests/kEpsilon/special-cases/lid-driven-curvature-corrected.i)
(test/tests/kEpsilon/lid-driven-cavity/lid-driven.i)
(test/tests/kEpsilon/special-cases/channel_ERCOFTAC_bulk_treatment.i)
(test/tests/kEpsilon/special-cases/channel_ERCOFTAC_low_Re.i)

References

Florian R Menter. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA journal, 32(8):1598–1605, 1994.

@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"
}

(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'
  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
    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}
  []
[]

[LinearFVKernels]
  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    advected_interp_method = ${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 = ${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 = LinearFVAnisotropicDiffusion
    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
  []
  [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
    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'
  []
[]

(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'
  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
    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
  []
[]

[LinearFVKernels]
  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    advected_interp_method = ${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 = ${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 = LinearFVAnisotropicDiffusion
    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
  []
  [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
    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'
  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
    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}
  []
[]

[LinearFVKernels]
  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    advected_interp_method = ${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 = ${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 = LinearFVAnisotropicDiffusion
    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
  []
  [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
    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'
  []
[]

(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'
  advected_interp_method = ${advected_interp_method}
[]

[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
  []
[]

[LinearFVKernels]
  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    advected_interp_method = ${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 = ${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 = LinearFVAnisotropicDiffusion
    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
  []
  [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
    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/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'
  advected_interp_method = 'upwind'
[]

[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}
  []
[]

[LinearFVKernels]
  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    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
    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 = LinearFVAnisotropicDiffusion
    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
  []
  [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
    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'
  []
[]

(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'
  advected_interp_method = 'upwind'
[]

[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}
  []
[]

[LinearFVKernels]
  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    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
    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 = LinearFVAnisotropicDiffusion
    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
  []
  [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
    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/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'
  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
    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}
  []
[]

[LinearFVKernels]
  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    advected_interp_method = ${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 = ${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 = LinearFVAnisotropicDiffusion
    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
  []
  [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
    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'
  []
[]

(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'
  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
    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}
  []
[]

[LinearFVKernels]
  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    advected_interp_method = ${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 = ${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 = LinearFVAnisotropicDiffusion
    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
  []
  [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
    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'
  []
[]

Model variants
Bulk formulation
Wall formulation
Interaction with - equation and viscosity var element = document.getElementById("moose-equation-9abfb639-887f-4d25-84d6-250894a9d73d");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});
Input Parameters
Input Files
References

kEpsilonTKESourceSink

Model variants

Bulk formulation

Strain-rate tensor and invariants

Shear production var element = document.getElementById("moose-equation-26504802-ec99-42fb-bdc0-ae04228526fd");katex.render("G_k", element, {displayMode:false,throwOnError:false});

Curvature / rotation correction var element = document.getElementById("moose-equation-8ff96247-5c6a-46ae-ba43-20184f91776b");katex.render("f_c", element, {displayMode:false,throwOnError:false});

Buoyancy production var element = document.getElementById("moose-equation-4785b38c-edca-44d3-8911-8a1dafebb261");katex.render("G_b", element, {displayMode:false,throwOnError:false});

nonlinear production var element = document.getElementById("moose-equation-ea140f4a-a26b-438b-926f-877c55a6f8c1");katex.render("G_{\\text{nl}}", element, {displayMode:false,throwOnError:false});

Compressibility correction var element = document.getElementById("moose-equation-d10bdada-b632-4336-a844-5609bc75a1e4");katex.render("\\gamma_M", element, {displayMode:false,throwOnError:false});