# Relevance of capillary interfaces simulation with the Shan-Chen multiphase LB model

Z. Benseghier, O. Millet, P. Philippe, A. Wautier, N. Younes, Erwan Liberge

► To cite this version:

Z. Benseghier, O. Millet, P. Philippe, A. Wautier, N. Younes, et al. Relevance of capillary interfaces simulation with the Shan-Chen multiphase LB model. Granular Matter, 2022, 24 (3), pp.82. 10.1007/s10035-022-01243-5. http://hal-03837948

HAL Id: hal-03837948

https://hal.science/hal-03837948

Submitted on 5 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Relevance of capillary interfaces simulation with the Shan-Chen multiphase LB model

Z. Benseghier · O. Millet · P. Philippe · A. Wautier · N. Younes · E. Liberge.

Received: date / Accepted: date

**Abstract** This paper presents a numerical study of capillary interfaces using the Single Component Multi-phase Shan-Chen model, which is based on the lattice Boltzmann method. Despite the simplicity of the model, it has been shown to be effective and the present study aims to test its ability to correctly reproduce the physics of multiphase systems. To this end, several benchmark simulations were carried out in the configurations of a drop on a flat wall and then on a spherical surface to characterize the wetting behavior and relate explicitly the contact angle to model parameters. In addition, the capillary forces induced by a liquid bridge between two spherical particles were numerically calculated. We show that the results obtained are in agreement with experimental and theoretical results from the literature. The model is thus accurate in addressing the wetting behavior and capillary interfaces in unsaturated granular soils despite the fact that surface tension and contact angles are not explicit parameters of the model. To this respect, explicit relationships with Shan-Chen parameters are provided.

**Keywords** capillary · contact angle · multi-phase · Lattice Boltzmann Method.

Z. Benseghier

LaSIE, UMR CNRS 7356, University of La Rochelle, 17042, La Rochelle Cedex 1, France.

E-mail: zeyd.benseghier@univ-lr.fr

O. Millet

LaSIE, UMR CNRS 7356, University of La Rochelle, 17042, La Rochelle Cedex 1, France.

E-mail: olivier.millet@univ-lr.fr

P. Philippe

INRAE, Aix Marseille Univ, RECOVER, 3275 route de Cézanne, 13100 Aix-en-Provence, France.

E-mail: pierre.philippe@inrae.fr

A. Wautier

INRAE, Aix Marseille Univ, RECOVER, 3275 route de Cézanne, 13100 Aix-en-Provence, France.

E-mail: antoine.wautier@inrae.fr

N. Younes

LaSIE, UMR CNRS 7356, University of La Rochelle, 17042, La Rochelle Cedex 1, France.

E-mail: nabil.younes@univ-lr.fr

E. Liberge

LaSIE, UMR CNRS 7356, University of La Rochelle, 17042, La Rochelle Cedex 1, France.

E-mail: erwan.liberge@univ-lr.fr
21

# 1 Introduction

22

23

24

25

26

27

28

29

30

31

32

Capillary effects play a fundamental role in the behavior of wet granular soils. When soil particles are subjected to humidity, capillary bridges occur and the induced adhesion forces lead to cohesion at the soil scale. The shape of such bridges is well described by Young-Laplace equation and ruled by the surface tension of the interfaces at stake. Consequently, a pressure difference will exist inside the capillary bridge, which in turn results in capillary forces between grains. In case capillary bridges merge or collapse (wetting or drying of the material), a progressive loss in cohesion may occur on a larger scale, such as in hydraulic earthworks (e.g. earth-dams, levees, and dikes), possibly with dramatic consequences on their mechanical stability. Therefore, it is crucial to first quantify accurately the capillary forces between solid grains in unsaturated media and next to perform reliable numerical computations at the scale of earthen structures.

33

34

35

36

37

38

39

40

41

Several numerical, experimental, and analytical studies have been carried out to quantify these forces. Mielniczuk et al. [19] have quantified capillary forces using an inverse problem technique for solving Young-Laplace equation when the capillary pressure is unknown. Their method was compared to experimental results and showed great similarities [19]. Recently, this approach has been extended to polydisperse particles, to capillary bridges between a sphere and a plane, or between two parallel planes [22, 23, 24, 25, 26, 27]. Duriez & Wan [8] have solved numerically the Young-Laplace equation in order to compute the attractive force between two grains when the capillary pressure is controlled and validated their strategy by comparing with experimental results.

42

43

44

45

46

47

48

49

50

51

52

Another way is to model all the phases before computing the capillary bridge directly from the simulation. Many works deal with gas-liquid-systems [6, 35], but only a few have focused on the calculation of the forces induced by capillary bridges. Miot et al. [20] used an energy minimization approach (with the open-source code Surface Evolver [3]) combined with the use of the virtual work principle to compute the shape of the capillary bridge and the corresponding forces transmitted to the grains. Even if the method can be applied to complex geometry and dynamic evolutions, it is limited to steady problems and the method relies on the use of well-chosen parameters to achieve the convergence of the results. Sun et al. [37] proposed a VOF-IB-DNS method, which combines Volume Of Fluid (VOF) to calculate the shape of the gas-liquid interface, Immersed Boundary (IB) to implement solid particles, and Finite Volume method to solve directly the Navier-Stokes (N-S) equations.

53

54

55

56

57

58

59

60

61

62

In this work, we propose to use the Lattice Boltzmann Method (LBM) to model wetting condition (contact angle) and capillary bridges between particles following the numerous advantages of this method in fluid dynamics, mainly in dealing with complex geometry. The LBM has been successfully developed for computational fluid dynamics since the 90's [2] and appears to be a suitable alternative numerical method. Based on the Boltzmann equation, the LBM considers the transport of the probability to find a particle according to time, space and velocity. Then, the macroscopic variables are obtained using momenta of the adapted distribution functions. The power of the LBM also lies in its simplicity and easy implementation on Graphic Processor Units (GPU) that considerably reduces the computation time [10].

63

64

65

66

67

There exists a considerable body of literature on modeling interface-tracking models with LBM, such as Cahn-Hilliard [38, 42, 44] or Allen-Cahn [9, 18] equations. These models have proved to be accurate when dealing with high-density ratios by solving, at the same time, the Navier-Stokes (N-S) equations as well as one of the interface-tracking equations indicated above. However, a primary problem of this strategy is the large computation and memory
costs, since two equations need to be solved simultaneously. Therefore, accurate computations of capillary forces at the scale of a REV seem out of reach with these approaches.

A possible solution to address this issue would be to simulate capillary bridges and compute the corresponding capillary forces based on a less accurate but more efficient model, especially with the Shan-Chen model that accounts for multi-phase flows without tracking the interface. This model relies on molecular interactions (pseudo-potential) between phases, which leads to phase separation, and can be divided into two categories. The first approach is Multi-Component Multi-Phase (MCMP), where two LB equations need to be solved when the different fluids are made of different materials (two immiscible liquids for instance). However, the maximum density ratio is of the order of $O(1)$, so that the method cannot handle water/air interfaces in particular [15]. This approach was used for partially saturated media [21, 28] and for studying the capillary forces between two particles [40]. The second family of models is Single-Component Multi-Phase (SCMP), where only one LB equation is required to solve N-S equations and a supplementary Equation Of State (EOS) is introduced to simulate the coexistence between liquid and gaseous phase for the same substance (for instance water liquid and water vapor). Some applications concern flows in unsaturated porous materials [7, 14, 30, 36]. The density ratio is controlled by the EOS and the maximum ratio that can be achieved is in the order $O(10^2)$ [15]. The latter is indeed more efficient and will be described in more detail and implemented in the remainder of this paper. While using this model for water, one has to keep in mind that the maximum density ratio of $O(10^2)$ imposes to work at high temperature and high pressure.

The overall aim of this paper is to study the computational efficiency as well as the accuracy of the SCMP model in dealing with contact angles for flat and convex solid frontiers and in computing capillary forces. The plan of the paper is organized as follows. In the first section, a brief introduction of LBM and the SCMP model is given in addition to fluid-solid interaction with wetting conditions (contact angle). Afterward, several validation benchmarks are performed to validate the model, including stationary spherical drop to verify the Laplace equation and deduce the surface tension value, and drop on a flat and spherical particle, to validate the contact angle implementation. Next, a section is devoted to capillary bridges between two spherical grains and induced capillary forces with a comparison between our results and experimental and theoretical ones from the literature [19, 29, 30]. Finally, the last section provides a general conclusion opening up future prospects.

# 2 Implementation of LBM for multi-phase flows and wetting condition

## 2.1 Standard Lattice Boltzmann method

The Lattice Boltzmann method with Bhatnagar-Gross-Krook (BGK) collision operator and a source term $F_\alpha$ can be written as follows:

