kEpsilonTKEDSourceSink

kEpsilonTKEDSourceSink is a finite volume elemental kernel that computes the turbulent source and sink terms for the $var element = document.getElementById("moose-equation-8b104c49-c60a-470c-8a6e-dbd2120f96f4");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ -equation in the k– $var element = document.getElementById("moose-equation-86756707-58ee-402b-b5b0-baf0762fe204");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ family of models. Together with kEpsilonTKESourceSink and kEpsilonViscosity, it forms the full set of transport and closure relations for the k– $var element = document.getElementById("moose-equation-5b599491-360c-41b0-87d4-34654cf64069");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ turbulence models implemented in OpenPronghorn.

note

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

The $var element = document.getElementById("moose-equation-ec818aa6-04fc-4c2c-88ac-057d9fe9bb40");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ -equation is written in conservative form as

$var element = document.getElementById("moose-equation-f2543224-852b-41d2-ba53-731eae64018c");katex.render("\\frac{\\partial \\rho \\epsilon}{\\partial t} + \\nabla \\cdot (\\rho \\mathbf{u} \\epsilon) = P_\\epsilon - D_\\epsilon,", element, {displayMode:true,throwOnError:false});$

where

$var element = document.getElementById("moose-equation-d971183a-c5db-4f10-a776-c55925ca984b");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ is the turbulent dissipation rate,
$var element = document.getElementById("moose-equation-7a02c308-9ade-469a-978e-dc657e243727");katex.render("\\mathbf{u}", element, {displayMode:false,throwOnError:false});$ is the mean velocity,
$var element = document.getElementById("moose-equation-e1af5f64-282c-4ffc-8d47-84ed475d7f66");katex.render("P_\\epsilon", element, {displayMode:false,throwOnError:false});$ is the assembled production-side term for the $var element = document.getElementById("moose-equation-a63ad7b9-c90f-4cbf-b2d4-d7f1c3e09bd7");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ -equation (shear, buoyancy, nonlinear production, and any variant-specific corrections such as Yap or the low-Re $var element = document.getElementById("moose-equation-241f8505-eeb7-4fe0-89bd-27c0c1bdc36e");katex.render("G'", element, {displayMode:false,throwOnError:false});$ term),
$var element = document.getElementById("moose-equation-b4206104-cd06-421a-9f24-27b01012e4c6");katex.render("D_\\epsilon", element, {displayMode:false,throwOnError:false});$ is the destruction term proportional to $var element = document.getElementById("moose-equation-5da6106c-61c4-4ada-8be0-a121c56e9b3c");katex.render("\\rho\\,\\epsilon^2/k", element, {displayMode:false,throwOnError:false});$ .

Some references split the right-hand side as $var element = document.getElementById("moose-equation-95ec080c-a8e7-41d0-a9cb-2e0896437161");katex.render("\\rho\\,(P_\\epsilon - D_\\epsilon + S_\\epsilon)", element, {displayMode:false,throwOnError:false});$ . In OpenPronghorn we fold any additional sources/corrections that would fall under $var element = document.getElementById("moose-equation-29066a98-311b-40ba-b39f-10bf2dd2f92a");katex.render("S_\\epsilon", element, {displayMode:false,throwOnError:false});$ directly into $var element = document.getElementById("moose-equation-a9475978-6a2b-43c5-8f08-eecebf86760d");katex.render("P_\\epsilon", element, {displayMode:false,throwOnError:false});$ , so the kernel assembles a single production term.

In the implementation, kEpsilonTKEDSourceSink forms the $var element = document.getElementById("moose-equation-6edd79dd-9ac2-406f-b0d1-75590401239d");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ -source term as $var element = document.getElementById("moose-equation-b7a8bce7-9a6c-42e9-b888-da34d7d92f05");katex.render("\\text{source}_\\epsilon = P_\\epsilon - D_\\epsilon = P_\\epsilon - \\rho C_{\\epsilon2} f_2 \\frac{\\epsilon^2}{k},", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-7c4c5f86-d09c-4955-b40f-58c1da883b82");katex.render("C_{\\epsilon2}", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-ece761f1-ad0d-42de-aed4-be3cd6923609");katex.render("f_2", element, {displayMode:false,throwOnError:false});$ depend on the selected k– $var element = document.getElementById("moose-equation-0ca0a8d0-80b6-400a-ba85-c4ed8325bee3");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ variant and damping model.

Model variants

The kernel supports the same k– $var element = document.getElementById("moose-equation-f46e16dd-ac3b-4b54-832c-79f99443299b");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ variants as kEpsilonTKESourceSink via the "k_epsilon_variant" parameter:

Standard — classical high-Reynolds-number k– $var element = document.getElementById("moose-equation-a9b8c2ae-b6a9-4a69-8c3e-8bc6a87a5354");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ model.
StandardLowRe — low-Reynolds-number k– $var element = document.getElementById("moose-equation-3cbf007c-7dda-457a-aa90-e0c4a46145ad");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ model with damping functions.
StandardTwoLayer — two-layer k– $var element = document.getElementById("moose-equation-8d87218a-8210-4678-979a-0439cb22ee85");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ formulation using a blending between outer k– $var element = document.getElementById("moose-equation-9a358cd5-e281-45a3-a11f-0e5e43b47087");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ and near-wall length scales.
Realizable — realizable k– $var element = document.getElementById("moose-equation-1728fe95-479f-49be-a4d2-d6373c705898");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ model with variable $var element = document.getElementById("moose-equation-ba06121c-1d5d-40ca-9c9b-e5e1a9efc9b0");katex.render("C_\\mu", element, {displayMode:false,throwOnError:false});$ .
RealizableTwoLayer — realizable model with a two-layer near-wall treatment.

