A parallel fluid solid coupling model using LAMMPS and Palabos based on the immersed boundary method

TL;DR

This paper presents a parallel, partitioned fluid-solid coupling model combining LAMMPS and Palabos via the immersed boundary method, validated through biological and engineering flow simulations, demonstrating good scalability and application potential.

Contribution

It introduces a novel partitioned coupling approach using LAMMPS and Palabos for fluid-solid interactions based on IBM, capable of simulating rigid and deformable solids in flow.

Findings

Linear scaling from 512 to 8192 cores for coupling simulations

Validated with blood flow and particle tests

Effective for simulating flexible drug carriers in blood flow

Abstract

The study of viscous fluid flow coupled with rigid or deformable solids has many applications in biological and engineering problems, e.g., blood cell transport, drug delivery, and particulate flow. We developed a partitioned approach to solve this coupled Multiphysics problem. The fluid motion was solved by Palabos (Parallel Lattice Boltzmann Solver), while the solid displacement and deformation was simulated by LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator). The coupling was achieved through the immersed boundary method (IBM). The code modeled both rigid and deformable solids exposed to flow. The code was validated with the Jeffery orbits of an ellipsoid particle in shear flow, red blood cell stretching test, and effective blood viscosity flowing in tubes. It demonstrated essentially linear scaling from 512 to 8192 cores for both strong and weak scaling cases. The computing time for the coupling increased with the solid fraction. An example of the fluid-solid coupling was given for flexible filaments (drug carriers) transport in a flowing blood cell suspensions, highlighting the advantages and capabilities of the developed code.

Tables

Table 1. Table 1 : The percentage of CPU time for the coupling functions.

ID	Functions	CPU time (%)
		p1	p2	p4	p8	p16
$t_1$	interpolateVelocity3D	10.12	10.86	10.07	8.72	7.12
$t_2$	spreadForce3D	28.25	27.27	26.45	23.75	18.23
$t_1+t_2$		38.37	38.12	36.51	32.46	25.35

Table 2. Table 2 : Simulation parameters used for the particle transport in blood flow. The shear modulus for the cell μ 0 = 6 ​ μ ​ N / m subscript 𝜇 0 6 𝜇 𝑁 𝑚 \mu_{0}=6\mu N/m . The stretching and bending parameters for filaments were selected based on [ 102 , 103 ] . The filaments were assumed to be straight in equilibrium, and the equilibrium length was chosen due to the resolution of the fluid lattice. The Morse potential was used to avoid the filament-cell and cell-cell overlapping.

Red blood cells
$k_{B}T/p\left(N\right)$	$k_b\left(J\right)$	$k_g\left(N/m\right)$	$k_l\left(N/m\right)$	$k_v\left(N/m^2\right)$
$1\times {10}^{-12}$	$2.3\times {10}^{-19}$	$2.1\times {10}^{-4}$	$2.1\times {10}^{-4}$	$2.2$
Filaments
$k_{B}T/p\left(N\right)$	$k_p\left(J\right)$	$k_b\left(J\right)$	$l_0\left(\mu m\right)$	${\theta }_0{(}^{\circ })$
$1.36\times {10}^{-12}$	$1.7\times {10}^{-12}$	$[4.1,410]\times {10}^{-20}$	0.33	$180$
Morse potential
$D_0\left(J\right)$	$\alpha$	$r_0\left(\mu m\right)$	$r_c\left(\mu m\right)$
$1.2\times {10}^{-18}$	$0.5$	$0.66$	$0.66$

Table 3. Table 3 : The apparent size and dispersion rate of the filaments in the blood flow with cells. The apparent size of the filaments were smaller than the actual contour length due to the bending and coiling. The dispersion rate was in the order of 10 − 11 ​ m 2 / s superscript 10 11 superscript 𝑚 2 𝑠 10^{-11}m^{2}/s , which was about 2 orders-of-magnitude larger than thermal diffusion for microparticles.

contour size($\mu m$)	$1.66$		$7.92$
stiffness($k_b/k_{B}T$)	$10$	$1000$	$10$	$1000$
apparent size($\mu m$)	$1.57\pm 0.31$	$1.67\pm 0.01$	$5.76\pm 2.4$	$7.89\pm 1.92$
dispersion rate($\times {10}^{-11}{m}^2/s$)	$9.76$	$8.64$	$6.26$	$7.26$

Equations

$$\underset{\text{streaming}}{\underbrace{f_i(\mathbf{x}+{\vec{\mathbf{c}}}_i,t+1)-f_i(\mathbf{x},t)}}=\underset{\text{BGK collision}}{\underbrace{\frac{1}{\tau }(f_i^{eq}-f_i(\mathbf{x},t))}}+F_i,$$

$$f_i^{eq}=w_i\rho \left(1+\frac{{\vec{\mathbf{c}}}_i\cdot \vec{\mathbf{u}}}{c_s^2}+\frac{({\vec{\mathbf{c}}}_i\cdot \vec{\mathbf{u}}{)}^2}{2c_s^4}-\frac{{\vec{\mathbf{u}}}^2}{2c_s^2}\right),$$

$$\rho =\sum f_i,\rho \vec{\mathbf{u}}=\sum f_i{\vec{\mathbf{c}}}_i+\frac{1}{2}\vec{\mathbf{g}},$$