$$f_{\alpha}(\mathbf{x} + \mathbf{c}_{\alpha}\Delta t, t + \Delta t) = f_{\alpha}(\mathbf{x}, t) - \frac{1}{\tau} [f_{\alpha}(\mathbf{x}, t) - f_{\alpha}^{eq}(\mathbf{x}, t)] + F_{\alpha}(\mathbf{x}, t) \quad (1)$$

where $f_\alpha$ is the particle distribution function at a given lattice position $\mathbf{x}$, time $t$, and along a direction $\alpha$, restricted to a limited number of discrete velocities. The equilibrium distribution function $f_\alpha^{eq}$ is given by

$$f_{\alpha}^{eq}(\rho, \mathbf{u}) = \rho w_{\alpha} \left[ 1 + \frac{\mathbf{u} \cdot \mathbf{c}_{\alpha}}{c_s^2} + \frac{(\mathbf{u} \cdot \mathbf{c}_{\alpha})^2}{2c_s^4} - \frac{\mathbf{u} \cdot \mathbf{u}}{2c_s^2} \right] \quad (2)$$
where $c_s = c/\sqrt{3}$ is the lattice speed of sound with $c = \frac{\Delta x}{\Delta t}$ the lattice speed. $\Delta x$ and $\Delta t$ denote the lattice size and time steps, in general chosen equal to one in lattice units ($\Delta x = \Delta t = 1$ lu). The D3Q19 scheme (three dimensional configuration with nineteen velocities) is widely used in the literature for 3D simulation and implemented in this study, whose discrete velocities $c_\alpha$ read

$$ c_{\alpha} = \begin{cases} (0, 0, 0) & \alpha = 0 \\ (\pm 1, 0, 0), (0, \pm 1, 0), (0, 0, \pm 1), & \alpha = 1, 2, \dots, 6 \\ (\pm 1, \pm 1, 0), (\pm 1, 0, \pm 1), (0, \pm 1, \pm 1) & \alpha = 7, 8, \dots, 18 \end{cases} \quad (3) $$

with the corresponding weights $w_0 = 1/3$, $w_{1-6} = 1/18$, and $w_{7-18} = 1/36$.

The source term $F_\alpha$ in Eq. (1) can include external body force, fluid-fluid interaction force $F_{int}$, and adhesive force $F_{ads}$ (i.e fluid-solid interaction force) which will be introduced in Sec. 2.3 <fnref n="1" />.

The density $\rho$ and velocity $u$ of the fluid, in lattice units, can be directly derived from the distribution functions as:

$$ \rho = \sum_{\alpha} f_{\alpha} \quad (4) $$

$$ u = \frac{1}{\rho} \sum_{\alpha} f_{\alpha} c_{\alpha}. \quad (5) $$

The relation between the relaxation time $\tau$ and the kinematic viscosity $\nu$ of the fluid classically reads <fnref n="2" />:

$$ \tau = \frac{\nu}{\Delta t c_s^2} + \frac{1}{2}. \quad (6) $$

## 2.2 SCMP Shan-Chen model

For single-component multi-phase system (SCMP) (e.g water and its vapor), we use the model proposed by Shan and Chen [32, 33], in which a fluid interaction force is introduced to induce phase separation. This interaction term reads

$$ F_{\text{int}}(\mathbf{x}) = -G\psi(\mathbf{x}) \sum_{\alpha} w_{\alpha} \psi(\mathbf{x} + \mathbf{c}_{\alpha} \Delta t) \mathbf{c}_{\alpha} \quad (7) $$

where $\psi(\mathbf{x})$ is the effective mass, which depends on the local density $\rho(\mathbf{x})$, and $G$ is the interaction strength. $w_\alpha$ are the lattice weights as previously defined. Note that the sum runs over all neighboring fluid nodes.

The effective mass (also called the pseudopotential) proposed in the original work of Shan and Chen [32, 33] takes the following form:

$$ \psi(\mathbf{x}) = \rho_0 \left[ 1 - \exp \left( -\frac{\rho}{\rho_0} \right) \right] \quad (8) $$

$$ F_{\alpha} = w_{\alpha} \left( \frac{c_{\alpha} - u}{c_{s}^{2}} + \frac{(c_{\alpha} - u) c_{\alpha}}{c_{s}^{4}} \right) \cdot F. $$

<footnote n="1">¹ Note that $F_\alpha$ can be linked to $F$ for instance following Guo et al. [12] through:</footnote>

<footnote n="2">² To recover Navier-Stokes equation with BGK collision operator.</footnote>
129

130

131

132

133

134

135

136

where $\rho_0$ is a normalization constant which is usually chosen as 1 in lattice units. It was found however that using this form does not satisfy the thermodynamic consistency and may lead to high spurious currents [5, 41]. Therefore the expression of the effective density given in (8) can only be used to describe mixtures of liquid and gas with ratio of a maximum density of the order $O(10)$ [5].

In order to increase accessible density ratios and reduce spurious currents, we chose to implement a different Equation Of State (EOS), namely the Carnahan–Starling (C–S) one using the method developed in [41], where the effective mass reads

$$ \psi(\rho) = \sqrt{\frac{2(P - \rho c_s^2)}{G c_s^2}}. \qquad (9) $$

137

138

The pressure $P$ applied on a given node is expressed according to the following Equation Of State (EOS) <fnref n="3" />

$$ P = \rho R T \frac{1 + b \rho / 4 + (b \rho / 4)^2 - (b \rho / 4)^3}{(1 - b \rho / 4)^3} - a \rho^2 \qquad (10) $$

139

140

141

142

143

144

where $a = 0.4963(RT_c)^2/P_c$, $b = 0.1873RT_c/P_c$, $R$ is the perfect gas constant and $T$ is the temperature, all expressed in lattice units. $T_c$ and $p_c$ denote the critical temperature and pressure, respectively (also expressed in lattice units). The expressions given for $a$ and $b$ consistently satisfy the conditions $(\partial P / \partial \rho)_T = (\partial^2 P / \partial \rho^2)_T = 0$ at the critical point. Note that the parameter $G$ no longer controls the interaction strength (as in Eq. (7)) but rather ensures that the term inside the square root in Eq. (9) is positive.

145

146

147

148

149

150

According to Yuan and Schaefer [41], these parameters are fixed at $a = 1$ lu, $b = 4$ lu and $R = 1$ lu. With these values, the interface thickness is approximately 5 lu for most $T$ and $\tau$ values, and the simulation remains stable [15]. For water, the corresponding critical properties in lattice units are $T_c = 0.0943$ lu, $P_c = 0.00442$ lu, and $\rho_c = 0.1136$ lu [41]. The coexistence densities for C-S EOS can be deduced from the results of [41] and are shown in Figure 1.

151

152

153

154

155

In our simulations, the temperature is fixed to $T = 0.75T_c$, with the corresponding gas and liquid densities $\rho_g \approx 0.015$ lu and $\rho_l \approx 0.33$ lu, which are deduced from Fig. 1. In this work, it was found that by decreasing the temperature, the simulation becomes unstable. However, according to [13, 15], in this temperature range (for $T \ge 0.7T_c$), the thermodynamic consistency is satisfied with the C-S EOS.

156

157

158

159

160

161

162

163

164

165

166

The conversion from lattice units to physical units of the fluid properties can be done using the reduced properties concept [4, 41] based on the equivalence of some dimensional numbers between lattice and physical units (see the Appendix). Thus, the corresponding physical properties are temperature $T=212.17$ °C = 485.32 °K, water liquid density $\rho_l = 814.62$ kg/m<fnref n="3" />, and water vapor density<fnref n="4" /> $\rho_g = 37.02$ kg/m<fnref n="3" />. It is important to quote that this choice of physical parameters $T$, $P$, and $\rho$ (similar to other studies found in the literature) does not correspond to those used in experiments performed for comparison, but is mandatory to ensure the stability of the LBM simulations. However, as we will see in the application examples, this choice of parameters enables to recover the experimental values obtained for multi-phase systems, that proves its relevance which is not justified in literature to our knowledge.

<footnote n="3">³ Note that for an "ideal" or "perfect" gas laws, the pressure is $P(\rho) = \rho c_s^2$</footnote>

<footnote n="4">⁴ The physical thermodynamic properties for water at the critical point are critical temperature $T_c=373.946$ °C = 647.096 °K, critical pressure $P_c=217.7$ atm = 220.6 bar = 22.06 MPa, and critical density $\rho_c=322$ kg/m³.</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1: Simulation (solid circles) and theoretical coexistence curves (solid line) for C-S EOS. The results and the figure are extracted from [41].</figcaption></figure>

2.3 Fluid-solid interaction and contact angle

Modeling capillary bridges requires the introduction of a solid phase in addition to the two fluid phases presented in the previous section. The adhesion force between the fluids and the solid nodes takes a similar form as in Eq. (7) and reads according to [1]

