Consider a wall-stress model that is an ODE in the wall-normal direction. This would result, e.g., by assuming equilibrium (i.e., parallel flow without acceleration or any pressure gradient) and an eddy-viscosity model but without solving this ODE analytically. Our example in this section is the compressible equilibrium wall-stress model (cf. Kawai and Larsson, 2012; Larsson et al, 2016) given by the momentum equation
\begin{equation}
\label{eqn:wm_momentum}
\tag{1}
\frac{d}{dy} \left[\left( \mu_{\rm wm} + \mu_{t,{\rm wm}} \right)\frac{dU_{\rm wm}}{dy}\right] = 0 \,,
\end{equation} the energy equation
\begin{equation}
\label{eqn:wm_energy}
\tag{2}
\frac{d}{dy}\left[c_p\left( \frac{\mu_{\rm wm}}{Pr} + \frac{\mu_{t,{\rm wm}}}{Pr_{t,{\rm wm}}} \right)\frac{dT_{\rm wm}}{dy}\right]=-\frac{d}{dy} \left[\left( \mu_{\rm wm} + \mu_{t,{\rm wm}} \right)U_{\rm wm} \frac{dU_{\rm wm}}{dy}\right]\,,
\end{equation} and the eddy-viscosity
\begin{equation}
\label{eqn:mut}
\tag{3}
\mu_{t,{\rm wm}} = \kappa \rho_{\rm wm} \sqrt{\frac{\tau_{w, {\rm wm}}}{\rho_{\rm wm}}} \, y\left[ 1 – \exp \left( -y^*/A^+ \right) \right]^2\,,\\y^* = \frac{y \sqrt{\rho \tau_w}}{\mu}\,,
\end{equation} where the scaled wall-distance in the Van Driest damping function is in semi-local scaling in order to work in strongly non-adiabatic conditions (Patel et al, 2015; Trettel and Larsson, 2016).
The model parameters are taken as \(\kappa=0.41\), \(A^+=17\) and \(Pr_{t,{\rm wm}}=0.9\). The coordinate \(y\) is wall-distance, and thus the equations should be solved between \(y=0\) (the wall) and \(y=h_{\rm wm}\). The solution to the wall-stress model is labeled \(U_{\rm wm}\) and \(T_{\rm wm}\) to emphasize that it is different from the velocity and temperature in the LES code itself. Similarly, the quantities \(\mu_{\rm wm}\), \(\rho_{\rm wm}\) and \(\tau_{w, {\rm wm}}\) are computed from the solution to the wall-stress model, not from the LES solution.
The boundary conditions on the wall-stress model solution at the top are \begin{align*} U_{\rm wm} (h_{\rm wm}) &= \left| \vec{u}_{{\rm wm-top}, \parallel} \right| \,, \\ T(h_{\rm wm}) &= T_{\rm wm-top} \,, \end{align*} (i.e., Dirichlet conditions at the velocity and temperature interpolated from the LES). At the wall, the boundary conditions are exactly what they would be in a RANS solver, i.e., \(U_{\rm wm}(0) = 0\) and either a Dirichlet or Neumann condition on \(T_{\rm wm}(0)\). This completes the specification of the wall-stress model on paper; what remains is to solve this numerically.
Finite-volume grid
The wall-stress model described in Eqns. (\ref{eqn:wm_momentum})-(\ref{eqn:mut}) describes the conservation of momentum and total energy, and thus it is quite naturally solved using a finite-volume method. This is true even if the LES code itself uses a different approach! The approach used in our own work has been to use a lightly stretched finite-volume grid defined by
\begin{align*} y_{{\rm f},j} &= \Delta y_w \, \frac{ r^j – 1 }{r-1} \,, \ \ j=0,1,\ldots,n \\ y_{{\rm c},j} &= \frac{ y_{{\rm f},j} + y_{{\rm f},j+1}}{2} \,, \ \ j=0,1,\ldots,n-1
\end{align*} where \(y_{{\rm f}, j}\) and \(y_{{\rm c},j}\) are the coordinates of the \(j\)th face and cell-centroid, respectively, and \(\Delta y_w\) is the spacing at the wall. Note that the size of the \(j\)th cell is \(\Delta y_j = y_{{\rm f}, j+1} – y_{{\rm f}, j} = \Delta y_w \, r^j\) (so \(r\) is the stretching parameter).
The parameters \(\Delta y_w\) and \(r\) should be chosen to produce an acceptable error (in the wall-stress model solution itself) at minimum cost (i.e., the number of wall-stress model cells \(n\)). For the second-order numerics discussed below, we have found that the optimal values of \((\Delta y_w^+, r)\) are approximately \((0.6, 1.016)\) for 0.005% error, \((0.8, 1.025)\) for 0.01% error, \((1.2, 1.066)\) for 0.05% error, and \((1.2, 1.10)\) for 0.1% error, where the error is the estimated error in the friction coefficient (note: the actual error depends on the Reynolds number and the flow in question, so these are just guidelines).
In addition, the centroid of the final cell should be at height \(h_{\rm wm}\) above the wall in order to impose the boundary condition at the right height (i.e., \(y_{{\rm c}, n-1} = h_{\rm wm}\)). In practice, this can be achieved by adjusting the value of \(\Delta y_w\) once the grid has been built (i.e., to build the grid, and then change all coordinates by the same factor).
Note that the ODE-grid must be created in outer units (i.e., not viscous units) since we need to match up the last cell center to \(h_{\rm wm}\), which is necessarily given in outer units. We therefore need to keep track of the maximum value of \(\Delta y_w^+\) during the LES run, since this needs to be sufficiently small for our desired error level. This is not a major issue in general: the grid-stretching means that one can err on the side of too small \(\Delta y_w\) at a very small additional cost.
Solving the finite-volume equations
A finite-volume discretization of the momentum equation (\ref{eqn:wm_momentum}) is
\begin{equation}\label{eqn:wm_momentum_discrete}\tag{4}\left. \left( \mu_{\rm wm} + \mu_{t,{\rm wm}} \right) \frac{dU_{\rm wm}}{dy} \right|_{{\rm f}, j+1}-\left. \left( \mu_{\rm wm} + \mu_{t,{\rm wm}} \right) \frac{dU_{\rm wm}}{dy} \right|_{{\rm f}, j}= 0 \,, \ \ j=0,1,\ldots,n-2\,,\end{equation} which should be solved for cells \(j=0,1,\ldots,n-2\) (i.e., not for the last cell \(n-1\) where we impose the boundary condition). The derivative of the velocity is then approximated using a second-order central finite-difference, i.e., as \[\left. \frac{dU_{\rm wm}}{dy} \right|_{{\rm f}, j} \approx \frac{U_{{\rm wm},j}-U_{{\rm wm},j-1}}{y_{{\rm c},j}-y_{{\rm c},j-1}}\,, \ \ j=1,2,\ldots,n-1\,.\]
A standard one-sided approximation is used at the boundary, i.e.,\[\left. \frac{dU_{\rm wm}}{dy} \right|_{{\rm f}, w} \approx\frac{U_{{\rm wm}, 0} – U_{{\rm wm}, w} }{y_{{\rm c},0} – 0}\,.\]
The discretization of the energy equation (\ref{eqn:wm_energy}) is analogous.
The initial guesses for \(U_{\rm wm}\) and \(T_{\rm wm}\) can be either linear profiles or, if they are stored, the solutions from the prior time step. From these initial guesses one computes \(\tau_{w, {\rm wm}}\) (by a standard finite-difference in the first off-wall cell), \(\rho_{\rm wm}\) at the face locations (from the equation-of-state, assuming constant \(p\) taken from the LES), \(\mu_{\rm wm}\) at the face locations (from the viscosity law), and \(\mu_{t, {\rm wm}}\) at the face locations (from Eqn. (\ref{eqn:mut})).
In our own implementations we then solve the discretized equations in a segregated fashion, i.e., we do: (i) solve the momentum equation (\ref{eqn:wm_momentum_discrete}) treating \(\mu_{\rm wm}\) and \(\mu_{t, {\rm wm}}\) as fixed; (ii) update \(\tau_{w, {\rm wm}}\) and \(\mu_{t, {\rm wm}}\); (iii) solve the energy equation while treating everything except \(T_{\rm wm}\) as fixed; (iv) update \(\tau_{w, {\rm wm}}\), \(\mu_{\rm wm}\), \(\rho_{\rm wm}\), and \(\mu_{t, {\rm wm}}\); and (v) repeat until convergence.
When treating the ODE-system in this segregated manner, each equation is a tri-diagonal system that can be solved efficiently by the tri-diagonal matrix algorithm (TDMA). Convergence is best measured by the values of \(\tau_{w, {\rm wm}}\) and \(q_{w, {\rm wm}}\) (or \(T_{w, {\rm wm}}\)), since these are the actual outputs of the wall-model. For example, one can iterate until these values change from iteration to iteration by less than some number (we use 0.01%).
After convergence, the wall-stress model would send \(\tau_{w, {\rm wm}}\) back to the LES solver where it is used in the coupling procedure as the magnitude of the wall-parallel stress \(\left| \vec{\tau}_{w, \parallel} \right|\). In addition, either or both of the heat flux \(q_{w, {\rm wm}}\) or the wall temperature \(T_{w, {\rm wm}}\) are sent back as well.