$$F_i=(1-\frac{1}{2\tau })w_i\left(\frac{{\vec{\mathbf{c}}}_i-\vec{\mathbf{u}}}{c_s^2}+\frac{{\vec{\mathbf{c}}}_i\cdot \vec{\mathbf{u}}}{c_s^4}{\vec{\mathbf{c}}}_i\right)\cdot \vec{\mathbf{g}}.$$

$$\nu =c_s^2(\tau -0.5)=\frac{\tau -0.5}{3}.$$

$$U({\mathbf{X}}_i)=U_{stretch}+U_{bending}+U_{area}+U_{volume},$$

$$U_{stretch}=\sum _{j\in \mathrm{1...}{N}_s}\left[\frac{k_{B}T{l}_m}{4p}\frac{3x_j^2-2x_j^3}{1-x_j}+\frac{k_p}{l_j}\right],$$

$$U_{bending}=\sum _{j=\mathrm{1...}{N}_s}{k}_b(1-cos({\theta }_j-{\theta }_0)),$$

$$U_{area}=\frac{k_g(A-A_0{)}^2}{2A_0}+\sum _{j=\mathrm{1...}{N}_t}\frac{k_l(A_j-A_{j0}{)}^2}{2A_{j0}},$$

$$U_{volume}=\frac{k_v(V-V_0{)}^2}{2V_0},$$

$${\mu }_0=\frac{\sqrt{3}{k}_{B}T}{4p{l}_m{x}_0}\left[\frac{x_0}{2(1-x_0{)}^3}-\frac{1}{4(1-x_0{)}^2}+\frac{1}{4}\right]+\frac{3\sqrt{3}{k}_p}{4l_0^3}.$$

$$U_{bending}^p=k_b^p({\theta }_j-{\theta }_0{)}^2,$$

$$U_{morse}=D_0[e^{-2\alpha (r-r_0)}-2e^{-\alpha (r-r_0)}],r<r_c.$$

$$\mathbf{U}(\mathbf{X},t)=\int \mathbf{u}(\mathbf{x},t)\delta (\mathbf{x}-\mathbf{X})d\mathbf{x},$$

$$\mathbf{f}(\mathbf{x},t)=\int \mathbf{F}(\mathbf{X},t)\delta (\mathbf{x}-\mathbf{X})d\mathbf{X},$$