$$F_{ads}(\mathbf{x}) = -G\psi(\mathbf{x}) \sum_{\alpha} w_{\alpha}\psi(\rho_{w})s(\mathbf{x} + \mathbf{c}_{\alpha}\Delta t)\mathbf{c}_{\alpha} \quad (11)$$

where s is the indicator function for the solid phase, equal to 1 for solid nodes and 0 for fluid nodes. The effective density, that is different from the real physical density, is denoted $\rho_w$. The total force acting on a fluid node $\mathbf{x}$ can be now written as:

$$F(\mathbf{x}) = F_{int}(\mathbf{x}) + F_{ads}(\mathbf{x}) \\ = -\psi(\mathbf{x}) \left( -G \sum_{\alpha}^{\text{fluid}} w_{\alpha} \psi(\mathbf{x} + \mathbf{c}_{\alpha} \Delta t) \mathbf{c}_{\alpha} + G \sum_{\alpha}^{\text{solid}} w_{\alpha} \psi(\rho_w) \mathbf{c}_{\alpha} \right). \quad (12)$$

To evaluate the effective density $\rho_w$, we use the method proposed by Li et al. [17], which is an improved version of the virtual-density scheme. $\rho_w$ is a local quantity, instead of a constant value in the whole solid domain, defined as

$$\rho_w = \begin{cases} \varphi\rho_{ave}(\mathbf{x}), & \varphi \ge 1, \text{ For } \theta \le 90^{\circ} \\ \rho_{ave}(\mathbf{x}) - \Delta\rho & \Delta\rho \ge 0, \text{ For } \theta \ge 90^{\circ}. \end{cases} \quad (13)$$

where $\varphi$ and $\Delta\rho$ are constants to be tuned to achieve different contact angles $\theta$. Moreover, the average fluid density $\rho_{ave}$ reads
$$ \rho_{ave}(\mathbf{x}) = \frac{\sum_{\alpha} w_{\alpha} \rho(\mathbf{x} + \mathbf{c}_{\alpha} \Delta t)(1 - s(\mathbf{x} + \mathbf{c}_{\alpha} \Delta t))}{\sum_{\alpha} w_{\alpha}(1 - s(\mathbf{x} + \mathbf{c}_{\alpha} \Delta t))}. \quad (14) $$

Note that in the definition given by (13), $\rho_w(\mathbf{x})$ is restricted to the range $\rho_g \le \rho_w(\mathbf{x}) \le \rho_l$.

To summarize, once the parameter $\varphi$ or $\Delta\rho$ is set and $\rho_w$ calculated from Eq. (13), the adhesive force between solid and fluid is evaluated by Eq. (11). As a consequence, a contact angle is implicitly obtained at equilibrium. It is important to notice that with this approach, the contact angle is a result of the LBM simulation at equilibrium, but is not a straightforward input data. As a result, the implicit relationship between the model parameters and the corresponding contact angle needs to be established to allow for an effective use of the model, which is one of the aims of the present study.

## 2.4 Numerical scheme

To specify the total force acting on the fluid $F$, we chose to implement the velocity shift force scheme [32, 33]. In this scheme, the source term in Eq. (1) is set to 0 ($F_\alpha = 0$) and the force is incorporated in the velocity. More precisely, when calculating the equilibrium distribution function $f^{eq}$ in Eq. (2), the following equilibrium velocity $u_{eq}$ is used instead of $u$ in Eq. (5):

$$ u_{eq} = u + \frac{\tau}{\rho} F \Delta t. \quad (15) $$

where $F(\mathbf{x})$ is given by Eq. (12).

The actual physical fluid velocity, which will be used in section 4 for stress calculation, is defined as [33]:

$$ u_p = u + \frac{1}{2\rho} F \Delta t. \quad (16) $$

The algorithm adopted to solve the Lattice Boltzmann equation with the BGK collision model Eq. (1) is presented in Fig. 2. The collision step is expressed as follows

$$ f_{\alpha}^{*}(\mathbf{x}, t) = f_{\alpha}(\mathbf{x}, t) - \frac{1}{\tau} [f_{\alpha}(\mathbf{x}, t) - f_{\alpha}^{eq}(\mathbf{x}, t)] \quad (17) $$

whereas the streaming step reads

$$ f_{\alpha}(\mathbf{x} + \mathbf{c}_{\alpha} \Delta t, t + \Delta t) = f_{\alpha}^{*}(\mathbf{x}, t). \quad (18) $$

The non-slip boundary condition at the solid frontier relies on the halfway bounce-back scheme [43], in which the unknown distribution functions at the boundary fluid node $\mathbf{x}_b$, next to the solid, are replaced by the ones in the opposite directions

$$ f_{\tilde{\alpha}}(\mathbf{x}_b, t + \Delta t) = f_{\alpha}^{*}(\mathbf{x}_b, t) \quad (19) $$

where $\tilde{\alpha}$ stands for the opposite direction of $\alpha$ and $f_{\alpha}^*$ is the post collision distribution function. Figure 3 shows an illustration of the position of fluid, fluid boundary, and solid nodes.
<figure><img src="image_3.png" /><figcaption>Fig. 2: Algorithm for LBM resolution with GPU implementation. Kernels are GPU functions that are executed in parallel by different GPU threads.</figcaption></figure>

## 3 Benchmark simulations

In this section three benchmark configurations are considered: a drop of water i) with no contacts, ii) lying on a solid plane or iii) lying on a solid sphere. These benchmark simulations will be used to derive the explicit relationship between the model parameters and the corresponding surface tension and contact angle.

### 3.1 Surface tension deduced from the simulation of a stationary spherical drop

In order to bridge the gap between the model parameters describing the shape of the solid gas interface and more physically relevant quantities, as surface tension and contact angle, a series of numerical simulations are performed to compute the equilibrium of an initially spherical water drop of radius $r_0$. In this simple configuration, the Laplace equation must be satisfied. For a spherical drop of radius $r$, this equation reads

$$ \Delta P = P_{in} - P_{out} = \frac{2\gamma}{r} \qquad (20) $$
<figure><img src="image_4.png" /><figcaption>**Fig. 3:** Representation of fluid, boundary fluid, and solid nodes.</figcaption></figure>

where $P_{in}$ and $P_{out}$ are the pressure inside and outside the drop far away from the interface, $\gamma$ is the surface tension, and $r$ is the droplet radius at equilibrium. More precisely, $r$ is calculated as the distance between the center of the droplet and the middle position of the interface: $r = (r_{in}+r_{out})/2$, where $r_{in}$ and $r_{out}$ are the inner and outer radius of the fluid/gas interface. Far from the interface, $\rho = \rho_l$ in the fluid and $\rho = \rho_g$ in the gas. Note that in LBM multi-phase simulations, the interface has a thickness that may not be neglected, and depends on the mesh resolution. Moreover, with the approach used here, the interface thickness cannot be imposed but it is a result of the thermodynamic equilibrium that actually follows the EOS (10).

A series of simulations was conducted by varying the initial drop radius $r_0$. The drop is placed in the middle of the domain of size $N_x \times N_y \times N_z = 256 \times 256 \times 256$, with $N_x, N_y,$ and $N_z$ the domain size in x, y, and z direction, respectively. Periodic boundary conditions are applied in all directions. The temperature ratio is fixed at $T/T_c = 0.75$, the relaxation time $\tau = 1$ lu, and gas and liquid densities are fixed at $\rho_g = 0.015$ lu and $\rho_l = 0.33$ lu, respectively.

Figure 4b shows pressure distribution along the x-axis for an initial drop radius $r_0 = 30$ lu. It can be seen that the pressure is almost constant inside and outside the drop. There is a positive and negative tips inside the interface, which is due to the strong non-linearity of C-S EOS (Eq. 10).

Figure 5 shows the pressure drop $\Delta P$ versus $2/r$ for two different resolutions. The fact that the relationship is linear is fully consistent with the existence of a surface tension corresponding to the slope of the curve. The comparison between the two mesh resolutions highlights that $\gamma$ is not mesh dependent. The slopes represent indeed the surface tension with a value $\gamma = 0.0111$ lu ($R^2 = 0.99968$) for the resolution $256 \times 256 \times 256$ and $\gamma = 0.0114$ lu ($R^2 = 0.99953$) for the resolution $384 \times 384 \times 384$. Note that a very similar value of 0.01025 was reported in a previous work using the same model parameters [13].
<figure><img src="image_5.png" /><figcaption>**Fig. 4:** (a) Density contours at equilibrium state for initial drop radius $r_0$ = 30 lu. (b) Pressure distribution along the horizontal axis that passes through the center of the domain.</figcaption></figure>

## 3.2 Contact angle of a drop on a flat solid surface