The $var element = document.getElementById("moose-equation-0de71c15-b46a-40be-b9ee-d1c51428be9c");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ -production term $var element = document.getElementById("moose-equation-e9aa06ad-bc0a-4c9b-85d5-8854fd2d51c5");katex.render("P_\\epsilon", element, {displayMode:false,throwOnError:false});$ is assembled differently for each variant, combining shear production, buoyancy contributions, nonlinear production, Yap correction, and low-Re extra production $var element = document.getElementById("moose-equation-2a36af18-5022-4a35-9474-40baa7b3484e");katex.render("G'", element, {displayMode:false,throwOnError:false});$ .

Strain-rate invariants and basic quantities

As in the k-equation kernel, kEpsilonTKEDSourceSink uses the invariants of the strain and rotation tensors:

$var element = document.getElementById("moose-equation-bde557e0-eb4c-409e-afec-c599223129b5");katex.render("S^2 = 2 S_{ij} S_{ij}", element, {displayMode:false,throwOnError:false});$ ,
$var element = document.getElementById("moose-equation-51517cef-1066-4349-a58d-639bf0782f26");katex.render("W^2 = 2 W_{ij} W_{ij}", element, {displayMode:false,throwOnError:false});$ ,
$var element = document.getElementById("moose-equation-1eaaa231-46df-44bb-9ee4-9e209fa03795");katex.render("\\nabla \\cdot \\mathbf{u}", element, {displayMode:false,throwOnError:false});$ ,

See the theory section for the definition of these ones. These invariants are used to construct:

the shear production of k, $var element = document.getElementById("moose-equation-722cf6b1-97a8-4375-b812-e64fa5adfcb7");katex.render("G_k = \\mu_t S^2 \\;\\; (\\text{plus compressible terms if enabled}),", element, {displayMode:true,throwOnError:false});$
the shear-based $var element = document.getElementById("moose-equation-ec4e0c4c-4eb7-48f1-8e92-3c4837452844");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ -production term $var element = document.getElementById("moose-equation-6befa72d-9e2e-4a32-8e66-d71cfdedbb7b");katex.render("S_k = \\mu_t S^2", element, {displayMode:false,throwOnError:false});$ for realizable variants (note that $var element = document.getElementById("moose-equation-997d48a3-e8d3-4e95-9f6d-962a3ff81263");katex.render("G_k", element, {displayMode:false,throwOnError:false});$ may be limited so a disctinction is made with the non-limited $var element = document.getElementById("moose-equation-23ae6447-2d1d-4ca2-b21d-3a2e60332056");katex.render("S_k", element, {displayMode:false,throwOnError:false});$ ),
the curvature correction factor $var element = document.getElementById("moose-equation-91b71c34-1cc5-4d52-a0fd-ccddc6fbfd02");katex.render("f_c", element, {displayMode:false,throwOnError:false});$ (for realizable models when enabled),
the normalized strain/rotation quantities used by the nonlinear constitutive relations.

$var element = document.getElementById("moose-equation-a2058b5a-84b9-48d6-8a72-726e8cf3e865");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ -production term $var element = document.getElementById("moose-equation-bc528d09-1729-461f-816b-1f096c86c916");katex.render("P_\\epsilon", element, {displayMode:false,throwOnError:false});$

The $var element = document.getElementById("moose-equation-2fa3353c-8ff1-4a01-8d5c-1f0ef2efe0b1");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ -production term $var element = document.getElementById("moose-equation-e3ef7212-0fe8-4d18-8f0b-c48399863374");katex.render("P_\\epsilon", element, {displayMode:false,throwOnError:false});$ is built from several contributions. Here, $var element = document.getElementById("moose-equation-4dd24161-82c1-4371-8c45-82c48ce6f20d");katex.render("P_\\epsilon", element, {displayMode:false,throwOnError:false});$ denotes the full non-destruction part of the $var element = document.getElementById("moose-equation-5b9eaac8-a24f-4b8f-a130-639bdcd360f4");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ source term (i.e., it already includes terms that are sometimes written separately as $var element = document.getElementById("moose-equation-234eba05-cf42-470f-bc55-ede5342f2a3a");katex.render("S_\\epsilon", element, {displayMode:false,throwOnError:false});$ , such as Yap or other near-wall/low-Re corrections):

