CFD Online Logo CFD Online URL
Home > Forums

Closing on wall functions - part 2: the iterative procedure

Register Blogs Members List Search Today's Posts Mark Forums Read

Rate this Entry

Closing on wall functions - part 2: the iterative procedure

Posted April 23, 2022 at 05:39 by sbaffini
Updated May 14, 2022 at 16:30 by sbaffini

The first step in the wall function approach delineated in the first post of this series requires determining \tau_w from

\tau_w = \frac{\left[\left(U_p - U_w\right) \left(\frac{\mu}{y_p}\right) - y_p\sum_{i=0}^{N}\frac{F_U^i}{i+1}\left(\frac{{s_U^i}^+}{{y^+}^{i+2}}\right)\bigg\rvert_{y_p^+} \right]}{\left(\frac{{s_U^{-1}}^+}{y^+}\right)\bigg\rvert_{y_p^+}}

with analogous steps also required for the temperature and scalars. If one has a mean to univocally/externally determine y_p^+, then velocity, temperature and scalars are all equivalent and the relation above (properly adjusted for temperature and scalars) can be used in one shot to determine the wall flux (or the wall value if the flux is known). For k (turbulent kinetic energy) based turbulence models one can, for example, assume:

y_p^+ = \frac{\rho y_p \sqrt{\sqrt{C_{\mu}}k_p}}{\mu}

For the Spalart-Allmaras model instead, one can assume (it will be shown in a following post how this is, indeed, relevant):

y_p^+ = \frac{\rho \widetilde{\nu}_p}{\mu \kappa}

with \kappa the von Karman constant used in the model. These relations are strictly valid in equilibrium cases (F_U^i=0 for all i) but are usually used in non equilibrium cases as well (there is a second advantage, besides non requiring iterations, that will be clear in a moment). If such variables/definitions are not available (say, in LES) one tipically defines:

y_p^+ = \frac{\rho y_p u_{\tau}}{\mu}


u_{\tau} = \sqrt{\frac{\left|\tau_w\right|}{\rho}}

thus one ends up with an implcit non linear relation for \tau_w that needs iterations to be resolved. Before going to the algorithmic part of these iterations, however, it is worth mentioning that the above choice for u_{\tau} is largely inconvenient for non equilibrium flows, as the zero flux outcome is only possible when y^+=0 instead of, generally, when the numerator becomes 0 (which is the other advantage of the y^+ formulations based on turbulence variables). Following the relevant literature on the topic, we thus redefine the y^+ as follows*:

y_p^+ = \frac{\rho y_p}{\mu} \sqrt{\frac{\left|\tau_w\right|}{\rho} + \sum_{i=0}^{N}\left[\frac{\mu \left|F_U^i\right|}{\rho^2\left(i+1\right)}\right]^{\frac{1}{3}}} = \sqrt{\left|\theta_p\right|+\sum_{i=0}^{N}\left|F_p^i\right|^{\frac{2}{3}}}

where we have introduced the nondimensional parameters:

\theta_p = \frac{\rho \tau_w y_p^2}{\mu^2}

F_p^i = \frac{\rho F_U^i y_p^3}{\mu^2\left(i+1\right)}

As both appear in the y_p^+ definition, it seems natural to transform the \tau_w equation as well, by multiplying it with the factor \rho y_p^2/\mu^2. The result is:

\theta_p = \frac{Re_p - \sum_{i=0}^{N}F_p^i\left(\frac{{s_U^i}^+}{{y^+}^{i+2}}\right)\bigg\rvert_{y_p^+}}{\left(\frac{{s_U^{-1}}^+}{y^+}\right)\bigg\rvert_{y_p^+}} = f\left(Re_p,F_p^i,y_p^+\left(\theta_p,F_p^i\right)\right) = f\left(\theta_p\right)

where we have introduced the additional nondimensional number:

Re_p = \frac{\rho \left(U_p-U_w\right) y_p}{\mu}

