# Parameter optimization of phase-field-based LBM model for calculating capillary forces

R. Bouchard, N. Younes, Olivier Millet, A. Wautier

<figure><img src="image_1.png" /></figure>

## To cite this version:

R. Bouchard, N. Younes, Olivier Millet, A. Wautier. Parameter optimization of phase-field-based LBM model for calculating capillary forces. Computers and Geotechnics, 2024, 172, pp.106391. 10.1016/j.compgeo.2024.106391. hal-04689846

HAL Id: hal-04689846

https://hal.inrae.fr/hal-04689846v1

Submitted on 6 Jan 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Parameter optimization of phase-field-based LBM model for calculating capillary forces

R. Bouchard<sup>a</sup>, N. Younes<sup>a,b,c</sup>, O. Millet<sup>a,*</sup>, A. Wautier<sup>b</sup>

<sup>a</sup>LaSIE, UMR CNRS 7356, University of La Rochelle, La Rochelle, France

<sup>b</sup>INRAE, Aix Marseille University, RECOVER research unit, Aix-en-Provence, France

<sup>c</sup>Schulich School of Engineering, University of Calgary, Calgary, Canada

## Abstract

Recently, Younes et al. (2022) [1] developed a phase-field-based Lattice Boltzmann Method (LBM) model capable of capturing the merging of capillary bridges in a triplet configuration without relying on geometrical criteria, thus enabling an easy transition from pendular to funicular regimes. However, the computational time of the formation of capillary bridges is quite significant particularly when it comes to coupling a macroscopic method (e.g. FEM) with a DEM-LBM coupling model to simulate engineering structures.
In this study, we present a parameter optimization strategy aimed at accelerating the computational efficiency of the LBM alone. We show that even by not accounting for transient evolutions, we can still accurately obtain capillary forces with an error margin limited to 5% at equilibrium.

**Keywords:** Optimization, LBM, GPU, capillary force, 3D

## 1. Introduction

In many geotechnical applications, earthen structures (e.g. earthen dikes) can be seen as granular material under unsaturated conditions. These shallow unsaturated areas may be consistently exposed to drying and wetting cycles due to changing conditions, such as intense rainfall and heatwave events associated with the phenomenon of climate change, making dikes vulnerable to failure.

<footnote>*Correspondence:</footnote>

<footnote>Email address: olivier.millet@univ-lr.fr (O. Millet)</footnote>
For decades, researchers have devised analytical and charts-based meth-
ods to assess slope stability across various soil conditions, including dry, un-
saturated, and fully saturated soils [2, 3, 4, 5, 6, 7, 8]. However, these methods
usually rely on simplifications and assumptions that limit their applicability
in certain engineering contexts, especially under unsaturated conditions that
are not well understood by either researchers or engineers. Therefore, there is
the need for numerical models able to account for microscale and the underlying
physics of the granular assemblies. Although, for the last decade, numerous
studies have been carried out to examine granular soils behavior using the
Discrete Element Method (DEM) [9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20],
they mostly focus on cases of dry and fully saturated granular assemblies.
Regarding unsaturated granular assemblies capillary forces are frequently cal-
culated based on the analytical solution of the Young-Laplace equation and
then imposed as additional forces in DEM computations. This approach has
brought insights into longstanding questions, including but not limited to the
existence of effective stress under unsaturated conditions at the plastic limit
[21, 22, 23, 24, 25, 26, 27]. This approach is, unfortunately, restricted only to
the pendular regime, where capillary bridges exist only between pairs of par-
ticles. To model the funicular and capillary regimes, various methods have
been proposed for small granular assemblies. One such method involves using
the energy minimization technique [28]. For example, Miot et al. (2021) and
Di Renzo et al. (2020) [29, 30] employed this method along with the virtual
work principle to calculate capillary forces for doublets, triplets, and quadrup-
lets of spherical particles. To model the funicular and capillary regimes, var-
ious methods have been proposed for smaller granular assemblies. One such
method involves using the energy minimization technique [28]. For example,
Miot et al. (2021) and Di Renzo et al. (2020) [29, 30] employed this method
along with the virtual work principle to calculate capillary forces for doublets,
triplets, and quadruplets of spherical particles. However, despite providing
accurate results, this method has two main drawbacks: (1) it is limited to
analyzing steady-state problems and (2) requires specific geometrical criteria
to address capillary bridge rupture or coalescence. This complexity makes
it challenging to apply this method to larger granular assemblies. More re-
cently, Younes et al. (2022) [1] used a phase-field-based method involving the
Navier-Stokes and Allen-Cahn equations. Their findings showed that isolated
capillary bridges in configurations with three spherical particles merge upon
contact. Not only the coalescence was readily addressed, but the authors also
formulated a capillary force expression applicable all over regimes: pendular,
funicular, and capillary regimes. The following year, Younes et al. (2023) [31] applied the aforementioned method and coupled it with DEM to simulate large granular assemblies under unsaturated conditions for a large range of degrees of saturation $S_r$. This study provided promising results, particularly regarding the Soil-Water Characteristic Curve and capillary stress.

While accurate simulations of wetting and drying cycles were reported, it is important to note that simulations are time-consuming, with a single simulation requiring several hours or even days to complete. Therefore, this work focuses on the parametric optimization of a 3D phase-field Lattice Boltzmann Method (LBM) model used in [1, 31]. This paper presents a parametric optimization of the LBM and the DEM-LBM coupling for small granular assemblies, aiming to significantly reduce computational time without altering the obtained results for capillary forces. The primary assumption of this work is that our interest lies only in the converged value of the capillary forces, not the overall behavior of the air-water interface in transient evolutions. The optimized parameters can be used in any engineering scale computations relying on continuum mechanics methods (e.g., MPM/FEM/FDM) where DEM-LBM will be used as a constitutive model accounting for the physics of the underlying air-water-solid interface behavior.

The paper is organized as follows: a brief introduction to the phase-field-based LBM model is presented, along with the parameters involved in the optimization process and the conversion coefficients linking lattice units to physical ones. The convergence of capillary bridges is next examined by varying several parameters to study their impact on the final result. Finally, the selection of the best set of optimized parameters is based on the investigation of capillary forces while moving a particle in a doublet configuration.

## 2. Description of the numerical model

### 2.1. *Description of the governing equations*

The modeling of partially saturated media as a multiphasic system requires simulating capillary interfaces separating gas from the liquid. The Navier-Stokes equations (NSE) modeling the behavior of the fluid were coupled to Allen-Cahn equations (ACE) to simulate the moving air/water interface [32, 1, 31].