### 3.2.1 Method for contact angle evaluation

Here, the simulation of a liquid drop spread on a flat surface is addressed. We study in particular the influence of the model parameter $\varphi$ on the contact angle $\theta$ without gravity. As a starting configuration of the simulation, the domain is initialized with the gas density $\rho_g$, except a drop region of radius $r_0$ (see Fig. 7) where prevails the liquid density $\rho_l$. The center of the drop is initially located at $(N_x/2, r_0, N_z/2)$. Periodic boundary conditions (BC) are imposed in both horizontal and front-back directions, whereas wall BC are applied at the top and at the bottom of the domain. Note that for the wall, we implemented the improved
<figure><img src="image_6.png" /><figcaption>Fig. 5: Pressure drop between inside and outside the droplet versus 2/r. The solid and dashed lines represent a linear fit with zero intercept consistent with the Laplace equation (20).</figcaption></figure>

virtual density scheme (see Sec. 2.3), which ensures the wetting condition implicitly (contact angle). Then, the simulation is launched until a final equilibrium situation is reached.

The contact angle is evaluated geometrically using Eq. (21) given below, with the assumption that the shape of the drop is spherical, which is the case in the absence of gravity:

$$ \theta = \begin{cases} 2 \tan^{-1}(\frac{2h}{b}) & \text{For } \theta \le 90^{\circ} \\ \arcsin(\frac{h-r}{r} + \frac{\pi}{2}) & \text{For } \theta > 90^{\circ}. \end{cases} \qquad (21) $$

Figure 6 shows the geometric parameters, namely h, r, and b, required to perform the contact angle calculation in Equation (21).

Figure 7 illustrates a drop on a plane surface with the boundary conditions used for the simulation and the shape of the drop at initial state and final state, once equilibrium is reached.

### 3.2.2 Sensitivity to the mesh resolution

We first analyze the mesh sensitivity on the value of the contact angle. For all meshes, the ratio $r_0/(N_x - 1)$ is kept fixed, where $N_x$ is the number of lattices nodes in x-direction ($r_0$ is thus constant in physical units).

The values of parameter $\varphi$ and contact angle $\theta$ are shown in Table 1. By comparison to a reference value chosen as the one given by the finest resolution available, we first observe that the relative error decreases as the mesh resolution increases. It is to note also that the error tends to roughly decrease with $\varphi$.

In the following, all the LBM simulations will be performed with the second mesh resolution, since the latter gives accurate results within reasonable calculation time.
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 6:** Definitions of several geometric parameters needed for the contact angle $\theta$ calculation for: (a) $\theta < 90^{\circ}$ and (b) $\theta > 90^{\circ}$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 7:** Diagram illustrating a drop on a flat surface and the boundary conditions used for the simulation.</figcaption></figure>
<table><caption>**Table 1:** Values obtained for the contact angle θ versus φ for different mesh resolutions. The errors are calculated based on the finest resolution values.</caption><thead><tr><th>Resolution:<br>(N<sub>x</sub> × N<sub>y</sub> × N<sub>z</sub>)</th><th colspan="2">128 × 128 × 128</th><th colspan="2">256 × 128 × 256</th><th colspan="2">384 × 128 × 384</th></tr><tr><th>r<sub>0</sub> (lu)</th><th colspan="2">10</th><th colspan="2">20</th><th colspan="2">30</th></tr><tr><th>φ</th><th>θ (°)</th><th>error (%)</th><th>θ (°)</th><th>error (%)</th><th>θ (°)</th><th>error (%)</th></tr></thead><tbody><tr><td>0.8</td><td>91.87</td><td>0.66</td><td>91.11</td><td>0.17</td><td>91.27</td><td>-</td></tr><tr><td>1.0</td><td>72.43</td><td>1.98</td><td>71.22</td><td>0.27</td><td>71.03</td><td>-</td></tr><tr><td>1.2</td><td>60.44</td><td>1.85</td><td>59.57</td><td>0.38</td><td>59.34</td><td>-</td></tr><tr><td>1.4</td><td>50.17</td><td>2.87</td><td>50.92</td><td>1.41</td><td>51.65</td><td>-</td></tr><tr><td>1.6</td><td>42.95</td><td>1.76</td><td>43.51</td><td>0.49</td><td>43.72</td><td>-</td></tr><tr><td>1.8</td><td>35.15</td><td>4.2</td><td>36.29</td><td>1.08</td><td>36.69</td><td>-</td></tr><tr><td>2.0</td><td>26.33</td><td>8.63</td><td>27.82</td><td>3.48</td><td>28.82</td><td>-</td></tr></tbody></table>
### 3.2.3 Influence of the initial drop volume on the contact angle

In line with the previous finding on the influence of the parameter $\varphi$ on the contact angle, a parametric study was performed by varying now the initial drop radius $r_0$ (corresponding to the initial drop volume) to verify that the observed contact angle (which is a response parameter of the LBM scheme) is not problem dependent.

Figure 8 shows the evolution of the contact angle $\theta$ versus $\varphi$ for different initial drop radius $r_0$. It can be concluded that, with the virtual density scheme implemented here, the initial drop volume does not influence significantly the contact angle in the absence of gravity. This is consistent with the fact that the contact angle is indeed a material property. This implicit definition of the contact angle in the numerical scheme is consequently physically sounded, at least when considering a drop on a plane.

<figure><img src="image_10.png" /><figcaption>**Fig. 8:** Values obtained for the contact angle $\theta$ versus the parameter $\varphi$ for different initial drop radii $r_0$ for a drop on a flat surface.</figcaption></figure>

## 3.3 Contact angle of a drop on a convex solid surface

### 3.3.1 Method for contact angle and drop volume evaluation

In this second configuration, we now consider a liquid drop on a convex solid surface, namely a spherical cap. Yet again, the drop is placed exactly at the top of the solid surface as sketched in Figure 9(a). Here, the objective is to confirm that the contact angle does not depend on the curvature of the solid surface. Periodic BC are applied in both horizontal and front-back directions as well as top-bottom limits. We apply the improved virtual density scheme at the solid boundary to ensure the wetting condition (contact angle) while the bounce-back scheme is implemented for the probability distribution function $f$.
To evaluate the contact angle, we use image processing tools to identify the drop and solid disk in its vertical median plane at the end of the simulation, as sketched in Figure 9(b). Once both discs have been identified, we measure the geometric parameters R, r, and d. The contact angle is then calculated using the generalized Pythagorean theorem

$$ \cos \theta = \frac{R^2 + r^2 - d^2}{2rR} \qquad (22) $$

We can also determine the so-called apparent contact angle $\theta' = \theta + \alpha$ (see Fig. 9).

<figure><img src="image_11.png" /><figcaption>**Fig. 9:** Scheme illustrating the geometric parameters required to calculate the contact angle $\theta$ on a sphere: (a) initial condition, (b) equilibrium condition.</figcaption></figure>

Additionally, the volume of the drop can be calculated. Indeed, in all our simulations, the final volume of the drop at equilibrium is not equal to its initial value due to the exchange of matter between the two fluid phases of the same component (water and its vapor). The calculation of the final drop volume requires the introduction of the two heights $h_1$ and $h_2$ as defined in Figure 10. From there, the volume of the drop is given by the difference between the volume of the sphere of radius $r$ and the intersection volume $V_{int}$ as colored in green in Fig. 9. This volume $V_{int}$ is simply the sum of two spherical caps of radii $R$ and $r$, and heights $h_1$ and $h_2$, respectively. From the formula of the volume of a spherical cap<fnref n="5" />, one gets:

$$ \begin{aligned} V_{int} &= V(R, h_1) + V(R, h_2) \\ &= \frac{\pi(R+r-d)^2(d^2+2dr-3r^2+2dR+6rR-3R^2)}{12d} \end{aligned} \qquad (23) $$

with

$$ h_1 = \frac{(r - R + d)(r + R - d)}{2d} \qquad (24) $$

<footnote n="5">⁵ The volume of a spherical cap of height $h$ and radius $R$ is $V_{cap}(R, h) = \frac{1}{3}\pi h^2(3R - h)$</footnote>
<figure><img src="image_12.png" /><figcaption>Fig. 10: Schematic defining the heights $h_1$ and $h_2$ used to calculate.</figcaption></figure>

$$h_2 = \frac{(R - r + d)(R + r - d)}{2d}. \qquad (25)$$

The volume of the drop is consequently equal to:

$$V = \frac{4}{3}\pi r^3 - V_{\text{int}}. \qquad (26)$$

### 3.3.2 Sensitivity to the mesh resolution

Figure 11 shows the equilibrium state by the end of the simulation for different mesh resolutions. The data obtained for the contact angle as a function of $\varphi$ are shown in Table 2.

