Direct numerical simulation of wave propagation in saturated random granular packings using coupled LBM-DEM

. Poroelasticity theory predicts wave velocities in a saturated porous medium through a coupling between the bulk deformation of the solid skeleton and porous ﬂuid ﬂow. The challenge emerges below the characteristic wavelengths at which hydrodynamic interactions between grains and pore ﬂuid become important. We investigate the pressure and volume fraction dependence of compressional- and shear-wave velocities in ﬂuid-saturated, random, isotropic, frictional granular packings. The lattice Boltzmann method (LBM) and discrete element method (DEM) are two-way coupled to capture the particle-pore ﬂuid interactions; an acoustic source is implemented to insert a traveling wave from the ﬂuid reservoir to the saturated medium. We extract wave velocities from the acoustic branches in the wavenumber-frequency space, for a range of conﬁning pressures and volume fractions. For random isotropic granular media the pressure-wave velocity data collapse on a single curve when scaled properly by the volume fraction.


Introduction
Understanding wave propagation in saturated granular media is vital for non-destructive soil testing, seismic inversion and oil exploration, among others. Extensive work has been done to investigate the dispersion and attenuation properties of dry granular media in both laboratory experiments and computer simulations [1,2,3]. The discrete element method (DEM) has been used to capture the microscopic kinematics associated with wave propagation in dry granular media. To simulate wave propagation in saturated granular media, however, the momentum and energy exchange between fluid, solid and/or gas phases must be taken into account in a fully coupled framework.
While conventional computational fluid dynamics methods are sufficient for modeling dilute suspensions of particles, direct numerical simulation techniques, such as the coupled lattice Boltzmann-discrete element method (LBM-DEM), are better suited for wave modeling because it captures the fluid flow and particle-fluid interactions at the (sub)pore scale. In the past decade, coupled LBM-DEM simulations have been applied to solve various geotechnical and geophysical problems, including sedimentation [3], hydraulic fracture [5], piping erosion [6], etc.
In the present work, the two-way coupled LBM-DEM is utilized to model the fluid-solid interaction at the pore scale, relevant to the propagation of pressure and shear waves in saturated poroelastic media. LBM resolves the pore fluid flow due to (oscillating) pressure gradients in the pores, while DEM captures the change of particle motion * e-mail: h.cheng@utwente.nl and interparticle forces, induced by wave perturbation. The fluid-solid coupling is handled with the momentum exchange method (MEM) [7] which computes the hydrodynamic forces on solid particles and updates local flow at the no-slip fluid-solid interfaces. To propagate elastic waves in the fluid-solid system, an oscillating fluid pressure boundary condition that emits pressure waves from the fluid is implemented.
In previous work [8], the hydro-micromechanical modeling framework was compared with the poroelasticity theory [9], for saturated, ordered, isotropic, packings of frictionless spheres. In this work, saturated, random, cubic packings of frictional spheres are considered. The ability of the hydro-micromechanical model is tested by computing the pressure-wave velocity relationship at various initial volume fractions. The space-time evolution of in-phase particle/fluid velocities are measured during wave propagation and processed in the wavenumber-frequency space. Along the acoustic branches, the wave velocities are extracted for different initial volume fractions.

Coupled LBM-DEM for wave modeling
The local hydrodynamics in the pore space of a granular medium is solved by the discrete Boltzmann equation with a Bhatnagar-Gross-Krook (BGK) collision operator. The BGK operator uses a single relaxation time for the probability distribution of local fluid velocities towards an equilibrium distribution. Interactions between contacting solid particles are modeled with DEM.
DEM represents granular materials as packings of solid particles with simplified geometries (e.g., spheres) and tiny interparticle overlaps [10]. The particles in contact interact with their neighbors via repulsive springs and viscous dashpots, resulting in momentum exchange between particles. After all forces acting on each particle, either from the neighboring particle or the surrounding fluid, are known, the kinematics of each particle is updated by Newton's second law via explicit time integration.
The Hertz-Mindlin contact laws are used to relate the interparticle forces between two touching solid particles and the overlaps in normal and tangential displacements. From the contact force network in a given microstructural configuration, the effective static stress tensor is given by the average over the particle assembly is the total number of contacts contained within the volume V of the granular assembly including the pore space, and d c is the branch vector connecting the centers of each pair of contacting particles in the force network.
The fluid-solid coupling scheme is implemented following the link-based momentum exchange method (MEM) [11]. The main idea is to project solid obstacles onto the lattice, and converting a node between "solid" and "fluid" dynamically, depending on the motion of the solid obstacles. Each solid node has at least one link with neighboring fluid nodes. Along the links, momentum is exchanged and no-slip boundary conditions are applied on the solid surface. From the exchanged momentum, the hydrodynamic force arises as a consequence of the no-slip condition applied to fluid advecting towards solid nodes. Interested readers are referred to [8] for further details.

