A Fourier/Modal-Spectral-Element Method for the Simulation of High-Reynolds Number Incompressible Stratified Flows in Domains with a Single Non-Periodic Direction

TL;DR

This work addresses the challenge of simulating high-$Re$ incompressible stratified flows by developing a high-order solver that couples a Fourier pseudo-spectral method in the horizontal with a Legendre-Galerkin modal spectral-element method in the vertical. It employs a third-order implicit–explicit time stepping (KIO) and a Helmholtz/Poisson solve with static condensation and Schur-complement domain decomposition to achieve an efficient, scalable solver that yields banded system matrices. The method is validated through a sequence of 2D/3D benchmarks, including vortex-dipole collisions and stratified wakes, and is extended to realistic high-$Re$ stratified wakes behind a sphere, capturing thin shear regions and internal-wave radiation with good agreement to literature. This framework enables high-$Re$ stratified turbulence studies in long, high-aspect-ratio domains, with potential applications to oceanic and atmospheric flows and to problems involving boundary layers, jets, and mixing driven by internal waves.

Abstract

We present the components of a high-order accurate Navier-Stokes solver designed to simulate high-Reynolds-number stratified flows. The proposed numerical model addresses some of the numerical and computational challenges that high-Reynolds-number simulations pose, facilitating the reproduction of stratified turbulent fluid dynamics typically observed in oceanic and atmospheric flows, namely the development of thin regions of high vertical shear, strongly layered turbulence at high Reynolds numbers and internal wave radiation. This Navier-Stokes solver utilizes a Fourier pseudo-spectral method in the horizontal direction and a modal spectral element discretization in the vertical. We adopt an implicit-explicit time discretization scheme that involves solving several one-dimensional Helmholtz problems at each time step. Static condensation and modal boundary-adapted basis functions result in an inexpensive algorithm based on solving many small tridiagonal systems. A series of benchmark studies is presented to demonstrate the robustness of the flow solver. These include two-dimensional and three-dimensional problems, concluding with a turbulent stratified wake generated by a sphere in linear stratification.

Figures (24)

• Figure 1: First seven modes of the modal boundary-adapted basis described by Eq.(\ref{['Eq:lobattoBasis']}). Linear or vertex modes $\psi_0(r)$ and $\psi_1(r)$ are colored in red. First two bubble modes $\psi_2(r)$ and $\psi_3(r)$ are colored in green, and higher order bubble modes $\psi_k(r)$, $4\leq k \leq 6$ in blue. Color code is introduced to facilitate the interpretation of the matrix sparsity patterns in Fig. \ref{['sparsityPatterns']}.
• Figure 2: Sparsity patterns of the global second order stiffness matrix ${K}$ (top) and global mass matrix ${M}$ (bottom) using the local modal boundary-adapted basis defined in Eq. (\ref{['Eq:lobattoBasis']}), and four elements with nine modes per element. Left: Lexicographic order. Right: separating interior from boundary modes. Dots indicate terms corresponding to edge modes (red), interior modes (blue), and the coupling terms between them (green).
• Figure 3: Vorticity field of translating vortex dipole interaction with bottom boundary: initial condition $a)~t=0$, before first collision $(b)~t=0.25$, right after sheet with high-amplitude vorticity detaches from the boundary $(c)~t=0.36$, and before the second collision. The Reynolds number, $Re \sim 1/\nu$, equals to $Re=2500$. Element interfaces indicated for the $256 \times 256$ case with grid parameters as indicated in Table \ref{['grid_parameters_vortex_xz']}.
• Figure 4: Vorticity contour plots of the dipole-wall collision at $t=0.51$ for $Re=1250$ (a), $2500$ (b), $2500$ (b), $5000$ (c), and $10000$ (d). Only a fraction of the computational domain is presented, with $-0.5 \leq x \leq 1.0$ and $-1.0 \leq z \leq 0.0$. Spatial resolution corresponds to $1024 \times 1024$ case from Table \ref{['grid_parameters_vortex_xz']}. Contour levels of non-dimensional vorticity in the range $\{-260, 260 \}$ values drawn every $50$, where pink colors indicate negative values, and green colors positive ones. As $Re$ grows, the thickness of the vorticity filament at the bottom boundary decreases. At $x=-0.2$, the estimated thickness of the filament is: (a) $25.6\times 10^{-3}$, (b) $19.5\times 10^{-3}$, (c) $14.0\times 10^{-3}$, and (d) $10.10\times 10^{-3}$. In the latter, the spatial resolution guarantees at least $14$ grid points, in the vertical direction, to represent this filament (specifically at $x=-0.2$).
• Figure 5: Total enstrophy vs time for the vortex dipole test case different $Re$ and spatial resolutions: $256 \times 256$ (dotted line), $512 \times 512$ (dashed line) and $1024 \times 1024$ (solid line). The collisions of the primary vortex with the no-slip wall generate peaks in the evolution of the Enstrophy.