Using as a reference value the one obtained with the finest resolution, we observe that the error decreases as the mesh size decreases. Here again, we will use the second mesh resolution to speed-up the calculation time as the error compared with the finest case remains limited to a few percents.

<table><caption>**Table 2:** Values obtained for the contact angle $\theta$ versus $\varphi$ for different mesh resolutions of a drop on a sphere. The errors are calculated based on the finest resolution.</caption><thead><tr><th>Resolution:<br>(N<sub>x</sub> × N<sub>y</sub> × N<sub>z</sub>)</th><th>128 × 128 × 128</th><th>256 × 256 × 256</th><th>320 × 320 × 320</th></tr></thead><tbody><tr><th scope="row">R (lu)</th><td>24</td><td>48</td><td>60</td></tr><tr><th scope="row">r<sub>0</sub> (lu)</th><td>14</td><td>30</td><td>37</td></tr><tr><th scope="row">φ</th><th>θ (°)</th><th>error (%)</th><th>θ (°)</th><th>error (%)</th><th>θ (°)</th><th>error (%)</th></tr><tr><th scope="row">0.8</th><td>87.96</td><td>2.98</td><td>91.6</td><td>1.05</td><td>90.65</td><td>-</td></tr><tr><th scope="row">1.0</th><td>67.7</td><td>3.4</td><td>69.65</td><td>0.62</td><td>70.09</td><td>-</td></tr><tr><th scope="row">1.2</th><td>54.25</td><td>7.21</td><td>55.71</td><td>4.71</td><td>58.46</td><td>-</td></tr><tr><th scope="row">1.4</th><td>46.04</td><td>2.55</td><td>48.13</td><td>1.88</td><td>47.24</td><td>-</td></tr><tr><th scope="row">1.6</th><td>38.95</td><td>3.85</td><td>39.25</td><td>3.11</td><td>40.51</td><td>-</td></tr><tr><th scope="row">1.8</th><td>33.47</td><td>2.98</td><td>32.16</td><td>1.06</td><td>32.5</td><td>-</td></tr><tr><th scope="row">2.0</th><td>19.85</td><td>14.72</td><td>22.25</td><td>4.38</td><td>23.27</td><td>-</td></tr></tbody></table>
<figure><img src="image_13.png" /><figcaption>**Fig. 11:** Equilibrium state obtained for different mesh resolutions with $\varphi = 1.2$: a) $N_x \times N_y \times N_z = 128 \times 128 \times 128$, b) $N_x \times N_y \times N_z = 256 \times 256 \times 256$, c) $N_x \times N_y \times N_z = 320 \times 320 \times 320$.</figcaption></figure>

### 3.3.3 Influence of the solid surface curvature radius

The influence of the solid curvature radius $R$ on the contact angle calculation is first studied. The initial drop radius is fixed at $r_0 = 30$ lu, while $R = 48$, 64, and 72 lu.

The corresponding data of the contact angle versus $\varphi$ is plotted in Figure 12. It can be seen that the results almost collapse on a unique curve, showing that the curvature radius of the solid surface has merely no influence on the contact angle, which is, here again, consistent with the fact that the contact angle is a material property independent from the problem geometry.
<figure><img src="image_14.png" /><figcaption>Fig. 12: Values obtained for the contact angle versus $\varphi$ with a fixed initial drop radius $r_0 = 30$ lu and three different solid radii for a drop on a sphere. The mesh resolution is 256×256×256.</figcaption></figure>

### 3.3.4 Influence of the initial drop volume

Similarly to section 3.3.3, we now examine the potential impact of the initial drop radius. To this end, Figure 13 shows the contact angle as a function of $\varphi$ for different initial drop radii $r_0$. The previous data of a drop on a flat surface are added for comparison. It can be seen that the data are almost all identical, especially for small values of $\varphi$ ($\varphi < 1.4$ or $\theta > 50^\circ$)$^6$. In the same way, we can conclude that the impact of the drop volume is negligible, confirming again that the contact angle can be definitely considered as a material property depending only on the model parameter $\varphi$.

It can be noted that based on the simulation data, the parameter $\varphi$ can be approximated by a polynomial of degree four in $\theta$ as shown in Fig. 14 and given in Eq. (27). Then, for any value of the contact angle $\theta$ in radians, this empirical interpolation enables to prescribe the parameter $\varphi$ to be used in the LBM simulation to recover the target value of $\theta$. Note finally that for a fixed value of $\theta$ corresponds a unique value of $\varphi$ for $\theta \le 90^\circ$.

$$ \varphi = 1.655 + 3.216\theta - 8.048\theta^2 + 5.697\theta^3 - 1.330\theta^4. \quad (27) $$

Finally, Figure 15 shows the dimensionless drop volume at equilibrium $V/R^3$ versus $\varphi$. We see that the method for calculating the volume, from (26), gives the same results and collapse with the LBM simulations, except for the larger volume where a small gap appears (of the order of 4%). Moreover, we can conclude that the final volume of the drop does not depend on $\varphi$ and therefore on the contact angle $\theta$. Equivalently, for an imposed volume $V_0$, the final volume at equilibrium does not depend on $\theta$, that is physically consistent.

<footnote>$^6$ For small values of $\theta$, the improved virtual-density scheme [17] appears less accurate for a drop on convex solid surfaces than on flat surfaces due to the discretization of the surface.</footnote>
<figure><img src="image_15.png" /><figcaption>**Fig. 13:** Values obtained for the contact angle versus $\varphi$ for three different initial drop radii $r_0$ and a fixed solid radius $R = 72$ lu. The data of a drop on a flat surface are added.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 14:** Empirical calibration of $\varphi$ versus $\theta$ with a 4th order polynomial approximation. The data correspond to the situation of an initial drop radius $r_0 = 15$ lu and a convex surface with a fixed radius of curvature $R = 72$ lu.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Fig. 15:** Values obtained for the dimensionless drop volume at equilibrium $V/R^3$ versus $\phi$ for different initial drop volumes $V_0$ and a fixed radius of curvature of the solid surface $R = 72$ lu.</figcaption></figure>
### 3.3.5 Comparison to experiments

To validate our results, we simulate the same experiments as in [39], where water drops of different volumes were deposed on a spherical PC (polycarbonate) surface with a radius of curvature $R = 6.35$ mm. The measured value of the contact angle was $\theta = 89^{\circ}$. Based on the proposed empirical relation (27), the wetting parameter is fixed at $\phi = 0.769$ giving a contact angle equal to the experimental value.

Figure 16 shows a comparison of the apparent contact angle $\theta' = \theta + \alpha$ (Fig. 9) between the experimental data of [39] and the values obtained with our simulation as a function of the reduced volume $V/R^3$, with $R$ the radius of curvature of the solid surface. The agreement is rather satisfactory, even if we observe a small deviation in the simulation results. This may be due to the solid surface which is not perfectly spherical in the simulation as it is represented through a stair-case approximation which may cause small errors in the measurement of the apparent contact angle.

<figure><img src="image_18.png" /><figcaption>**Fig. 16:** Values obtained for the apparent contact angle $\theta'$ versus dimensionless drop volume $V/R^3$ for a radius of curvature of the solid surface $R = 6.35$ mm. The experimental results of [39] are included for water drops on PC surface ($\theta = 89^{\circ}$). The error bar in the simulation is induced by the mesh resolution with a one-pixel uncertainty in the calculation of the different geometric quantities.</figcaption></figure>

## 4 Capillary forces induced by a liquid bridge between two spherical particles

As already explained in the introduction, the aim of this work is to develop an efficient LBM method, fast and accurate enough, to simulate capillary interfaces between several solid particles and to compute the resulting capillary forces.
The last study being addressed in the present work therefore concerns the calculation of the capillary force between two solid spheres connected by a liquid bridge. To this end, two solid particles with the same radius $R = 102$ lu ($R = 8$ mm) are placed in the domain and separated with a dimensionless separation distance $D^* = D/R = 0.2$. The parameter $\varphi$ was fixed at 2.0 which corresponds to a contact angle $\theta = 21.4^{\circ}$ from Eq. (27).

The lattice size in physical units is fixed at $\Delta x_{phy} = L_c/(N_x - 1)$, where $L_c = 0.03$ m is the domain length in x-direction. Snapshots of the liquid bridge shape during the simulation is shown in Figure 17, corresponding to condensation phenomenon.

<figure><img src="image_19.png" /><figcaption>Fig. 17: The numerical (top) and theoretical (bottom) shapes of the liquid capillary bridge profiles during the simulation for different dimensionless volumes ($V^* = V/R$): (a) $V^* = 0.016$, (b) $V^* = 0.041$, (c) $V^* = 0.085$, and (d) $V^* = 0.15$.</figcaption></figure>

## 4.1 Brief reminder on capillary bridges modelling and associated forces