Oscillating pressure boundary condition
An oscillating fluid pressure boundary is employed to send acoustic waves from the fluid boundary into the saturated porous material. Our approach consists in enforcing LBM boundary nodes to emit the correct numbers of fluid particles, which macroscopically leads to oscillations in the velocity component along the propagation direction. To model fluid flow in porous media, pressure boundary conditions are typically set with constant local fluid densities. The flow is then driven by an imposed pressure gradient between two opposite boundaries [12].    To propagate plane waves from the pressure boundaries, the densities at the boundary nodes are locked to periodic functions that oscillate around ρ f for finite time t n , where ρ f is the fluid density. For the pressure boundary nodes that agitate a single wavefront, the local fluid densities varies according to where ω is the input angular frequency. It is important to ensure the maximum perturbation ρ ρ f and the computational Mach number M = |u max |/c 1, so that nonlinear wave effects are avoided and LBM remains a valid approximation of the Navier-Stokes equations. Readers are referred to [8] for details.

Elastic wave propagation in fluid-saturated, random, isotropic granular packings
Five granular packings of monodisperse particles (r = 5 mm) are randomly generated and compacted to a confining pressure σ c = 10 MPa, using different initial interparticle friction coefficients (see Table 1), resulting in different initial volume fractions φ 10MPa . Similar to [8], the properties of glass particles (Young's modulus E = 70 GPa and Poisson's ratio ν = 0.2) are used and the packings are much longer in the longitudinal than the transverse direction, as shown in Fig. 1. Before the simulation of wave propagation, each packing is loaded to obtain isotropic stress states in a pressure range of σ c ∈ (10, 200) MPa. After equilibrium is achieved at each stress state (with a sufficiently long relaxation), a pressure wave is generated from one of the fluid boundaries with an input angular frequency, e.g., ω/2π = 5.2 kHz for Figs. 4c and 4d. Fig. 1 shows a snapshot of the simulation of elastic wave propagation in the saturated granular packing. The first and last layers of particles (perpendicular to the x-axis) close to the pressure boundaries are fixed in space, while the other four side boundaries are periodic.
During the simulation, the space-time evolution of the particle and fluid velocities are recorded (see Figs. 3 and 2), and further processed in the wavenumber-frequency space to obtain the acoustical branches (see Figs. 4a and 4b). . Fig. 2 shows that the wavefront increasingly widens and the energy dissipates as the P-wave travels in space and time. Note that an extremely slow and diffusive wave exists near the source, suggesting the presence of a secondary slow P-wave, which has been observed in regular granular crystals [8] and predicted by [9].
Figs. 4a and 4b show the frequency domain response obtained from the discrete Fourier transform (DFT). The acoustic branches, highlighted red, suggest a linear relationship between wavenumber k and frequency ω, as expected for the in-phase motion of the fluid and particles within the primary wave front. For each combination of confining pressure and solid volume fraction, we fit these acoustic branches with straight lines, thereby obtaining the pressure-wave velocity relationship. Furthermore, another input frequency (ω/2π = 52 kHz) is used to agitate P-and S-waves in the granular packing at a confining pressure σ c = 10 MPa, with φ 10MPa = 0.6322 (see Figs. 4c and 4d), to explore the frequency dependence of wave velocities.

