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