Let us recall that the shape of a capillary bridges between two spherical particles is axisymmetric of revolution along the x-axis connecting the centers of the two particles. It is given by a solution of the Young-Laplace equation:

$$ \frac{y''(x)}{(1 + y'^2(x))^{3/2}} - \frac{1}{y(x)\sqrt{1 + y'^2(x)}} = -\frac{\Delta P}{\gamma} = H \quad (28) $$

where $H = -\frac{\Delta P}{\gamma}$ denotes the mean curvature which is constant in all the capillary bridges, while $\Delta P = P_{in} - P_{out}$ stands for the pressure difference between inside and outside of the
<figure><img src="image_20.png" /><figcaption>Fig. 18: Illustration of a capillary bridge between two identical spherical particles. δ and θ are the filling and contact angle, respectively. Figure extracted from [11].</figcaption></figure>

liquid bridge, and γ is the surface tension of the fluid. In Figure 18, δ denotes the filling angle, y* the gorge radius, and θ the wetting angle. D is the separation distance between the particles. At static equilibrium, the capillary force is a first integral of the Young–Laplace equation [11, 19] and can be calculated at the gorge radius y* as

$$F_{cap} = \pi\gamma H y^{*2} + 2\pi\gamma y^{*} \qquad (29)$$

or at the contact line following

$$F_{cap} = \pi\gamma H R^2 \sin^2\delta + 2\pi\gamma R \sin\delta \sin(\delta + \theta). \qquad (30)$$

Therefore, the capillary force classically includes two contributions: the first one results from the pressure difference inside and outside of the bridge (Laplace pressure) and the second contribution is due to surface tension along the wetted interfaces.

Note that equation (28) can be also solved analytically using a cylindrical approximation [29] leading to a direct and convenient relationship between the dimensionless capillary forces $F_{cap}^* = F_{cap}/(2\pi\gamma R)$, the dimensionless inter-particles distance $D^* = D/R$, and the dimensionless volume $V^* = V/R^3$:

$$F_{cap}^* = \cos\theta \left( 1 - \frac{1}{\sqrt{1 + \frac{2V^*}{\pi D^{*2}}}} \right). \qquad (31)$$

This relationship is however only valid for small volumes, typically for $V^* < 0.01$. Richefeu et al. [31] further proposed the following empirical relation for the capillary force based on fitting the numerical solution of Laplace-Young equation for two particles of same diameter:

$$F_{cap}^* = \cos\theta \exp\left(\frac{-D^*}{0.9\sqrt{V^*}}\right). \qquad (32)$$

These two approximate expressions of capillary force will be used in section 4.3 to validate our LBM numerical simulations of capillary bridges.
## 4.2 Numerical calculation of the capillary forces

The total force exerted by the fluid on the solid can be calculated based on the integration of the fluid stress tensor $\sigma$ on a surface $\Omega$ close to the solid particle as illustrated in Figure 19:

$$F_t = \int_{\Omega} \sigma \cdot n dA \qquad (33)$$

where $n$ is the outer unit normal vector to the surface $\Omega$ and $dA$ stands for an elementary area element. For viscous Newtonian fluids, the fluid stress tensor is given by

$$\sigma_{ij} = -P\delta_{ij} + \tau_{ij}. \qquad (34)$$

In expression (34), $P$ is the fluid pressure given by Eq. (10) whereas the viscous stress tensor $\tau_{ij}$, reads

$$\tau_{ij} = \mu(\partial_j u_i + \partial_i u_j) \qquad (35)$$

where $\mu$ is the dynamic fluid viscosity ($\mu = \rho v$).

Note that we follow almost the same integration technique as in [34] but with different stress tensor formula, since they use a free-energy method to model two-phase liquid-vapor flows which are different than the shan-chen model.

The discretized form of Eq. (33), that will be used for LBM numerical calculations reads

$$F_t = \sum_{x \in \Omega} \sigma \cdot n dA \qquad (36)$$

where $x \in \Omega$ accounts for all lattice points located on the voxelized surface $\Omega$, as shown in the zoom of Fig. 19. The area element $dA$ is equal to 1 in lattice units. Note that the center of the cubes of the voxelized surface $\Omega$ are considered to be located at the lattice points to avoid interpolation of fluid quantities. The viscous stress tensor $\tau_{ij}$ is calculated from the LBM simulation using finite difference method applied to Eq. (35).

It is important to stress that Eqs. (33) or (36) do not include the surface tension term that is present in expression (29) or (30) of capillary forces. However since the interface has a finite thickness in the LBM computation (Fig. 4b), we will prove through the LBM simulations that the surface tension term is implicitly included in the contribution of the stress in the interface.

Before concluding on this point, we will check the mesh dependency of the capillary force measurement. To do so, the capillary forces will be computed at two different distances $\epsilon$ away from the solid particle for three different resolutions, as indicated in Table. 3 and Fig. 19. Note that $\epsilon$ is constant in physical units (mm) but varies in lattice units (lu) according to the chosen mesh resolution. The reason for using a strictly positive epsilon is to get rid of non physical phases junction of the fluid pressure on the solid/fluid interface.

The capillary forces calculated with the three mesh resolutions and the two values of $\epsilon$ are displayed in Figure 20 versus the dimensionless volume $V^*$ of the water bridge and for $D^* = 0.2$. As can be seen, the resolution has a limited effect on the capillary force, especially for the smallest epsilon value. Therefore, the simulations and calculations of the associated capillary forces will be done with the finest mesh 384×450×384 in the following to minimize the error on the wetting angle. The GPU simulation time was about 100 min for this resolution using a GPU card Quadro RTX 5000.

<footnote n="7">7 The capillary force at this position was linearly interpolated between 7th and 8th lattice.</footnote>
<figure><img src="image_21.png" /><figcaption>**Fig. 19:** Snapshot of a simulation of a liquid bridge between two spherical particles. $\Omega$ is the integration surface. In this example, the dimensionless separation distance is $D^* = D/R = 0.2$ and the reduced liquid volume is $V^* = V/R^3 = 0.18$.</figcaption></figure>

<table><caption>**Table 3:** Distance of integration $\epsilon$ in physical (mm) and lattice (lu) units for three different mesh resolutions.</caption><thead><tr><th>Resolution:<br>(N<sub>x</sub> × N<sub>y</sub> × N<sub>z</sub>)</th><th>256 × 384 × 256</th><th>320 × 380 × 320</th><th>384 × 450 × 384</th></tr><tr><th>ε (mm)</th><th colspan="3">ε (lu)</th></tr></thead><tbody><tr><td>0.471</td><td>4</td><td>5</td><td>6</td></tr><tr><td>0.706</td><td>6</td><td>7.5<sup>7</sup></td><td>9</td></tr></tbody></table>

The stress integration results are given in Figure 21. As a general trend, it can be observed that the capillary forces increase as the volume increases in the range $0 \le V^* \le 0.4$. Also, the capillary forces increase with the distance $\epsilon$ from the solid surface where they are calculated. It thus appears that the choice of the integration surface $\Omega$ is particularly crucial. The reason why the integration surface $\Omega$ is shifted at some distance $\epsilon$ from the solid particles is that a thin liquid film is created near the outer solid surface due to the implementation of the wetting condition as shown in Figure 22 where the density profile is plotted versus the distance from the solid surface. Therefore, $\epsilon$ should be chosen higher than or equal to 5 lu ($\epsilon \ge 5$ lu). In the sequel, for the comparison between the capillary force calculated with LBM simulation and the existing results from the literature, we will consider $\epsilon = 5$ lu.
<figure><img src="image_22.png" /><figcaption>**Fig. 20:** Normalized capillary forces versus normalized volume for three different resolutions at two distances of integration $\epsilon = 0.471$ mm and $\epsilon = 0.706$ mm.</figcaption></figure>

## 4.3 Comparison with the literature

In Figure 23, we compare the results obtained with our LBM simulations using $\epsilon = 5$ lu during the condensation of the liquid bridge with the two theoretical expressions of capillary force given by Eq. (29) and (30), and with the approximations given by Eq. (31) and (32). The experimental data of [19] are also added. It can be observed that the LBM simulations with $\epsilon = 5$ lu give values of the capillary force very close to the gorge radius expression (29) and in good agreement with the contact line expression (30), the cylindrical approximation of Eq. (31), and the experiments of [19]. Based on the founding, we can conclude that the capillary force is a result of the integration of the pressure and the viscous part, and the fluid velocity does not tend toward zero at equilibrium, which is specific to the shan-chen model with the velocity shift force scheme. In other words the viscous part contribute to the capillary force, even when the capillary bridge is at equilibrium.