Results and discussion
As a first test for LBM-DEM modeling framework, we explore the dependence of wave velocities on volume fraction. Figs. 5a-a shows the relationships between the pressure and the P-wave velocities for different initial volume fractions, at an input frequency of 5.2 kHz. As expected, wave velocities increase monotonically, as the packing fraction and the confining pressure increases. For the compressional wave, it appears that all curves are parallel to each other. By extracting the P-wave velocities at σ c = 10 MPa and relating the values with the initial volume fraction, a linear scaling function is found for the compressional wave velocities, F p (φ) = φ 10MPa + 0.217. As shown in Figs. 5a-b, all curves collapse universally on the same trend, regardless of the initial solid volume fractions. The trend suggests a potential approach to improve the theoretical prediction of P-wave velocities considering the effect of the microstructure [13] in random, isotropic, monodisperse, granular materials.
Cheng et al. [8] showed that the S-wave velocities in saturated, isotropic, frictionless granular crystals are not much affected by the hydrodynamic coupling due to the pore fluid flow. Similar observations have been reported both experimentally and theoretically. Compared to the hydrodynamic coupling, the role of the porosity or solid volume fraction is more significant in modifying the Swave velocities, as shown in Fig. 5b-a. For the shear wave velocity, we use a different normalization function, that is F s (e) = (2.17−e) 2 /(1+e) with the void ratio e = (1−φ)/φ. This empirical function is widely used in geotechnical engineering [13] to remove the dependence of the shear modulus of a dry soil sample on the void ratio. Fig. 5b-b shows the curves normalized by F s (e). It appears that all curves collapse at confining pressures lower than 100 MPa. The marked discrepancy at the high confining pressures may be attributed to the very high pressure and/or very high solid fractions for which the empirical correlation from soil mechanics fails to work.
The (normalized) pressure-wave velocity relationships are obtained with an input frequency of 5.2 kHz (see 4).
Although not shown here, we have observed that further reducing this frequency does not lead to a significant increase of the wave velocities. However, with an input frequency of 52 kHz, the wave velocities obtained from the wavenumber-frequency domain increase by 5% for the Pwave and 1% for the S-wave, as shown in Figs. 4c and 4d. The decreasing trend, albeit with a small difference, is inline with the poroelasticity theory and many experimental observations, that is the wave velocities of saturated porous media increases with an increasing input frequency of the traveling wave. Such behavior is fundamentally different from that of dry granular materials, in which a decreasing wave velocity is typically observed with an increasing input frequency.

Conclusions
A hydro-micromechanically coupled LBM-DEM framework has been used to simulate wave propagation in saturated poroelastic granular media. The hydrodynamics in the pore fluid is resolved with LBM, and the translational and rotational motion of solid particles with DEM. The novelty of the work lies in the application of the two-way coupled LBM-DEM method to study the fluidsolid interaction at local level, induced by small-amplitude waves traveling in saturated media. An oscillating pressure boundary is used to propagate planar pressure waves from the fluid phase into the saturated mixture. We consider fluid-saturated, random, isotropic granular packings of frictional spheres, extending our previous investigations on regular granular crystals [8]. Granular packings obtained at the same initial confining pressure but with different volume fractions are loaded under isotropic compression. After stabilizing each isotropic stress state, a pressure wave is agitated by the fluid nodes with a varying den-sity. The in-phase particle/fluid velocity fields and their discrete Fourier transform in the wavenumber-frequency space allow to obtain the wave velocities. The frequency dependence is briefly investigated with two input frequencies of the traveling waves.
From the coupled LBM-DEM simulations, pressurewave velocity relationships are obtained for different initial volume fractions. For random isotropic granular media the pressure-wave velocity data collapse on a single curve when scaled properly by the initial volume fraction, suggesting a potential approach to improve theoretical poroelasticity. Because the pore fluid flow does not contribute significantly to the shear waves, a normalization function widely used in geotechnical engineering is sufficient to remove the effect of volume fraction on the shear waves. The increasing trend of the wave velocities with higher input frequency is also well captured by the coupled LBM-DEM simulations. Further work involves improving existing theoretical predictions of wave velocities for fluid-saturated granular materials, incorporating the dependence on input wave frequency, solid volume fraction, and anisotropy.