Analysis of Cerebrospinal Fluid Dynamics - II. LBM

Cerebrospinal fluid (CSF) dynamics in the human central nerve system (CNS) is an extremely interesting topic in bio-fluid mechanics. The disorders of the CNS components disrupt the CSF pulsatility which is involved in a number of vital functions. The CFD analysis can be used to quantify the functional interactions of the CNS member. The numerical simulations of CSF filled in brain ventricles and subarachnoid spaces (SAS) are performed to build up a basis for testing and optimizing the CNS therapeutics via understanding of the CSF pathophysiology. In this post a recent collaboration of the SimLab - Highly Scalable Fluids and Solids Engineering at the Jülich Supercomputing Centre with Korea University Medical Center is presented.

Numerical Methods and Algorithms

In this project, several software components will be utilized to develop HPC- and CFD services for the medical community. That is, for CFD computations the modular framework ZFS, which is primarily developed at AIA, will be used. Interactive supercomputing, i.e., black-boxed access to HPC resources and simulation results will be provided by an intuitive Jupyter Hub environment. Virtual surgeries and hence in-situ computational steering will be implemented by interfacing ZFS with a web-based application to modify the anatomical shape. In the context of ML, the package Tensorflow will be employed. Since the lattice-Boltzmann method in ZFS will consume most of the computational resources, the following discussion concentrates on the explanation of the involved numerical algorithms.

ZFS consists of a parallel Cartesian grid generator and several solver modules to account for different physics. The code is written in C++, hybrid MPI/OpenMP parallelized and has in parts been ported to GPGPUs (Kraus et al. 2014, Nicolini et al. 2016) and Intel Xeon Phi Knights Landing platforms. Parallel I/O makes use of HPC I/O libraries (Li et al. 2003, Folk & Pourmal 2010) .

The computational domain is discretized by a hierarchical Cartesian octree mesh, which is generated by the method described in Lintermann et al. 2014. Initially, a single cubic cell containing the whole computational domain is generated on every process, which is uniformly refined up to a level $l_\alpha$ by isotropic subdivision into eight cubic child cells. Cells outside the geometry are removed from the tree. On level $l_\alpha$ all cells on $l<l_\alpha$ are deleted. A Hilbert space-filling curve Sagan 1994 decomposes the remaining cells on the number of processors. In parallel, each process continues to refine the mesh uniformly up to level $l_\beta$. Individual regions of the mesh are further refined to level $l_\gamma$ to meet resolution requirements. The algorithm ensures neighboring cells on different levels to have a level difference of maximum 1. Uneven work-load distributions are accounted for by a load-balancing method. The mesh is the written to disk in parallel. A list, consisting in principle of the cells on level $l_\alpha$, is used in the simulation to optimally subdivide the problem on the computational cores. Note that the number of processes employed for the simulation can be different from that of the grid generation. Furthermore, the geometry can be parallelized using the approach from Lintermann 2016. For more information, the reader is referred to Lintermann et al. 2014, Lintermann 2016.

The LBM has been shown to be well-suited for the simulation of respiratory flows in the low Mach number and low to moderate Reynolds number regime (Freitas et al. 2008, Lintermann & Schröder 2017). It solves the lattice-Bhatnagar-Gross-Krook (LBGK) equation (Bhatnagar et al. 1954, Hänel 2004) for the particle probability distribution functions $f_i$ (PPDFs) for a direction $i$.

The PPDF $f_i$ with $i=27$ corresponds to the rest particle distribution. The D3Q15 and D3Q19 models are subsets of the D3Q27 model.

\begin{equation} f_i\left(\mathbf{x}+\mathbf{\Xi}_i\delta t,t+\delta t\right) = f_i\left(\mathbf{x},t\right) + \omega_F\delta t\cdot\left(f_i^{eq}\left(\mathbf{x},t\right)-f_i\left(\mathbf{x},t\right)\right). \end{equation}

In this equation, $\mathbf{\Xi}_i$ represents the discrete velocity vector, $\mathbf{x}$ is the spatial location, $f_i^{eq}$ is the Maxwellian equilibrium distribution function, $\omega_F$ is the particle collision frequency, $t$ is the time, and $\delta t$ is the time increment. The collision frequency is determined by

\begin{equation} \omega_F = \frac{c_s^2}{\nu+\delta t c_s^2/2} \end{equation}

with the speed of sound $c_s$ and the kinematic viscosity $\nu$ of the fluid. The Maxwellian equilibrium function is given by

\begin{equation} f_i^{eq}(r,t)=\rho \underbrace{t_p \bigg[1+\frac{v_a\xi_{i,a}}{c_s^2}+\frac{v_a v_b}{2c_s^2}\cdot\left(\frac{\xi_{i,a}\xi_{i,b}}{c_s^2}-\delta_{ab}\right)\bigg]}_{\chi}, \end{equation}

where $t_p$ is a direction dependent weight factor, $\rho$ and $v_{a,b}$ are the fluid density and velocities, and $\delta_{ab}$ is the Kronecker delta with $a,b\in{1,2,3}$. The most frequently applied discretization models to solve the LBGK equation in three dimensions are the D3Q19 and D3Q27 models (Qian et al. 1992) with 19 and 27 discrete PPDF directions, respectively. Since former model is rotational not invariant (White & Chong 2011, Kuwata & Suga 2015), the D3Q27 model is preferred for the computations in this project. The figure exemplarily shows the corresponding directions of the D3Q27 model, with the D3Q15 and D3Q19 being subsets of the D3Q27.

Since for quasi-incompressible flow, a decoupling of the energy equation from the Navier-Stokes equation is obtained, the LBGK equation cannot be used to solve for the temperature. Therefore, a multi-distribution function (MDF) approach is used to solve a passive scalar equation for the temperature. This means, it is solved for an additional set of PPDFs via the TLBGK equation

\begin{equation} g_i\left(\mathbf{x}+\mathbf{\Xi}_i\delta t,t+\delta t\right) = g_i\left(\mathbf{x},t\right) + \omega_T\delta t\cdot\left(g_i^{eq}\left(\mathbf{x},t\right)-g_i\left(\mathbf{x},t\right)\right), \end{equation}

where $\omega_T$ is a function of the Prandtl number and $g_i^{eq}=T\cdot\chi$ is the equilibrium distribution function for the temperature.

The no-slip boundary condition employs the interpolated bounce-back formulation by Bouzidi et al., which is second-order accurate. The wall temperature is set equal to body temperature $T_w=309.15$K. At the nostrils and at the pharynx the combined pressure and Saint-Vernant/Wanzel boundary condition presented in Lintermann et al. 2014 is used, i.e., the pressure at the pharynx is continuously reduced until a target Reynolds number is reached at the pharynx, and the density and velocity at the nostrils are recalculated by an extrapolation of the momentum from the next inner cell layer using the LBM-reformulated equation of Saint-Vernant/Wanzel. The temperature at the nostrils is set to $T_n=293.15$K.