To finish, Figure 24 shows the dimensionless mean curvature $H^* = H \times R$, with $H = -\Delta P/\gamma$, versus $V^* = V/R^3$, where $\Delta P$ is deduced from the LBM simulation results. Two different shapes are encountered: nodoid shape, which corresponds to $H > 0$, and unduloid shape for $H < 0$, as described in detail in [11]. The theoretical $H$ values are determined using inverse problem technique to solve Young-Laplace equation, the same image processing in-house code as in [11] is used.

The effect of the wetting angle on the capillary force is shown in Fig. 25. The stress integration method is used with $\epsilon = 5$ lu. We see that the forces decreases as the wetting angle increases.
<figure><img src="image_23.png" /><figcaption>**Fig. 21:** Normalized capillary force versus normalized volume during the condensation of a liquid bridge between two spherical particles of same radius *R* at different distances *ε*.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Fig. 22:** Density profile versus the distance *ε* from the solid surface out of the capillary bridge region, where ρ_g = 0.015 lu corresponds to the gas phase. The density of the liquid phase ρ_I = 0.333 lu.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Fig. 23:** Comparison of normalized capillary force versus normalized volume calculated using stress integration method with $\varepsilon = 5$ lu, gorge radius (Eq. (29)), contact line (Eq. (30)), cylindrical approximation (Eq. (31)), approximation given by Richefeu (Eq. (32)), and experimental data of [19].</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Fig. 24:** Variation of the dimensionless mean curvature $H^* = H \times R$ of the capillary bridge between two spheres versus dimensionless volume. $H > 0$ corresponds to a nodoid shape and $H < 0$ corresponds to an unduloid shape.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>**Fig. 25: Normalized capillary force versus normalized volume for different contact angles.**</figcaption></figure>
## 5 Conclusion

In this paper, we have studied the Single Component Multi-Phase (SCMP) Shan-Chen model for simulating capillary interfaces in various configurations and for calculating the induced capillary forces. The SCMP model was implemented in 3D and solved using GPU. The relevance of the model was verified by performing several benchmarks, such as drop on flat and spherical walls and stationary spherical drops. Regarding the latter, results found for the surface tension and the contact angle are consistent. Moreover, we have highlighted a direct relation between the parameters of the virtual density scheme and the induced wetting angle $\theta$ that enables to easily adjust the parameters to directly determine the desired value of $\theta$. Next, we have considered capillary bridges between two spherical particles and analyzed in details the capillary interface, and then calculated the associated capillary force. Although we have revealed some limitations of the SCMP Shan-Chen model (high values of temperature and pressure to guarantee the numerical stability which do not correspond to classical experimental values, difficulty to control the interface thickness), we have proved however that the computation of the capillary forces at some distance from the solid grain, provides results quite close to analytical expressions existing in literature and to experiments performed in similar conditions while avoiding spurious terms in the calculation. This way, we can justify the values of the parameters usually used in LBM SCMP model in literature, which are given in lattice units and never compared to usual values encountered in lab experiments or in natural conditions. Therefore, accounting for the very small computation time needed to stabilize the capillary interface, this approach may be considered promising to simulate the behavior of a partially saturated granular assembly containing a large number of particles. Also, this approach is suited for time evolution problems, such as partially saturated granular media subjected to condensation/evaporation cycles for instance.

## Acknowledgments

The authors would like to express their sincere gratitude to the French National Centre for Space Studies (CNES) and to the NEEDS program for having supported this work. They also acknowledge the International Research Network GeoMech (IRN CNRS) for enabling intensive and productive interactions between all of them. In particular the authors would like to thank Marie Miot for insightful discussions on the different strategies to compute capillary forces in numerical simulations.

## Appendix: Conversion between physical and lattice units

The LBM simulations with Shan-Chen model involve four dimensional quantities: time, length, mass and temperature. The conversion between lattice units and physical units is done with use of the reduced properties concept [4, 41] for the fluid properties (namely density $\rho$, pressure $P$, and temperature $T$) and with one additional conversion factors for length.

$$ \rho_R = \frac{\rho}{\rho_c}, \quad P_R = \frac{P}{P_c}, \quad T_R = \frac{T}{T_c}, \qquad (37) $$

where the subscript "R" and "c" are the reduced and critical properties, respectively. According to this concept, the reduced properties in lattice and physical units should be equal. As an example, we have $\rho_R^{lu} = \rho_R^{phy}$, leading to $\rho^{phy} = \rho^{lu}\rho_c^{phy}/\rho_c^{lu}$. The other properties can be obtained in a similar manner.

The conversion factors for length, time and density are respectively $C_l = \Delta x_{phy}/\Delta x_{lu}$, $C_t = \Delta t_{phy}/\Delta t_{lu}$, and $C_\rho = \rho_{phy}/\rho_{lu} = \rho_c^{phy}/\rho_c^{phy}$. The conversion factor for force can be deduced as $C_f = C_\rho C_l^4/C_t^2$.
494 To get the value in physical units, the corresponding value in lattice units is multiplied by the conversion factor that has the same physical units. For example, for length we have $L_{phy} = L_{lu} \times C_l$.

496 For capillary interfaces with no flow, where the Reynolds number is no longer relevant, the time scale is not given by viscosity <fnref n="8" />. The conversion between surface tension $\gamma$ in lattice and physical units reads

$$ \gamma^{phy} = \gamma^{lu} C_{\rho} \frac{C_l^3}{C_l^2} \qquad (38) $$

498 In the presence of gravity, the conversion between gravity in lattice and physical units is given as

$$ g^{phy} = g^{lu} \frac{C_l}{C_l^2} \qquad (39) $$

499 Alternatively, physical units can be related to lattice units through dimensionless numbers (e.g. the Bond number (Bo)) instead of using the conversion factors. As an example, the dimensionless Bond number is defined as

$$ Bo = \frac{(\rho_l - \rho_g)gr^2}{\gamma} \qquad (40) $$

502 where $g$ is the gravity, $\gamma$ is the surface tension and $r$ is the length scale (e.g. drop radius in the case studied in section 3.2). By setting $Bo_{lattice} = Bo_{physical}$, the conversion leads to

$$ \frac{(\rho_l^{phy} - \rho_g^{phy})g^{phy}(r^{phy})^2}{\gamma^{phy}} = \frac{(\rho_l^{lu} - \rho_g^{lu})g^{lu}(r^{lu})^2}{\gamma^{lu}}. \qquad (41) $$

504 Expression (41) enables to fix the gravity in lattice units $g^{lu}$ that can be used in LBM simulation when taking it into account.

## References

507 1. Benzi, R., Biferale, L., Sbragaglia, M., Succi, S., Toschi, F.: Mesoscopic modeling of a two-phase flow in the presence of boundaries: The contact angle. Phys. Rev. E 74, 021509 (2006). DOI 10.1103/PhysRevE.74.021509

508 2. Benzi, R., Succi, S., Vergassola, M.: The lattice Boltzmann equation: theory and applications. Physics Reports 222(3), 145–197 (1992). DOI 10.1016/0370-1573(92)90090-M

510 3. Brakke, K.A.: The surface evolver. Experimental Mathematics 1(2), 141–165 (1992). DOI 10.1080/10586458.1992.10504253

511 4. Cartwright, H.: Molecular thermodynamics. by donald a. mcquarrie and john d. simon.(1999) university science books, 55d gate five road, causalito ca 94965, usa. 672 pp $78.00, isbn 1-891389-05-x. The Chemical Educator 4(3), 120–121 (1999)

517 5. Chen, L., Kang, Q., Mu, Y., He, Y.L., Tao, W.Q.: A critical review of the pseudopotential multiphase lattice boltzmann model: Methods and applications. International Journal of Heat and Mass Transfer 76, 210–236 (2014). DOI https://doi.org/10.1016/j.ijheatmasstransfer.2014.04.032

518 6. Connington, K.W., Lee, T., Morris, J.F.: Interaction of fluid interfaces with immersed solid particles using the lattice boltzmann method for liquid–gas–particle systems. Journal of Computational Physics 283, 453–477 (2015). DOI https://doi.org/10.1016/j.jcp.2014.11.044

523 7. Delenne, J.Y., Richefeu, V., Radjai, F.: Liquid clustering and capillary pressure in granular media. Journal of Fluid Mechanics 762, R5 (2015). DOI 10.1017/jfm.2014.676

525 8. Duriez, J., Wan, R.: Contact angle mechanical influence in wet granular soils. Acta Geotechnica 12(1), 67–83 (2017)

527 9. Fakhari, A., Mitchell, T., Leonard, C., Bolster, D.: Improved locality of the phase-field lattice-boltzmann model for immiscible fluids at high density ratios. Physical Review E 96(5), 053301 (2017)

529 10. Fan, Z., Qiu, F., Kaufman, A., Yoakum-Stover, S.: GPU cluster for high performance computing. IEEE/ACM SC2004 Conference, Proceedings pp. 297–308 (2004)