$var element = document.getElementById("moose-equation-000baa76-039d-48a0-b24a-c18cdb2357b3");katex.render("G_k", element, {displayMode:false,throwOnError:false});$ : shear production of k,
$var element = document.getElementById("moose-equation-30b6f6a1-caf9-48d7-acb5-590c627114f7");katex.render("G_{\\text{nl}}", element, {displayMode:false,throwOnError:false});$ : nonlinear production (optional quadratic/cubic constitutive relations),
$var element = document.getElementById("moose-equation-21593821-23ad-4a72-b301-7bd07a2532f5");katex.render("G'", element, {displayMode:false,throwOnError:false});$ : low-Re extra production (StandardLowRe only),
$var element = document.getElementById("moose-equation-4b977cf1-0f5d-47bd-90e0-a9163995baa2");katex.render("G_b", element, {displayMode:false,throwOnError:false});$ : buoyancy production,
$var element = document.getElementById("moose-equation-937c8708-0161-45e6-a2fb-028e3890b476");katex.render("Y_y", element, {displayMode:false,throwOnError:false});$ : Yap correction (two-layer and low-Re variants),
$var element = document.getElementById("moose-equation-f5b5cb1e-8248-4d71-9d4d-b353ebfc36df");katex.render("f_c", element, {displayMode:false,throwOnError:false});$ : curvature correction factor (realizable variants).

The exact combination depends on the chosen k– $var element = document.getElementById("moose-equation-7d680d0e-3f16-41bf-ac9b-cfda2dd6eb87");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ variant.

Standard k– $var element = document.getElementById("moose-equation-160ff5af-9823-48de-bc73-5e45baa208e2");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$

For the Standard high-Re k– $var element = document.getElementById("moose-equation-74931ea1-7117-40aa-87ae-f707b624b859");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ model, the $var element = document.getElementById("moose-equation-d76c3762-fb33-4550-9573-5da798558c1e");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ -production is

$var element = document.getElementById("moose-equation-fe1fe3a9-7a5f-4a1f-aeb4-9b12bf1d2879");katex.render("G_k^{\\text{lim}} + G_{\\text{nl}} + C_{\\epsilon3} G_b,", element, {displayMode:true,throwOnError:false});$

where

$var element = document.getElementById("moose-equation-779975b9-12d3-4016-86b4-1bbacecc878b");katex.render("G_k^{\\text{lim}}", element, {displayMode:false,throwOnError:false});$ is the limited shear production $var element = document.getElementById("moose-equation-4d7c5980-261d-4abe-ab83-fcc0c0ac59e7");katex.render("G_k^{\\text{lim}} = \\min(G_k, C_{PL}\\,\\rho\\,\\epsilon),", element, {displayMode:true,throwOnError:false});$
$var element = document.getElementById("moose-equation-b16a79a2-90e5-49f6-809d-cd3fcf6603bc");katex.render("G_{\\text{nl}}", element, {displayMode:false,throwOnError:false});$ is the nonlinear production term (if nonlinear_model != none),
$var element = document.getElementById("moose-equation-e84df53f-9d3e-4a98-87c9-bb0980633105");katex.render("C_{\\epsilon3}", element, {displayMode:false,throwOnError:false});$ is the buoyancy coefficient (usually model-dependent).

The production limiter is consistent with the one used in the k-equation kernel.

StandardTwoLayer k– $var element = document.getElementById("moose-equation-fc817866-eddb-4fc5-afd9-166825e5fb3e");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$

For the StandardTwoLayer model, the $var element = document.getElementById("moose-equation-94667820-4728-439a-ba82-e11efc620a5c");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ -production is

$var element = document.getElementById("moose-equation-82c28946-5873-4ab5-85c4-94a6256cbb3b");katex.render("G_k + G_{\\text{nl}} + C_{\\epsilon3} G_b + \\rho C_{\\epsilon1} Y_y,", element, {displayMode:true,throwOnError:false});$

where

the shear production $var element = document.getElementById("moose-equation-0a412f1d-7cd8-4bc7-b84d-317bcacb097d");katex.render("G_k", element, {displayMode:false,throwOnError:false});$ is taken without limiter,
the Yap correction $var element = document.getElementById("moose-equation-eb1e1fa7-0060-4b68-bfd7-81dfc4ba30b6");katex.render("Y_y", element, {displayMode:false,throwOnError:false});$ provides an additional near-wall sink/source and is multiplied by $var element = document.getElementById("moose-equation-abe39b81-4e60-4d95-bfa4-e55cdd1c5b33");katex.render("\\rho C_{\\epsilon1}", element, {displayMode:false,throwOnError:false});$ .

The Yap term is especially important near walls and is detailed in a separate section below.

StandardLowRe k– $var element = document.getElementById("moose-equation-15c20fac-4dec-421f-8c0f-40087f751773");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ (low-Re model)

For the StandardLowRe model, low-Re corrections are included via damping functions and the extra production term $var element = document.getElementById("moose-equation-2a5ebd23-7268-4210-baaa-d7784a637f75");katex.render("G'", element, {displayMode:false,throwOnError:false});$ . The $var element = document.getElementById("moose-equation-d4ab259c-7ac8-4b00-8f21-57631738e2fe");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ -production is

$var element = document.getElementById("moose-equation-713a79ca-1f3a-4cff-ad98-b189f28884d3");katex.render("G_k + G_{\\text{nl}} + G' + C_{\\epsilon3} G_b + \\rho C_{\\epsilon1} Y_y,", element, {displayMode:true,throwOnError:false});$

where

