# Wall-stress model: ODE

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

\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 \,, the energy 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]\,, and the eddy-viscosity

\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}\,, 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
$$\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\,,$$ 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.