(wall_model_tbl_solve)=

# Solving the TBL equation

The {ref}`thin boundary layer model <wall_model_wall_functions>` leaves us with a second-order ordinary differential equation for the tangential velocity along the wall normal.
This page works through its discretization on the wall-normal stencil and its solution by the tridiagonal matrix algorithm (TDMA), ending with the recovery of the friction velocity that the wall model feeds back to the solver.

## Finite-difference stencil

The distance $h_{\mathrm{wm}}$ between $\mathrm{P_{0}}$ and $\mathrm{P_{N}}$ is divided in $N$ points to which the simplified TBL equation is solved through a finite difference scheme:

$$
\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)}
$$ (wm_tbl_finite_diff)

where $\nu_{\mathrm{T}}=\nu+\nu_{\mathrm{e}}$

$$
\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]
$$ (wm_derivative_eddy)

The $\nu$ factor in $\nu_{\mathrm{e}}$ cancels because the derivative is taken with respect to the physical wall-normal coordinate $x_{n}$ while $\nu_{\mathrm{e}}$ is expressed in wall units $x_{n}^{+}=x_{n}u^{*}/\nu$, which is why the prefactor reduces to $\kappa u^{*}$.

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

```{figure} /_static/img/theory/wall_model/tbl_fd_stencil.svg
---
align: center
width: 50 %
---
The wall-normal grid for the thin-boundary-layer solve. The TBL equation is discretised on N points from the wall $\mathrm{P_0}$ to the matching point $\mathrm{P_N}$; a three-point finite-difference molecule at an interior point couples it to its neighbours, producing the tridiagonal system solved by the TDMA.
```

```{admonition} Equilibrium vs non-equilibrium: the pressure-gradient term
---
class: important
---
This one term decides which model you are running. Setting $\partial p / \partial x_{t} = 0$ gives the **equilibrium** model, valid where the streamwise pressure gradient is small. Retaining it gives the **non-equilibrium** model, the one to use wherever that gradient is not negligible: separation, reattachment, and adverse-pressure-gradient regions.
```

The equilibrium lattice-BC solve described on this page sets $\partial p / \partial x_{t} = 0$, so the source terms $d_{\left(i\right)}$ vanish. The non-equilibrium model evaluates a single tangential pressure gradient per node by central-differencing the LBM density along the local tangential direction, and applies it uniformly across the interior rows. Since $p = \rho c_{s}^{2}$,

$$
\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}}
$$ (wm_pressure_grad_approx)

## Tridiagonal system and TDMA

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

$$
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)}
$$ (wm_tdma_standard)

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

$$
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)}
$$ (wm_tdma_a)

$$
b_{\left(i\right)} = -2\nu_{\mathrm{T}}^{\left(i\right)}
$$ (wm_tdma_b)

$$
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)}
$$ (wm_tdma_c)

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

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

$$
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}}}
$$ (wm_tdma_dnminus)

The matrix of a linear system is then written as:

$$
\begin{bmatrix}   b_{\left(1\right)} & c_{\left(1\right)} & 0 & 0 & \cdots & 0 & 0 \\   a_{\left(2\right)} & b_{\left(2\right)} & c_{\left(2\right)} & 0 & \cdots & 0 & 0 \\   0 & a_{\left(3\right)} & b_{\left(3\right)} & c_{\left(3\right)} & \cdots & 0 & 0 \\   0 & 0 & a_{\left(4\right)} & b_{\left(4\right)} & \cdots & 0 & 0 \\   0 & 0 & 0 & \cdots & a_{\left(N-2\right)} & b_{\left(N-2\right)} & c_{\left(N-2\right)} \\   0 & 0 & 0 & \cdots & 0 & a_{\left(N-1\right)} & b_{\left(N-1\right)} \end{bmatrix} \begin{bmatrix}   u^{\left(1\right)} \\   u^{\left(2\right)} \\   u^{\left(3\right)} \\   u^{\left(4\right)} \\   \vdots \\   u^{\left(N-1\right)} \end{bmatrix} = \begin{bmatrix}   d_{\left(1\right)} \\   d_{\left(2\right)} \\   d_{\left(3\right)} \\   d_{\left(4\right)} \\   \vdots \\   d_{\left(N-1\right)} \end{bmatrix}
$$

For the first term, $i=1$:

$$
c_{\left(1\right)}^{'}=\frac{c_{\left(1\right)}}{b_{\left(1\right)}}
$$ (wm_tdma_c1)

$$
d_{\left(1\right)}^{'}=\frac{d_{\left(1\right)}}{b_{\left(1\right)}}
$$ (wm_tdma_d1)

After which, the following coefficients up to $N-1$ are calculated:

$$
c_{\left(i\right)}^{'}=\frac{c_{\left(i\right)}}{b_{\left(i\right)}-a_{\left(i\right)}c_{\left(i-1\right)}^{'}}
$$ (wm_tdma_cn)

$$
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)}^{'}}
$$ (wm_tdma_dn)

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 finite-difference scheme:

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

## Wall-clustered grid

The stencil above places the $N$ points at a uniform physical spacing $h_{\mathrm{wm}}/N$. Because the wall model is driven by the wall-normal gradient at $\mathrm{P_0}$, resolution is best spent near the wall. Nassu therefore allows the grid to be clustered toward the wall while keeping the tridiagonal (Thomas) solver untouched: the spacing variation is absorbed into a coordinate map, not a wider stencil.

The solve is carried out on a uniform **computational coordinate** $\xi \in \left[0, 1\right]$ with constant spacing $\Delta\xi = 1/N$, related to the physical wall-normal distance $x_n$ by an exponential map clustered at the wall ($\xi = 0$):

$$
x_{n}\left(\xi\right) = h_{\mathrm{wm}}\,\frac{e^{\beta\xi}-1}{e^{\beta}-1}, \qquad J\left(\xi\right) \equiv \frac{\mathrm{d}x_{n}}{\mathrm{d}\xi} = h_{\mathrm{wm}}\,\frac{\beta\, e^{\beta\xi}}{e^{\beta}-1}
$$ (wm_stretch_map)

The clustering strength $\beta \ge 0$ is the `TDMA_stretch_beta` parameter (default $\beta = 2$, which places the first node inside the viscous sublayer for a typical near-wall cell of 10-50 wall units); $\beta \to 0$ recovers the uniform grid $x_n = h_{\mathrm{wm}}\,\xi$ with $J = h_{\mathrm{wm}}$, and the wall-unit position becomes $x_{n}^{+}=x_{n}\left(\xi_i\right)\left(u^{*}/\nu\right)$ in {eq}`wm_derivative_eddy`.

Rather than expand the variable-coefficient operator, the TBL equation is discretised in its **conservative form**. Mapping $\partial_{x_n}\!\left(\nu_{\mathrm{T}}\,\partial_{x_n} u_t\right)$ into $\xi$ gives

$$
\frac{\partial}{\partial \xi}\!\left(\frac{\nu_{\mathrm{T}}}{J}\frac{\partial u_t}{\partial \xi}\right) = J\,\frac{1}{\rho}\frac{\partial p}{\partial x_{t}}
$$ (wm_stretch_conservative)

A finite-volume balance on the uniform $\xi$-grid, with the flux coefficient $\Gamma = \nu_{\mathrm{T}}/J$ evaluated at the cell faces $i\pm\tfrac{1}{2}$, keeps the system tridiagonal:

$$
a_{\left(i\right)} = \Gamma_{i-\frac{1}{2}}, \qquad
b_{\left(i\right)} = -\left(\Gamma_{i-\frac{1}{2}}+\Gamma_{i+\frac{1}{2}}\right), \qquad
c_{\left(i\right)} = \Gamma_{i+\frac{1}{2}}, \qquad
d_{\left(i\right)} = \Delta\xi^{2}\,J^{\left(i\right)}\,\frac{1}{\rho}\frac{\partial p}{\partial x_{t}}^{\left(i\right)}
$$ (wm_stretch_coeffs)

The face viscosity $\nu_{\mathrm{T}}$ in $\Gamma_{i\pm\frac{1}{2}}$ is the van Driest value at the face wall-distance. The conservative form absorbs the $\partial\nu_{\mathrm{T}}/\partial x_n$ term that appears explicitly in {eq}`wm_tdma_a`-{eq}`wm_tdma_c`, so no separate eddy-viscosity-derivative coefficient is needed. The forward/backward Thomas sweeps {eq}`wm_tdma_cn`-{eq}`wm_tdma_dn` are identical; only the coefficient assembly changes. The friction velocity {eq}`wm_friction_velocity_tdma` carries the metric at the wall,

$$
u^{*}=\sqrt{\nu\,\frac{\left(-3u_{t}^{(0)}+4 u_{t}^{(1)}-u_{t}^{(2)}\right)}{2\,\Delta\xi\,J\left(0\right)}}
$$ (wm_friction_velocity_stretch)

which reduces to {eq}`wm_friction_velocity_tdma` for $\beta \to 0$ since $\Delta\xi\,J\left(0\right) = h_{\mathrm{wm}}/N$.

```{eval-rst}
.. footbibliography::
```