531 11. Gagneux, G., Millet, O.: Analytic calculation of capillary bridge properties deduced as an inverse problem from experimental data. Transport in Porous Media 105, 117–139 (2014)

<footnote n="8">⁸ The time scale for flows that are controlled by Reynolds number is classically given by: $\Delta t_{phy} = c_s^2(\tau - \frac{1}{2})\Delta x_{phy}^2/\nu_{phy}$, where $c_s$ and $\tau$ in lattice units [16].</footnote>
533

534

535

536

537

538

539

540

541

542

543

544

545

546

547

548

549

550

551

552

553

554

555

556

557

558

559

560

561

562

563

564

565

566

567

568

569

570

571

572

573

574

575

576

577

578

579

580

581

582

583

584

585

586

587

588

589

590

591

592

12. Guo, Z., Zheng, C., Shi, B.: Discrete lattice effects on the forcing term in the lattice boltzmann method. Phys. Rev. E **65**, 046308 (2002). DOI 10.1103/PhysRevE.65.046308

13. Huang, H., Krafczyk, M., Lu, X.: Forcing term in single-phase and shan-chen-type multiphase lattice boltzmann models. Phys. Rev. E **84**, 046710 (2011). DOI 10.1103/PhysRevE.84.046710

14. Huang, H., Li, Z., Liu, S., Lu, X.y.: Shan-and-chen-type multiphase lattice boltzmann study of viscous coupling effects for two-phase flow in porous media. International Journal for Numerical Methods in Fluids **61**(3), 341–354 (2009). DOI https://doi.org/10.1002/fld.1972

15. Huang, H., Sukop, M., Lu, X.: Multiphase lattice Boltzmann methods: Theory and application. John Wiley & Sons (2015)

16. Krueger, T., Kusumaatmaja, H., Kuzmin, A., Shardt, O., Silva, G., Viggen, E.: The Lattice Boltzmann Method: Principles and Practice. Graduate Texts in Physics. Springer (2016)

17. Li, Q., Yu, Y., Luo, K.H.: Implementation of contact angles in pseudopotential lattice boltzmann simulations with curved boundaries. Phys. Rev. E **100**, 053313 (2019). DOI 10.1103/PhysRevE.100.053313

18. Liang, H., Xu, J., Chen, J., Wang, H., Chai, Z., Shi, B.: Phase-field-based lattice boltzmann modeling of large-density-ratio two-phase flows. Physical Review E **97**(3), 033309 (2018)

19. Mielniczuk, B., Millet, O., Gagneux, G., El Youssoufi, M.S.: Characterisation of pendular capillary bridges derived from experimental data using inverse problem method. Granular Matter **20**(1), 14 (2018). DOI 10.1007/s10035-017-0784-8

20. Miot, M., Wautier, A., Veylon, G., Philippe, P., Nicot, F.: Numerical modeling of capillary forces in mesoscale assemblies of grains: from pendular to funicular regimes. accepted for publication in Granular Matter (2021)

21. Montellá, E.P., Yuan, C., Chareyre, B., Gens, A.: Hybrid multi-scale model for partially saturated media based on a pore network approach and lattice boltzmann method. Advances in Water Resources **144**, 103709 (2020). DOI https://doi.org/10.1016/j.advwatres.2020.103709

22. Nguyen, H.N.G., Millet, O., Gagneux, G.: Exact calculation of axisymmetric capillary bridge properties between two unequal-sized spherical particles. Mathematics and Mechanics of Solids **24**(9), 2767–2784 (2018). DOI 10.1177/1081286518787842

23. Nguyen, H.N.G., Millet, O., Gagneux, G.: On the capillary bridge between spherical particles of unequal size: analytical and experimental approaches. Continuum Mechanics and Thermodynamics **31**(1), 225–237 (2018). DOI 10.1007/s00161-018-0658-2

24. Nguyen, H.N.G., Millet, O., Gagneux, G.: Liquid bridges between a sphere and a plane - classification of meniscus profiles for unknown capillary pressure. Mathematics and Mechanics of Solids **24**(10), 3042–3060 (2019). DOI 10.1177/1081286519831047

25. Nguyen, H.N.G., Millet, O., Zhao, C.F., Gagneux, G.: Theoretical and experimental study of capillary bridges between two parallel planes. European Journal of Environmental and Civil Engineering **0**(0), 1–11 (2020). DOI 10.1080/19648189.2019.1706055

26. Nguyen, H.N.G., Zhao, C.F., Millet, O., Gagneux, G.: An original method for measuring liquid surface tension from capillary bridges between two equal-sized spherical particles. Powder Technology **363**, 349–359 (2020). DOI https://doi.org/10.1016/j.powtec.2019.12.049

27. Nguyen, H.N.G., Zhao, C.F., Millet, O., Selvadurai, A.: Effects of surface roughness on liquid bridge capillarity and droplet wetting. Powder Technology **378**, 487–496 (2021). DOI https://doi.org/10.1016/j.powtec.2020.10.016

28. Pan, C., Hilpert, M., Miller, C.T.: Lattice-boltzmann simulation of two-phase flow in porous media. Water Resources Research **40**(1) (2004). DOI https://doi.org/10.1029/2003WR002120

29. Pitois, O., Moucheront, P., Chateau, X.: Liquid bridge between two moving spheres: An experimental study of viscosity effects. Journal of Colloid and Interface Science **231**(1), 26–31 (2000). DOI https://doi.org/10.1006/jcis.2000.7096

30. Richefeu, V., Radjai, F., Delenne, J.Y.: Lattice boltzmann modelling of liquid distribution in unsaturated granular media. Computers and Geotechnics **80**, 353–359 (2016). DOI https://doi.org/10.1016/j.compgeo.2016.02.017

31. Richefeu, V., Radjai, F., El Youssoufi, M.: Stress transmission in wet granular materials. European Physical Journal E **21**(4), 359–369 (2006). DOI 10.1140/epje/i2006-10077-1

32. Shan, X., Chen, H.: Lattice boltzmann model for simulating flows with multiple phases and components. Physical Review E **47**(3), 1815–1819 (1993). DOI 10.1103/PhysRevE.47.1815. Cited By 2297

33. Shan, X., Chen, H.: Simulation of nonideal gases and liquid-gas phase transitions by the lattice boltzmann equation. Physical Review E **49**(4), 2941–2948 (1994). DOI 10.1103/PhysRevE.49.2941. Cited By 954

34. Shinto, H., Komiyama, D., Higashitani, K.: Lattice boltzmann study of capillary forces between cylindrical particles. Advanced Powder Technology **18**(6), 643–662 (2007). DOI https://doi.org/10.1163/156855207782514950. URL https://www.sciencedirect.com/science/article/pii/S0921883108608602
35. Singh, P., Joseph, D.D.: Fluid dynamics of floating particles. Journal of Fluid Mechanics **530**, 31–80 (2005). DOI 10.1017/S0022112005003575

36. Sukop, M.C., Or, D.: Lattice boltzmann method for modeling liquid-vapor interface configurations in porous media. Water Resources Research **40**(1) (2004). DOI https://doi.org/10.1029/2003WR002333

37. Sun, X., Sakai, M.: Direct numerical simulation of gas-solid-liquid flows with capillary effects: An application to liquid bridge forces between spherical particles. Phys. Rev. E **94**, 063301 (2016). DOI 10.1103/PhysRevE.94.063301

38. Swift, M.R., Osborn, W., Yeomans, J.: Lattice boltzmann simulation of nonideal fluids. Physical review letters **75**(5), 830 (1995)

39. Wu, D., Wang, P., Wu, P., Yang, Q., Liu, F., Han, Y., Xu, F., Wang, L.: Determination of contact angle of droplet on convex and concave spherical surfaces. Chemical Physics **457**, 63 – 69 (2015). DOI https://doi.org/10.1016/j.chemphys.2015.05.020

40. Yang, L., Sega, M., Harting, J.: Capillary-bridge forces between solid particles: Insights from lattice boltzmann simulations (2020)

41. Yuan, P., Schaefer, L.: Equations of state in a lattice boltzmann model. Physics of Fluids **18**(4), 042101 (2006). DOI 10.1063/1.2187070

42. Zheng, H., Shu, C., Chew, Y.T.: A lattice boltzmann model for multiphase flows with large density ratio. Journal of computational physics **218**(1), 353–371 (2006)

43. Zou, Q., He, X.: On pressure and velocity boundary conditions for the lattice Boltzmann BGK model. Physics of Fluids **9**(6), 1591–1598 (1997). DOI 10.1063/1.869307

44. Zu, Y., He, S.: Phase-field-based lattice boltzmann model for incompressible binary fluid systems with density and viscosity contrasts. Physical Review E **87**(4), 043301 (2013)