Simulation of Capillary-Driven Kinetics with Multi-Phase-Field and Lattice-Boltzmann Method

TL;DR

This paper introduces a combined multi-phase-field and lattice Boltzmann computational approach to simulate capillary-driven motion of solid particles, enabling detailed study of liquid-phase sintering and particle dynamics.

Contribution

The paper presents a novel integrated method for simulating capillary-driven solid particle motion, validated against analytic solutions and applied to complex multi-particle systems.

Findings

Method accurately reproduces analytic solutions for simple cases.

Liquid morphology significantly influences particle dynamics.

Suitable for studying liquid-phase sintering processes.

Abstract

We propose a combined computational approach based on the multi-phase-field and the lattice Boltzmann method for the motion of solid particles under the action of capillary forces. The accuracy of the method is analyzed by comparison with the analytic solutions for the motion of two parallel plates of finite extension connected by a capillary bridge. The method is then used to investigate the dynamics of multiple spherical solid bodies connected via capillary bridges. The amount of liquid connecting the spheres is varied, and the influence of the resulting liquid-morphology on their dynamics is investigated. It is shown that the method is suitable for a study of liquid-phase sintering which includes both phase transformation and capillary driven rigid body motion.

Tables

${\rho }_0$	reference density	1.000
${\rho }_L$	liquid density	1.888
$\tau$	relaxation time	1.100
$𝒢$	coupling constant	-5.000
${\rho }_{\alpha }$	mass density of ${\varphi }_{\alpha }$	various
${\varrho }_{\alpha }$	wetting parameter of ${\varphi }_{\alpha }$	2.000
$M_{\alpha \beta }$	phase-field mobility	${10}^{-5}$
$\eta$	phase-field interface width	5.000

Equations

$$dU={\sigma }_{SL}d{A}_{SL}+{\sigma }_{SV}d{A}_{SV}+{\sigma }_{LV}d{A}_{LV}-p_{V}d{V}_V-p_{L}d{V}_L$$

$$dU=-\left({\sigma }_{SV}-{\sigma }_{SL}\right)d{A}_{SL}+{\sigma }_{LV}d{A}_{LV}-\Delta pd{V}_L$$

$$\Delta p={\sigma }_{LV}\kappa$$

$$dU={\sigma }_{LV}d{A}_{LV}$$

$$F_p=2{\sigma }_{LV}L$$

$$F_c=\pi {\sigma }_{LV}R$$

$$h_p=h_0-\frac{1}{2}\frac{F_p}{m}{t}^2-\frac{1}{2}\frac{F_p}{m}{t}^2$$

$$h_p=h_0-\frac{2\sigma L}{m}{t}^2$$

$$h_c\left(t\right)\approx h_0-\frac{\pi \sigma R_0}{m}{t}^2$$

$${\Theta }_S=\arccos \left(\frac{{\chi }_s^2-2}{2}\right)-\arccos \left(\frac{{\chi }_s}{2}\right)$$

$${\chi }_s=\frac{2}{\sqrt[4]{3}}\frac{\sqrt{\left(1-c\right)\sqrt{c\left(1-c\right)}}}{1-c}$$

$${\Theta }_c=\frac{\pi }{2}-\arccos \left(1-{\chi }_c\right)$$

$${\chi }_c=-\cos \left(\frac{1}{3}\arccos \left(\frac{1-9c}{1-c}\right)+\frac{\pi }{3}\right)+\frac{1}{2}$$

$${\ddot{\vec{X}}}_B=M_B^{-1}{\vec{F}}_B$$

$${\dot{\vec{\omega }}}_B={\underline{\underline{I}}}_B^{-1}{\vec{T}}_B$$

$${\vec{u}}_B\left(\vec{x}\right)={\dot{\vec{X}}}_B+{\vec{\omega }}_B\times \left(\vec{x}-{\vec{X}}_B\right)$$

$${\dot{\varphi }}_{\alpha }={\vec{u}}_B\cdot \nabla {\varphi }_{\alpha }$$