$var element = document.getElementById("moose-equation-33a733fd-d85b-4e71-87c3-fc7a10d79ab6");katex.render("G'", element, {displayMode:false,throwOnError:false});$ is the low-Re extra production (see below),
$var element = document.getElementById("moose-equation-82ac194c-c316-4fe0-80cb-3a07c5b5fbb0");katex.render("f_2", element, {displayMode:false,throwOnError:false});$ damping appears in the destruction term,
Yap correction $var element = document.getElementById("moose-equation-31156ee8-1841-42d1-b9b9-d4aa79294e0e");katex.render("Y_y", element, {displayMode:false,throwOnError:false});$ is applied similarly to the two-layer model.

Realizable k– $var element = document.getElementById("moose-equation-8e5d9205-ad66-4f27-9d2c-c463dda3a14c");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$

For the Realizable high-Re model, the $var element = document.getElementById("moose-equation-e8c94e68-216e-48af-bee2-02defdc11ab7");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ -production is written in terms of the shear-based

production $var element = document.getElementById("moose-equation-b0e29a1d-8e3d-4b71-a606-2670dbd56a2e");katex.render("S_k", element, {displayMode:false,throwOnError:false});$ : $var element = document.getElementById("moose-equation-cfeae242-bb44-4096-b417-66bdf30a6b65");katex.render("f_c S_k + C_{\\epsilon3} G_b,", element, {displayMode:true,throwOnError:false});$

where

$var element = document.getElementById("moose-equation-d7fc810e-1745-48f4-9f88-59ec90d4dd7f");katex.render("S_k = \\mu_t S^2", element, {displayMode:false,throwOnError:false});$ is the shear contribution,
$var element = document.getElementById("moose-equation-6d35e585-4dc2-4b24-ba28-5572cf80c6be");katex.render("f_c", element, {displayMode:false,throwOnError:false});$ is the curvature correction factor (if curvature_model != none),
$var element = document.getElementById("moose-equation-56290535-f749-455a-bea8-cc098307d9c5");katex.render("C_{\\epsilon3} G_b", element, {displayMode:false,throwOnError:false});$ is the buoyancy contribution.

In realizable models the nonlinear constitutive relation affects the viscosity (via realizable $var element = document.getElementById("moose-equation-bb39a18c-9715-4838-9df7-7bc8f094b149");katex.render("C_\\mu", element, {displayMode:false,throwOnError:false});$ ) and the k-equation, but $var element = document.getElementById("moose-equation-c499c86c-b7c8-41db-83f1-c40afd1919f8");katex.render("P_\\epsilon", element, {displayMode:false,throwOnError:false});$ itself is expressed in terms of $var element = document.getElementById("moose-equation-c51d44de-e332-45c8-88b3-fc1f59dc846a");katex.render("S_k", element, {displayMode:false,throwOnError:false});$ rather than $var element = document.getElementById("moose-equation-9f1dce97-e386-49db-a991-13acf6cb78b4");katex.render("G_k", element, {displayMode:false,throwOnError:false});$ or $var element = document.getElementById("moose-equation-28d5266f-5e91-4f85-8454-05bbc23b04ec");katex.render("G_{\\text{nl}}", element, {displayMode:false,throwOnError:false});$ .

RealizableTwoLayer k– $var element = document.getElementById("moose-equation-4aab7f7d-b9d5-4a97-9162-e06126736133");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$

The RealizableTwoLayer model combines the realizable bulk behavior with a two-layer near-wall enhancement via Yap:

$var element = document.getElementById("moose-equation-c64f81f7-4dd0-4f9d-9fb9-6e0a8bb744f3");katex.render("f_c S_k + C_{\\epsilon3} G_b + \\rho C_{\\epsilon1} Y_y,", element, {displayMode:true,throwOnError:false});$

