Wall function

The wall function models the velocity profile between two points of computational domain. The coupling with LES is then performed through velocity adjustments near wall boundaries.

Log Law

Following the log wall of the wall, the average velocity profile for a determined terrain category, give its roughness can be written as (Davidson[1]):

(1)\[u\left(z\right)=\frac{u^{*}}{\kappa} \mathrm{ln}\left(\frac{z}{z_{0}}\right)\]

where \(\kappa=0.41\) is the Von Kármán constant and \(z_{0}\) the roughness length. In a numerical simulation, the friction velocity \(u^{*}\) can be calculated from the tangential velocity at the reference point estabilished \(\mathrm{P_{N}}\):

(2)\[u^{*} = \kappa u_{t,\mathrm{P_{N}}} \left[\mathrm{ln}\left(\frac{z_{\mathrm{P_{N}}}}{z_{0}}\right)\right]^{-1}\]

Note

One great use case for this approach is terrain category or roughness models. Because the roughness of the terrain can be modeled by varying the \(z_{0}\) value. Further details regarding how \(z_{0}\) models the terrain profile can be found in Atmospheric Boundary Layer

Thin boundary layer (TBL)

A different approach used for a smooth wall representation consists on modeling the region between \(\mathrm{P_{0}}\) and \(\mathrm{P_{N}}\) through the thin boundary layer equation (TBL) for normal \(n\) and tangential \(t\) directions (Suga et al.[2]):

(3)\[\frac{\partial}{\partial x_{n}}\left(\left(\nu+\nu_{\mathrm{e}}\right)\frac{\partial u_{t}}{\partial x_{n}}\right)=\frac{\partial u_{t}}{\partial t}+\frac{\partial}{\partial x_{n}}\left(u_{n} u_{t}\right)+\frac{1}{\rho}\frac{\partial p}{\partial x_{t}}\]

where \(\nu_{\mathrm{e}}\) is the eddy-viscosity, which we assume to follow the Johnson-King mixing-length model with Van Driest damping function (Johnson and King[3]):

(4)\[\nu_{\mathrm{e}}=\nu\kappa y^{+}\left(1 - e^{-\frac{y^{+}}{A^{+}}}\right)^{2}\]

where \(\kappa=0.41\) is the von Karman constant, and \(A^{+}=17\). For the implemented models, the TBL is simplified to:

(5)\[\frac{\partial}{\partial x_{n}}\left(\left(\nu+\nu_{\mathrm{e}}\right)\frac{\partial u_{t}}{\partial x_{n}}\right)=+\frac{1}{\rho}\frac{\partial p}{\partial x_{t}}\]

then:

(6)\[\frac{\partial \nu_{\mathrm{T}}}{\partial x_{n}}\frac{\partial u_{t}}{\partial x_{n}}+\nu_{\mathrm{T}}\frac{\partial^{2}u_{t}}{\partial x_{n}^{2}} = +\frac{1}{\rho}\frac{\partial p}{\partial x_{t}}\]

The distance \(h_{\mathrm{wm}}\) between \(\mathrm{P_{0}}\) and \(\mathrm{P_{N}}\) is divided in \(N\) points to which the EDP above is solved through a finite difference scheme:

(7)\[\left[\nu_{\mathrm{T}}-\frac{h_{\mathrm{wm}}}{2N}\frac{\partial \nu_{\mathrm{T}}}{\partial x_{n}}\right]^{\left(i\right)}u_{t}^{\left(i-1\right)} - 2\nu_{\mathrm{T}}^{\left(i\right)}u_{t}^{\left(i\right)}+\left[\nu_{\mathrm{T}}+\frac{h_{\mathrm{wm}}}{2N}\frac{\partial \nu_{\mathrm{T}}}{\partial x_{n}}\right]^{\left(i\right)}u_{t}^{\left(i+1\right)} = \left(\frac{h_{\mathrm{wm}}}{N}\right)^{2}\frac{1}{\rho}\frac{\partial p}{\partial x_{t}}^{(i)}\]

where \(\nu_{\mathrm{T}}=\nu+\nu_{\mathrm{e}}\)

(8)\[\frac{\partial \nu_{\mathrm{T}}}{\partial x_{n}}=\kappa u^{*}\left[1-\frac{2\left(A^{+}-x_{n}^{+}\right)}{A^{+}}e^{-\frac{x_{n}^{+}}{A^{+}}}+\frac{\left(A^{+}-2x_{n}^{+}\right)}{A^{+}}e^{-\frac{2x_{n}^{+}}{A^{+}}}\right]\]

The value for \(x_{n}^{+}\) in a \(i\) point is then \(x_{n}^{+}=\left(h_{\mathrm{wm}}i/N\right)\left(u^{*}/\nu\right)\).

Important

For equilibrium models, \(\partial p / \partial x_{t}\) is simply set to zero, while non-equilibrium models consider this on the equation.

Equilibrium models are suited for cases where the pressure’s gradients is small. For cases and places where this gradient is not negligible, a non-equilibrium model is advised.

The set of velocity equations generated can be solved with a TDMA (tridiagonal matrix algorithm) to which the velocity equation may be written as:

(9)\[a_{\left(i\right)}u^{\left(i-1\right)} + b_{\left(i\right)}u^{\left(i\right)} + c_{\left(i\right)}u^{\left(i+1\right)} = d_{\left(i\right)}\]

since \(u^{\left(0\right)}=0\), \(a_{1}=0\) and:

(10)\[a_{\left(i\right)} = \left[\nu_{\mathrm{T}}-\frac{h_{\mathrm{wm}}}{2N}\frac{\partial \nu_{\mathrm{T}}}{\partial x_{n}}\right]^{\left(i\right)}\]

(11)\[b_{\left(i\right)} = -2\nu_{\mathrm{T}}^{\left(i\right)}\]

(12)\[c_{\left(i\right)} = \left[\nu_{\mathrm{T}}+\frac{h_{\mathrm{wm}}}{2N}\frac{\partial \nu_{\mathrm{T}}}{\partial x_{n}}\right]^{\left(i\right)}\]

(13)\[d_{\left(i\right)} = \left(\frac{h_{\mathrm{wm}}}{N}\right)^{2}\frac{1}{\rho}\frac{\partial p}{\partial x_{t}}^{(i)}\]

For \(2\leq i \leq N-2\). We consider that \(u^{\left(N\right)}=u_{t,\mathrm{P_{N}}}\) and \(c_{\left(N-1\right)}=0\), hence:

(14)\[d_{\left(N-1\right)} = \left(\frac{h_{\mathrm{wm}}}{N}\right)^{2}\frac{1}{\rho}\frac{\partial p}{\partial x_{t}}^{\left(N-1\right)} - \left[\nu_{\mathrm{T}}+\frac{h_{\mathrm{wm}}}{2N}\frac{\partial \nu_{\mathrm{T}}}{\partial x_{n}}\right]^{\left(N-1\right)}u_{t,\mathrm{P_{N}}}\]

The matrix of of a linear system is then written as:

\[ \left[ \begin{array}{ @{}% no padding l@{\quad}% some padding r@{}% no padding >{{}}r@{}% no padding >{{}}l@{}% no padding } b_{\left(1\right)} & c_{\left(1\right)} & 0 & 0 & ... & 0 & 0 \\ a_{\left(2\right)} & b_{\left(2\right)} & c_{\left(2\right)} & 0 & ... & 0 & 0 \\ 0 & a_{\left(3\right)} & b_{\left(3\right)} & c_{\left(3\right)} & ... & 0 & 0 \\ 0 & 0 & a_{\left(4\right)} & b_{\left(4\right)} & ... & 0 & 0 \\ 0 & 0 & 0 & ... & a_{\left(N-2\right)} & b_{\left(N-2\right)} & c_{\left(N-2\right)} \\ 0 & 0 & 0 & ... & 0 & a_{\left(N-1\right)} & b_{\left(N-1\right)} \end{array} \right] \left[ \begin{array}{ @{}% no padding l@{\quad}% some padding r@{}% no padding >{{}}r@{}% no padding >{{}}l@{}% no padding } u^{\left(1\right)} \\ u^{\left(2\right)} \\ u^{\left(3\right)} \\ u^{\left(4\right)} \\ u^{\left(N-2\right)} \\ u^{\left(N-1\right)} \end{array} \right] = \left[ \begin{array}{ @{}% no padding l@{\quad}% some padding r@{}% no padding >{{}}r@{}% no padding >{{}}l@{}% no padding } d_{\left(1\right)} \\ d_{\left(2\right)} \\ d_{\left(3\right)} \\ d_{\left(4\right)} \\ d_{\left(N-2\right)} \\ d_{\left(N-1\right)} \end{array} \right] \]

For the first term, \(i=1\):

(15)\[c_{\left(1\right)}^{'}=\frac{c_{\left(1\right)}}{b_{\left(1\right)}}\]

(16)\[d_{\left(1\right)}^{'}=\frac{d_{\left(1\right)}}{b_{\left(1\right)}}\]

After which, the following coefficients up to \(N-1\) are calculated:

(17)\[c_{\left(i\right)}^{'}=\frac{c_{\left(i\right)}}{b_{\left(i\right)}-a_{\left(i\right)}c_{\left(i-1\right)}^{'}}\]

(18)\[d_{\left(i\right)}^{'}=\frac{d_{\left(i\right)}-a_{\left(i\right)}d_{\left(i-1\right)}^{'}}{b_{\left(i\right)}-a_{\left(i\right)}c_{\left(i-1\right)}^{'}}\]

Knowing that \(u_{t}^{\left(N\right)}=u_{t,\mathrm{P_{N}}}\), the remaining velocities are calculated backwards by:

\[u_{t}^{\left(i\right)}=d_{\left(i\right)}^{'}-c_{\left(i\right)}^{'}u_{t}^{\left(i+1\right)}\]

Along these lines, the friction velocity can be calculated with a second-order finnite difference scheme:

(19)\[u^{*}=\sqrt{\nu \frac{\left(-3u_{t}^{(0)}+4 u_{t}^{(1)}-u_{t}^{(2)}\right)}{2\rho_{\mathrm{w}}h_{\mathrm{wm}}/N}}\]

Important

The implementation of wall models within IBM is directed on correcting the the near-wall regions velocities as described in IBM Wall Model.

[1]

Peter Alan Davidson. Turbulence: an introduction for scientists and engineers. Oxford university press, 2015.

[2]

Kazuhiko Suga, Tomoki Sakamoto, and Yusuke Kuwata. Algebraic non-equilibrium wall-stress modeling for large eddy simulation based on analytical integration of the thin boundary-layer equation. Physics of Fluids, 2019.

[3]

Dennis A Johnson and LS King. A mathematically simple turbulence closure model for attached and separated turbulent boundary layers. AIAA journal, 23(11):1684–1692, 1985.