$$e\left(\vec{x}_{i},\vec{x}_{j}\right)=e_{0}\left\{\begin{array}[]{ll}\left|\frac{\left\|\vec{x}_{i}-\vec{x}_{j}\right\|-r_{c}}{r_{c}}\right|^{n}&\textrm{if }\left\|\vec{x}_{i}-\vec{x}_{j}\right\|<r_{c}\\ 0&\textrm{else}\end{array}\right.$$

$$\vec{\mathsf{f}}\left({\vec{x}}_i,{\vec{x}}_j\right)=\frac{\partial e\left({\vec{x}}_i,{\vec{x}}_j\right)}{\partial \left({\vec{x}}_i-{\vec{x}}_j\right)}$$

$$=e_0\frac{n}{r_c}\frac{{\vec{x}}_i-{\vec{x}}_j}{\Vert {\vec{x}}_i-{\vec{x}}_j\Vert }{\left|\frac{\Vert {\vec{x}}_i-{\vec{x}}_j\Vert -r_c}{r_c}\right|}^{n-1}$$

$$\vec{c}_{i}=\left\{\begin{array}[]{ll}(0,0,0)&i=0\\ (\pm c,0,0),(0,\pm c,0),(0,0,\pm c)&i=1,\ldots,6\\ (\pm c,\pm c,0),(0,\pm c,\pm c),(\pm c,0,\pm c)&i=7,\ldots,18\\ (\pm c,\pm c,\pm c)&i=19,\ldots 26\end{array}\right.$$

$$f_i^{\ast }\left(\vec{x},t\right)=f_i\left(\vec{x},t\right)-\frac{1}{\tau }\left[f_i\left(\vec{x},t\right)-f_i^{eq}\left(\rho \left(\vec{x},t\right),\vec{u}\left(\vec{x},t\right)\right)\right]+\Delta t{F}_i$$

$$f_i\left(\vec{x},t+\Delta t\right)=f_i^{\ast }\left(\vec{x}-\overrightarrow{c_i}\Delta t,t\right)$$

$$f_i^{eq}\left(\rho ,\vec{u}\right)=\rho \sum _i{w}_i\left[1-\frac{\vec{u}\cdot \vec{u}}{2c_s^2}+\frac{\overrightarrow{c_i}\cdot \vec{u}}{c_s^2}+\frac{{\left(\overrightarrow{c_i}\cdot \vec{u}\right)}^2}{2c_s^4}\right]$$

$$w_{i}=\left\{\begin{array}[]{ll}8/27&i=0\\ 2/27&i=1,\ldots,6\\ 1/54&i=7,\ldots,18\\ 1/216&i=19,\ldots 26\end{array}\right.$$

$$F_i=\left(1-\frac{1}{2\tau }\right)w_i\left[\frac{{\vec{c}}_i-\vec{v}}{c_s^2}+\frac{{\vec{c}}_i\cdot \vec{v}}{c_s^4}{\vec{c}}_i\right]\cdot \vec{\mathsf{f}}$$

$$\rho \vec{v}=\sum _i{\vec{c}}_i{f}_i+\frac{\Delta t}{2}\vec{F}$$

$$f_{i}\left(\vec{x},t+\Delta t\right)=\left\{\begin{array}[]{ll}f^{*}_{-i}\left(\vec{x},t\right)+g\left(\vec{x},\vec{c_{i}},\vec{u}_{B},t\right)&\textrm{if }\vec{x}-\vec{c_{i}}\Delta{t}\textrm{ is solid}\\ f^{*}_{i}\left(\vec{x}-\vec{c_{i}}\Delta{t},t\right)&\textrm{else}\end{array}\right..$$

$$g\left(\vec{x},\overrightarrow{c_i},{\vec{u}}_B\left(\vec{x}\right)\right)=\frac{2\rho \left(\vec{x}\right)}{c_s^2}{w}_i\overrightarrow{c_i}\cdot {\vec{u}}_B\left(\vec{x}-{\vec{c}}_i\Delta t/2\right)+\delta \rho \left(\vec{x}\right){\delta }_{i0}$$

$$\delta \rho \left(\vec{x}\right)=\frac{2\rho \left(\vec{x}\right)}{c_s^2}{w}_i\overrightarrow{c_i}\cdot {\vec{u}}_B\left(\vec{x}-{\vec{c}}_i\Delta t/2\right)$$

Full text

††: Modelling Simulation Mater. Sci. Eng.

Keywords: Liquid-phase sintering, Wetting, Capillarity, Phase-field, Lattice Boltzmann, Multi-phase fluids

1 Introduction

During the sintering process, liquids can form intricate structures and bridges between multiple solid particles which by themselves can have a complex topology and can form large structures with other solid grains. These capillary bridges between the solid particles lead to compaction of the sample. Various theoretical approaches consider the force exerted by capillary bridges on solids. An early analysis of capillary forces in liquid-phase sintering processes of spherical particles can be found in [1]. The elementary capillary bridge between two identical spheres offers already a wide range of possible investigations such as its shape, the capillary force exerted by the bridge on a spherical particle, or the wetting angle and the amounts of liquid for which the bridge exists before its point of rupture [2, 3, 4]. A more complex scenario is considered in recent works [5, 6] where also the effect of unequal sized spheres is investigated.

Villanueava et. al [7] proposed a combined approach of a multicomponent and multi-phase-field model with the Navier-Stokes equations for the simulation of liquid-phase sintering, Also within the lattice Boltzmann methods [8, 9], models have been proposed to simulate the dynamics of solid particles, multiple fluids [10, 11], and more recently with the integration of liquid-gas interfaces with solid bodies [12]. Lately, Sun and Sakai [13] have performed an intensive numerical study on capillary bridges between two, three, and four spherical particles, where they used the so-called direct numerical simulation method. In the spherical case of a capillary bridge between two spherical bodies, they also studied the motion of the bodies under the action of the capillary force.

However, we are not aware of any work addressing the full multiphysics problem of liquid phase sintering, simultaneously accounting for capillary forces, rigid body motion and phase transformation kinetics at solid-liquid and solid-solid interfaces. Therefore, we present a combined approach of the so-called multi-phase-field method and the lattice Boltzmann method (Sec. 3). We show the reliability of the model, by comparing the obtained results with analytic solutions for the force and motion of two finite parallel plates connected by a cylindrical liquid bridge. The analytic solutions are introduced in Sec. 2.In addition to the comparison with the analytic solution, we show an investigation of the resolution dependence of the capillary force acting on the two plates, their motion due to the action of the capillary force and the wetting angle (see Sec. 4.2). Beyond these benchmark scenarios, we study the dynamics of two, three, and four spherical bodies under the action of the capillary force for various amounts of liquid fractions in Sec. 4.3. By doing this investigation, we show the differences in dynamics but also the similarities. We summarize the results in the Sec. 5.

2 Theory

As a test for our model, we consider the force of a single capillary bridge on two parallel plates and the resulting motion of these plates in Sec. 2.1. An introductive example is the topology of a capillary bridge between two spherical bodies which depends on the amount of liquid and the wetting angle. We discuss this dependency in Sec. 2.2.

2.1 Capillary Bridge Between Two Plates

In principle, a capillary bridge between two plates can form a variety of different surface shapes with a constant curvature which depends on the wetting angle, the distance between the plates, and the amount of liquid between the plates. These shapes can be separated into two basic classes, convex and concave bridges. In addition to this classification, analytic descriptions of these shapes can be found, such as a sphere, a cylinder, or a catenoid [14].

Here, we consider only simple examples to test our model, a planar capillary bridge with a $90^{\circ}$ wetting angle and a cylindrical bridge with the same wetting angle. For simplicity, the planar capillary bridge can be considered as a 2D problem. Nevertheless, the planar capillary bridge example is simulated in a three-dimensional setup with periodic boundaries.

This approach not only offers the possibility to calculate the resulting force on the plates (see. Sec. 2.1.1) but also allows the integration of the equation of motion in order to obtain an approximate solution for the dynamics at early times as the plates start to move toward each other (see. Sec. 2.1.2).

2.1.1 Force on the Plates

We consider two scenarios as illustrated in Fig. 1. The depicted capillary bridges are determined by their radius $R$, the distance between the plates $h$ and the wetting angle $\Theta$.

To calculate the resulting force on the plates, we use the total differential of the internal energy $U$ of such a capillary bridge which is given by

[TABLE]

. $\sigma_{{SL}}$ denotes the surface tension of the solid-liquid interface, $\sigma_{{SV}}$ of the solid-vapor interface, and $\sigma_{LV}$ of the liquid-vapor interface with the associated areas $A_{{SL}}$, $A_{{SV}}$, and $A_{{LV}}$. $p_{{V}}$ and $p_{{L}}$ are pressures of the vapor and liquid phases with the respective volumes $V_{{V}}$ and $V_{{L}}$. Because of the considered geometry, one can see that $dA_{{SV}}=-dA_{{SL}}$ and $dV_{{V}}=-dV_{{L}}$ so that Eq. (1) can be simplified

[TABLE]

, where $\Delta p=p_{{L}}-p_{{V}}$ is the pressure difference between vapor and liquid.

The equilibrium shape of the liquid minimizes the internal energy $U$. In the absence of solids, in the case of a spherical liquid droplet ($dA_{{SL}}=0$), it can be shown by minimization of the total internal energy, i.e., by setting $dU=0$, that the pressure difference between vapor and liquid is proportional to the curvature $\kappa$ of the surface,

[TABLE]

, which is called the Young-Laplace equation. We consider here the curvature as the sum of the principal curvatures $\kappa=\sum^{D-1}_{i}\kappa_{i}$ where $D$ denotes the dimension. Because a pressure gradient would result in a motion of the fluid, Eq. (3) states that the equilibrium surface of the liquid must be of constant curvature. In the presence of solid-fluid interfaces, the energy minimization principle delivers also, the well-known Young’s equation $\cos{\Theta}=\left(\sigma_{{SV}}-\sigma_{{SL}}\right)/{\sigma_{{LV}}}$, where $\Theta$ is the wetting angle.

We consider a wetting of $90^{\circ}$ so that the first term in Eq. (2) vanishes. Further, we neglect the effect of evaporation due to a change of liquid-vapor interface curvature so that the liquid volume can be considered as constant and hence the last term in Eq. (2) vanishes also. The change of the internal thus determined by the change of liquid-vapor interface area

[TABLE]

.

In the planar case (Fig. 1a), the magnitude of the force is then given by

[TABLE]

, where $L$ is the length of the planar capillary wall. Later, in the simulation, we use periodic boundary conditions. $L$ would then be the length of the simulation domain perpendicular to the depicted planar capillary bridge in Fig. 1a.

For the cylindrical capillary bridge (see Fig. 1b), one obtains

[TABLE]

. In the following, we drop the index of $\sigma_{{LV}}$ to simplify the notation. A more detailed derivation can be found in A.

2.1.2 Motion of the Plates

Based on the analytic descriptions of the forces Eqs. (6) and (5), the distance $h$ of the plates as a function of time can be calculated. This dynamic solution provides a more sophisticated benchmark for the dynamics of the present method. To reduce the complexity of the analysis, we assume that the motion of the plates is slow compared to the motion of the fluid. In this case, the bridge can be assumed to remain in its cylindrical shape during the motion of the plates. Also, the kinetic energy of the liquid can be neglected. These assumptions are valid if the mass of the plates is much larger than the mass of the liquid.

In the planar case, the computation is simple because the integral of the constant force Eq. (5) delivers

{dgroup}

[TABLE]

[TABLE]

,

where $m$ is the mass of a single plate, $h_{0}$ the initial distance, and $t$ the time. The factor of two accounts for the fact that both plates are mobile and subject to the same magnitude of the force.

In the cylindrical case, we consider the inverse function $h_{{c}}\left(t\right)$ can be approximated with

[TABLE]

where the change of the capillary radius is neglected.

2.2 Liquid Bridge Between two Spherical Solids

As an introductive scenario, we consider a liquid bridge between two spherical solids of the same size in contact with each other. Even in this simplified case, the analytic description of the bridge can be complicated. The possible equilibrium configurations of such a system are determined by the wetting angle and the liquid volume fraction $c=V_{{L}}/\left(V_{{L}}+V_{{S}}\right)$, where $V_{{S}}$ is the solid volume. The choice of these two parameters offers the possibility to discuss the topology of this system independent of its scale. For two spherical solids in contact, two ideal solutions for the topology of the liquid bridge exist, a spherical bridge in the limit of a large liquid volume fraction $(c=0.75,\Theta=0)$ and a cylindrical bridge for the other extreme case of low liquid content $(c=0.2,\Theta=0)$. Both solutions exist because they form a surface of constant curvature. By calculating the wetting angle between a spherical liquid bridge as a function of the liquid volume fraction, one obtains the relation

{dgroup}

[TABLE]

[TABLE]

.

The same way, a relationship between the wetting angle and the liquid volume fraction can be found in the case of a cylindrical bridge,

{dgroup}

[TABLE]

[TABLE]

.

A derivation of Eq. (11) and Eq. (13) can be found in B. Equation. (11) and Eq. (13) are depicted in Fig. 2.

However, more interesting are the parameter regions delimited by the spherical and cylindrical solutions. For all pairs of wetting angles and liquid volume fraction below the cylindrical solution in Fig. 2, one of the radii of curvature of the bridge becomes negative. Between the cylindrical and the spherical solution lies the region of parameters where the bridges have a convex shape but not yet a spherical one. In the region above the spherical solution, the spherical solids are not brought into contact by the liquid capillary bridge. Instead, the solids are separated by the bridge. As long as the solid surface is not entirely hydrophobic, the spherical solids will not separate from the liquid bridge. Furthermore, the liquid bridge has a spherical shape for liquid volume fractions and wetting-angles above the plotted spherical solution.

Instead of determining the shape of the liquid bridge, the liquid fraction and the wetting angle determine the distance between the spherical solids. In the following, we focus on the investigation of the parameter range below the spherical solution. To illustrate these possible shapes, pictures obtained from simulation are added to Fig. 2 according to their contact angle and liquid volume fraction. The pictures are obtained from simulations with the method described in Sec. 3.

3 Model

In the following section, we give a brief compendium of the implemented model which consists basically of a combined approach of the multi-phase-field method for the modeling of solid phases, including solid-solid as well as solid-liquid phase transformation, and the Lattice Boltzmann method for modeling the fluid flow, liquid-vapor phase separation, and solid-fluid interaction. An overview of the modeling of solids and their dynamics is shown in Sec. 3.1. The Lattice Boltzmann method is discussed in Sec. 3.2 with its coupling to the dynamics of solid bodies. In Sec. 3.3, we explore some numerical aspects and parameters.

3.1 Modeling of Solids

The phase-field method is a method for solving interfacial problems, and it has been applied to various kinds of problems such as solidification [15], grain-growth [16], surface or phase-boundary diffusion [17], and elastic deformation of solid bodies due to surface tension [18]. Here, we use the well-established multi-phase-field method which is summed up in two review articles [19, 20].

The model allows us to distinguish and track each solid particle while for example being deformed or being in the process of a phase transition. This tracking of a rigid body and its topology is achieved by a continuous indicator function called the phase-field $\phi_{\alpha}$. Given a point $\vec{x}$ in space, $\phi_{\alpha}\left(\vec{x}\right)=1$ means that the phase $\alpha$ is present at that point and $\phi_{\alpha}\left(\vec{x}\right)=0$ means that it is not. All values of $\phi_{\alpha}$ between [math] and $1$ are considered as surface or interface of the particle or phase $\alpha$ with its surroundings. The considered system may consist of $N$ phases so that all phase-fields must satisfy the sum constraint $\sum_{\alpha}^{N}\phi_{\alpha}=1$. A rigid body $B$ must consist of at least one phase-field, but it may even consist of multiple phase-fields. We consider the dynamics of the solid particles as modeled by the advection of the phase-fields according to the velocity of the solid bodies. The velocity is calculated according to forces and torques acting on the bodies (see Sec. 3.1.1). Between two colliding bodies, a repulsive force is used (see Sec. 3.1.2).

As mentioned above, the phase-field method is commonly used to describe phenomena such as phase transformations or even surface/phase-boundary diffusion which may lead to growth, shrinkage or deformation of solid particles. These processes can be considered with the present model, but they are not considered in this work, and the reader is referred to [19, 20, 17].

3.1.1 Solid Dynamics

Consider a rigid body $B$ which is subject to the force $\vec{F}_{B}$ and the torque $\vec{T}_{B}$. The center of mass $\vec{X}_{B}$ of the solid body is accelerated with

[TABLE]

, where $M_{B}$ is the mass of the solid body. The angular velocity $\vec{\omega}_{B}$ changes with

[TABLE]

, where $\underline{\underline{I}}_{B}$ is the tensor of inertia. The resulting velocity $\vec{u}_{B}$ of $B$ at the point $\vec{x}$ can be calculated with

[TABLE]

. Each of the phase-fields $\phi_{\alpha}$ associated with the rigid body $B$ is advected with

[TABLE]

. For the advection Eq. (17) of the phase-fields, we use a high-order scheme with directional splitting and a monotonized central flux limiter [21].

3.1.2 Solid-Solid Interactions

Collisions are nearly inevitable when simulating rigid bodies inside a dynamic fluid flow. Thus, we utilize a repulsive solid-solid interaction to account for collisions. This repulsive force between two rigid bodies is derived from a potential energy density $e\left(\vec{x}_{i},\vec{x}_{j}\right)$ between a point $\vec{x}_{i}$ of the rigid body $B_{1}$ and a point $\vec{x}_{i}$ of the rigid body $B_{2}$

[TABLE]

where $r_{c}$ is a characteristic interaction distance and the energy density $e_{0}$ determines the interactions strength. Furthermore, a point $\vec{x}$ is considered to be a part of the solid body $B$ if one of its associated phase-fields is non zero. The force density $\vec{\mathsf{f}}_{ij}$ of a rigid body $B_{1}$ at $\vec{x}_{j}$ acting on a rigid body $B_{2}$ at $\vec{x}_{i}$ can obtain as

{dgroup}

[TABLE]

[TABLE]

.

By the summation over all pairs of $\vec{x}_{i}$ and $\vec{x}_{j}$ within the characteristic interaction length, the resulting total forces on the rigid bodies can be calculated. To keep the solid-solid interaction as short ranged as possible, we have chosen $r_{c}=3\,\Delta x$ and $n=4$ in our simulations where $\Delta x$ is the lattice spacing.

3.2 Modeling of Fluids

We use for the simulation of fluids the lattice Boltzmann method which is summed up in Sec. 3.2.1. The interaction of the fluid with the solid bodies is shown in Sec. 3.2.2. In Sec. 3.2.3, we introduce the implemented method for the liquid-vapor phase separation.

3.2.1 The Lattice Boltzmann Method

Consider the distribution function $f\left(\vec{x},\vec{v},t\right)$ which is the probability to find a pseudo particle at the position $\vec{x}$ with the physical velocity $\vec{v}$ at the time $t$ [22]. The continuous velocity space $\vec{v}$ is discretized with a finite number of physical velocities $\vec{c_{i}}$ of the pseudo particles. Although many schemes have been proposed for the discretization of the velocity space, we only consider here the so-called three-dimensional twenty-seven velocity model (D3Q27) with

[TABLE]

where $c=\Delta x/\Delta t$ is a shorthand and $\Delta t$ is the time step. The fluid density $\rho$ is calculated as the sum of the distribution function overall velocities $\rho=\sum_{i}f_{i}$ where the shorthand $f_{i}\left(\vec{x},t\right)\equiv f\left(\vec{x},\vec{c}_{i},t\right)$ has been used. Likewise, the fluid velocity $\vec{u}$ is calculated as average overall physical velocities $\vec{c_{i}}$ weighted by the amount of pseudo particles $f_{i}$ with $\vec{u}=1/{\rho}\sum_{i}\vec{c_{i}}f_{i}$. This way, the lattice Boltzmann equations can be written as

{dgroup}

[TABLE]

[TABLE]

which describe the collision step Eq. (22) and the streaming step Eq. (23) of the pseudoparticles, $f_{i}$. The right-hand side of Eq. (22) consists of two parts, the relaxation of the distribution function towards the local equilibrium distribution $f_{i}^{{eq}}$ and the forcing term $F_{i}$. The local equilibrium distribution function can be calculated with

[TABLE]

, which is a second-order expansion of the local Maxwell-Boltzmann distribution with the weights

[TABLE]

. $\tau$ is the relaxation time and $c_{s}$ the speed of sound in the lattice Boltzmann fluid with $c_{s}=c/\sqrt{3}$.

We use the forcing term as published in [23] with

[TABLE]

, where $\vec{\mathsf{f}}$ is the local body force density acting on the fluid and $\vec{v}$ is the fluid velocity which can be calculated with

[TABLE]

. We apply body force densities originating from the solid-fluid interaction (Sec. 3.2.2) and from non-ideal fluid contributions (Sec. 3.2.3) which cause a phase separation between the liquid and vapor phase.

3.2.2 Solid-Fluid Interaction

A rigid body and a fluid exert force on each other due to their motion and inertia. In the following, we introduce so-called “bounce-back” method [24] which has been implemented.

To use this solid-fluid interaction method a point in the simulation domain has to be either solid or liquid. Because the phase-field method uses a continuous transition between two phases, a point is considered to be solid if the phase-field representing the fluid phase content is less then $0.05$. If a point and its adjacent points are considered to be fluid, the lattice method is applied (see Sec.3.2.1). Near a solid surface, the bounce-back adds to the right-hand side of Eq. (23) and accounts for fluid which “bounces back” form the solid surface

[TABLE]

If no point $\vec{x}-\vec{c_{i}}\Delta{t}$ in the vicinity of $\vec{x}$ is solid, the bounce-back is not calculated so that Eq. (23) and Eq. (28) deliver the same result. However, if the point $\vec{x}-\vec{c_{i}}\Delta t$ lies within a solid body, there is no fluid population that can be streamed to point $\vec{x}$. Hence, the pseudo particles coming from the point $f_{i}\left(\vec{x}\right)$ with the velocity $\vec{c_{i}}$ are assumed to “bounce-back” from the solid surface in the middle between the two points. Consequently, $f^{*}_{-i}\left(\vec{x},t\right)$ denotes the population of the pseudo particles with the velocity $-\vec{c}_{i}$. If a solid surface is present in the vicinity, the corresponding part of the distribution function is reflected. The solid surface is here assumed to be in the middle between the fluid and the rigid body node at $\vec{x}-\vec{c_{i}}\Delta t/2$. This adds to the consequence the remaining fluid has to be distributed to the neighboring points if a certain point changes from liquid to solid.

Further, we consider the movement of the solid body with the velocity $\vec{u}_{B}$. The movement of the solid body leads to a change of momentum of the adjacent fluid which is accounted for by

[TABLE]

in Eq. (28). Whereas $\delta\rho\left(\vec{x}\right)$ is an auxiliary term that ensures the local mass conservation with

[TABLE]

, where the sum over $i$ is executed for each point $\vec{x}-\vec{c}_{i}\Delta t$ that is solid so that the total change of mass caused by the first term in Eq. (29) is compensated. This mass correction $\delta\rho\left(\vec{x}\right)$ is then added to the resting fluid population with the term $\delta_{i0}$ which is only non-zero for $i=0$ so that no slip of the fluid along the solid surface is introduced.

The collision and bounce back of the fluid at the solid surface results in a momentum exchange between solid and fluid. The resulting local change of the momentum density $\Delta\vec{m}_{B}$ of the rigid body at a solid point $\vec{x}$ is given by

[TABLE]

where the sum is executed for each point $\vec{x}+\vec{c}_{i}\Delta t$ that is fluid. By integration over the entire solid body $B$, corresponding total force $\vec{F}_{B}$ and torque $\vec{T}_{B}$ can then be obtained with

{dgroup}

[TABLE]

[TABLE]

,

where $\vec{X}_{B}$ denotes the center of mass of the solid body (see Sec. 3.1.1).

In addition to the introduced bounce-back, our simulations reveal that yet another solid-fluid interaction is necessary to conserve the fluid mass and momentum with a rigid body moving through it. Consider a moving solid body in the discrete simulation space, a former fluid node at the point $\vec{x}_{n}$ may become a solid node because a rigid body is passing by. The remaining fluid distribution function $f_{i}\left(\vec{x}_{n}\right)$ has to be distributed to the neighboring fluid points so that the fluid mass and momentum are conserved. This conservation can be achieved with

[TABLE]

, where $N_{F}\left(\vec{x}_{n}\right)$ is the set of fluid neighbor nodes in a sphere of the radius $R_{A}$ around $\vec{x}_{n}$ with $N_{F}\left(\vec{x_{n}}\right)={\left.\left\{\vec{x}\,\right|\left\|\vec{x_{n}}-\vec{x}\right\|\leq R_{A}\textrm{~{}and }\vec{x}\neq\vec{x_{n}}\right\}}.$ In our simulations, the choice of $R_{A}=\Delta x$ has shown to be sufficient.

In the case that a former solid node becomes fluid, it is sufficient to set the fluid density to zero.

3.2.3 Non-Ideal Fluids

In the following, we briefly present the implemented method to simulate liquid and vapor phase separation and their interaction with solid surfaces. For more detailed information, the reader is referred to [25, 26, 22, 27]. The liquid and vapor phase separation is introduced by a force density $\vec{\mathsf{f}}$ acting on the fluid (see Eq. (23) and (26)) with

[TABLE]

, where $\mathcal{G}$ is a coupling constant and $\psi$ a so-called lattice version of the mean-field potential,

[TABLE]

, where $\rho_{0}$ is a reference density. The force density Eq. (35) can be interpreted as the sum of force densities which act between the pseudo particles at $\vec{x}$ and the surrounding pseudo particles at $\vec{x}+\vec{c_{i}}\Delta{}t$. Within this model, a convenient way to account for wetting effects is to identify a solid node as an immobile “fluid” with a fictions density $\varrho_{B}$. By the integration over all force densities between a fluid point $\vec{x}$ and an adjacent point $\vec{x}+\vec{c_{i}}\Delta{}t$ of the rigid body $B$, the resulting force $\vec{F}_{B}$ and torque $\vec{T}_{B}$ acting on it can then be obtained

{dgroup}

[TABLE]

[TABLE]

where again the sum over $i$ is executed for each point $\vec{x}+\vec{c}_{i}\Delta t$ that is solid.

With $\vec{\mathsf{f}}=-\nabla\cdot\underline{\underline{p}}$ it is possible to identify the corresponding pressure tensor

[TABLE]

. For the bulk pressure or the equation of states, one thus has

[TABLE]

. The pressure tensor can be used to calculate the surface tension. According to [28], the interface tension $\sigma$ is given by the difference of the normal component of the pressure tensor $p_{\bot}$ and its tangential component $p_{\parallel}$ across the interface with

[TABLE]

, where $n$ is the direction normal to the interface. If one considers only a planar interface, so that $\psi$ changes only along the normal direction of the interface, the interface tension can be calculated with

[TABLE]

, which is obtained by inserting Eq. (39) into Eq. (41).

Considering Eq. (42) and Eq. (36), one can see that the interface tension depends on the respective density profile. The density profile across the liquid-vapor interface is determined by the coupling constant $\mathcal{G}$. The density profile across the solid-vapor and the solid-liquid can be adjusted tuning the above introduced “solid density” $\varrho_{B}$. This way, the contact angle of the solid body can be adjusted, and since the term “solid density” is misleading in this context, we refer to $\varrho_{B}$ as the wetting parameter. For a more detailed explanation, we refer the reader to [27].

3.3 Simulation Procedure

As mention in Sec. 3.1.1, we use a high-resolution scheme for the advection, Eq. (17), of the phase-fields with directional splitting and a monotonized central flux limiter [21]. This higher-order scheme serves to reduce the well-known numerical problem of diffuse interface spreading under advection. An advantage of the multi-phase-field method proposed here is that accounting for phase-field kinetics stabilizes the interface profile further, thus improving numerical stability. As known from the standard literature on multi-phase-field [19, 20], the kinetics of $\phi$ read,

{dgroup}

[TABLE]

[TABLE]

where $M_{\alpha\beta}$ is the mobility of the interface between phases $\alpha$ and $\beta$, $\sigma_{\alpha\beta}$ is the interface energy between the phase-fields and $\eta$ is the interface width (see [20] for more details).

The restoration of the phase-field profile comes with a trade-off. Eq. (43) restores the phase-field profile but also introduces a shrinkage of the phase-field

[TABLE]

according to the local curvature. Because of this shrinkage, the last term $\Delta g_{\alpha\beta}$ is necessary which is dynamically determined to conserve the volume of a phase-field $\phi_{\alpha}$ with

[TABLE]

. Similar as in Sec. 3.1.2, $\left|\delta\phi_{\alpha}\right|$ denotes the number of nodes in the interface region of $\phi_{\alpha}$ with $\delta\phi_{\alpha}=\left.\left\{\vec{x}\right|0.0<{\phi\left(\vec{x}\right)}<1.0\right\}$. To minimize this solid-phase sintering, we chose the phase-field mobility as small as reasonable.

If not stated otherwise, the following values of the introduced parameters have been used for the present simulations:

4 Results

We first investigate the accuracy of the surface tension calculation in Sec. 4.1. This investigation is the basis for the testing of the present model. In Sec. 4.2, we test our model by the simulation of capillary bridges between two finite parallel plates. The results obtained from these simulations are compared with analytic solutions given in section 2.1.1. To demonstrate the capability to model the process of liquid-phase-sintering, we also investigate a more complicated situation in Sec. 4.3, the dynamics of multiple particles connected by a capillary bridge. This investigation provides a first step towards the modeling of larger and more complex systems with more solid particles.

4.1 Surface Tension

In the following, we test the consistency of the surface tension calculations Eq. (42) with the Laplace pressure Eq. (3). By evaluating Eq. (42) for droplets of different radii, we obtained the value $3.297\cdot 10^{-2}$ in lattice units for the surface tension. The value is used to predict the pressure gradient, by inserting it into Eq. (3). We show this prediction in Fig. 3. Also, the differences for different radii, obtained with Eq. (40), are depicted. One can see that the obtained pressure gradients are in agreement with the ones predicted by the Laplace pressure, Eq. (3). A deviation from the predicted value can only be seen for large curvatures, which is expected because the spatial resolution decreases.

Another possibility to obtain the surface tension is to calculate the pressure gradient with Eq. (40) and use it together with the radius as input for Eq. (3) to calculate the resulting surface tension. This way, we yield $3.583\cdot 10^{-2}$ for the value of surface tension which differs from the value obtained with Eq. (42) by approximately $8\%$. In the following chapter, we show that this value for the surface tension provides a better estimate for the acting capillary force on a solid body. Hence, $\sigma=3.583\cdot 10^{-2}$ is used in the following.

4.2 Capillary Bridge Between Two Plates

We continue the testing of the present model by investigating a capillary bridge between two finite parallel plates as schematically depicted in Fig. 1. The choice of finite plates ensures that the ambient pressure on both sides of the plate is the same. We use the analytic solution to approximate the resulting motion of the plates and compare them with the simulation results.

In Sec. 2.1.2, we introduced the approximative solution Eq. (9) for the motion of the plates connected by a cylindrical capillary bridge and the analytic solution Eq. (8) for the simple planar case. These solutions are a good approximation if the solid motion is slow and the fluid is near its equilibrium state. To realize this assumption in the simulation, we assign to the solid a high mass density $\rho_{S}$ as compared to the liquid mass density $\rho_{L}$ with $\rho_{S}/\rho_{L}=1.72\cdot 10^{4}$. The initial simulation setup is schematically depicted for the planar setup in Fig. 1a and for the cylindrical setup in Fig. 1b. Figure 4 shows the decreasing distance of the plates over time due to the action of the capillary force. The distance of the plates has been rescaled by the length $R$ which is the initial radius of the cylindrical capillary bridge. Furthermore, we have rescaled the time by a characteristic time $t_{c}$ with $t_{c}=\sqrt{m/\sigma}$ where $m$ is the mass of the solid plates. By using $R$ and $t$, one can rewrite Eqs. (8) and (9)

{dgroup}

[TABLE]

[TABLE]

.

To obtain the rescaled equations of motion Eq. (48), we have assumed for the length of the planar capillary bridge that $L=R$. The simulation results are depicted in Fig. 4 together with the rescaled equations of motion Eq. (48). We also show as insets in Fig. 4 the particular simulation setups at $t=0$ and at the end of the displayed time interval. A deviation of the simulations from the predictions is expected because the fluid motion has been neglected in the derivation of Eq. (48). Consequently, we show only the very beginning of the motion in Fig. 4 so that the theoretical prediction is accurate and can be used to test the simulations.

One can see that in both cases, the planar one in Fig. 4a and the cylindrical one Fig. 4b, the results agree with the predictions of Eqs. (47) and (48) in the beginning. Nevertheless, for later times the simulation results deviate from the predictions. This effect is more significant in the cylindrical case than in the planar one.

4.3 Dynamics of Multiple Particles

4.3.1 Two Particles Connected by a Capillary Bridge

We continue our investigation with two spherical particles connected by a capillary bridge with a particle radius of $R=30$, mass density of $\rho=4$, and liquid fraction of $c=0.5$ (see Fig. 5). The closest distance between the particles $h$ and the velocity $\dot{h}$ are shown in Fig. 5. Because we want to consider solids with a very high wettability, a wetting-parameter of $\varrho_{B}=2.0$ has been chosen. Above in section 3.3, we have introduced the phase-field dynamics Eq. (43) which models the physical process of sintering. To show the influence of the phase-field dynamics, the simulation has been performed twice, once with the phase-field dynamics and once without it where the phase-fields are only advected.

In Fig. 5, the closest distance between the surfaces of the particles $h$ and the corresponding velocity $\dot{h}$ are shown. The dashed lines in Fig 5 indicate the simulations where the phase-fields are only advected. Without the phase-field dynamics, it takes more time for the two particles to come into contact with each other. The reason for this can be seen when looking at the velocity profile without the phase-field dynamics. In the beginning, the two particles are strongly accelerated towards each other, until about two thousand time steps when the velocity reaches its peak value (velocity minimum of the dashed line in Fig. 5). After two-thousand time-steps, the two particles are continuously decelerated until the particles come into contact. By looking at Fig. 2, one can see that there exists an equilibrium solution for $c=0.5$ and low wetting angles where the particles are in contact. The existence of this equilibrium solution means that the capillary force must be attractive during the whole simulation time. Furthermore, the phase-fields are purely advected, and the solid-solid interaction has only a range of 3. Hence, the observed deceleration in Fig. 5 can only be caused by the dynamics of the fluid. We do not consider an explicit drag force here. Consequently, the bounce-back effect Eq. (28) can only cause this deceleration.

When the phase-fields are not only advected (solid lines in Fig 5), one can see that the time, until the two particles come into contact, is reduced. Even more visible is the influence of the phase-field dynamics in the velocity profile. At about $t=3000$, the magnitude of the velocity increases strongly, and a short time after that the velocity drops to nearly zero. This strong acceleration occurs at about a distance of $10$ and can be understood in the context of the phase-field method. The value of the present phase-fields varies continuously from [math] to $1$ within a distance $\eta=5$. This defines the interface thickness. Hence, the phase-fields begin to touch each other about a distance of $2\eta=10$. At this distance, the phase-fields start to merge and form a common interface between the two solid particles.

When $\eta$ is comparable to the width of a grain-boundary, the length-scale can be seen as physically meaningful, and the described process of merging can be regarded as the first step in the solid-phase sintering process. In contrast, if $\eta$ is larger than a physical grain boundary width, the dynamics of the phase-field and the resulting variation of the distance below $2\eta$ shall be viewed as a numerical feature to model solid-phase sintering of larger scales.

At a distance below $2\eta$, the phase-fields representing the particles start to deform. This occurs due to the solid-solid sintering process, where the particles reduce the surface area at the expense of a larger grain-boundary area. (The grain-boundary energy is assumed to be smaller than the surface energy). During this process, the centers of mass come closer together, as compared to hard and inert bodies at contact. This feature shows itself as a small negative distance in Fig. 5.

For more details about the modeling of surface diffusion with phase-field and the evolution of stress in two spherical particles undergoing the process of solid-phase sintering, the reader is referred to [17, 29, 18].

4.3.2 Effect of Liquid Fraction

Above, we investigated the dynamics of a capillary bridge between two spherical particles. A more intriguing investigation, however, is the influence of the liquid fraction on the dynamics of the system. The influence of the liquid fraction on the topology of a liquid bridge between two spheres has been shown above in Sec. 2.2.

We start with two spherical particles with a radius of $R=30$ and a mass density of $\rho=4$ which are initialized with different amounts of liquid between them. We continue with the investigation of three and four particles, as a forecast of a system with many particles. Similar to section 4.3.1, before the particles are allowed to move, the simulations run for fifteen thousand simulation time steps so that the liquid reaches its equilibrium shape. The time $t=0$ refers to the point in time when the particles are allowed to move. Furthermore, we investigate here three different liquid fraction $c=0.2$, $c=0.3$, and $c=0.5$. The simulation setup for $c=0.2$ at $t=0$ and $t=5000$ are depicted as insets in the corresponding plots in Fig. 6. For a liquid fraction of $c=0.5$, we show again with a dashed line for a simulation where the phase-field dynamics, Eq. (43), is switched off so that the phase-fields are purely advected.

We start the discussion with the concentration dependency of the distance and velocity profiles of two particles. In the case of a small amount of liquid ($c=0.2$), the velocity of the particles increases linearly with time until they come into contact. For a liquid fraction of $c=0.3$, one can see an increase in the velocity of the particles at the beginning, compared to the case of $c=0.2$. By looking at the velocity profile of $c=0.5$, the particles are more strongly accelerated at the beginning than for $c=0.2$, but for later times one can see a clear deceleration. To show this more clearly, we also plot the velocity profile corresponding to the case where the phase-fields are purely advected (dashed-line). By comparing the dashed and the solid lines of the velocity profile, one can see that the strong acceleration followed by an even stronger deceleration is caused by the dynamics of the phase-field. This process of deceleration for larger amounts of liquids, $c=0.2$ and $c=0.3$, can be explained as the effect of drag (see Sec. 4.3.1). In Fig. 6, one can see that the effect of drag is more significant for larger amounts of liquid, as to be expected.

By looking at the distance or velocity profiles in Fig. 6, one can easily identify the time when the particles are in contact. In the velocity profile, this point in time is indicated by the steep drop of the velocity to zero. The contact time is the shortest for a liquid fraction of $c=0.3$ and the largest for $c=0.5$. This observation is also true for the three and four-particle systems.

In the second setup, we placed three spherical particles, with the same radius and distance as in the first case, in such a way that they form an equilateral triangle. Different amounts of liquid are initialized in the center of the triangle so that the resulting liquid volume fractions are the same as in the first case. Fig. 6 shows the decreasing average distance $\tilde{h}=1/3\sum_{i}h_{i}$ between the three particles. Because the particles form an equilateral triangle, the force on each particle is the same, and hence the particles form an equilateral triangle during the entire simulation. This means that the average distance is equal to the distance between the particle $\tilde{h}=h_{2}=h_{3}=h_{3}$.

The results we obtain for three particles are qualitatively similar to the results obtained for two particles. However, one can see that the time, until the particles come into contact, is reduced for $c=0.2$ and $c=0.3$. In the case of $c=0.5$, the effect of drag is stronger than in the case of two particles.

In the last setup, we placed four spherical particles, with the same radius and distance as in the first case, in such a way that they form a tetrahedron. Again, we have initialized spherical liquid drops of different sizes in the center of the tetrahedron and waited for ten thousand simulation time steps so that the liquid reaches its equilibrium before the particles are allowed to move ($t=0$). Likewise, Fig. 6 shows the decreasing average distance $\tilde{h}=1/4\sum_{i}h_{i}$ between the particles for different liquid volume fractions.

The results we obtain for four particles are qualitatively similar to the results obtained for two and three particles. The time until the particles come into contact is reduced for $c=0.2$ and $c=0.3$ compared to two particles but is comparable to the case of three particles. Furthermore, the above-discussed effect of drag is even stronger than in the case of three particles. The fact that the effect of drag increases with the number of particles can be explained by the higher amount of solid surface in contact with the liquid which leads to a higher drag force.

5 Summary and Conclusion

We present a combined phase-field-lattice Boltzmann model for the simulation of liquid state sintering. The accuracy of the present method is carefully investigated by considering a liquid bridge between two parallel plates of finite size. We show that the obtained results for their motion are in agreement with the present theoretical predictions. The capability of the present model to simulate liquid state sintering is demonstrated by the investigation of the dynamics of two, three, and four solid particles under the action of a capillary attraction between them. In all the cases investigated, increasing the amount of liquid first accelerates the compaction process. A liquid fraction higher than a certain threshold, however, slows down the motion of spheres and leads to an increase of the time necessary for the particles to come into contact. We show that this is caused by the dynamics of the fluid, i.e. viscous drag. There is thus an optimum choice for liquid content concerning the compaction process. The exact value of this optimum parameter will depend in general on the powder packing fraction, grain shape, and size distribution. The present method provides a versatile tool to explore this important issue.

Acknowledgments

Financial support by the German Research Foundation DFG under the grant VA205/17-1 is gratefully acknowledged.

Appendix A Liquid Bridge Between two Plates — Force

We consider a cylindrical capillary between two finite plates. As shown above, the change of internal energy is given by $\mathrm{d}U=\sigma\,\mathrm{d}A$ with $A=2\pi R$ one obtains $\mathrm{d}U=2\pi\sigma\,\mathrm{d}R$ Further, we consider the energy of the liquid phase $U_{{L}}$ as a function of the distance $h_{c}$ between the two plates. Because the liquid volume is constant $V_{L}=\pi R_{0}^{2}h_{0}=\pi R^{2}h_{c}$, the radius $R$ can be expressed as function of the initial radius $R_{0}$, initial distance $h_{0}$ and distance $h_{c}$ with $R=R_{0}\sqrt{h_{0}/h_{c}}$. Hence, the force is obtained as

[TABLE]

.

Appendix B Liquid Bridge Between two Solid Spheres

B.1 Spherical Bridge

Equation (11) describes the necessary wetting angle to form spherical capillary bridge as a function of the liquid fraction. To derive this relation, we consider two spherical solid bodies in contact, as depicted in Fig. 7.

B.1.1 Small Spherical Capillary Bridge

The contact angle $\Theta_{{S}}$ is defined as the angle between the tangent lines of the liquid and solid surface at the triple point. It is thereby equal to the angle between the surface normals vectors of the liquid and solid surface at the triple point and can be calculated with

[TABLE]

, where $\theta_{1}$ and $\theta_{2}$ are the angles between the surface normals and the depicted dashed line (see. Fig. 7). These angles can be calculated by

{dgroup}

[TABLE]

[TABLE]

,

where $\chi_{{S}}\in\left[0,2\right]$ is a dimensionless scaling factor. We introduced this scaling factor in order to relate the radius of the capillary bridge $R_{{L}}$ to the radius of the solid spheres $R_{{S}}$ with $R_{{L}}=\chi_{{S}}R_{{S}}$. Again, one can obtain a system of equations for the side lengths of the depicted triangles by using the Pythagorean theorem,

{dgroup}

[TABLE]

[TABLE]

[TABLE]

,

with the solution

{dgroup}

[TABLE]

[TABLE]

[TABLE]

.

In Eq. (58), we expressed the geometric parameters depicted in Fig. 7 as function of the radius $R_{{S}}$ and the scaling parameter $\chi_{{S}}$. As to be expected, the contact angle can be expressed as a function of $\chi_{{S}}$ with

[TABLE]

.

Furthermore, we replace the scaling factor $\chi_{{S}}$ by the liquid volume fraction $c$. Therefore, the volume of the liquid has to be calculated by subtracting the volume of the spherical caps $V_{{scA}}$ and $V_{{scB}}$ (see Fig. 7) form the volume of a liquid sphere $V_{{Ls}}$ with

[TABLE]

and

{dgroup}

[TABLE]

[TABLE]

[TABLE]

and Eq. (58). The volumes of the spherical caps, $V_{{scA}}$ and $V_{{scB}}$, can be parametrized with $R_{{S}}$ and $\chi_{{S}}$ by inserting Eq. (56) into Eq. (62) and Eq. (57) into Eq. (63) so that

{dgroup}

[TABLE]

[TABLE]

.

Inserting Eqs. (64), (65), and (61) into Eqs. (60) delivers a simple parametrization for the liquid volume

[TABLE]

. The same way as for the cylindrical bridge, we can calculate the liquid volume fraction with

{dgroup}

[TABLE]

[TABLE]

.

By rearranging the terms of Eq. (68), we obtain a fourth order polynomial equation for $\chi_{{S}}$

[TABLE]

, which has a real and positive solution

[TABLE]

. Inserting Eq. (70) into Eq. (59) delivers the solution shown in Fig. 2.

B.1.2 Large Spherical Capillary Bridge

In the case of relatively large amount of liquid, a liquid bridge as depicted in Fig. 8 may form.

The angles $\theta_{1}$ and $\theta_{2}$ are obtained in the same manner as above with

{dgroup}

[TABLE]

[TABLE]

[TABLE]

,

where $\theta_{1}$ can be simplified to

[TABLE]

. The height $h_{{S}}$ can be calculated in the same manner as in B.1.1, and Eq. (59) can be obtained again.

B.2 Cylindrical bridge

Equation (13) describes the wetting angle necessary to form a cylindrical capillary bridge as a function of the liquid fraction. To derive this relation, we consider two spherical solid bodies in contact, as depicted in Fig. 9.

To calculate the contact angle $\Theta_{{c}}$, we consider the angle $\theta$ between two lines. The first line starts at the center of the sphere $A$ ends in the triple point $B$. The second line also starts at the center of the sphere $A$ but ends in the contact point $C$ between the spheres so that

$\theta=\arccos{\left(\frac{R-h}{R}\right)}$, where $h$ is half of the height of the cylindrical capillary bridge. By considering the intercept theorem, we obtain the contact angle

[TABLE]

, where we have introduced the dimensionless scaling parameter $\chi_{{c}}\in\left[0,1\right]$ with $h=\chi_{{c}}R$.

The volume of the liquid can be determined by calculating the volume of a cylinder and subtraction of the solid sphere caps, which cover part of the cylinder’s volume,

{dgroup}

[TABLE]

[TABLE]

[TABLE]

.

By using the Pythagorean theorem $R^{2}=R_{{S}}^{2}+\left(R-h\right)^{2}$, the radius of the capillary bridge $R_{{S}}$ in Eq. (78) can be replaced by the radius of the solid sphere $R_{{S}}$. With this replacement, we obtain for the fluid volume

[TABLE]

. Furthermore, we replace height $h$ again with $h=\chi_{{c}}R$ so that the expression for the fluid volume

[TABLE]

can be further simplified. With Eq. (80) and the volume of the solid, $V_{{S}}=\frac{8}{3}\pi R^{3}$, the liquid volume fraction can be calculated

{dgroup}

[TABLE]

[TABLE]

.

Equation (82) can be transformed in a third order polynomial

[TABLE]

. The polynomial equation (83) has three real roots for $c\in\left(0,0.2\right)$. By substituting $\chi_{{c}}=t+\frac{1}{2}$, we can write Eq. (83) in the form,

[TABLE]

with $p=-\frac{1}{4}$ and $q=\frac{2c}{1-c}-\frac{1}{4}$. With the use of the cosine power rule

[TABLE]

, we find three real solutions of Eq. (84)

[TABLE]

with $n\in\left\{0,1,2\right\}$. By back substitution, we also obtain three real solutions for the scaling factor

[TABLE]

. The physically meaningful solution of $\chi_{{c}}$ has to be limited to the interval $\left[0,1\right]$, which is only fulfilled by $\chi_{{c}2}$. Inserting Eq. (87) into Eq. (75) delivers the result shown in Fig. 2.