with the same definitions of $var element = document.getElementById("moose-equation-835ab82f-5aaf-4e25-8dc1-fd0120fb9c5b");katex.render("S_k", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-cdc80782-8e72-492f-a724-6bafa2c0fc01");katex.render("f_c", element, {displayMode:false,throwOnError:false});$ , and $var element = document.getElementById("moose-equation-060e0315-e25d-4736-b8c6-88ee6372ab66");katex.render("G_b", element, {displayMode:false,throwOnError:false});$ as in the realizable model and the Yap term as described below.

nonlinear $var element = document.getElementById("moose-equation-95dbef60-4a37-42e3-a9b3-eb3aad8885ec");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ -production $var element = document.getElementById("moose-equation-2c3e98d1-f478-40a0-b145-22d581f122e2");katex.render("G_{\\text{nl}}", element, {displayMode:false,throwOnError:false});$

When "nonlinear_model" is set to quadratic or cubic, kEpsilonTKEDSourceSink includes the contribution $var element = document.getElementById("moose-equation-7a417d13-809c-4731-b385-9e15b36a77c4");katex.render("G_{\\text{nl}}", element, {displayMode:false,throwOnError:false});$ in the $var element = document.getElementById("moose-equation-4aadec2f-e0e2-45cd-b5e8-8143e1deb021");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ -production for Standard-family variants:

Standard: $var element = document.getElementById("moose-equation-6ca021f5-c764-4f9f-aa05-c92df8110ef9");katex.render("\\text{Production}= G_k^{\\text{lim}} + G_{\\text{nl}} + C_{\\epsilon3} G_b,", element, {displayMode:true,throwOnError:false});$
StandardTwoLayer: $var element = document.getElementById("moose-equation-8658d54d-736b-46ad-957d-14a29372b10e");katex.render("\\text{Production} = G_k + G_{\\text{nl}} + C_{\\epsilon3} G_b + \\rho C_{\\epsilon1} Y_y,", element, {displayMode:true,throwOnError:false});$
StandardLowRe: $var element = document.getElementById("moose-equation-9062864e-8d4b-45e2-a67a-97f51af229e0");katex.render("\\text{Production} = G_k + G_{\\text{nl}} + G' + C_{\\epsilon3} G_b + \\rho C_{\\epsilon1} Y_y.", element, {displayMode:true,throwOnError:false});$

The term $var element = document.getElementById("moose-equation-8d42c039-f8d1-439e-9e7d-b5d8113c5d57");katex.render("G_{\\text{nl}}", element, {displayMode:false,throwOnError:false});$ is based on the quadratic or cubic nonlinear Reynolds stress models and contracted with $var element = document.getElementById("moose-equation-8f7de42a-eaaa-4238-ab64-7c15949ddf38");katex.render("\\nabla \\mathbf{u}", element, {displayMode:false,throwOnError:false});$ .

Curvature correction in $var element = document.getElementById("moose-equation-9e2c6fab-d566-4993-a408-9811ed81d004");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ -production

For realizable variants, the curvature correction factor $var element = document.getElementById("moose-equation-3bf7a6df-42e1-4cc0-b438-bdc12775d59c");katex.render("f_c", element, {displayMode:false,throwOnError:false});$ is applied multiplicatively to $var element = document.getElementById("moose-equation-19ea6c2c-e71d-4a72-9adf-185342965703");katex.render("S_k", element, {displayMode:false,throwOnError:false});$ in the $var element = document.getElementById("moose-equation-8cbf013b-3c69-43e5-91d2-dfd314882182");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ -production:

Realizable: $var element = document.getElementById("moose-equation-e09d8b45-183f-4679-b731-6337a74b2aef");katex.render("\\text{Production} = f_c S_k + C_{\\epsilon3} G_b,", element, {displayMode:true,throwOnError:false});$
RealizableTwoLayer: $var element = document.getElementById("moose-equation-91deca2e-79a7-46a9-9fb7-d000f67cd986");katex.render("\\text{Production} = f_c S_k + C_{\\epsilon3} G_b + \\rho C_{\\epsilon1} Y_y.", element, {displayMode:true,throwOnError:false});$

The same curvature model (e.g. curvature_model = standard) is used in both the k-equation and the $var element = document.getElementById("moose-equation-241df5c1-5bdd-4352-b7af-e7f7713c6f7f");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ -equation kernels for consistency.

Yap correction

The Yap correction is an additional near-wall term that modifies the $var element = document.getElementById("moose-equation-91310a57-af46-4845-ae40-6846433c7b06");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ -equation source in two-layer and low-Re variants. It is enabled when "use_yap" is true and a wall-distance functor is available.

Yap length scales

The Yap correction is built from two length scales:

a turbulence length scale $var element = document.getElementById("moose-equation-86ea733b-fd36-464f-b195-9d27c852d655");katex.render("l", element, {displayMode:false,throwOnError:false});$ ,
an $var element = document.getElementById("moose-equation-da2a6a3d-ca92-4353-a72e-8d7fad501152");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ length scale $var element = document.getElementById("moose-equation-966cdd5c-4f7a-4825-a688-0f3ff9d6e838");katex.render("l_\\epsilon", element, {displayMode:false,throwOnError:false});$ .

The implementation uses a limited turbulence time scale to obtain $var element = document.getElementById("moose-equation-a397f5ec-e671-4cd3-8d17-55fdce57ab6f");katex.render("l", element, {displayMode:false,throwOnError:false});$ . The effective time scale $var element = document.getElementById("moose-equation-70ca7af8-ef5d-495d-9d7d-4cb91647795e");katex.render("T_{\\text{eff}}", element, {displayMode:false,throwOnError:false});$ is computed as

$var element = document.getElementById("moose-equation-b25e9b2b-1318-4b9a-8842-c08e2c704145");katex.render("T_1 = \\frac{k}{\\epsilon},", element, {displayMode:true,throwOnError:false});$

$var element = document.getElementById("moose-equation-43a25d27-aca5-4bcc-9506-4f50e6647069");katex.render("T_2 = \\frac{6\\nu}{\\epsilon},", element, {displayMode:true,throwOnError:false});$

$var element = document.getElementById("moose-equation-87dfe8b7-1f25-4ae5-86d0-76c67d3e1adc");katex.render("T_3 = T_1^{1.625} T_2,", element, {displayMode:true,throwOnError:false});$

$var element = document.getElementById("moose-equation-47066f0a-ebd3-4740-9414-c7922c50246d");katex.render("T_{\\text{eff}} = T_3^{1/(1.625 + 1)},", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-5f12305f-9797-4ab3-b9b9-408be710b6e8");katex.render("\\nu = \\mu / \\rho", element, {displayMode:false,throwOnError:false});$ is the kinematic viscosity. The Yap turbulence length scale is then

$var element = document.getElementById("moose-equation-8443ef98-caf8-44a6-88f6-255096da9942");katex.render("l = \\sqrt{T_{\\text{eff}} \\epsilon}.", element, {displayMode:true,throwOnError:false});$

The $var element = document.getElementById("moose-equation-678f72ba-b9ae-475a-ac2f-fa72e50b6ea2");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ length scale is defined as

$var element = document.getElementById("moose-equation-fd706c97-34c7-4521-b85a-150a232aa730");katex.render("l_\\epsilon = C_l d,", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-3d742d66-1749-48ed-98b8-2843e14f9a59");katex.render("d", element, {displayMode:false,throwOnError:false});$ is the wall distance and

$var element = document.getElementById("moose-equation-7b5964f1-7540-4f99-a53d-9fd2ed3bcd49");katex.render("C_l = 0.42 C_\\mu^{-3/4}", element, {displayMode:true,throwOnError:false});$

with $var element = document.getElementById("moose-equation-66fa2ed2-1e52-4321-b8b1-bbdaab08c48e");katex.render("C_\\mu", element, {displayMode:false,throwOnError:false});$ the model’s closure constant (or a realizable value when applicable).

Yap correction term $var element = document.getElementById("moose-equation-e6567155-7d90-4ffc-af32-0072ef601f62");katex.render("Y_y", element, {displayMode:false,throwOnError:false});$

The Yap correction term $var element = document.getElementById("moose-equation-f9fbaf7a-3534-486f-b7ed-15e48853c5c5");katex.render("\\gamma_y", element, {displayMode:false,throwOnError:false});$ is computed as

$var element = document.getElementById("moose-equation-b2c3449f-b935-4472-a307-af9b5295c1b1");katex.render("\\gamma_y = C_w \\frac{\\epsilon^2}{k} \\max\\left[\\left(\\frac{l}{l_\\epsilon} - 1\\right)\\left(\\frac{l}{l_\\epsilon}\\right)^2, 0\\right],", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-48753035-55c7-4446-b129-a7e8d0c6fbb8");katex.render("C_w", element, {displayMode:false,throwOnError:false});$ is a model coefficient. In the $var element = document.getElementById("moose-equation-88bb0ec8-1790-40f9-9dfb-c2ff40f7c462");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ -equation, this is translated into a source term

$var element = document.getElementById("moose-equation-bebeae3a-619d-4d51-a3a4-41cb1ff203d8");katex.render("\\rho C_{\\epsilon1} Y_y = \\rho C_{\\epsilon1} \\gamma_y,", element, {displayMode:true,throwOnError:false});$

and added to $var element = document.getElementById("moose-equation-1682b37a-8ae6-41f9-a959-6ae17a5f5a31");katex.render("P_\\epsilon", element, {displayMode:false,throwOnError:false});$ for the variants that include Yap.

The Yap contribution appears in:

StandardTwoLayer,
StandardLowRe,
RealizableTwoLayer.

Low-Re extra production $var element = document.getElementById("moose-equation-61b3b7f5-88fb-454e-b5ad-937c6a2cfa9d");katex.render("G'", element, {displayMode:false,throwOnError:false});$

For the StandardLowRe model, an additional low-Re production term $var element = document.getElementById("moose-equation-359c089b-2fdf-4b0f-8bfc-8d6014df88bc");katex.render("G'", element, {displayMode:false,throwOnError:false});$ is included in $var element = document.getElementById("moose-equation-e6ee206b-5d8a-4d3e-819c-4f2782661c17");katex.render("P_\\epsilon", element, {displayMode:false,throwOnError:false});$ when "use_low_re_Gprime" is true and a wall-distance functor is provided.

The term $var element = document.getElementById("moose-equation-28e164b2-c1fc-4e0c-be56-e8da93abfc04");katex.render("G'", element, {displayMode:false,throwOnError:false});$ is defined as

$var element = document.getElementById("moose-equation-f32d7d55-ed4e-4ebf-9318-7827f75e6e5a");katex.render("G' = D f_2 \\left(G_k + 2 \\mu_t \\frac{k}{d^2}\\right) \\exp(-E \\mathrm{Re}_d^2),", element, {displayMode:true,throwOnError:false});$

where

$var element = document.getElementById("moose-equation-a07058be-251f-475b-a68d-670f906d8c1c");katex.render("D", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-f99a810a-5938-41cb-869d-d4a30e3d84ab");katex.render("E", element, {displayMode:false,throwOnError:false});$ are model coefficients,
$var element = document.getElementById("moose-equation-7eb6d00c-3c62-4fc7-8d38-e048e15b27da");katex.render("f_2", element, {displayMode:false,throwOnError:false});$ is the low-Re damping function, $var element = document.getElementById("moose-equation-941e4633-a117-4ccd-9ff2-2d1a08b9f606");katex.render("f_2 = 1 - C \\exp(-\\mathrm{Re}_t^2),", element, {displayMode:true,throwOnError:false});$ with $var element = document.getElementById("moose-equation-8063dfda-5613-4970-a648-171dba562734");katex.render("C", element, {displayMode:false,throwOnError:false});$ a low-Re coefficient,
$var element = document.getElementById("moose-equation-1f7d5d63-cdaf-4604-8460-7dcb3711f3c3");katex.render("\\mathrm{Re}_d = \\sqrt{k} d / \\nu", element, {displayMode:false,throwOnError:false});$ is the wall-distance Reynolds number,
$var element = document.getElementById("moose-equation-79941c5b-43bc-4ff5-ac49-2c309f082d87");katex.render("\\mathrm{Re}_t = k^2 / (\\nu \\epsilon)", element, {displayMode:false,throwOnError:false});$ is the turbulence Reynolds number.

In the $var element = document.getElementById("moose-equation-5cdb2735-ea00-42ad-a26a-71f11e3c3a58");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ -production for the StandardLowRe variant the production term is

$var element = document.getElementById("moose-equation-9a5570bd-5507-47af-9fed-43d71b44dc2c");katex.render("G_k + G_{\\text{nl}} + G' + C_{\\epsilon3} G_b + \\rho C_{\\epsilon1} Y_y.", element, {displayMode:true,throwOnError:false});$

The term $var element = document.getElementById("moose-equation-b88aa32c-a8bc-46d7-9ffc-1dff2f523315");katex.render("G'", element, {displayMode:false,throwOnError:false});$ vanishes as $var element = document.getElementById("moose-equation-94aed024-13fd-4816-9057-ac7f083bb730");katex.render("\\mathrm{Re}_d \\to \\infty", element, {displayMode:false,throwOnError:false});$ due to the exponential damping and is most active in the near-wall low-Re region.

Destruction term and damping

The destruction term in the $var element = document.getElementById("moose-equation-c0dbc7b1-3bbb-4bcd-9edc-0cf9a29349d5");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ -equation is modeled as

$var element = document.getElementById("moose-equation-75d27551-8311-4355-ba1b-9c3f947a6ccc");katex.render("D_\\epsilon = \\rho C_{\\epsilon2} f_2 \\frac{\\epsilon^2}{k},", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-90f02aff-efa4-4ef8-af49-cbc86f43fe79");katex.render("C_{\\epsilon2}", element, {displayMode:false,throwOnError:false});$ is the model coefficient and $var element = document.getElementById("moose-equation-27d528a5-88d6-4c76-9cc4-e8aca25ebf5d");katex.render("f_2", element, {displayMode:false,throwOnError:false});$ is a damping function that depends on the k– $var element = document.getElementById("moose-equation-9985d8b3-417e-4e4c-87bd-dd6db0f9ea22");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ variant:

Standard, StandardTwoLayer, Realizable, RealizableTwoLayer: $var element = document.getElementById("moose-equation-01e1c35b-7be7-4dfe-a139-e5795db53650");katex.render("f_2 = 1,", element, {displayMode:true,throwOnError:false});$
StandardLowRe: $var element = document.getElementById("moose-equation-17750a22-23a8-442e-b9f9-1eb8af569a8e");katex.render("f_2 = 1 - C \\exp(-\\mathrm{Re}_t^2),", element, {displayMode:true,throwOnError:false});$ consistent with the low-Re formulation used for $var element = document.getElementById("moose-equation-78fef753-c276-4b7f-a1be-52d61df1fad0");katex.render("G'", element, {displayMode:false,throwOnError:false});$ .

The same function is used for both the extra production term and the destruction damping in the low-Re model.

Input Parameters

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.
tkeCoupled turbulent kinetic energy (k). 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 (k). 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

C1_eps1.44C1 epsilon production coefficient.
Default:1.44
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:C1 epsilon production coefficient.
C2_eps1.92C2 epsilon destruction coefficient.
Default:1.92
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:C2 epsilon destruction coefficient.
C3_eps1C3 epsilon buoyancy coefficient (may be modified by Gb sign).
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:C3 epsilon buoyancy coefficient (may be modified by Gb sign).
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_lowRe1Low-Re f2 constant C used in the standard low-Re k-epsilon model.
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Low-Re f2 constant C used in the standard low-Re k-epsilon model.
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 constant multiplier.
Default:10
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Production limiter constant multiplier.
Ct6Ct coefficient used in Yap / ambient epsilon source terms.
Default:6
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Ct coefficient used in Yap / ambient epsilon source terms.
Cw0.83Yap correction constant Cw for epsilon equation.
Default:0.83
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Yap correction constant Cw for epsilon equation.
D_lowRe1Low-Re extra production coefficient D for G'.
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Low-Re extra production coefficient D for G'.
E_lowRe0.00375Low-Re extra production coefficient E for G'.
Default:0.00375
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Low-Re extra production coefficient E for G'.
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 epsilon.
Default:Standard
C++ Type:MooseEnum
Options:Standard, StandardLowRe, StandardTwoLayer, Realizable, RealizableTwoLayer
Controllable:No
Description:k-epsilon model variant used for epsilon.
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 epsilon equation.
Default:False
C++ Type:bool
Controllable:No
Description:Include buoyancy production Gb in the epsilon equation.
use_compressibilityFalseInclude compressibility correction gamma_M in the epsilon formulation (mainly in the k-equation, kept for parity).
Default:False
C++ Type:bool
Controllable:No
Description:Include compressibility correction gamma_M in the epsilon formulation (mainly in the k-equation, kept for parity).
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.
use_low_re_GprimeFalseInclude additional low-Re production G' in epsilon.
Default:False
C++ Type:bool
Controllable:No
Description:Include additional low-Re production G' in epsilon.
use_yapFalseInclude Yap correction gamma_y in epsilon.
Default:False
C++ Type:bool
Controllable:No
Description:Include Yap correction gamma_y in epsilon.
vThe velocity in the y direction. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
C++ Type:MooseFunctorName
Unit:(no unit assumed)
Controllable:No
Description:The velocity in the y direction. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
wThe velocity in the z direction. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
C++ Type:MooseFunctorName
Unit:(no unit assumed)
Controllable:No
Description:The velocity in the z direction. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
wall_distanceDistance to the closest wall; used for Yap and low-Re corrections. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
C++ Type:MooseFunctorName
Unit:(no unit assumed)
Controllable:No
Description:Distance to the closest wall; used for Yap and low-Re corrections. 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)

(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
Strain - rate invariants and basic quantities
- production term var element = document.getElementById("moose-equation-a2058b5a-84b9-48d6-8a72-726e8cf3e865");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false}); var element = document.getElementById("moose-equation-bc528d09-1729-461f-816b-1f096c86c916");katex.render("P_\\epsilon", element, {displayMode:false,throwOnError:false});
nonlinear - production var element = document.getElementById("moose-equation-95dbef60-4a37-42e3-a9b3-eb3aad8885ec");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false}); var element = document.getElementById("moose-equation-2c3e98d1-f478-40a0-b145-22d581f122e2");katex.render("G_{\\text{nl}}", element, {displayMode:false,throwOnError:false});
Curvature correction in - production var element = document.getElementById("moose-equation-9e2c6fab-d566-4993-a408-9811ed81d004");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});
Yap correction
Low - Re extra production var element = document.getElementById("moose-equation-61b3b7f5-88fb-454e-b5ad-937c6a2cfa9d");katex.render("G'", element, {displayMode:false,throwOnError:false});
Destruction term and damping
Input Parameters
Input Files