$$\varphi (r)=\{\begin{array}{ll}& 0,\text{otherwise}\\ & \frac{1}{4}\left(1+cos\left(\frac{\pi r}{2}\right)\right),|r|\le 2,\end{array}$$

$$f^{FSI}(\mathbf{X},t)={\rho }_f(\mathbf{u}(\mathbf{X},t)-\mathbf{U}(\mathbf{X},t))/\delta t,$$

$${\mathbf{F}}_{ij}^{lub}=-\frac{3\pi \nu r}{s}{\hat{\mathbf{x}}}_{ij}{\hat{\mathbf{x}}}_{ij}\cdot ({\mathbf{u}}_i-{\mathbf{u}}_j)$$

$$tan\theta =\frac{b}{a}tan\frac{ab\gamma t}{a^2+b^2},$$

Full text

††journal: Journal of Computational Science

1 Introduction

Fluid flows containing solid particles are common in engineering and medicine, e.g., suspension flows [1], sedimentation[2], cell transport in blood flow [3, 4, 5, 6, 7], and platelet deposition on blood vessel walls [8, 9, 10]. The dynamic behavior in such phenomena is complex due to interactions between individual particles as well as interactions between the particles and the surrounding fluid and bounding walls. Moreover, the presence of highly deformable particles, such as blood cells, vesicles and polymers, make it particularly challenging to accurately describe the dynamics in such systems. Understanding the interactions between the particulate components and the fluid is essential for optimized particulate design or detailed particulate flow behavior prediction, e.g., enhance particle mixing[11], fluid coking[12], drug carrier design[13, 4, 14], cell separation[15], and blood clotting[10, 16, 17].

The present paper is focused on blood flows containing a high concentration of deformable cells that are similar in size to the vessel diameter, thus requiring explicit consideration of the particle mechanics. Readers interested in the modeling of other classes of particulate flows, e.g., those relevant to industrial applications, are referred to other studies, e.g., Ref. [18, 19, 20]. Numerous methods have been developed to model blood flows containing cells such as red blood cells and platelets. The boundary element method is an example of a very efficient technique for these types of flows[21, 22, 23] as it formulates boundary value problems as boundary integral equations. Thus, it only requires discretization of the surface rather than the volume. However, the boundary element method requires explicit knowledge of a fundamental or analytical solution of the differential equations, e.g., linear partial differential equations. Thus, it is limited to Stokes flow conditions. The Arbitrary-Lagrangian-Eulerian method (ALE) is another approach that has been widely used to model fluid-solid interactions[24, 25]. In ALE methods the fluid mesh boundary conforms to the solid boundaries on the interface, i.e., the nodes at the fluid-solid interfaces are shared. However, the ALE technique is computationally very expensive because repeated mesh generation is necessary for flows where particles experience large deformations (e.g., for red blood cells).

By contrast, non-conforming mesh methods eliminate the mesh regeneration step. In these methods, an Eulerian mesh is used for the fluid and a Lagrangian mesh is used for the solid. These two meshes are independent, i.e., they do not share nodes across the interfaces. In this case, boundary conditions are applied through imposed constraints on the interfaces. The Immersed Boundary Method (IBM) is a popular example of a non-conforming mesh technique[26, 27]. Here, the velocity and force boundary conditions at the interfaces are imposed though interpolation functions that transfer velocities and forces across the domains. Consequently, the IBM approach may be regarded as an interface between two essentially separate simulators–one for evolving the particulate phase and the other for evolving the fluid. As we will show in this paper, the decoupled nature of the two components leads to a high degree of versatility in the context of software development. More generally, the present work represents an example of an emerging paradigm in multi-physics/multiscale modeling in which multiple existing (and independent) software packages are connected by a relatively simple interface to generate new functionality. Examples in the literature of such approaches include software packages for general fluid-solid coupling[28], sedimentation[29], atomic-continuum coupling[30], and fluid flow coupled with the discrete element method[19].

Another important aspect related to the implementation of methods for solving particulate flows is the portability to high performance computing (HPC) platforms so that large system sizes and long simulation times relevant to the phenomena of interest are accessible. An example of a reported HPC-enabled blood flow simulation is the work of [31], in which blood flow simulations in patient-specific coronary arteries at spatial resolutions ranging from the centimeter scale down to 10 $\mu m$ were performed on 294,912 cores. In Ref.[32], 2,500 deformable red blood cells in suspension under flow were simulated on IBM Blue Gene/P supercomputers. Additional examples include large-scale simulations of blood flow in the heart [33, 34, 35], cell separation in microfluidic flows[36, 15], and blood flow in the brain[37]. However, there are only a few parallel open sourced fluid solid coupling codes, such as the immersed boundary (IB) method with support for adaptive mesh refinement (AMR) IBAMR[38], the vascular flow simulation tool SimVascular[39, 40], the CFDEM project using computational fluid dynamics and discrete element methods[41, 19]. A general open source fluid solid coupling tool with versatile functionalities (e.g., deformable solids) and independent fluid and solid simulator is still missing.

In order to overcome challenges related to programming, model standardization, dissemination and sharing, and efficient implementation on parallel computing architectures, here we introduce a simple but effective implementation of the immersed boundary method for simulating deformable particulate flows using popular open source software. This project was also inspired by the previous work from LBDEM[41, 42]. The fluid solver is based on the lattice Boltzmann method(LBM), chosen because of its efficient parallelization across multiple processor environments[43, 32]. LBM is a versatile fluid flow solver engine and has been used to model diverse situations, such as flows in porous media with complex geometries[44, 45, 46] and multiphase flows[47, 48]. The solid particles in the fluids, particularly, the deformable blood cells are modeled as a coarse grained cell membrane model using particle based solid solver. Many LBM based fluid solid coupling work can be found in literature[49, 50, 15, 51, 52], some of them haven’t demonstrated the large scale parallel performance[51, 52]. Further, none of them are open sourced yet. In the present paper we employ a general LBM fluid solver Palabos[53, 54] and couple it with the LAMMPS software package [55] for describing particle dynamics. Even LAMMPS provides a LBM fluid solver through USER-LB extension[56, 57], the functionality of the fluid solver is limited, e.g., it can only model fluid with simple geometries and apply boundary velocities on the z direction. Thus, it is still necessary to build a large-scale parallel simulation tool for fluid-solid coupling based on widely used, general purpose open source codes. The implementation of a stable, efficient large-scale parallel solver for either fluid flow or particle dynamics is nontrivial, requiring training in scientific computing and software engineering. Both Palabos and LAMMPS are efficiently parallelized and clearly demonstrated by ample documentation and actively supported from large on-line communities:LAMMPS has an active mailing list that has hundreds questions and answers posted daily. These codes are also extensively validated by numerous examples and publications. Moreover, user-developers may easily extend the functionality of either package by implementing additional features, e.g., interaction potentials, integrators, fluid models, etc.

To the best of our knowledge, this is the first time Palabos has been coupled with LAMMPS in an immersed boundary method framework. The remainder of the paper is structured as follows. First, short introductions are provided in Section 2 describing the fluid solver (2.1), the solid solver (2.2), the immersed boundary method (2.3), and the spatial decomposition for the coupling (2.4). Next in Section 3.1, a validation test and convergence study are presented for a single ellipsoid in a shear flow. The parallel performance of the IBM solver is studied in Section 3.4. Finally, an example of flexible filament transport in blood cell suspensions is described in Section 3.5, highlighting the advantages and capabilities of the new solver. Conclusions and discussions are provided in Section 4.

2 Methods

2.1 Lattice Boltzmann fluid solver: Palabos

Palabos is an open source computational fluid dynamics solver based on the lattice Boltzmann method (LBM). It is designed in C++ with parallel features using Message Passing Interface (MPI). It has been employed widely in both academic and industrial settings. The LBM has been used extensively in blood flow modeling [58, 59, 60, 10, 61, 46]. Reviews of the underlying theory for the LBM can be found in the literature [62, 63, 64, 65]. LBM is usually considered as a second-order accurate method in space and time [66]. The fundamental quantity underpinning the LBM is the density distribution function $f_{i}(\boldsymbol{x},t)$ in phase space $(\boldsymbol{x},\vec{\boldsymbol{c}}_{i})$, where $t$ denotes the time and $\vec{\boldsymbol{c}}_{i}$ denotes the lattice velocity. The evolution of the density distribution function involves streaming and collision processes,

[TABLE]

where the simplest Bhatnagar$-$Gross$-$Krook (BGK) scheme [67] is used for the collision term and $F_{i}$ is the body force term that will be used to represent immersed cell boundaries[68]. The equilibrium distribution$f_{i}^{eq}$ is given by

[TABLE]

where $w_{i}$ are the weight coefficients and $c_{s}$ is the speed of sound. The density $\rho$ and velocity $\vec{\boldsymbol{u}}$ may be calculated as

[TABLE]

where $\vec{\boldsymbol{g}}$ is the external force density vector that is related to $F_{i}$ as

[TABLE]

see Ref.[68] for more details on force terms. The fluid viscosity $\nu$ is related to the relaxation parameter $\tau$ as

[TABLE]

In all LBM simulations reported in this paper, the fluid domain is discretized using a uniform D3Q19 lattice; see ref.[64]. The fluid density $\rho$, $\vec{\boldsymbol{u}}$ and density distribution function $f_{i}$ are initialized at the equilibrium distribution calculated from Eqn.4 based on the initial fluid velocity. During each time step, streaming and collision steps are performed on $f_{i}$ according to Eqn.1. Specifically, the $f_{i}$ is translated in the direction of the discretized velocity vector $\vec{\boldsymbol{c}}_{i}$ during the streaming step; then, $f_{i}$ is updated based on the equilibrium distribution, $f_{i}^{eq}$, the relaxation parameter, $\tau$, and the force density, $F_{i}$, which is passed to the LBM from the immersed solid objects, e.g., blood cells.

2.2 Particle based solid solver: LAMMPS for deformable cells and particles

LAMMPS was originally designed as a molecular dynamics simulation tool[55]. In molecular dynamics, a potential function is defined to model the interactions between atoms. The force on each atom is calculated as the derivative of the potential with respect to the atomic coordinates and the atomic system motion is updated based on numerical integrations of Newton’s 2nd law of motion. The LAMMPS package has now been extended to include a variety of additional dynamical engines, including peridynamics[69, 70], smooth particle hydrodynamics[71, 72], dissipative particle dynamics[73, 74], and stochastic rotation dynamics[75]. LIGGGHTS(LAMMPS improved for general granular and granular heat transfer simulations) is also an extension of LAMMPS for discrete element method particle simulation[19]. Many predefined potentials, functions, and ODE integrators in LAMMPS make it extremely powerful for modeling atomic, soft matter[76], and biological systems[77, 8].

A coarse-grained membrane model consisting of many interacting particles[78, 49, 79, 80, 46] is used to simulate red blood cells, as shown in Fig.1(a). The membrane model can bear stretching and bending. Constraints to maintain the constant membrane surface area and enclosed cell volume are imposed through harmonic potentials. The viscosity ratio of cytoplasm over blood plasma is about 5. Here we treat the blood cell internal and external viscosity to be the same for saving computational cost. Readers interested in different viscosity models are referred to Ref.[5, 81, 82]. The potential function for a red blood cell (RBC) used in the current work is given by

[TABLE]

where the stretching energy $U_{stretch}$ is used to represent the cytoskeleton’s resistance to deformation. The bending energy $U_{bending}$ represents the rigidity of the membrane bilayer imparted to the cytoskeleton. The last two terms are the constraints for maintaining constant membrane surface area and cell volume. The stretching potential is given by:

[TABLE]

where $l_{m}$ is the maximum bond length, the $j$th bond length ratio is $x_{j}=l_{j}/l_{m}$. $l_{m}$ was set to be 2 times the equilibrium bond length. $N_{s}$ is the number of springs, $p$ is the persistence length, $k_{B}$ is the Boltzmann constant, $T$ is the temperature, and $k_{p}$ is the repulsive potential constant. Once $p$ is specified, $k_{p}$ may be found using the value of $x_{j}=0.5$ at the equilibrium where the net force is zero.

The bending energy is defined as

[TABLE]

where $k_{b}$ is the bending constant and $\theta_{j}$ is the instantaneous angle formed by the two outward surface norms of two adjacent triangular meshes that share the same edge $j$. $\theta_{0}$ is the corresponding equilibrium, or spontaneous, angle.

Constraints for constant membrane surface area and cell volume are imposed though area/volume dependent harmonic potentials,

[TABLE]

[TABLE]

where $k_{g},k_{l}$ are the global and local area constraint constants, $N_{t}$ is the number of triangular surfaces, $A,A_{0}$ are the instantaneous and spontaneous total surface area of the cell membrane, $A_{j},A_{j0}$ are the instantaneous and spontaneous surface area for the $j$th triangle surface, $k_{v}$ is the volume constraint constant, and $V,V_{0}$ are the instantaneous and the equilibrium cell volume.

The parameters used in the coarse-grained membrane model can be related to membrane properties used in continuum model[79, 77], e.g., the shear modulus, $\mu_{0}$, which is given by

[TABLE]

where $x_{0}=l_{0}/l_{m}$ and $l_{0}$ is the bond length at equilibrium. Eqn.11 holds assuming the membrane surface area is preserved during simulation, thus the contribution from $U_{area}$ can be safely ignored. Interested readers can refer to Ref.[79, 77] for details.

For polymer particles, as shown in Fig.1(b), the stretching energy is the same as Eqn.(7), while the bending energy is a harmonic function of the angle deviation,

[TABLE]

where the superscript $p$ refers to polymer and the other variables are defined as in Eqn.(8).

To avoid the overlapping of the particles from different solid objects, e.g, cells, polymers, etc. a Morse potential was used for inter-particle interaction,

[TABLE]

where $D_{0}$ is the energy scale, $\alpha$ controls the width of the potential, $r$ is the distance between particles from different solid objects, $r_{0}$ is the equilibrium distance, $r_{c}$ is the cutoff distance.

All the parameters used in the simulation are listed in Table 2.

2.3 Fluid-solid coupling: the immersed boundary method

The immersed boundary method (IBM) was used to model the coupling between fluid and solid. The combination of LBM and IBM was first used to model fluid-particle interaction problems in Ref. [83]. Details of the immersed boundary method formulation may be found in Refs. [26, 27, 83]. Briefly, the solid velocity at each particle position is obtained through velocity interpolation from local fluid nodes, while fluid forces are obtained by spreading the local solid forces. Specifically, for an immersed solid with coordinates $\boldsymbol{X}$, the velocity $\boldsymbol{U}(\boldsymbol{X},t)$ is interpolated from the local fluid velocity $\boldsymbol{u}(\boldsymbol{x},t)$, while the solid force, $\boldsymbol{F}(\boldsymbol{X},t)$, calculated from Eqn.6 is spread out into the local fluid grid points as a force density $\boldsymbol{f}(\boldsymbol{x},t)$:

[TABLE]

[TABLE]

where $\delta(\boldsymbol{x})$ is the delta function. In a typical numerical implementation in three dimensions, $\delta(\boldsymbol{x})$ is constructed as the product of one-dimensional functions, i.e., $\delta(\boldsymbol{x})=\phi(x)\phi(y)\phi(z)$ with $\phi(r)$ defined as

[TABLE]

where $r$ is the distance between solid particles and fluid nodes. Different choices of the interpolation function influence the coupling accuracy, the influence range of the solid particles on the fluid, and the computational cost [84, 85]. Schematic representations of the velocity interpolation and force spreading process used in the present implementation are shown in Fig.2.

The coupling strategy described above becomes numerically unstable for rigid objects. In such situations, a different fluid-solid interaction (FSI) approach must be applied [86, 87]. In the present IBM formulation for rigid particles, the FSI force on each particle is used to assemble the total force and torque on the whole rigid body. The translation and rotation of the rigid object are then updated based on Newton’s 2nd law of motion. The FSI force can be expressed as:

[TABLE]

where $\rho_{f}$ is the fluid density, $\boldsymbol{u}(\boldsymbol{X},t)$ is the fluid velocity at the solid boundary node $\boldsymbol{X}$, which has to be obtained through interpolation by Eqn.(14). $\boldsymbol{U}(\boldsymbol{X},t)$ is the solid nodal velocity. The idea behind Eqn.(17) is that the local fluid particles with incoming velocity $\boldsymbol{u}(\boldsymbol{X},t)$ will collide with the solid boundary with outgoing velocity $\boldsymbol{U}(\boldsymbol{X},t)$. The FSI force is the change of the momentum divided by the collision time $\delta t$. Other approaches have been proposed to model rigid objects in IBM, e.g., a virtual boundary formulation was used in Ref.[88].

Another challenge in modeling fluid mediated solid transport is to correctly capture the fluid flow in the lubrication layer between two approaching solids or between one solid and a wall. For example, when the gap between the cell membrane and the wall is very small, e.g., the gap is smaller than a lattice space, the LBM fluid solver can not resolve the fluid flow in the thin layer. One approach is to use a finer mesh for the whole fluid domain or refine the mesh near the boundary layers[89, 90, 91]. This approach would increase the LBM simulation cost as more lattice space is used. It also requires some efforts to handle the density distribution passage over the interface if a multi-grid is used in LBM. Another approach is to introduce an lubrication force to repel the cell membrane so that there are enough fluid within the gap. The lubrication repelling force was introduced to the lattice Boltzmann method in Ref[92, 93]. Following their work, the lubrication force derived from the lubrication theory between two identical spheres is

[TABLE]

where $r$ is the spherical radius, $\nu$ is the fluid dynamic viscosity, $s$ is dimensionless gap $s=R/r-2$ where $R$ is the central distance between two spheres. $\boldsymbol{x}_{ij}$ is the position vector difference between sphere $i$ and $j$, defined as $\boldsymbol{x}_{ij}=\boldsymbol{x}_{i}-\boldsymbol{x}_{j}$, $\hat{\boldsymbol{x}}_{ij}$ is the unit vector. $\boldsymbol{u}$ is the spherical velocity. Eqn.18 can also be extended to the case where a sphere approach a stationary wall by setting $\boldsymbol{u}_{j}=0$. Interested readers can find more details there. In the present paper, we resolve the lubrication effect using smaller fluid mesh size, e.g., the fluid mesh size was set to 100 nm when a blood cell squeezing though a 5 $\mu$m tube, as shown in section 3.3.

2.4 Spatial decomposition for fluid-solid coupling

Spatial decomposition is adopted by both Palabos and LAMMPS for parallel computing. In our fluid-solid coupling approach, the same spatial decomposition was applied to both the fluid and the solid so that the coupling can be effectively handled by the same compute core for a given region. Consistent spatial decomposition for the two domains ensures that individual processors have access to both the fluid grid points and solid particles within the same sub-domain. An illustration of the partitioning process is shown in Fig.3, where the whole fluid-solid system is partitioned onto 8 compute cores. Ghost layers near the boundaries of each sub-domain are used for communication between neighboring cores. This approach may not be optimal for a system where solid particles are highly heterogeneously distributed across the entire domain. However, for most cases of interest the distribution of solid particles is quite homogeneous, and in any case, the majority of the computation is dedicated to the fluid solver. The issue of load optimization for heterogeneous systems is deferred to future work.

The FSI coupling algorithm consists almost entirely of two routines to carry out velocity interpolation and force spreading functions. As mentioned above the fluid solver usually represents the bulk of the computational demand because the solid fraction is usually quite small. Based on this intrisic asymmetry, we employed Palabos as the driving code, while LAMMPS was called as an external library. To fully access all the members(e.g., particle positions, velocities, forces, etc.) from LAMMPS, a pointer to LAMMPS was used and passed to the two interpolation functions. An outline of the functions implemented in our IBM algorithm is shown in List1.

where $lattice$ is the structure used in Palabos to store the population distribution functions for the LBM, while $wrapper$ is a pointer to an instance of LAMMPS. The advantage of the IBM in terms of software development is that only two functions are needed to couple a fluid solver to a solid solver. Our IBM implementation therefore only requires a few hundred lines of code beyond what is implicitly contained in the Palabos and LAMMPS packages. Detailed implementations can be found at https://github.com/TJFord/palabos-lammps.

3 Results

3.1 Validation I: Ellipsoid in shear flow

To validate our IBM implementation, we considered the trajectory of a rigid ellipsoid in a shear flow in Stokes flow regime, commonly referred to as the Jeffery’s orbit[94], where the rotation angle of the ellipsoid, $\theta$, satisfies

[TABLE]

where $\theta$ is the angle formed by the major axis of the ellipsoid and the shear flow direction, $a,b$ are the major and minor semi-axis lengths, $\gamma$ is the shear rate, and $t$ is the time. A 2D illustration of the ellipsoid is shown in Fig.4–the third semi-axis length is assumed to be the same as $b$.

An ellipsoid with semi-axis lengths $a=6,b=c=4.5$ was placed at the center of a channel with height $H=60$. The top and bottom velocities were set to $U=0.01$ and a viscosity of $\nu=1/6$ was prescribed. These parameters give a shear rate $\gamma=0.00033$ and a Reynolds number of $Re=aU/\nu=0.36$ to approximate Stokes flow. The IBM numerical result for the orientation angle $\theta$ and the analytical solution given by Eqn.(19) are plotted in Fig.5. The agreement is good, validating the FSI coupling. Shown in Fig.6 is the relative error, defined as $\dfrac{|\theta_{sim}-\theta_{the}|}{\theta_{the}}$ where $\theta_{sim}$ is the simulation data and $\theta_{the}$ is the theoretical data from Eqn.19. The relative error is seen to oscillate with the period of the ellipsoid and is, on average, constant at about 2.62%.

As expected, the relative error $\epsilon$ is a power function of the fluid grid resolution $N$–Fig. 7 shows that the relative error scales as the grid spacing following $\epsilon\propto N^{-0.96}$.

3.2 Validation II: Red blood cell stretching test

The RBC model was also validated and compared with optical tweezer stretching results under different loadings. The cell diameter without deformation was 8 $\mu m$, bending stiffness $k_{b}=2.3\times 10^{-19}J$, $k_{g}=k_{l}=2.1\times 10^{-4}N/m$, $k_{v}=2.2N/m^{2}$. The stretching constant $k_{B}T/p$ depended on the number of discretized membrane particles, e.g., $k_{B}T/p=1\times 10^{-12}N$ for 1320 surface particles, $k_{B}T/p=7.14\times 10^{-13}N$ for 2562 surface particles, $k_{B}T/p=3.57\times 10^{-13}N$ for 10242 surface particles. However, the shear modulus $\mu_{0}$ was kept constant with $\mu_{0}=6$ $\mu N/m$ for all the meshes. A pair of forces was applied to two membrane patches on each side, occupying about 5% of the total number of particles. The applied force was shared uniformly among these particles. Velocity verlet algorithm was used to integrate the Newton’s equation. Viscous damping was used to stabilize the system. The deformation of RBCs under stretching (Fig.8(a)), particularly, the axial and transversal diameter (Fig.8(b)), were recorded and shown in Fig.8. In Fig.8(b), the top curve was the axial diameter, and the bottom curve was the transversal diameter. Both curves agreed well with experimental optical tweezer stretching data[79] for three different mesh sizes, validating the coarse grained molecular dynamics cell model.

3.3 Validation III: Effective blood viscosity

Blood viscosity changes as it flows through different tube diameters and hematocrits. The viscosity decreases as the tube diameter in the range of 10$\sim$ 300 $\mu m$ due to the Fahraus-Lindqvist effect[95]. The effective viscosity of blood in tube flow was studied and compared with experimental results[96]. The hematocrit defined as the actual volume ratio between cells and the fluid was 30%. The tube diameter was ranging from 5 $\mu m$ to 30 $\mu m$. The fluid mesh size and the number of cell membrane nodes were 100 nm and 10242 for 5 and 6 $\mu m$ tube, 200 nm and 2562 for 8, 10, 15 $\mu m$, 333 nm and 2562 for 30 $\mu m$. The fluid mesh size was selected in a way such that it can resolve the fluid flow between cells and the tube wall. The mechanical properties were exact the same, as shown in section 3.2. A pressure gradient was applied to drive the flow. The Reynold’s number defined with cell diameter and average flow speed was around 0.017 $\ll$ 1. Periodic boundary conditions were applied at the inlet and outlet. The simulation ran for enough time to reach quasi-steady state, e.g., the volume rate was steady. The effective viscosity was normalized with respect to the plasma viscosity, and shown in Fig.9. The simulation data clearly showed a decreasing viscosity as the tube diameter decreased from 30 to 8 $\mu m$. The viscosity increased as the tube diameter decreased further, which was resulted from the physical membrane contact with the wall. The viscosity curve agreed well with the empirical fitting from Ref.[96], which validated that the model can capture the collective blood cell behavior under flow. The effective viscosity simulated by boundary integral method[97] was also provided for comparison. Both methods can predict the blood viscosity reasonably well.

3.4 Scalability: Parallel performance

The parallel performance of the code was examined using a benchmark simulation of red blood cells flowing in a rectangular box. The fluid domain was discretized into a $360\times 360\times 720$ grid in x,y,z directions with one inlet and one outlet and $4$ nonslip side walls. The fluid flow was driven by a body force of $1.34\times 10^{-7}$ with periodic inlet and outlet boundary conditions in the z direction. The solid phase consisted of $15,680$ red blood cells, discretized into $20,697,600$ particles connected by $61,998,720$ bonds, $41,332,480$ angles, and $61,998,720$ dihedrals in total. The hematocrit was about 40.8%. The cells were initially uniformly distributed in the flow domain. The simulations were executed on the IBM Blue Gene/Q system from Argonne National Laboratory. A ghost layer with thickness of three fluid lattice spacings was employed for inter-processor communication within the fluid solver. The ghost cutoff distance for LAMMPS was set to 1.5 times the LBM lattice spacing. The CPU time for 100 time steps was recorded. Both strong scaling and weak scaling were considered. For weak scaling cases, the fluid domain size increased from $360\times 360\times 45$ for 512 cores to $360\times 360\times 720$ for 8192 cores. The parallelization performance was analyzed using two parameters: speed up ($SU$) and scaling efficiency ($\gamma$). The speed up and scaling efficiency were calculated as $SU(N)=\dfrac{t_{0}^{s}}{t_{N}^{s}},\gamma(N)=\dfrac{t_{0}^{s}N_{0}}{t_{N}^{s}N}$ for strong scaling cases, and $SU(N)=\dfrac{t_{0}^{w}N_{0}}{t_{N}^{w}N},\gamma(N)=\dfrac{t_{0}^{w}}{t_{N}^{w}}$ for weak scaling cases, where $t_{0}$ was the reference simulation time for $N_{0}$ cores, $t_{N}$ was the simulation time for $N$ cores. The subscript $s,w$ denoted strong and weak scaling cases, respectively. Here $N_{0}$ was taken as 512. The speed up of the simulation time as well as scaling efficiency were shown in Fig.10. It showed that, as the number of cores increased from 512 to 8192, the speed up increased from 1 to 10.5 for strong scaling case, from 1 to 13.2 for weak scaling cases, see Fig.10(a). Similarly, the scaling efficiency dropped from 1 to 0.656 for strong scaling case, and from 1 to 0.824 for weak scaling case, see Fig.10(b). The dashed lines showed the ideal cases with no overhead and communication cost. It demonstrated that the combination of Palabos and LAMMPS using the immersed boundary method can achieve high performance among thousands of processors.

The scaling efficiency was found to depend strongly on the computer hardware. For example, less parallel efficiency was obtained on another computer with Intel Xeon E5630 2.53 GHz processors (16 cores in total) compared with 2 AMD Opteron 8-core 6128 processors with a clock speed of 2.0 GHz. While we did not study architecture dependence in more detail here, these results suggest that further work is needed to assess architecture-dependent performance issues. Finally, we note that the ghost layer thickness also influenced the performance of the IBM code. In the present study we used 3 layers of lattice points, while other studies have employed linear interpolation kernels that only require a single layer of lattice points for the communication layer[98].

The extra computational time induced by the coupling can be analyzed by GNU profiler tools. We considered the blood flow in a microfluidic channel with size of $60\times 60\times 120\mu m$. $1,470$ cells were placed in the channel, resulting a hematocrit of 31% in the flow. In total there were $1,940,400$ cell membrane particles and $11,664,000$ fluid lattices, as shown in Fig.11. On average, the extra computational time required for the information exchange between the fluid and solid solvers was about 34%, e.g., time spent on the two coupling functions: interpolateVelocity3D and spreadForce3D. The computational time for the two coupling functions expressed in terms of the percentage over the total CPU time is shown in Table 1. The function interpolateVelocity3D cost about 9% on average, while function spreadForce3D cost about 25%. The percentage for both functions decreased as more cores were used. This was due to the increased communication cost among many processors. As expected, the coupling time would increase as the hematocrit level, as the interpolation and spreading was done for each cell membrane particle.

3.5 Case Study: Transport of flexible filaments in flowing red blood cell suspensions

Next, the validated IBM code was applied to a problem relevant to drug carrier delivery in microcirculation flows to demonstrate its capabilities in a more complex setting. This is a multi-physics modeling problem as it involves fluid flow, large cell deformations, and polymer transport[6, 4, 99, 4, 7]. Of particular interest is the observation that particles with different sizes, shapes, or flexibility can exhibit distinct transport properties in the blood stream. For example, long flexible filaments persist in the circulation up to one week after intravenous injection in rodent models. This is about ten times longer than for spherical counterparts[100]. The flexibility and filament length may influence transport and contribute to the anomalously long circulation time.

We simulated filament transport in tube flow of red blood cell suspensions and analyzed filament migration, apparent size, spatial distribution, and dispersion rate. The fluid domain was represented by a cylindrical tube with diameter of $30\mu m$ and length $50\mu m$, as shown in Fig.12. No-slip conditions were applied at the cylinder wall while periodic boundary conditions were applied at the inlet and outlet for both fluid and cells. The relaxation time was set to $\tau=1$ giving a dimensionless kinetic viscosity of $1/6$. The LBM lattice spacing, $dx$, was fixed at $0.33\mu m$, the time step was set to $dt=1.82\times 10^{-8}s$. One hundred red blood cells were placed in the domain corresponding to a cell concentration of about 27%, which is within the physiological range[101]. Each cell was described by 1,320 surface nodes, 3,954 stretching bonds, 2,636 area segments, and 3,954 bending angles. A force density was applied to drive the flow such that the shear rate at the wall was $1000s^{-1}$ without the red blood cells. The actual wall shear rate was less as the blood cell component reduces the flow rate. A total of 156 filaments were then introduced at one end of the cylinder as shown in Fig.12. Each filament was a polymer chain modeled as beads connected by stretching and bending springs. Different lengths ($2\mu m$ and $8\mu m$) and different bending stiffnesses ($4.1\times 10^{-20}$ J, or about $10k_{B}T$, and $4.1\times 10^{-18}$ J, corresponding to about $1000k_{B}T$) were used to study the size and stiffness effect on filament transport in red blood cell suspensions.

The parameters used in the simulations are listed in Table.2. The corresponding shear modulus for red blood cell membranes was $\mu_{0}=6\mu N/m$ based on Eqn.11.

Simulation snapshots at $8\times 10^{6}$ time steps are shown in Fig.13 for cases with different filament size and flexibility. Different configurations are observed for filaments with different stiffnesses. The stiffer filaments were generally straight when the length was $2\mu m$, and exhibited a small amount of curvature for $8\mu m$, see Fig.13(b) and 13(d). By contrast, the more flexible filaments exhibited highly bent or coiled shapes, particularly for the $8\mu m$ filaments, see Fig.13(a) and 13(c). The apparent filament size, defined as the maximum extent for the filaments, is shown in Table.3. The apparent size for stiffer filaments was very close to the actual contour length, while the apparent size for flexible filaments was smaller than the actual contour length, due to the bending and coiling of the filaments.

The average radial position $\left<r\right>$ normalized to the vessel radius $R$ and the mean square displacement (MSD) for the filaments are plotted in Fig.14(a) and Fig.14(b). Filament migration towards the vessel wall was observed for all filament types, as indicated by the increasing values of $\left<r\right>/R$. The averaged mean radial positions for flexible filaments ($10k_{B}T$) was stabilized around $0.7R$ for both sizes. For stiffer filaments ($1000k_{B}T$), the $\left<r\right>/R$ kept increasing within the simulation time, showing better margination properties than flexible ones. This agrees well with previous findings for deformable particles[13]. Less margination for flexible ones lead to low binding to the vessel wall and high concentration in the blood stream, which can also explain the long circulation time for long flexible filaments in rodents[100]. The dispersion rate, defined as the slope of the MSD curve, is shown in Table 3. The dispersion rates are all of the order of $10^{-11}m^{2}/s$, which is about 2 orders-of-magnitude larger than the background thermal diffusion expected for particles of this size. This finding is also consistent with other observations for small particles in blood suspensions, such as platelets[23, 104, 105], particles[4, 13] and experimental measurements for microparticles[106].

The fraction of filaments present in the $3\mu m$ cell-free layer (CFL) near the vessel wall is shown in Fig.15. Only the data in the last 0.15 s were selected for analysis due to the transient effect in the beginning. The fraction of the 2$\mu m$ filaments in the CFL increases approximately from 28% (the baseline value assuming a uniform distribution) to $43.3\%\pm 1.7\%$ for soft filaments, and to $46.1\%\pm 3.4\%$ for stiff filaments. However, the fraction of 8$\mu m$ filaments in the CFL reached $47.5\%\pm 2.4\%$ for the flexible ones, and $51.1\%\pm 1.9\%$ for stiffer ones. These results suggest that longer and stiffer filaments in blood cell suspensions marginated quickly than the shorter and softer ones. The margination for rigid filaments also required less traveling distance in the flow direction, as the effective viscosity was increased for the mixture of rigid filaments in blood cell suspensions.

4 Conclusions

We introduced an efficient implementation of the Immersed Boundary Method by coupling a Lattice Boltzmann fluid solver with a particle-based solid solver using open source codes, namely Palabos and LAMMPS, respectively. The coupling was achieved by a simple interface that only requires a few hundred lines of code, dramatically reducing software development time, increasing robustness and facilitating transferability. The coupling was demonstrated to be computationally efficient. The extra computing time depends on the solid fraction, e.g., for a hematocrit of 30%, it required only about a 34% increase in the computational cost beyond what was required for executing fluid and solid solvers. The IBM code was also shown to scale linearly from 512 to 8192 cores for both strong and weak scaling cases, demonstrating great parallelization efficiency in massively parallel environments.

The validated IBM code was used to analyze polymer filament transport in red blood cell flows; such filaments are of interest for drug delivery applications. Simulations demonstrated how filament flexibility reduce the effective size of filament, so that flexible filaments can traverse through the blood cell suspensions easier than stiffer ones. The stiffer filaments showed better margination properties than flexible ones, consistent with previous findings[13]. These results highlight the complex physics that may underlie the long circulation time found for such particles in rodents [100].

5 Acknowledgment

This work was supported by NIH 1U01HL131053-01A1. We thank the high performance computing support from Gaea at Northern Illinois University and the Mira at the Argonne National Laboratory.