However, even if in reality an interface has a very thin (negligible) thickness, it is currently impossible to simulate such sharp behavior numerically.
To overcome this, we are forced to impose a certain thickness on our interface<fnref n="1" />, noted here W. This leads to the evaluation of a characteristic phase function $\phi(\mathbf{r}, t)$ in the fluid domain, which, depending on its value, tells us whether we are in air or in water. Conventionally $\phi = 0$ in the gas, $\phi = 1$ in the liquid; $\phi$ represents the interface for $\phi \in ]0; 1[$ as a diffuse interface zone where the two phases coexist. However, it has been proved that when $\phi = \frac{1}{2}$, corresponds to the interface separating gas from liquid [32, 1, 31].

<figure><img src="image_3.png" /><figcaption>Figure 1: Scheme of the interface. Taken from [1]</figcaption></figure>

At the equilibrium, the phase-field function is defined as

$$ \phi^{eq}(\xi) = \frac{1}{2} \left[ 1 + \tanh \left( \frac{2\xi}{W} \right) \right] \qquad (1) $$

where W is the interface thickness. From Eq. (1), we can define the parameter $\lambda$ as the derivative of $\phi^{eq}$ with respect to $\xi$, i.e.:

$$ \lambda = |\nabla \phi|^{eq} = \frac{4\phi(1-\phi)}{W}. \qquad (2) $$

We can then write the ACE, which is derived from the classical advection equation and allows us to compute the evolution of the phase field $\phi(\mathbf{r}, t)$, as

$$ \frac{\partial \phi(\mathbf{r}, t)}{\partial t} + \nabla \cdot [\phi(\mathbf{r}, t) \mathbf{u}] = \nabla \cdot [M(\nabla \phi(\mathbf{r}, t) - \lambda \mathbf{n})] \qquad (3) $$

<footnote n="1">¹Here it will be equal to 5 lattices.</footnote>
where $\mathbf{u}$ is the fluid velocity governed by Navier-Stokes equation, $M$ the interface mobility, $\mathbf{n} = \frac{\nabla\phi}{|\nabla\phi|}$ is the interface unit normal vector and $\mathbf{r}$ the space variable. In Eq. (3), the velocity is decomposed into two components:

(i) advected external fluid velocity $\mathbf{u}$

(ii) normal interface speed $\mathbf{u}_n$

The first component is usually determined from the Navier-Stokes equation, whereas the second, driven by the thermodynamics of the problem, is taken to be proportional to the interface curvature through mobility such that $\mathbf{u}_n = M\kappa\mathbf{n}$ where $\kappa$ is the curvature [33, 34].

As such, the ACE provides an implicit tracking of the interface which makes it suitable for problems in which the interface motion depends on external conditions, in particular, wetting condition [32, 1, 35, 36].

<figure><img src="image_4.png" /><figcaption>Figure 2: Explanatory diagram of the phase evolution through the interface.</figcaption></figure>

As previously mentioned, to consider viscous, immiscible, incompressible two-phase fluid flow, the Navier-Stokes equations are considered, with a supplementary source term $\mathbf{F}_s$ characterizing the surface tension force
$$
\left\{
\begin{aligned}
\rho \left( \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) &= -\nabla p + \nabla \cdot \left[ \rho \nu (\nabla \mathbf{u} + \nabla \mathbf{u}^T) \right] + \mathbf{F}_s && (4a) \\
\nabla \cdot \mathbf{u} &= 0 && (4b)
\end{aligned}
\right.
$$

110

111

112

113

where **u** is the fluid velocity, **v** is the kinematic viscosity, **p** is the hydrodynamic pressure. The fluid density field ρ(**r**,*t*) is calculated with linear interpolation between liquid and gas densities ρ*<sub>l</sub>* and ρ*<sub>g</sub>* , i.e.

$$
\rho(\mathbf{r}, t) = \phi(\rho_l - \rho_g) + \rho_g \quad (5)
$$

114

115

116

The link between the ACE and NSE equations is made at two levels. Firstly,
in a direct way via the density which depends explicitly on the phase function.
But also through the surface tension term F_s given by (see [32, 35, 36, 1, 37])

117

:

$$
\mathbf{F}_s = \mu_\phi \nabla \phi \qquad (6)
$$

118

119

where $\mu_\phi$ is the chemical potential defined as the minimum of the mixing energy ([32, 1, 38])

$$
\mu_{\phi} = 4\beta\phi(\phi - 1)\left(\phi - \frac{1}{2}\right) - k\nabla^2\phi_{\xi} \quad (7)
$$

120

121

122

In (7), $\nabla^2\phi_\xi$ denotes the second derivative according to the normal coor-
dinate $\xi$, $k$ and $\beta$ are constants related to the surface tension $\gamma$ and the
interface thickness $W$, i.e.

$$
k = \frac{3}{2}\gamma W \qquad \beta = \frac{12\gamma}{W}. \tag{8}
$$

123

124

125

126

It is worth noting that equilibrium is established whenever the chemical po-
tential µφ vanishes [32, 1, 38]. As such, φeq can be deduced as the solution of
the differential equation (7) with µφ equal to zero with the parametrization
of Fig. 1.

127

2.2. *Capillary force calculation*

128

129

130

The capillary forces $\hat{F}_f$ exerted by the fluid on the solid particles are computed using the method proposed in [1]

$$
\hat{F}_f = \int_{\Omega} \sigma_f n_{\Omega} dS + \gamma \int_{\Gamma} m dl \quad (9)
$$
where $\Gamma$ is the intersection of $\Omega$ with the capillary interface as shown in Fig. 3 and where $\boldsymbol{m} = (\boldsymbol{n} \times \boldsymbol{n}_{\Omega}) \times \boldsymbol{n}$. They can be decomposed into two parts, a pressure force $F_p$ and a surface tension force $F_{\gamma}$, respectively the first and the second term in the right hand side of (9). Expression (9) is then descritized on the nodes (lattices) leading to expression (72) of [1], which allows to handle all capillary bridges configuration. The method proposed is very general for computing capillary forces in any configuration (capillary doublets or triplets, before or after coalescence, etc...)

<figure><img src="image_5.png" /><figcaption>Figure 3: Capillary bridge and integration domain within a sharp interface (Fig. from [1]).</figcaption></figure>

## 2.3. Lattice Boltzmann discretization

The kinetic theory describes the evolution of a system of molecular particles defined by their positions and velocities with a probabilistic framework. This theory was successfully developed in the 19th century, mainly by Maxwell and Boltzmann, to model gases as a set of molecules. In recent
years, kinetic theory has been applied to a large number of problems in bi- ology, medicine, and the study of human behavior. We are interested here in the Lattice Boltzmann method, based on the Boltzmann equation, which was born from this theory:

$$ \underbrace{(\partial_t f + \mathbf{v} \cdot \nabla f)}_{\text{propagation}}(t, \mathbf{r}, \mathbf{v}) = \underbrace{Q(f(t, \mathbf{r}))(\mathbf{v})}_{\text{collision operator}} \qquad (10) $$

where $f$ is the Probability Density Functions (PDFs) of the presence of a particle, $t$ is the time variable, $\mathbf{r}$ the space variable, $\mathbf{v}$ the velocity of the particles and $Q$ the collision operator which models the variation in the velocities of the molecules due to collisions. In this work, the Bhatnagar-Gross-Krook (BGK) collision operator is used [35, 36, 32, 1] as well as two PDFs, $g(\mathbf{r}, t)$ and $h(\mathbf{r}, t)$ which correspond to ACE and NSE, respectively. The associated discretizations used are given by

$$ g_i(\mathbf{r} + \mathbf{v}_i \Delta t, t + \Delta t) - g_i(\mathbf{r}, t) = -\frac{\Delta t}{\Gamma_g} [g_i(\mathbf{r}, t) - g_i^{eq}(\mathbf{r}, t)] + \Delta t G_i(\mathbf{r}, t) \quad (11a) $$

$$ \underbrace{h_i(\mathbf{r} + \mathbf{v}_i \Delta t, t + \Delta t) - h_i(\mathbf{r}, t)}_{\text{propagation}} = \underbrace{-\frac{\Delta t}{\Gamma_h} [h_i(\mathbf{r}, t) - h_i^{eq}(\mathbf{r}, t)]}_{\text{BGK operator}} + \underbrace{\Delta t H_i(\mathbf{r}, t)}_{\text{source term}} \quad (11b) $$

where $\Delta t$ is the LBM time step, $\Gamma_g$ and $\Gamma_h$ are the dimensionless relaxation times for $g(\mathbf{r}, t)$ and $h(\mathbf{r}, t)$, respectively. In this paper, unlike $\Gamma_g$ which is a constant, $\Gamma_h$ is computed a linear interpolation of the gas ($\tau^g$) and liquid ($\tau^l$) dimensionless relaxation times:

$$ \Gamma_h = \phi (\tau_h^l - \tau_h^g) + \tau_h^g \qquad (12) $$

The discrete velocities ($\mathbf{v}_i$) depend on the LBM discretization scheme (D2Q9, D3Q15, D3Q19 etc..).

At the equilibrium the PFDs are expressed as :

$$ g_i^{eq}(\mathbf{r}, t) = \omega_i \phi \left( 1 + \frac{\mathbf{v}_i \cdot \mathbf{u}}{c_s^2} \right) \qquad (13a) $$

$$ h_i^{eq} = \begin{cases} \frac{p}{c_s^2}(\omega_i - 1) + \rho s_i(\mathbf{u}), & \text{if } i = 0 \\ \frac{p}{c_s^2}\omega_i + \rho s_i(\mathbf{u}), & \text{if } i \neq 0 \end{cases} \qquad (13b) $$
164 with

$$s_i(\mathbf{u}) = \omega_i \left[ \frac{\mathbf{v}_i \cdot \mathbf{u}}{c_s^2} + \frac{(\mathbf{v}_i \cdot \mathbf{u})^2}{2c_s^4} - \frac{\mathbf{u} \cdot \mathbf{u}}{2c_s^2} \right] \quad (14)$$

where $\mathbf{u} = \mathbf{u}(\mathbf{r}, t)$ is the fluid velocity, $\rho$ is the local density defined in Eq. (5), $c_s$ the lattice speed of the sound velocity (dependent of the discretization) and $p$ is the pressure. The source terms $H_i$ and $G_i$ can be written as:

$$H_i = \left(1 - \frac{\Delta t}{2\tau_h}\right) \omega_i \left( \frac{\mathbf{v}_i \cdot \mu_\phi \nabla \phi}{c_s^2} + \frac{(\rho_l - \rho_g) \mathbf{u} \nabla \phi : \mathbf{v}_i \mathbf{v}_i}{c_s^2} \right) \quad (15a)$$

$$G_i = \left(1 - \frac{\Delta t}{2\tau_g}\right) \frac{\omega_i \mathbf{v}_i \cdot [\partial_t(\phi \mathbf{u}) + c_s^2 \lambda \mathbf{n}]}{c_s^2} \quad (15b)$$

Readers can refer to [32, 1] for more details about the phase-field-based model as well as the boundary conditions (wetting and bounce-back conditions) used in this paper. Furthermore, it has been choosen to use the scheme D3Q109 solve the Eqs. (11a) and (11b), which imposed the speed of sound in lattice units $c_s = \frac{1}{\sqrt{3}}$.

## 2.4. LBM parameters

As all the variables used in the BE resolution are dimensionless<fnref n="2" />, it is necessary to create a system of conversions from dimensional to dimensionless variables. It turns out that with only three units such as the meter, time, and mass it is possible to describe all the variables in our system of equations. Consequently you have to set 3 independent coefficient of conversions. The link between the physical and scaled variables<fnref n="3" /> is made using conversion coefficients<fnref n="4" /> which, for the three main variables of the problem, are presented in the table 1.

<footnote n="2">²We use the notation $lu$ for "lattice unit", which describes dimensionless parameters.</footnote>

<footnote n="3">³All scaled variables, linked to the LBM, will be surmounted by the notation $\sim$.</footnote>

<footnote n="4">⁴All the coefficient of conversion are noted as $C_{physical\ parameter}$.</footnote>
<table><caption>Table 1: List of parameters and associated units.</caption><thead><tr><th>Parameters</th><th>Symbols</th><th>Units</th><th>Symbols (lu)</th><th>Coefficient of conversion</th></tr></thead><tbody><tr><td>Density</td><td>ρ</td><td>kg.m<sup>-3</sup></td><td>ρ̃</td><td>C<sub>ρ</sub> (kg.m<sup>-3</sup>)</td></tr><tr><td>Size mesh</td><td>Δx</td><td>m</td><td>Δx̃</td><td>C<sub>L</sub> (m)</td></tr><tr><td>Time step</td><td>Δt</td><td>s</td><td>Δt̃</td><td>C<sub>t</sub> (s)</td></tr></tbody></table>

183 The link between parameters is given by

$$ \Delta x = C_L \Delta x \qquad \Delta t = C_t \Delta t \qquad \rho = C_\rho \rho \qquad (16) $$

184 We can now express all the other coefficients of conversions from $C_L$, $C_t$ and $C_\rho$. As the unit of the Newton can be expressed as $kg.m.s^{-2}$, we can define the forces conversion coefficient $C_F$ as:

$$ C_F = C_\rho \frac{C_L^4}{C_t^2} \qquad (17) $$

187 The other conversion coefficients are shown in Table 2.

<table><caption>Table 2: List of the physical parameters and the associated coefficients of conversion</caption><thead><tr><th>Associated physical parameter</th><th>γ</th><th>P</th><th>ν</th></tr></thead><tbody><tr><td>Coefficient of conversion</td><td>C<sub>γ</sub> = C<sub>ρ</sub> C<sub>L</sub><sup>3</sup> / C<sub>t</sub><sup>2</sup></td><td>C<sub>P</sub> = C<sub>ρ</sub> C<sub>L</sub><sup>2</sup> / C<sub>t</sub><sup>2</sup></td><td>C<sub>ν</sub> = C<sub>L</sub><sup>2</sup> / C<sub>t</sub></td></tr><tr><td>Unit</td><td>N.m<sup>-1</sup></td><td>N.m<sup>-2</sup></td><td>m.s<sup>-1</sup></td></tr></tbody></table>

188 It is important to note that in a given numerical simulation, the physical parameters (kinematic viscosity and densities) are given from by the physics of the problem, but the mesh size and the time step of the simulation must be adapted. This is why, in order to respect certain mathematical and physical considerations, if the mesh size is defined in an arbitrary way by the user, the time step has to be calculated directly from the parameters of the simulation. 189 190 191 192 193 It has been chosen to express the LBM time step conversion coefficient $C_t$, as a function of the surface tension conversion coefficient $C_\gamma$, which gives:

$$ C_t = \sqrt{\frac{C_\rho C_L^3}{C_\gamma}} \qquad (18) $$
where $C_γ$ becomes a parameter that has to be set at the beginning of the simulation. The choice of calculating the time step using (18), which is linked to the surface tension is arbitrary. It could have been calculated according to the coefficient of conversion of the kinematic viscosity $C_ν$. Usually, $\tilde{\Delta}x$ and $\tilde{\Delta}t$ are both taken equal to 1 lu, furthermore, the value of $C_ρ$ has been selected to be equal to 1 as well. This sets of the 3 independent scaling coefficients of conversion between physical units and lattice units

Table 3 lists the different LBM parameters pointing out which ones can be freely defined by the user.

<table><caption>Table 3: List of LBM parameters and conversion coefficients</caption><thead><tr><th>Parameter</th><th>Allocation</th><th>Value</th></tr></thead><tbody><tr><td>ρ</td><td>Imposed by the physics</td><td></td></tr><tr><td>C<sub>ρ</sub></td><td>User choice for lu conversion</td><td>1 kg.m<sup>-3</sup>.lu<sup>-1</sup></td></tr><tr><td>ρ̃</td><td>Calculated</td><td>ρ</td></tr><tr><td>Δ<sub>x</sub></td><td>User choice for space discretization</td><td></td></tr><tr><td>C<sub>L</sub></td><td>Calculated</td><td>Δ<sub>x</sub></td></tr><tr><td>Δ̃<sub>x</sub></td><td>User choice for lu conversion</td><td>1 lu</td></tr><tr><td>γ</td><td>Imposed by the physics</td><td></td></tr><tr><td>C<sub>γ</sub></td><td>Calculated</td><td>γ/γ̃</td></tr><tr><td>γ̃</td><td>User choice for lu conversion</td><td></td></tr><tr><td>Δ<sub>t</sub></td><td>Calculated</td><td>C<sub>t</sub></td></tr><tr><td>C<sub>t</sub></td><td>Calculated</td><td>√(C<sub>ρ</sub>C<sub>L</sub><sup>3</sup>/C<sub>γ</sub>)</td></tr><tr><td>Δ̃<sub>t</sub></td><td>User choice for time discretization</td><td>1 lu</td></tr></tbody></table>

## 3. Optimization of LBM parameters

In this section, we aim to study the influence of the LBM parameters mentioned previously (mainly, $\tilde{\gamma}$, $\nu^{l,g}$ and $C_ρ$) on the value of capillary force at convergence, which is the most important calculated parameters that need to be updated in a coupled DEM-LBM simulation for partially saturated granular assembly.

To optimize our DEM/LBM coupled code for simulating the behavior of
thousands of particles under varying relative humidity, we aim to significantly improve computing speed. We are not interested in the behavior of the fluid during the simulation, but only in the value of the capillary forces at the equilibrium. Therefore, we will modify some parameters (DEM and/or LBM) in order to speed up the calculation time as much as possible. Consequently transient evolutions will no longer be meaningful

The kinematic viscosity $\nu$ in NSE (4a) is linked to the relaxation time $\Gamma_h$ through the relation:

$$ \nu = c_s^2 \left( \Gamma_h - \frac{\Delta t}{2} \right) \qquad (19) $$

It is possible to rewrite Eq. (19) to determine an expression of the relaxation time in the lattice units. Indeed, we have:

$$ \nu^{\ell,g} = c_s^2 \left( \Gamma_h^{\ell,g} - \frac{\Delta t}{2} \right) = \tilde{c}_s^2 \frac{C_L^2}{C_t^2} \left( \tilde{\tau}_h^{\ell,g} C_t - \frac{\tilde{\Delta}t C_t}{2} \right) = \underbrace{\tilde{c}_s^2 \left( \tilde{\tau}_h^{\ell,g} - \frac{\tilde{\Delta}t}{2} \right)}_{\tilde{\nu}^{\ell,g}} \frac{\tilde{C}_L}{C_t} \qquad (20) $$

where $\tilde{\tau}_h^{\ell,g}$ is the dimensionless value of $\Gamma_h^{\ell,g}$, which leads to

$$ \tilde{\tau}_h^{\ell,g} = \frac{C_t \nu^{\ell,g}}{\tilde{c}_s^2 C_L^2} + \frac{\tilde{\Delta}t}{2} \qquad (21) $$

This expression of the dimensionless relaxation time is essential for the stability of the code. Indeed, according to [39], the relaxation time must be greater than $\frac{1}{2}$ for the LBM scheme to be stable. Nevertheless, still according to [39], for values of $\tilde{\tau}_h^{\ell,g}$ greater than 1, oscillations in the results may appear. Therefore, to ensure stability, it is preferable to always keep this parameter in the interval $]1/2, 1[$.

By analyzing Eq. (21), we can note that by construction the stability of the code is always insured thanks to the parameter $\tilde{\Delta}t/2$.

In the second term of the right-hand side of (21), the parameters $c_s$ and $C_L$ are fixed upstream of the simulation and are not modifiable. Moreover, as the kinematic viscosities of the fluids $\nu^{\ell,g}$ (liquid and gas) are only involved
237 in the calculation of the relaxation time in the LBM method and not in
238 the resolution of the static equilibrium of capillary bridges, by adjusting this
239 parameter, it is possible to artificially modify the value of the relaxation time
240 without influencing the final result of the capillary forces.

241 On the other hand, it is possible to modify the LBM time step $C_t$ given by
242 Eq. (18) to play with the time discretization of the problem. Once $C_L$ and $C_\rho$
243 are fixed, the only parameter that controls the LBM time step is $C_\gamma$ directly
244 linked to the LBM surface tension $\tilde{\gamma}$ in lattice units. Ideally, it is preferable to
245 maximize $\Delta t_{\text{LBM}}$ for numerical efficiency, while also guaranteeing numerical
246 stability. The pros and cons of large values $\tilde{\gamma}$ are discussed in the following.

### 3.1. Case study description

247 A single capillary bridge between two spherical particles is considered.
248 It has been already proved in [1] that with the LBM simulations performed
249 with the code described in section 2, an initialization of a droplet at the
250 mid-distance between particles will ultimately converge to a capillary bridge
251 as shown in Fig. 4.

<figure><img src="image_6.png" /><figcaption>Figure 4: From the initialization of a droplet at mid-distance to a converged capillary bridge.</figcaption></figure>

253 For this simulation, we chose a discretization of the domain $N_x \times N_y \times N_z =$
254 173 × 173 × 356, the particles having a radius of $R = 4$ mm. This discretization
255 leads to $\Delta x \approx 5 \times 10^{-5}$ m fixed once for all in the simulation.
<table><caption>Table 4: Simulation's physical parameters.</caption><thead><tr><th>Parameters</th><th>Symbols</th><th>Values</th><th>Units</th></tr></thead><tbody><tr><td>Surface tension (water)</td><td>γ</td><td>0.0720</td><td>N.m<sup>-1</sup></td></tr><tr><td>Kinetic viscosity (water)</td><td>ν<sup>l</sup></td><td>1.3 × 10<sup>-6</sup></td><td>m<sup>2</sup>.s<sup>-1</sup></td></tr><tr><td>Kinetic viscosity (air)</td><td>ν<sup>g</sup></td><td>1.8 × 10<sup>-5</sup></td><td>m<sup>2</sup>.s<sup>-1</sup></td></tr><tr><td>Radii</td><td>R</td><td>4</td><td>mm</td></tr><tr><td>Initial distance</td><td>d</td><td>0.5</td><td>mm</td></tr><tr><td>Temperature</td><td>T</td><td>10</td><td>°C</td></tr><tr><td>Contact angle</td><td>θ</td><td>30</td><td>°</td></tr></tbody></table>

## 3.2. From a droplet to a capillary bridge

As the shape of a capillary bridge can be complex, it's easier and more natural to let a water droplet converge to a stable capillary bridge, than to directly initialize a bridge between two grains. The problem with this method is that the first convergence step induces a non-negligible additional computation time. In this section, we will study the influence of the parameter $\tilde{\gamma}$ on the convergence of the drop through four values: 15, 75, 150, and 200 $lu$.

We can observe in Fig. 5 that during the first iterations after the creation of the water droplet, the calculated capillary force involves rather large fluctuations before convergence. Furthermore, Fig. 5 emphasizes that the value of the surface tension in lattice units affects very little the value of the capillary force at convergence. The difference between the lowest and highest values of capillary forces $F^{cap}$ is $40\mu N$, equivalent to 2%, making it negligible and therefore can be disregarded in the rest of the analysis. On the other hand, increasing the surface tension in lattice units enables the system to reach longer physical times, as the LBM time step of the simulation is linked to the surface tension through (18). Table 5, summarizes the value of the LBM time step $\Delta t_{LBM}$ and the relaxation times associated to the value of $\tilde{\gamma}$. The other parameters remain the same that those presented in table 4.
<table><caption>Table 5: Value of relaxation times and LBM time step for different $\tilde{\gamma}$ for water ($\gamma = 0.0720 \text{ N.m}^{-1}$).</caption><thead><tr><th>γ̃ (lu)</th><th>C<sub>γ</sub> (N.m<sup>-1</sup>)</th><th>(τ̃<sub>h</sub><sup>l</sup>; τ̃<sub>h</sub><sup>g</sup>) (s)</th><th>Δt<sub>LBM</sub> (s)</th></tr></thead><tbody><tr><td>15</td><td>4.8 × 10<sup>-3</sup></td><td>(0.508 ; 0.610)</td><td>5.103 × 10<sup>-6</sup></td></tr><tr><td>75</td><td>9.6 × 10<sup>-4</sup></td><td>(0.518 ; 0.746)</td><td>1.141 × 10<sup>-5</sup></td></tr><tr><td>150</td><td>4.8 × 10<sup>-4</sup></td><td>(0.525 ; 0.848)</td><td>1.614 × 10<sup>-5</sup></td></tr><tr><td>200</td><td>3.6 × 10<sup>-4</sup></td><td>(0.529 ; 0.902)</td><td>1.863 × 10<sup>-5</sup></td></tr></tbody></table>

We can also compare the values of the capillary forces $F^{cap}$ obtained by the simulation to experimental results detailed in [40] for a volume of 10 µL between two spherical glass beads with a separation distance of 0.5 mm, marked as the continuous purple line on Fig. 5. An error of 4% for the converged value of $F^{cap}$ is observed which is quite acceptable and reliable.
<figure><img src="image_7.png" /><figcaption>Figure 5: Evolution of the capillary force for different values of gamma in terms of (a) the number of iterations and (b) time.</figcaption></figure>
Finally, we compare in Fig. 6 the relative errors $\epsilon$ which are calculated

as

$$ \epsilon = 100 \times \frac{F_i^{cap} - F_{converged}^{cap}}{F_{converged}^{cap}} \qquad (22) $$

where $F_{converged}^{cap}$ is the value of the capillary force at the last iteration (here 9000) and $F_i^{cap}$ is the value of the capillary force at iteration i. We observe that for $\tilde{\gamma} = 75$ lu fewer oscillations occur and allow a precision under 1% of error in only 5500 iterations, which is the minimum between the four tested values. Thus, during the next study, the value used the the surface tension in lattice unit ($\tilde{\gamma}$) will be fixed to 75 lu and we will let at least 5500 iterations for the convergence of the droplet.

<figure><img src="image_8.png" /><figcaption>Figure 6: Evolution of the relative error for the capillary force of a water doublet for different surface tension during the convergence versus the number of iterations.</figcaption></figure>

### 3.3. Convergence of the capillary forces for a moving capillary bridge

Unlike in [1], where each simulation corresponded to a separation distance, here we will continuously space the upper particle in order to accurately model the experiments presented in [40]. After each displacement of the top grain, we let the capillary bridge converge to a new stable shape. The following parametric study will be conducted to determine the best set
296 of parameters to minimize the number of LBM iterations necessary to get an accurate value of the capillary force after moving particles.
297

<figure><img src="image_9.png" /><figcaption>Figure 7: Scheme of the numerical experiment of moving particles of a displacement $\delta$.</figcaption></figure>

<table><caption>Table 6: LBM physical parameters for an air/water interface and moving particles.</caption><thead><tr><th>Parameters</th><th>Symbols</th><th>Values</th><th>Units</th></tr></thead><tbody><tr><td>Surface tension (water)</td><td>γ</td><td>0.0720</td><td>N.m<sup>-1</sup></td></tr><tr><td>Kinetic viscosity (water)</td><td>ν<sup>l</sup></td><td>1.3 × 10<sup>-6</sup></td><td>m<sup>2</sup>.s<sup>-1</sup></td></tr><tr><td>Kinetic viscosity (air)</td><td>ν<sup>g</sup></td><td>1.8 × 10<sup>-5</sup></td><td>m<sup>2</sup>.s<sup>-1</sup></td></tr><tr><td>Radii</td><td>R</td><td>4</td><td>mm</td></tr><tr><td>Distance</td><td>d</td><td>0.5</td><td>mm</td></tr><tr><td>Temperature</td><td>T</td><td>10</td><td>°C</td></tr><tr><td>Volume of water</td><td>V</td><td>10</td><td>μL</td></tr><tr><td>Contact angle</td><td>θ</td><td>30</td><td>°</td></tr></tbody></table>

We saw in the last sections that 5500 LBM iterations were necessary to obtain the first convergence of a spherical drop of water to a capillary bridge, so we will wait for 6000 LBM iterations before moving the top particle of a distance of 0.2mm along the z-axis, as in the experiments presented [41, 40]. In Fig. 8a we show the evolution of capillary forces for 18 000 iterations, corresponding to two consecutive displacements of particles.
### 3.3.1. Influence of $\tilde{\gamma}$ on the convergence of the capillary forces

We can observe in Fig. 8a the three plateaus correspond to the displacement of the particle. When the separation distance between particles increases, the associated capillary force decreases. Furthermore, and as previously seen in section 3.2, we notice that the value of $\tilde{\gamma}$ has not a significant influence on the converged value of the associated capillary force. Moreover, the difference on $F^{cap}$ is about 1% for the extremum values of $\tilde{\gamma}$. Finally, as previously quoted, the larger $\tilde{\gamma}$, the larger simulated time.

In this new simulation configuration, the calculation of the relative error is still done using (22) where $F_{converged}^{cap}$ is the converged value of the concerned plateau$^5$.

We can note that after each displacement, the capillary force exhibits fluctuations with a maximum amplitude of 5% before it converges to a plateau (Fig. 8b). For the future simulation, we will maintain the value of $\tilde{\gamma}$ to 75 $lu$, to minimize the computation time while conserving the stability of the code.

This means 3500 iterations if we want a relative error of less than 1%, or only 1000 iterations if we want an error of less than 5%.

<footnote>$^5$Thus, there is three different value of $F_{converged}^{cap}$, one for each plateau.</footnote>
<figure><img src="image_10.png" /><figcaption>(a) Evolution of the capillary force for different values of $\tilde{\gamma}$ during the particle displacement in terms of the number of iterations</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) Evolution of the relative error ($\epsilon$) on $F^{cap}$ ($N.m^{-1}$) after the first and the second displacement in terms of iterations.</figcaption></figure>

Figure 8
### 3.3.2. Influence of $C_ρ$ on the convergence of the capillary forces

According to (18), the LBM time step depends on the coefficient of conversion of the density $C_ρ$. Usually considered equal to 1 lu, it may be increased to increase the LBM time step. In some situations, this may constitute an advantage, because of the limitation of increasing $\tilde{\gamma}$, due to the stability condition on the interface velocity. Indeed, according to [39] the LBM stability is ensured if

$$|\mathbf{u}| < \sqrt{\frac{2}{3} \frac{\Delta x}{\Delta t}} \qquad (23)$$

where $\mathbf{u}$ denotes the fluid velocity. By definition of $\Delta x$ (16) and $\Delta t$ (18), we obtain

$$|\mathbf{u}| < \sqrt{\frac{2\gamma}{3C_\rho\tilde{\gamma}C_L}} \qquad (24)$$

This implies that when we increase the value of $C_ρ$ or $\tilde{\gamma}$, we reduce the maximal fluid velocity keeping the numerical stability. On the other hand, the factor $C_L$ (representing the mesh size), should be the smallest as possible to ensure better stability for the LBM simulations. In the numerical simulations of capillary bridges that are considered in this work, the fluid velocity is very small (we are close to the quasi-static equilibrium) and the condition (24) is always fulfilled. Indeed, Fig. 9 emphasizes, that the limitation value induced by the Eqs. (23) and (24), in the most limited case ($C_ρ = 4$ and $\tilde{\gamma} = 200$), is never reached by the maximal value of the velocity field during the simulation. So that $C_ρ$ can be increased without any real consequence on the numerical stability.
<figure><img src="image_12.png" /><figcaption>Figure 9: comparison of the evolution of the maximum velocity of the fluid phase for $\tilde{\gamma} = 200$ and $C_\rho = 4$ with the stability condition according to [39].</figcaption></figure>

During our simulations, we consider values of $C_\rho$ between 1 and 4, because higher values require a decrease in the value of the surface tension in lattice units ($\tilde{\gamma}$) to meet the condition $\tau < 1$. As expected, we can see in figure 10 that the final simulation time differs according to the value of $C_\rho$. Once again, for values ranging from 1 to 3, we observe only small variations of the final values of the capillary forces. On the other hand, we can see in Fig. 10b that the parameter $C_\rho$ does not seem to affect at all the convergence speed of the capillary forces, but limit a little bit the oscillations. To conclude, a safe strategy to accelerate LBM computation involves tripling the value of $C_\rho$ in our upcoming calculations. However, in order to test the influence of other parameters, we will maintain the value of $C_\rho$ to 1 during the next studies.
<figure><img src="image_13.png" /><figcaption>(a) Evolution of the capillary force for different values of $C_ρ$ during the particle displacement in terms of the number of iterations</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(b) Evolution of the relative error ($\epsilon$) on $F^{cap}$ versus the number of iterations during simulations with particle displacement.</figcaption></figure>

Figure 10
### 3.3.3. Influence of $\nu^l$ and $\nu^g$ on the convergence of the capillary forces

As previously stated, in order to maintain a stable simulation and prevent potential oscillations in the results, we need to keep the relaxation times $\tau^g$ and $\tau^l$ within the interval $[0.5;1[$. What's more, we know empirically that the further away we are from the bounds of this interval, the more stable our code is. However, kinematic viscosity is a physical parameter that only appears in the relaxation time calculation Eq. (21). Varying the value of this parameter allows us to modulate the stability of the code, without modifying the value of the capillary forces obtained at static equilibrium.

The protocol remains the same as in the previous sections, i.e. we let the initial drop converge towards a capillary bridge during 6000 iterations. Then, we increase the separation distance between the two particles of $\delta = 0.2\text{mm}$ and we let the capillary bridge converge again for 6000 new iterations. This procedure is repeated two times.

In this simulation, we use the values of parameters presented in the table 7 for $\nu^l$ and the associated relaxation time $\tilde{\tau}_h^l$ (the relaxation time used for the gas is constant as $\nu^g$ is constant).

<table><caption>Table 7: Evolution of the relaxation time in function of the liquid kinetic viscosity.</caption><thead><tr><th>ν<sup>l</sup>(m.s<sup>-2</sup>)</th><th>(τ̃<sub>h</sub><sup>l</sup>; τ̃<sub>h</sub><sup>g</sup>) (s)</th><th>C<sub>t</sub> (s)</th></tr></thead><tbody><tr><td>0.9 × 10<sup>-6</sup></td><td>(0.512 ; 0.746)</td><td>1.14 × 10<sup>-5</sup></td></tr><tr><td>1.1 × 10<sup>-6</sup></td><td>(0.515 ; 0.746)</td><td>1.14 × 10<sup>-5</sup></td></tr><tr><td>1.3 × 10<sup>-6</sup></td><td>(0.518 ; 0.746)</td><td>1.14 × 10<sup>-5</sup></td></tr><tr><td>1.5 × 10<sup>-6</sup></td><td>(0.520 ; 0.746)</td><td>1.14 × 10<sup>-5</sup></td></tr><tr><td>1.7 × 10<sup>-6</sup></td><td>(0.523 ; 0.746)</td><td>1.14 × 10<sup>-5</sup></td></tr></tbody></table>
<figure><img src="image_15.png" /><figcaption>Figure 11: Evolution of the capillary force for different values of $v^l$ during the particle displacement in terms of time.</figcaption></figure>

As expected, Fig. 11 confirms that the kinematic viscosity of the fluid has no influence on the resulting capillary force at convergence but seems to induce a higher amplitude of the oscillation for the first convergence and do not influence the value of time step. This result allows us to affirm that we can increase our time step $\Delta t$ by decreasing the viscosity in a proportional way while keeping a stable relaxation time, in order to decrease our calculation time (see Eq. (21)).

Furthermore, the same numerical study was performed by changing the value of $v^g$ (table 8) and keeping fixed the value of $v^l = 1.3 \times 10^{-6}$. In a similar way, changing the value of $v^g$ has no influence on the values of capillary forces at equilibrium (Fig. 12).
<table><caption>Table 8: Variation of the relaxation time in function of the liquid kinetic viscosity.</caption><thead><tr><th>ν<sup>g</sup>(m.s<sup>-2</sup>)</th><th>(τ̃<sub>h</sub><sup>l</sup>; τ̃<sub>h</sub><sup>g</sup>) (s)</th><th>C<sub>t</sub> (s)</th></tr></thead><tbody><tr><td>1.1 × 10<sup>-5</sup></td><td>(0.518 ; 0.691)</td><td>1.14 × 10<sup>-5</sup></td></tr><tr><td>1.5 × 10<sup>-5</sup></td><td>(0.518 ; 0.719)</td><td>1.14 × 10<sup>-5</sup></td></tr><tr><td>1.8 × 10<sup>-5</sup></td><td>(0.518 ; 0.746)</td><td>1.14 × 10<sup>-5</sup></td></tr><tr><td>2.1 × 10<sup>-5</sup></td><td>(0.518 ; 0.773)</td><td>1.14 × 10<sup>-5</sup></td></tr><tr><td>2.7 × 10<sup>-5</sup></td><td>(0.518 ; 0.801)</td><td>1.14 × 10<sup>-5</sup></td></tr></tbody></table>

<figure><img src="image_16.png" /><figcaption>Figure 12: Evolution of the capillary force during the particle displacement force for different values of $ν^g$ during the particle displacement in terms of time.</figcaption></figure>

These findings confirm that in simulations aimed at swiftly capturing capillary forces within capillary bridges, without accounting for the produced flow's physical aspects, adjusting the physical values of the kinematic viscosities (gas and liquid) can uphold the stability of the code.

## 4. Conclusion

In this work, we studied the influence of parameters linked to LBM simulations on the convergence speed of capillary forces. We have proposed a
strategy for optimizing parameters in LBM simulations.

First, a study on the convergence of a droplet to form a stable, immobile capillary bridge, for different values of surface tension in lattice units $\tilde{\gamma}$ revealed that 5500 iterations are required for this crucial step.

Additionally, examining a capillary bridge established between two spherical particles, where the upper sphere moves upwards, allowed us to identify the optimal parameter set. This minimizes calculation times while guaranteeing both code stability as well as accurate results.

As for the moving particle, our initial simulations, conducted with non-optimized parameters, took roughly 4 hours using a Quadro RTX 5000 16 GB. However, we demonstrated that it is possible to enhance and stabilize our simulations by adjusting the values of the lattice surface tension $\tilde{\gamma}$, the coefficient of conversion of the density $C_\rho$, and the kinetic viscosities of the fluids $\nu^{\ell,g}$. Specifically, employing 6000 iterations for the initial convergence, followed by 2000 iterations after each particle movement and using the optimal set of parameters, i.e., $\tilde{\gamma} = 80$, $C_\rho = 3 \text{ kg} \cdot \text{m}^{-3}$, $\nu^g = 1.4 \times 10^{-5} \text{ m} \cdot \text{s}^{-2}$, and $\nu^\ell = 1.5 \times 10^{-6} \text{ m} \cdot \text{s}^{-2}$, alongside decoupling time steps between the DEM and LBM, a simulation completes in just 1 hour while maintaining an accuracy level of 5%. Finally, in this study, it must be noted that we considered a capillary bridge volume of $10 \mu\text{L}$. However, in our upcoming triaxial tests on real-sized sand grains, the capillary bridge volume will be around $4 \times 10^{-5} \mu\text{L}$. This smaller volume will necessitate significantly fewer iterations to reach convergence, thereby reducing the computational time required for the simulations (see appendix A).

## 5. Acknowledgement

In memory of our colleague A.P.S. Selvadurai, with whom this work was initiated as part of the Stabdigue project.

The authors would like to express their sincere gratitude to the French National Centre for Space Studies (CNES) for having supported this work.

For the purpose of Open Access, a CC-BY public copyright licence has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission.
A. Appendix

To study the effect of the initial volume on the capillary force value as well as the time computation, three different volumes have been selected 4, 10 and 25 µL. The other physical parameters used are presented in table 6.

<table><caption>Table A.9: LBM physical parameters.</caption><thead><tr><th>Parameters</th><th>Symbols</th><th>Values</th><th>Units</th></tr></thead><tbody><tr><td>Surface tension (water)</td><td>γ</td><td>0.0720</td><td>N.m<sup>-1</sup></td></tr><tr><td>Kinetic viscosity (water)</td><td>ν<sup>l</sup></td><td>1.3 × 10<sup>-6</sup></td><td>m<sup>2</sup>.s<sup>-1</sup></td></tr><tr><td>Kinetic viscosity (air)</td><td>ν<sup>g</sup></td><td>1.8 × 10<sup>-5</sup></td><td>m<sup>2</sup>.s<sup>-1</sup></td></tr><tr><td>Radii</td><td>R</td><td>4</td><td>mm</td></tr><tr><td>Distance</td><td>d</td><td>0.5</td><td>mm</td></tr><tr><td>Temperature</td><td>T</td><td>10</td><td>°C</td></tr><tr><td>Contact angle</td><td>θ</td><td>30</td><td>°</td></tr></tbody></table>

Fig. A.13 illustrates the significant influence of water volume on both the value of the capillary force and the time to convergence. The figure shows that the larger the volume, the larger the value of the capillary force, and the longer the convergence is.
<figure><img src="image_17.png" /><figcaption>Figure A.13: Evolution of the capillary force $F^{cap}$ during the convergence of the droplet to the capillary bridge in terms of the number of iterations for different initial water volumes $V \in \{4, 10, 25\}$ µL.</figcaption></figure>
[1] N. Younes, Z. Benseghier, O. Millet, A. Wautier, F. Nicot, R. Wan, Phase-field lattice boltzmann model for liquid bridges and coalescence in wet granular media, Powder Technology (2022) 117942.

[2] K. Terzaghi, 45th James Forrest Lecture, 1939. Soil mechanics - a new chapter in engineering science, Journal of the Institution of Civil Engineers 12 (7) (1939) 106-142, publisher: ICE Publishing. doi:10.1680/ijoti.1939.14534.

[3] K. Terzaghi, Mechanism of Landslides, in: S. Paige (Ed.), Application of Geology to Engineering Practice, Geological Society of America, 1950, p. 0. doi:10.1130/Berkey.1950.83.

[4] L. R. Daryono, A. Hendratno, M. Nukman, E. Hartantyo, Paper Title for Landslide Identifications Based on Geology and Geoelectricity Mapping Analysis, Study Area: Road Province Pacitan-Ponorogo, East Java, Indonesia., Vol. 2018, European Association of Geoscientists & Engineers, 2018, pp. 1-6. doi:10.3997/2214-4609.201800371.

[5] G. Tang, J. Huang, D. Sheng, S. W. Sloan, Stability analysis of unsaturated soil slopes under random rainfall patterns, Engineering Geology 245 (2018) 322-332. doi:10.1016/j.enggeo.2018.09.013.

[6] A. W. Bishop, The use of the slip circle in the stability analysis of slopes, Geotechnique 5 (1) (1955) 7-17.

[7] A. W. Bishop, N. Morgenstern, Stability coefficients for earth slopes, Geotechnique 10 (4) (1960) 129-153.

[8] X. Yang, S. VANAPALLI, Slope stability analyses of outang landslide based on the peak and residual shear strength behavior, Adv. Eng. Sci. 51 (004) (2019) 55-68.

[9] C.-F. Zhao, N. P. Kruyt, O. Millet, Capillary bridges between unequal-sized spherical particles: Rupture distances and capillary forces, Powder Technology 346 (2019) 462-476.

[10] C. Zhao, N. P. Kruyt, M. Pouragha, R. Wan, Fabric response to stress probing in granular materials: Two-dimensional, anisotropic systems
465

466

467

146 104695. doi:10.1016/j.compgeo.2022.104695.

URL https://www.sciencedirect.com/science/article/pii/
S0266352X22000593

468

469

470

[11] G. Ma, W. Zhou, X.-L. Chang, Modeling the particle breakage of rockfill materials with the cohesive crack model, Computers and Geotechnics 61 (2014) 132-143. doi:10.1016/j.compgeo.2014.05.006.

471

472

473

474

[12] G. Ma, W. Zhou, R. A. Regueiro, Q. Wang, X. Chang, Modeling the fragmentation of rock grains using computed tomography and combined FDEM, Powder Technology 308 (2017) 388-397. doi:10.1016/j.powtec.2016.11.046.

475

476

477

478

[13] Z.-Y. Yin, P.-Y. Hicher, C. Dano, Y.-F. Jin, Modeling Mechanical Behavior of Very Coarse Granular Materials, Journal of Engineering Mechanics 143 (1) (2017) C4016006, publisher: American Society of Civil Engineers. doi:10.1061/(ASCE)EM.1943-7889.0001059.

479

480

481

482

[14] W. Zhou, L. Yang, G. Ma, X. Chang, Y. Cheng, D. Li, Macro-micro responses of crushable granular materials in simulated true triaxial tests, Granular Matter 17 (4) (2015) 497-509. doi:10.1007/s10035-015-0571-3.

483

484

485

486

[15] H. Zhu, F. Nicot, F. Darve, Meso-structure organization in two-dimensional granular materials along biaxial loading path, International Journal of Solids and Structures 96 (2016) 25-37. doi:10.1016/j.ijsolstr.2016.06.025.

487

488

489

490

[16] H. Zhu, H. N. G. Nguyen, F. Nicot, F. Darve, On a common critical state in localized and diffuse failure modes, Journal of the Mechanics and Physics of Solids 95 (2016) 112-131. doi:10.1016/j.jmps.2016.05.026.

491

492

493

[17] A. Tordesillas, Force chain buckling, unjamming transitions and shear banding in dense granular assemblies, Philosophical Magazine 87 (32) (2007) 4987-5016. doi:10.1080/14786430701594848.

494

495

496

[18] A. Tordesillas, D. M. Walker, Q. Lin, Force cycles and force chains, Physical Review E 81 (1) (2010) 011302, publisher: American Physical Society. doi:10.1103/PhysRevE.81.011302.
497

498

499

[19] A. Wautier, Analyse micro-inertielle des instabilités mécaniques dans les milieux granulaires, application à l'érosion interne, These de doctorat, Aix-Marseille (Sep. 2018).

500

501

502

[20] F. Nicot, F. Darve, The H-microdirectional model: Accounting for a mesoscopic scale, Mechanics of Materials 43 (12) (2011) 918-929. doi: 10.1016/j.mechmat.2011.07.006.

503

504

505

[21] J. Duriez, R. Wan, Stress in wet granular media with interfaces via homogenization and discrete element approaches, Journal of Engineering Mechanics-ASCE 142 (12) (2016).

506

507

[22] J. Duriez, R. Wan, Contact angle mechanical influence in wet granular soils, Acta Geotechnica 12 (1) (2017) 67-83.

508

509

510

[23] J. Duriez, M. Eghbalian, R. Wan, F. Darve, The micromechanical nature of stresses in triphasic granular media with interfaces, Journal of the Mechanics and Physics of Solids 99 (2017) 495-511.

511

512

513

[24] J. Duriez, R. Wan, A micromechanical $\mu$unsat effective stress expression for stress-strain behaviour of wet granular materials, Geomechanics for Energy and the Environment 15 (2018) 10-18.

514

515

516

517

[25] L. Scholtès, P.-Y. Hicher, F. Nicot, B. Chareyre, F. Darve, On the capillary stress tensor in wet granular materials, International journal for numerical and analytical methods in geomechanics 33 (10) (2009) 1289-1313.

518

519

520

[26] L. Scholtes, B. Chareyre, F. Nicot, F. Darve, Micromechanics of granular materials with capillary effects, International journal of engineering science 47 (1) (2009) 64-75.

521

522

523

[27] L. Scholtès, Modélisation micromécanique des milieux granulaires partiellement saturés, Ph.D. thesis, Institut National Polytechnique de Grenoble-INPG (2008).

524

525

[28] K. A. Brakke, The surface evolver, Experimental mathematics 1 (2) (1992) 141-165.

526

527

528

[29] M. Miot, G. Veylon, A. Wautier, P. Philippe, F. Nicot, F. Jamin, Numerical analysis of capillary bridges and coalescence in a triplet of spheres, Granular Matter 23 (3) (2021) 1-18.
529

530

531

[30] A. Di Renzo, G. Picarelli, F. P. Di Maio, Numerical investigation of funicular liquid bridge interactions between spherical particles, Chemical Engineering & Technology 43 (5) (2020) 830-837.

532

533

534

[31] N. Younes, A. Wautier, R. Wan, O. Millet, F. Nicot, R. Bouchard, DEM-lbm coupling for partially saturated granular assemblies, Computers and Geotechnics 162 (2023) 105677.

535

536

537

[32] Z. Benseghier, O. Millet, P. Philippe, A. Wautier, N. Younes, E. Liberge, Relevance of capillary interfaces simulation with the shan-chen multi-phase lb model, Granular Matter 24 (3) (2022) 1-18.

538

539

[33] Y. Sun, C. Beckermann, Sharp interface tracking using the phase-field equation, Journal of Computational Physics 220 (2) (2007) 626-653.

540

541

542

[34] P.-H. Chiu, Y.-T. Lin, A conservative phase field method for solving incompressible two-phase flows, Journal of Computational Physics 230 (1) (2011) 185-204.

543

544

545

[35] H. Liang, J. Xu, J. Chen, H. Wang, Z. Chai, B. Shi, Phase-field-based lattice boltzmann modeling of large-density-ratio two-phase flows, Phys. Rev. E 97 (2018) 033309.

546

547

548

[36] H. Liang, H. Liu, Z. Chai, B. Shi, Lattice boltzmann method for contact-line motion of binary fluids with high density ratio, Phys. Rev. E 99 (2019) 063306.

549

550

551

552

[37] A. Fakhari, D. Bolster, Diffuse interface modeling of three-phase contact line dynamics on curved boundaries: A lattice boltzmann model for large density and viscosity ratios, Journal of Computational Physics 334 (2017) 620-638.

553

554

555

[38] O. Penrose, P. C. Fife, Thermodynamically consistent models of phase-field type for the kinetic of phase transitions, Physica D: Nonlinear Phenomena 43 (1) (1990) 44-62.

556

557

558

[39] T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, E. M. Viggen, The lattice boltzmann method, Springer International Publishing 10 (978-3) (2017) 4-15.
[40] B. Mielniczuk, O. Millet, G. Gagneux, M. S. El Youssoufi, Characterisation of pendular capillary bridges derived from experimental data using inverse problem method, *Granular Matter* 20 (1) (2018) 14.

[41] B. Mielniczuk, T. Hueckel, M. S. El Youssoufi, Laplace pressure evolution and four instabilities in evaporating two-grain liquid bridges, *Powder Technology* 283 (2015) 137–151.