and highlighted the general dependence of the RHS of the \theta_p equation from \theta_p itself trough y_p^+. Despite certainly having an extra \mu^2 factor for all the terms (which could then be simplified to make all the numbers involved smaller), the advantage of this form is that being non dimensional it can be also solved once and for all for all the required nondimensional values and used at runtime as an interpolation table. Indeed, at this stage, it has to be reminded that the procedure is still applicable to generic \mu_t/\mu functions, and their integration might be costly as well. However, even in this case one would need to solve the equation above at least once for each set of relevant parameters. And having a nondimensional equation helps developing an iterative procedure.

I have found that with the constant non equilibrium term only, for different \mu_t/\mu parameterizations, the above equation can be easily solved with a Newton-Raphson method with no more than 10 iterations (starting always from the same laminar solution, which is typically very wrong for high Re):

\theta_p^{n+1} = \frac{f\left(\theta_p^n\right)-\theta_p^n f'\left(\theta_p^n\right)}{1-f'\left(\theta_p^n\right)}

provided that f'\left(\theta_p^n\right) is set to 0 whenever it is equal or greater than 1. This specific modification, which reverts the method to fixed point iterations, is needed not only to avoid a possible division by 0, but also because the fixed point iterations will push the solution out of a problematic branch of the function.

I provide, again without proof, the expression for the derivative f'\left(\theta_p\right) needed in the method:

f'\left(\theta_p\right) = -\frac{SIGN\left(\theta_p\right)}{2\left(\frac{{s_U^{-1}}^+}{y^+}\right)\bigg\rvert_{y_p^+}} \left[f\left(\theta_p\right) {q^{-1}}^+ + \sum_{i=0}^{N}F_p^i {q^i}^+\right]

where the {q^i}^+ functions have been defined in the first post of this series.

In equilibrium and favourable pressure gradient cases the above equation really has no surprises. However, for adverse pressure gradients the matter is more complicated and up to 3 solutions are typically possible. The procedure above has been found to always pick up the non separated solution closest to the equilibrium one as long as possible, while the separated solution is taken only when it is the only solution to the equation. The algorithm might still get trapped in a positive solution branch when the actual solution is negative, but this seems largely ininfluential in the overall picture (considering also the nature of the solution, which is an approximate companion problem for a bc) and only affects the specific point where the switch of the solution branches actually happens. Of course, this reasoning is based on a certain class of common turbulent viscosity profiles; changing them to properly adapt to the non equilibrium case is actually the key in properly solving the problem.

In the next two posts in the series I will provide two closed form analytical solutions for the integrals {s_{U,T}^i}^+ and {p^i}^+. The first will be the equivalent of the standard wall function, but for arbitrary N. The main advantage of its derivation in the present more rigorous framework is for the thermal case, where the resulting function just works for Pr/Pr_t \le 1, while for the Pr/Pr_t > 1 case involves the Jayatilleke function only for the determination of the switching position between the linear ang logarithmic part and not its actual value (where, however, it is still implicitly present).

The second solution will be a reworking of the already presented Musker-Monkewitz solution. I will only provide a template of the solution here, but that template has practically shown to be easily extendable to arbitrary N. In particular, aside from a polynomial term, the same solution template works for arbitrary N, the difference being only in the function coefficients. I recall here that the Musker-Monkewitz solution is an all y+ solution that by a math trick can be used for arbitrary Pr/Pr_t ratios (see, for example, here), the underlying turbulent voscosity has the correct C{y^+}^3 behavior near the wall and the correct \kappa y^+ one in the logarithmic region. Also, it can be made the exact near wall behavior of the Spalart-Allmaras model by a trivial modification of any correct implementation of the model (more on this in a later post).

* Note that the y_p^+ definition used above has a very specific context, which is the evaluation of the integrals involving the turbulent viscosity ratio, because it is the variable used to parameterize it. In practice, it helps regularizing the implicit function for \tau_w. However, the true y^+, only dependent on \tau_w, remains the only relevant one for any other mean (including postprocessing).
Posted in Uncategorized
Views 592 Comments 0 Edit Tags Email Blog Entry
« Prev     Main     Next »
Total Comments 0



All times are GMT -4. The time now is 06:35.