Vielbein Lattice Boltzmann approach for fluid flows on spherical surfaces

TL;DR

The paper tackles the challenge of simulating compressible fluid flows on curved manifolds, specifically spherical surfaces, where traditional Cartesian LBMs struggle to incorporate curvature. It introduces a vielbein-based lattice Boltzmann method (vLBM) that maps velocity space to local vielbein components $v^{\hat{\theta}}$ and $v^{\hat{\varphi}}$, and employs Gauss–Hermite quadrature to exactly recover hydrodynamic moments up to order $Q-1$. The method combines BGK collisions, Hermite-based velocity-space discretization, flux-based advection (including WENO-5 and Komissarov schemes), and pole boundary conditions, with axisymmetric analytical NS solutions used as benchmarks; it demonstrates accurate propagation of sound and shear waves, a Sod-like shock, and physically meaningful vortex dynamics on the sphere. The work also analyzes isotropy and provides performance insights on GPUs, outlining future IMEX/semi-Lagrangian extensions for large-scale turbulent simulations on curved surfaces and potential applications to geophysical and biomembrane systems.

Abstract

In this paper, we develop a lattice Boltzmann scheme based on the Vielbein formalism for the study of fluid flows on spherical surfaces. The Vielbein vector field encodes all details related to the geometry of the underlying spherical surface, allowing the velocity space to be treated as on the Cartesian space. The resulting Boltzmann equation exhibits inertial (geometric) forces that ensure that fluid particles follow paths that remain on the spherical manifold, which we compute by projection onto the space of Hermite polynomials. Due to the point-dependent nature of the advection velocity in the polar coordinate $θ$ , exact streaming is not feasible, and we instead employ finite-difference schemes. We provide a detailed formulation of the lattice Boltzmann algorithm, with particular attention to boundary conditions at the north and south poles. We validate our numerical implementation against two analytical solutions of the Navier-Stokes equations derived in this work: the propagation of sound and shear waves. Additionally, we assess the robustness of the scheme by simulating the compressible flow of an axisymmetric shock wave and analyzing vortex dynamics on the spherical surface.

Figures (13)

• Figure 1: Spherical surface parametrized by the angles $\theta\in[0,\pi]$ and $\varphi\in[0,2\pi)$, with radius $R$. The vector $\vb{r}$ marks the position vector for a generic point on the spherical surface.
• Figure 2: North pole view (with the letter N indicating the pole) of the boundaries defined by Eq. \ref{['eq:bcs_theta']} for $\theta=0$ in the case of positive velocity. The black arrows represent the advection velocities at different grid points. In red we highlight the components of a velocity vector traveling towards the north pole. Notice that the "translated" advection velocities have the opposite sign compared to those of the population residing at the same node (black arrows).
• Figure 3: Full solution for the velocity profile of a sound wave with initial velocity $u^{\hat{\theta}}_0(\theta) = U_0 \theta(\pi - \theta)$. The velocity is normalized by the initial amplitude $U_0$, while the coordinate is shown as $\theta/\pi$. The numerical data (colored symbols) is shown to match the analytical solution (black lines) at every time $t$ (in lattice units) considered.
• Figure 4: Evolution of different amplitudes for the sound wave solution. In the left panel, the numerical data (colored symbols) for the modes $n=1,2,3$ is plotted as a function of time, matching the analytical solution (black line) everywhere. In the right panel, the exponential decay of the fastest decaying mode considered ($n=3$) is made apparent by considering a larger time window. The solutions are normalized by the initial velocity amplitude $U_0$ and by the integration constant.
• Figure 5: Convergence test for the sound wave solution. The L2 error evaluated between the analytical and numerical solution for (from left to right) the $n=1$, $n=2$ and $n=3$ modes, using different advection schemes, is plotted against the grid size $N_\theta$. U1, U2 and U3 indicate the first, second and third-order upwind schemes, while W5 is the weno5 scheme. Dotted lines refer to simulations in which we have employed the Komissarov scheme.