kEpsilonTKEDSourceSink

Model variants

Strain-rate invariants and basic quantities

Standard k–var element = document.getElementById("moose-equation-160ff5af-9823-48de-bc73-5e45baa208e2");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});

StandardTwoLayer k–var element = document.getElementById("moose-equation-fc817866-eddb-4fc5-afd9-166825e5fb3e");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});

StandardLowRe k–var element = document.getElementById("moose-equation-15c20fac-4dec-421f-8c0f-40087f751773");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false}); (low-Re model)

Realizable k–var element = document.getElementById("moose-equation-8e5d9205-ad66-4f27-9d2c-c463dda3a14c");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});

RealizableTwoLayer k–var element = document.getElementById("moose-equation-4aab7f7d-b9d5-4a97-9162-e06126736133");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});

Curvature correction in var element = document.getElementById("moose-equation-9e2c6fab-d566-4993-a408-9811ed81d004");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});-production

Yap correction

Yap length scales

Yap correction term var element = document.getElementById("moose-equation-e6567155-7d90-4ffc-af32-0072ef601f62");katex.render("Y_y", element, {displayMode:false,throwOnError:false});

Low-Re extra production var element = document.getElementById("moose-equation-61b3b7f5-88fb-454e-b5ad-937c6a2cfa9d");katex.render("G'", element, {displayMode:false,throwOnError:false});

