ABL.jl: Turbulence-Resolving Simulations of Atmospheric Boundary-Layer Flows

ABL.jl can simulate a range of physical problems, but its primary purpose is to simulate the evolution of turbulent flows in simplified representations of the atmospheric boundary-layer.

The workflow generally consists of defining a Model by configuring the desired computational domain, resolution of the discretization, and physical processes. This model can then be evolved in time using evolve!, collecting output data along the way. The package also contains functionality to help process and analyze the results.

The workflow is described in more detail on the Setup & Workflow page.

Computational Domain

ABL.jl computes solutions in a semiperiodic Cartesian domain of size $L_1 \times L_2 \times L_3$, i.e.

\[q(x_1, x_2, x_3) = q(x_1+L_1, x_2, x_3) = q(x_1, x_2+L_2, x_3)\]

where $x_i$ are the coordinates and $q$ is any physical quantity defined in the domain. Usually $x_1$ is the primary horizontal (streamwise) coordinate, $x_2$ is the secondary horizontal (cross-stream) coordinate, and $x_3$ is the vertical coordinate.

To allow for grid transforms, the $x_3$-coordinates are defined as a function of $\zeta \in \left[ 0, 1 \right]$, e.g. $x_3(\zeta) = L_3 \zeta$.

The computational domain is therefore defined as

\[(x_1, x_2, x_3) \in (-\infty, +\infty) \times (-\infty, +\infty) \times \left[x_3(\zeta=0), x_3(\zeta=1)\right] \,.\]

Boundary conditions have to be specified at $x_3(\zeta=0)$ and $x_3(\zeta=1)$.

The computational domain including boundary conditions is configured with the SemiperiodicDomain type.

Discretization

In this semiperiodic domain, a quantity $q$ can be represented as the Fourier series

\[q(x_1, x_2, x_3) = \sum_{ \begin{subarray}{l} -\infty < \kappa_1 < \infty \\ -\infty < \kappa_2 < \infty \end{subarray}} \hat{q}^{\kappa_1\kappa_2}\!\left(x_3(\zeta)\right) \,\mathrm{e}^{i 2 \pi \kappa_1 x_1 / L_1} \,\mathrm{e}^{i 2 \pi \kappa_2 x_2 / L_2} \,.\]

To compute numerical solutions, the series is truncated to $\left| \kappa_1 \right| < N_1/2$ and $\left| \kappa_2 \right| < N_2/2$, where the integers $N_i$ are prescribed to define the resolution. Note that the number of resolved wavenumbers is always odd even when $N_i$ are even numbers, as the wavenumbers are symmetric around $\kappa_i=0$.

For $\zeta$, the range $[0, 1]$ is split into $N_3$ sections. Numerical solutions are then computed either at the center $\zeta_C$ of each section or at the interfaces $\zeta_I$ between the sections, with

\[\zeta_C \in \left\{ \frac{1/2}{N_3}, \frac{3/2}{N_3}, \dots, \frac{N_3-1/2}{N_3} \right\} \quad \text{and} \quad \zeta_I \in \left\{ \frac{1}{N_3}, \frac{2}{N_3}, \dots, \frac{N_3-1}{N_3} \right\} \,.\]

A staggered arrangement can reduce discretization errors and avoid instabilities due to odd–even decoupling.

The derivatives $\partial/\partial x_1$ and $\partial/\partial x_2$ can be computed exactly. For $\partial/\partial x_3$, the derivatives are approximated with central second-order finite differences, i.e.

\[\frac{\partial \hat{q}^{\kappa_1 \kappa_2}}{\partial x_3} \bigg|_\zeta = \frac{\mathrm{d}\zeta}{\mathrm{d}x_3} \bigg|_\zeta \frac{\mathrm{d}\hat{q}^{\kappa_1\kappa_2}}{\mathrm{d}\zeta} \bigg|_\zeta = \frac{\mathrm{d}\zeta}{\mathrm{d}x_3} \bigg|_\zeta \frac{\hat{q}^{\kappa_1\kappa_2}(\zeta+\delta\zeta/2) - \hat{q}^{\kappa_1\kappa_2}(\zeta-\delta\zeta/2)}{\delta\zeta}\]

with $\delta\zeta=1/N_3$. The derivatives of the grid transform $x_3(\zeta)$ are computed analytically.

The configuration options for the resolution are described in the documentation of the Model. Some processes have additional settings that influence the precision with which they are discretized, such as the number of collocation points.

Physical Processes

The model simulates the evolution of a number of quantities $q_i$ that are each represented by a vector $\bm{\hat{q}_i}$ of the $N_1 N_2 N_3$ modal/nodal values $\hat{q}_i^{\kappa_1\kappa_2\zeta}$ ($N_3-1$ for quantities defined at $\zeta_I$-nodes). Together they make up the state $\bm{\hat{s}} = \begin{pmatrix} \bm{\hat{q}_1}, \bm{\hat{q}_2}, \dots \end{pmatrix}^\intercal$. The state generally consists of the three components of the velocity field plus perhaps a number of scalar fields.

The evolution of this state is governed by an equation in the form of

\[\frac{\partial \bm{\hat{s}}}{\partial t} = \sum_i P_i(\bm{\hat{s}}, t) \,,\]

where $P_i$ are functions that describe the contribution of each physical process (molecular diffusion, advection, etc.) to the rate of change of the state.

The processes that can be passed to the Model and their configuration are described here.

Time Evolution