.. _wall_model_wall_functions: ************* 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 (:footcite:t:`davidson2015turbulence`): .. math:: u\left(z\right)=\frac{u^{*}}{\kappa} \mathrm{ln}\left(\frac{z}{z_{0}}\right) :label: wm_velocity_terrain where :math:`\kappa=0.41` is the Von Kármán constant and :math:`z_{0}` the roughness length. In a numerical simulation, the friction velocity :math:`u^{*}` can be calculated from the tangential velocity at the reference point estabilished :math:`\mathrm{P_{N}}`: .. math:: u^{*} = \kappa u_{t,\mathrm{P_{N}}} \left[\mathrm{ln}\left(\frac{z_{\mathrm{P_{N}}}}{z_{0}}\right)\right]^{-1} :label: wm_friction_velocity_log .. 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 :math:`z_{0}` value. Further details regarding how :math:`z_{0}` models the terrain profile can be found in :ref:`Atmospheric Boundary Layer ` ========================= Thin boundary layer (TBL) ========================= A different approach used for a smooth wall representation consists on modeling the region between :math:`\mathrm{P_{0}}` and :math:`\mathrm{P_{N}}` through the thin boundary layer equation (TBL) for normal :math:`n` and tangential :math:`t` directions (:footcite:t:`suga2019algebraic`): .. math:: \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}} :label: wm_tbl where :math:`\nu_{\mathrm{e}}` is the eddy-viscosity, which we assume to follow the Johnson-King mixing-length model with Van Driest damping function (:footcite:t:`johnson1985mathematically`): .. math:: \nu_{\mathrm{e}}=\nu\kappa y^{+}\left(1 - e^{-\frac{y^{+}}{A^{+}}}\right)^{2} :label: wm_eddy_viscosity where :math:`\kappa=0.41` is the von Karman constant, and :math:`A^{+}=17`. For the implemented models, the TBL is simplified to: .. math:: \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}} :label: wm_tbl_simpl then: .. math:: \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}} :label: wm_tbl_simpl2 The distance :math:`h_{\mathrm{wm}}` between :math:`\mathrm{P_{0}}` and :math:`\mathrm{P_{N}}` is divided in :math:`N` points to which the EDP above is solved through a finite difference scheme: .. math:: \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)} :label: wm_tbl_finite_diff where :math:`\nu_{\mathrm{T}}=\nu+\nu_{\mathrm{e}}` .. math:: \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] :label: wm_derivative_eddy The value for :math:`x_{n}^{+}` in a :math:`i` point is then :math:`x_{n}^{+}=\left(h_{\mathrm{wm}}i/N\right)\left(u^{*}/\nu\right)`. .. important:: For **equilibrium models**, :math:`\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 equation :math:numref:`wm_tbl_finite_diff` is modeled linearly between the points :math:`\mathrm{P_{0}}` and :math:`\mathrm{P_{N}}`: .. .. math:: .. \frac{\partial p}{\partial x_{t}}^{(i)}=\left(1-\frac{i}{N}\right)\frac{\partial p}{\partial x_{t}}^{(0)}+\frac{i}{N}\frac{\partial p}{\partial x_{t}}^{(N)} .. :label: wm_pressure_grad_lin .. where the pressure gradients at :math:`\mathrm{P_{0}}` and :math:`\mathrm{P_{N}}`, for the non-equilibrium model, are estimated with: .. .. math:: .. \frac{\partial p}{\partial x_{t}}= \frac{\left[\rho\left(x_{t}+\Delta x_{t}\right)-\rho\left(x_{t}-\Delta x_{t}\right)\right]c_{s}^{2}}{2\Delta x_{t}} .. :label: wm_pressure_grad_approx .. For simplification, we can assume :math:`\frac{\partial p}{\partial x_{n}} \approx 0`. Consequently :math:`\frac{\partial p}{\partial x_{t}} \approx \frac{\partial p}{\partial x_{t}}|_{\mathrm{P_{N}}}=\mathrm{cte}`. The set of velocity equations generated can be solved with a TDMA (tridiagonal matrix algorithm) to which the velocity equation may be written as: .. math:: 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)} :label: wm_tdma_standard since :math:`u^{\left(0\right)}=0`, :math:`a_{1}=0` and: .. math:: 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)} :label: wm_tdma_a .. math:: b_{\left(i\right)} = -2\nu_{\mathrm{T}}^{\left(i\right)} :label: wm_tdma_b .. math:: 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)} :label: wm_tdma_c .. math:: d_{\left(i\right)} = \left(\frac{h_{\mathrm{wm}}}{N}\right)^{2}\frac{1}{\rho}\frac{\partial p}{\partial x_{t}}^{(i)} :label: wm_tdma_d For :math:`2\leq i \leq N-2`. We consider that :math:`u^{\left(N\right)}=u_{t,\mathrm{P_{N}}}` and :math:`c_{\left(N-1\right)}=0`, hence: .. math:: 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}}} :label: wm_tdma_dnminus The matrix of of a linear system is then written as: .. math:: :nowrap: \[ \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, :math:`i=1`: .. math:: c_{\left(1\right)}^{'}=\frac{c_{\left(1\right)}}{b_{\left(1\right)}} :label: wm_tdma_c1 .. math:: d_{\left(1\right)}^{'}=\frac{d_{\left(1\right)}}{b_{\left(1\right)}} :label: wm_tdma_d1 After which, the following coefficients up to :math:`N-1` are calculated: .. math:: c_{\left(i\right)}^{'}=\frac{c_{\left(i\right)}}{b_{\left(i\right)}-a_{\left(i\right)}c_{\left(i-1\right)}^{'}} :label: wm_tdma_cn .. math:: 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)}^{'}} :label: wm_tdma_dn Knowing that :math:`u_{t}^{\left(N\right)}=u_{t,\mathrm{P_{N}}}`, the remaining velocities are calculated backwards by: .. math:: 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: .. math:: u^{*}=\sqrt{\nu \frac{\left(-3u_{t}^{(0)}+4 u_{t}^{(1)}-u_{t}^{(2)}\right)}{2\rho_{\mathrm{w}}h_{\mathrm{wm}}/N}} :label: wm_friction_velocity_tdma .. important:: The implementation of wall models within IBM is directed on correcting the the near-wall regions velocities as described in :ref:`IBM Wall Model `. .. footbibliography::