Destruction term and damping

Input Parameters

Required Parameters

Optional Parameters

Contribution To Tagged Field Data Parameters

Advanced Parameters

Parallel Ghosting Parameters

Input Files

k_epsilon_variant

nonlinear_model

use_yap

use_low_re_Gprime

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

Standard k– $var element = document.getElementById("moose-equation-160ff5af-9823-48de-bc73-5e45baa208e2");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$

StandardTwoLayer k– $var element = document.getElementById("moose-equation-fc817866-eddb-4fc5-afd9-166825e5fb3e");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$

StandardLowRe k– $var element = document.getElementById("moose-equation-15c20fac-4dec-421f-8c0f-40087f751773");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ (low-Re model)

Realizable k– $var element = document.getElementById("moose-equation-8e5d9205-ad66-4f27-9d2c-c463dda3a14c");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$

RealizableTwoLayer k– $var element = document.getElementById("moose-equation-4aab7f7d-b9d5-4a97-9162-e06126736133");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$

Curvature correction in $var element = document.getElementById("moose-equation-9e2c6fab-d566-4993-a408-9811ed81d004");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ -production

Yap correction term $var element = document.getElementById("moose-equation-e6567155-7d90-4ffc-af32-0072ef601f62");katex.render("Y_y", element, {displayMode:false,throwOnError:false});$

Low-Re extra production $var element = document.getElementById("moose-equation-61b3b7f5-88fb-454e-b5ad-937c6a2cfa9d");katex.render("G'", element, {displayMode:false,throwOnError:false});$