Simulation of Capillary-Driven Kinetics with Multi-Phase-Field and Lattice-Boltzmann Method
Raphael Schiedung Marvin Tegeler Dmitry Medvedev Fathollah Varnik


unsrt






 
 



compat=newest


dmath
dmath



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

^1Ruhr Universität Bochum, Interdisciplinary Center for Advanced Materials Simulation (ICAMS), Universitätsstr. 150, 44801 Bochum, Germany
^2Lavrentyev Institute of Hydrodynamics SB RAS, Lavrentyev prosp. 15, 630090, Novosibirsk, Russia
^2Novosibirsk State University, Pirogova str. 2, 630090, Novosibirsk, Russia
raphael.schiedung@rub.de, fathollah.varnik@rub.de

    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.




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







§ 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 <cit.>.
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 <cit.>.
A more complex scenario is considered in recent works <cit.> where also the effect of unequal sized spheres is investigated.

Villanueava et. al <cit.> 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 <cit.>, models have been proposed to simulate the dynamics of solid particles, multiple fluids <cit.>, and more recently with the integration of liquid-gas interfaces with solid bodies <cit.>.
Lately, Sun and Sakai <cit.> 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. <ref>).
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. <ref>.
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. <ref>).
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. <ref>.
By doing this investigation, we show the differences in dynamics but also the similarities.
We summarize the results in the Sec. <ref>.



§ 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. <ref>.
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. <ref>.



 §.§ 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 <cit.>.

Here, we consider only simple examples to test our model, a planar capillary bridge with a 90^∘ 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. <ref>) 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. <ref>).



  §.§.§ Force on the Plates


We consider two scenarios as illustrated in Fig. <ref>.
The depicted capillary bridges are determined by their radius R, the distance between the plates h and the wetting angle Θ.

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


    dU =
        σ_SL  dA_SL+
        σ_SV  dA_SV+
        σ_LV  dA_LV-
        p_V  dV_V-
        p_L  dV_L
.
σ_SL denotes the surface tension of the solid-liquid interface, σ_SV of the solid-vapor interface, and σ_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. (<ref>) can be simplified


    dU
        =
        - (σ_SV - σ_SL)   dA_SL+
        σ_LV  dA_LV -
        Δ p   dV_L
,
where Δ 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 d U = 0, that the pressure difference between vapor and liquid is proportional to the curvature κ of the surface,


    Δp = σ_LVκ
,
which is called the Young-Laplace equation.
We consider here the curvature as the sum of the principal curvatures κ = ∑^D-1_i κ_i where D denotes the dimension.
Because a pressure gradient would result in a motion of the fluid, Eq. (<ref>) 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Θ
    = (σ_SV-σ_SL)/σ_LV
,
where Θ is the wetting angle.

We consider a wetting of 90^∘ so that the first term in Eq. (<ref>) 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. (<ref>) vanishes also.
The change of the internal thus determined by the change of liquid-vapor interface area


    dU
        =
        σ_LV  dA_LV
.

In the planar case (Fig. <ref>), the magnitude of the force is then given by


    F_p = 2 σ_LV L
,
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. <ref>.

For the cylindrical capillary bridge (see Fig. <ref>), one obtains


    F_c = πσ_LV R
.
In the following, we drop the index of σ_LV to simplify the notation.
A more detailed derivation can be found in <ref>.



  §.§.§ Motion of the Plates


Based on the analytic descriptions of the forces Eqs. (<ref>) and (<ref>), 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. (<ref>) delivers



    h_p = h_0 - 1/2F_p/m t^2 - 1/2F_p/m t^2


    h_p= h_0 - 2σ L/m t^2
,

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(t) can be approximated with


    h_c(t) ≈ h_0 - πσ R_0/m t^2

where the change of the capillary radius is neglected.



 §.§ 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/(V_L + V_S)
,
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,Θ = 0) and a cylindrical bridge for the other extreme case of low liquid content (c = 0.2,Θ = 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


    
    Θ_S
            = arccos(χ_s^2-2/2)
            - arccos(χ_s/2)

    
    χ_s = 2/√(3)√((1-c)√(c(1-c)))/1-c
.
    


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


    
    Θ_c
            =
            π/2 - arccos(1 - χ_c)

    
    χ_c
            =
            - cos(
                1/3 arccos(1-9 c/1-c)
                +π/3)
            + 1/2
.
    

A derivation of Eq. (<ref>) and Eq. (<ref>) can be found in <ref>.
Equation. (<ref>) and Eq. (<ref>) are depicted in Fig. <ref>.



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. <ref>, 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. <ref> according to their contact angle and liquid volume fraction.
The pictures are obtained from simulations with the method described in Sec. <ref>.



§ 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. <ref>.
The Lattice Boltzmann method is discussed in Sec. <ref> with its coupling to the dynamics of solid bodies.
In Sec. <ref>, we explore some numerical aspects and parameters.



 §.§ 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 <cit.>, grain-growth <cit.>, surface or phase-boundary diffusion <cit.>, and elastic deformation of solid bodies due to surface tension <cit.>.
Here, we use the well-established multi-phase-field method which is summed up in two review articles <cit.>.

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 ϕ_α.
Given a point x⃗ in space,

    ϕ_α (x⃗) = 1

means that the phase α is present at that point and

    ϕ_α (x⃗) = 0

means that it is not.
All values of ϕ_α between 0 and 1 are considered as surface or interface of the particle or phase α with its surroundings.
The considered system may consist of N phases so that all phase-fields must satisfy the sum constraint

    ∑_α^N ϕ_α = 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. <ref>).
Between two colliding bodies, a repulsive force is used (see Sec. <ref>).

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 <cit.>.



  §.§.§ Solid Dynamics


Consider a rigid body B which is subject to the force F⃗_B and the torque T⃗_B.
The center of mass X⃗_B of the solid body is accelerated with


    Ẍ⃗̈_B = M^-1_BF⃗_B
,
where M_B is the mass of the solid body.
The angular velocity ω⃗_B changes with


    ω̇⃗̇_B
        =
        I^-1_BT⃗_B
,
where I_B is the tensor of inertia.
The resulting velocity u⃗_B of B at the point x⃗ can be calculated with


    u⃗_B(x⃗)
        =
        Ẋ⃗̇_B + ω⃗_B×(x⃗ - X⃗_B)
.
Each of the phase-fields ϕ_α associated with the rigid body B is advected with


    ϕ̇_α
        =
        u⃗_B·∇ϕ_α
.
For the advection Eq. (<ref>) of the phase-fields, we use a high-order scheme with directional splitting and a monotonized central flux limiter <cit.>.



  §.§.§ 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(x⃗_i,x⃗_j) between a point x⃗_i of the rigid body B_1 and a point x⃗_i of the rigid body B_2


    e(x⃗_i,x⃗_j)
        =
        e_0
        {[ |‖x⃗_i-x⃗_j‖-r_c/r_c|^n if ‖x⃗_i - x⃗_j ‖ < r_c;                       0                    else ].

where r_c is a characteristic interaction distance and the energy density e_0 determines the interactions strength.
Furthermore, a point 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 𝖿⃗_ij of a rigid body B_1 at x⃗_j acting on a rigid body B_2 at x⃗_i can obtain as



    𝖿⃗(x⃗_i,x⃗_j )
        = ∂ e(x⃗_i,x⃗_j )/∂(x⃗_i-x⃗_j)


    =
        e_0
        n/r_cx⃗_i-x⃗_j/‖x⃗_i-x⃗_j‖|‖x⃗_i-x⃗_j‖-r_c/r_c|^n-1
.

By the summation over all pairs of x⃗_i and 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 Δ x and n=4 in our simulations where Δ x is the lattice spacing.



 §.§ Modeling of Fluids


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



  §.§.§ The Lattice Boltzmann Method


Consider the distribution function f ( x⃗, v⃗, t ) which is the probability to find a pseudo particle at the position x⃗ with the physical velocity v⃗ at the time t <cit.>.
The continuous velocity space v⃗ is discretized with a finite number of physical velocities 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


    c⃗_i
        =
        {[                               (0,0,0)                                 i = 0;       (± c,0,0), (0,± c,0), (0,0,± c)                             i = 1,…,6; (± c,± c,0), (0,± c,± c), (± c,0,± c)                            i = 7,…,18;                         (± c,± c,± c)                            i = 19,…26 ].

where c=Δ x/Δ t is a shorthand and Δ t is the time step.
The fluid density ρ is calculated as the sum of the distribution function overall velocities

    ρ= ∑_i f_i
 where the shorthand

    f_i ( x⃗, t )
    ≡f ( x⃗, c⃗_i, t )

has been used.
Likewise, the fluid velocity u⃗ is calculated as average overall physical velocities c⃗_⃗i⃗ weighted by the amount of pseudo particles f_i with

    u⃗ = 1/ρ ∑_i c⃗_⃗i⃗ f_i
.
This way, the lattice Boltzmann equations can be written as


    
    f_i^* (x⃗,t)
            =
            f_i(x⃗,t)
            -1/τ[ f_i(x⃗,t) -f_i^eq(ρ(x⃗,t), u⃗(x⃗,t))]
            +Δ t F_i

    
    f_i(x⃗,t + Δ t )
            =
            f_i^*(x⃗- c⃗_⃗i⃗Δt,t)

    

which describe the collision step Eq. (<ref>) and the streaming step Eq. (<ref>) of the pseudoparticles, f_i.
The right-hand side of Eq. (<ref>) 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


    f_i^eq(ρ, u⃗)
        =
        ρ∑_i w_i
        [1-
            u⃗·u⃗/2 c_s^2 +
            c⃗_⃗i⃗·u⃗/c_s^2 +
            (c⃗_⃗i⃗·u⃗)^2/2 c_s^4]
,
which is a second-order expansion of the local Maxwell-Boltzmann distribution with the weights


    w_i
        =
        {[       8/27      i = 0;       2/27  i = 1,…,6;       1/54 i = 7,…,18;      1/216 i = 19,…26 ].
.
τ is the relaxation time and c_s the speed of sound in the lattice Boltzmann fluid with

    c_s =c/√(3)
.

We use the forcing term as published in <cit.> with


    F_i
        =
        (1-1/2τ) w_i [c⃗_i-v⃗/c_s^2 + c⃗_i ·v⃗/c_s^4c⃗_i]
        ·𝖿⃗
,
where 𝖿⃗ is the local body force density acting on the fluid and v⃗ is the fluid velocity which can be calculated with


    ρv⃗
        =
        ∑_i c⃗_i f_i + Δ t/2F⃗
.
We apply body force densities originating from the solid-fluid interaction (Sec. <ref>) and from non-ideal fluid contributions (Sec. <ref>) which cause a phase separation between the liquid and vapor phase.



  §.§.§ 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 <cit.> 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.<ref>).
Near a solid surface, the bounce-back adds to the right-hand side of Eq. (<ref>) and accounts for fluid which “bounces back” form the solid surface


    f_i (x⃗, t + Δ t )
        =
        {[ f^*_-i(x⃗, t) +
               g(x⃗, c⃗_⃗i⃗, u⃗_B,t )                           if x⃗-c⃗_⃗i⃗Δt is solid;                            f^*_i (x⃗-c⃗_⃗i⃗Δt, t)                                              else ]..

If no point x⃗-c⃗_⃗i⃗Δt in the vicinity of x⃗ is solid, the bounce-back is not calculated so that Eq. (<ref>) and Eq. (<ref>) deliver the same result.
However, if the point x⃗-c⃗_⃗i⃗Δ t lies within a solid body, there is no fluid population that can be streamed to point x⃗.
Hence, the pseudo particles coming from the point f_i(x⃗) with the velocity c⃗_⃗i⃗ are assumed to “bounce-back” from the solid surface in the middle between the two points.
Consequently, f^*_-i(x⃗, t) denotes the population of the pseudo particles with the velocity -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 x⃗-c⃗_⃗i⃗Δ 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 u⃗_B.
The movement of the solid body leads to a change of momentum of the adjacent fluid which is accounted for by


    g (x⃗, c⃗_⃗i⃗, u⃗_B(x⃗) )
        =
        2ρ(x⃗)/c_s^2
        w_i.c⃗_⃗i⃗·u⃗_B(x⃗-c⃗_iΔ t /2).
        + δρ(x⃗) δ_i0

in Eq. (<ref>).
Whereas δρ(x⃗) is an auxiliary term that ensures the local mass conservation with


    δρ(x⃗) =
        2ρ(x⃗)/c_s^2
        w_ic⃗_⃗i⃗·u⃗_B(x⃗-c⃗_iΔ t /2)
,
where the sum over i is executed for each point x⃗-c⃗_iΔ t that is solid so that the total change of mass caused by the first term in Eq. (<ref>) is compensated.
This mass correction δρ(x⃗) is then added to the resting fluid population with the term δ_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 Δm⃗_B of the rigid body at a solid point x⃗ is given by


    Δm⃗_B (x⃗)
        =
        -
        c⃗_⃗i⃗[
            2 f^*_i(x⃗+c⃗_⃗i⃗Δt) +
            g(x⃗+c⃗_⃗i⃗Δt, c⃗_⃗i⃗, u⃗_B)
        ]

where the sum is executed for each point x⃗+c⃗_iΔ t that is fluid.
By integration over the entire solid body B, corresponding total force F⃗_B and torque T⃗_B can then be obtained with



    F⃗_B = 1/Δt∫_BΔm⃗_B (x⃗)d V


    T⃗_B
        =
        1/Δt∫_B(x⃗-c⃗_⃗i⃗Δt/2-X⃗_B) ×Δm⃗_B (x⃗) d V
,
    

where X⃗_B denotes the center of mass of the solid body (see Sec. <ref>).

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 x⃗_n may become a solid node because a rigid body is passing by.
The remaining fluid distribution function f_i(x⃗_n) has to be distributed to the neighboring fluid points so that the fluid mass and momentum are conserved.
This conservation can be achieved with

    f^new_i(x⃗) = f_i(x⃗) + f^eq_i( ρ(x⃗_⃗m⃗), u⃗(x⃗_⃗m⃗) )/|N_F(x⃗_⃗m⃗)|
,
where N_F(x⃗_n) is the set of fluid neighbor nodes in a sphere of the radius R_A around x⃗_n with

    N_F(x⃗_⃗n⃗) = .{x⃗ | ‖x⃗_⃗n⃗ - x⃗ ‖≤R_A  and x⃗ ≠x⃗_⃗n⃗}.

In our simulations, the choice of R_A = Δ 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.



  §.§.§ 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 <cit.>.
The liquid and vapor phase separation is introduced by a force density 𝖿⃗ acting on the fluid (see Eq. (<ref>) and (<ref>)) with


    𝖿⃗(x⃗) = -
        𝒢ψ(x⃗)∑_i w_i ψ(x⃗ + c⃗_⃗i⃗Δt )c⃗_⃗i⃗
,
where 𝒢 is a coupling constant and ψ a so-called lattice version of the mean-field potential,


    ψ(x⃗) = 1 - e^-ρ(x⃗)/ρ_0
,
where ρ_0 is a reference density.
The force density Eq. (<ref>) can be interpreted as the sum of force densities which act between the pseudo particles at x⃗ and the surrounding pseudo particles at x⃗ + c⃗_⃗i⃗Δ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 ϱ_B.
By the integration over all force densities between a fluid point x⃗ and an adjacent point x⃗ + c⃗_⃗i⃗Δt of the rigid body B, the resulting force F⃗_B and torque T⃗_B acting on it can then be obtained



    F⃗_B
        = -
        𝒢∫_F dx⃗ ψ(x⃗)  
        
        w_i ψ(x⃗ + c⃗_⃗i⃗Δt )c⃗_⃗i⃗


    T⃗_B =
        -𝒢∫_F dx⃗ ψ(x⃗)  
        
        ( x⃗ + c⃗_⃗i⃗ - X⃗_B ) ×
        w_i ψ(ϱ_B )c⃗_⃗i⃗


where again the sum over i is executed for each point x⃗+c⃗_iΔ t that is solid.

With

    𝖿⃗ = - ∇·p

it is possible to identify the corresponding pressure tensor


    p_ij =
        .[
            c_s^2 ρ +
            1/2 c_s^2 𝒢ψ^2 +
            1/2 c_s^4 𝒢ψΔψ +
            𝒢 c_s^4/4|∇ψ|^2
        ] δ_ij.
        - 1/2 c_s^4 𝒢∂_i ψ∂_j ψ
.
For the bulk pressure or the equation of states, one thus has


    p = c_s^2ρ + 1/2 c_s^2 𝒢ψ^2
.
The pressure tensor can be used to calculate the surface tension.
According to <cit.>, the interface tension σ is given by the difference of the normal component of the pressure tensor p_ and its tangential component p_∥ across the interface with


    σ= ∫_0^ηdn
        (p_- p_∥)
,
where n is the direction normal to the interface.
If one considers only a planar interface, so that ψ changes only along the normal direction of the interface, the interface tension can be calculated with


    σ=-
        1/2 c_x^4 𝒢∫|∂ψ/∂n|^2  d n
,
which is obtained by inserting Eq. (<ref>) into Eq. (<ref>).

Considering Eq. (<ref>) and Eq. (<ref>), 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 𝒢.
The density profile across the solid-vapor and the solid-liquid can be adjusted tuning the above introduced “solid density” ϱ_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 ϱ_B as the wetting parameter.
For a more detailed explanation, we refer the reader to <cit.>.



 §.§ Simulation Procedure


As mention in Sec. <ref>, we use a high-resolution scheme for the advection, Eq. (<ref>), of the phase-fields with directional splitting and a monotonized central flux limiter <cit.>.
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 <cit.>, the kinetics of ϕ read,



    ϕ̇_̇α̇
        =
        ∑_β≠α^N
        M_αβ/N(
        [
            σ_αβ( I_α - I_β) +
            ∑_γ≠α,β(σ_βγ - σ_αγ) I_γ] +
        π^2/4ηΔ g_αβ)


    I_α = ∇^2ϕ_α + π^2/η^2ϕ_α,


where M_αβ is the mobility of the interface between phases α and β, σ_αβ is the interface energy between the phase-fields and η is the interface width (see <cit.> for more details).

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


    Δϕ_α = ∫[ϕ_α(t_ref) - ϕ_α(t)]  dV

according to the local curvature.
Because of this shrinkage, the last term Δ g_αβ is necessary which is dynamically determined to conserve the volume of a phase-field ϕ_α with


    Δ g_αβ = - Δ g_βα = 4η/π^2Δϕ_α/|δϕ_α|
.
Similar as in Sec. <ref>, |δϕ_α| denotes the number of nodes in the interface region of ϕ_α with δϕ_α = .{x⃗| 0.0 < ϕ(x⃗) < 1.0}.
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:




§ RESULTS


We first investigate the accuracy of the surface tension calculation in Sec. <ref>.
This investigation is the basis for the testing of the present model.
In Sec. <ref>, 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 <ref>.
To demonstrate the capability to model the process of liquid-phase-sintering,
we also investigate a more complicated situation in Sec. <ref>, 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.



 §.§ Surface Tension


In the following, we test the consistency of the surface tension calculations Eq. (<ref>) with the Laplace pressure Eq. (<ref>).
By evaluating Eq. (<ref>) for droplets of different radii, we obtained the value 3.297 · 10^-2 in lattice units for the surface tension.
The value is used to predict the pressure gradient, by inserting it into Eq. (<ref>).
We show this prediction in Fig. <ref>.
Also, the differences for different radii, obtained with Eq. (<ref>), are depicted.
One can see that the obtained pressure gradients are in agreement with the ones predicted by the Laplace pressure, Eq. (<ref>).
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. (<ref>) and use it together with the radius as input for Eq. (<ref>) to calculate the resulting surface tension.
This way, we yield 3.583 · 10^-2 for the value of surface tension which differs from the value obtained with Eq. (<ref>) 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, σ = 3.583 · 10^-2 is used in the following.



 §.§ 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. <ref>.
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. <ref>, we introduced the approximative solution Eq. (<ref>) for the motion of the plates connected by a cylindrical capillary bridge and the analytic solution Eq. (<ref>) 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 ρ_S as compared to the liquid mass density ρ_L with ρ_S/ρ_L = 1.72 · 10^4.
The initial simulation setup is schematically depicted for the planar setup in Fig. <ref> and for the cylindrical setup in Fig. <ref>.
Figure <ref> 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= √(m/σ) where m is the mass of the solid plates.
By using R and t, one can rewrite Eqs. (<ref>) and (<ref>)



    h_p/R = h_0/R - 2(t/t_c)^2


    h_c/R≈h_0/R - π(t/t_c)^2
.
    

To obtain the rescaled equations of motion Eq. (<ref>), we have assumed for the length of the planar capillary bridge that L=R.
The simulation results are depicted in Fig. <ref> together with the rescaled equations of motion Eq. (<ref>).
We also show as insets in Fig. <ref> 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. (<ref>).
Consequently, we show only the very beginning of the motion in Fig. <ref> 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. <ref> and the cylindrical one Fig. <ref>, the results agree with the predictions of Eqs. (<ref>) and (<ref>) 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.



 §.§ Dynamics of Multiple Particles





  §.§.§ 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 ρ = 4, and liquid fraction of c=0.5 (see Fig. <ref>).
The closest distance between the particles h and the velocity ḣ are shown in Fig. <ref>.
Because we want to consider solids with a very high wettability, a wetting-parameter of ϱ_B = 2.0 has been chosen.
Above in section <ref>, we have introduced the phase-field dynamics Eq. (<ref>) 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. <ref>, the closest distance between the surfaces of the particles h and the corresponding velocity ḣ are shown.
The dashed lines in Fig <ref> 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. <ref>).
After two-thousand time-steps, the two particles are continuously decelerated until the particles come into contact.
By looking at Fig. <ref>, 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. <ref> 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. (<ref>) can only cause this deceleration.


When the phase-fields are not only advected (solid lines in Fig <ref>), 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 0 to 1 within a distance η = 5.
This defines the interface thickness.
Hence, the phase-fields begin to touch each other about a distance of 2 η = 10.
At this distance, the phase-fields start to merge and form a common interface between the two solid particles.

When η 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 η is larger than a physical grain boundary width, the dynamics of the phase-field and the resulting variation of the distance below 2η shall be viewed as a numerical feature to model solid-phase sintering of larger scales.

At a distance below 2η, 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. <ref>.

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 <cit.>.



  §.§.§ 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. <ref>.


We start with two spherical particles with a radius of R = 30 and a mass density of ρ = 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 <ref>, 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. <ref>.
For a liquid fraction of c=0.5, we show again with a dashed line for a simulation where the phase-field dynamics, Eq. (<ref>), 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. <ref>).
In Fig. <ref>, 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. <ref>, 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. <ref> shows the decreasing average distance h̃ = 1/3 ∑_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 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. <ref> shows the decreasing average distance h̃ = 1/4 ∑_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.



§ 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.






§ 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

    d U =  σ d A

with A=2π R one obtains

    d U =  2πσ 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 = π R_0^2 h_0 = π 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√(h_0/h_c).
Hence, the force is obtained as

    F= (d U/d h_c)_V_L = πσ R_0√(h_0/h_c)
.




§ LIQUID BRIDGE BETWEEN TWO SOLID SPHERES





 §.§ Spherical Bridge

Equation (<ref>) 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. <ref>.



  §.§.§ Small Spherical Capillary Bridge





The contact angle Θ_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


    Θ_S = θ_2 - θ_1
,
where θ_1 and θ_2 are the angles between the surface normals and the depicted dashed line (see. Fig. <ref>).
These angles can be calculated by


    
    θ_1 = arccos(R_S - h_S/R_S)

    
    θ_2 = π - arccos(h_S/χ_S R_S)
,
    

where χ_S∈[0,2] 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 = χ_S R_S.
Again, one can obtain a system of equations for the side lengths of the depicted triangles by using the Pythagorean theorem,


    
    χ_S R_S = h_L + h_S

    
    R_S^2 = a^2 + (R_S - h_S)^2

    
    χ_S^2 R_S^2 = a^2 + h_S^2
,
    

with the solution


    
    h_S = χ_S^2/2 R_S

    
    h_L = 2χ_S-χ_S^2/2 R_S

    
    a = 1/2√(4χ_S^2-χ_S^4) R_S
.
    

In Eq. (<ref>), we expressed the geometric parameters depicted in Fig. <ref> as function of the radius R_S and the scaling parameter χ_S.
As to be expected, the contact angle can be expressed as a function of χ_S with


    Θ_S
        = arccos(χ_S^2/2 - 1) - arccos(χ_S/2)
.

Furthermore, we replace the scaling factor χ_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. <ref>) form the volume of a liquid sphere V_Ls with


    V_L = V_Ls  - 2 V_scA - 2 V_scB

and


    
    V_Ls = 4π/3χ_S^3 R^3

    
    V_scA = π/3 h_S^2 (3R_S - h_S)

    
    V_scB = π/3 h_L^2 (3χ_S R_S-h_L)


and Eq. (<ref>).
The volumes of the spherical caps, V_scA and V_scB, can be parametrized with R_S and χ_S by inserting Eq. (<ref>) into Eq. (<ref>) and Eq. (<ref>) into Eq. (<ref>) so that


    
    V_scA = π R_S^3 (-1/24χ_S^6 + 1/4χ_S^4 )

    
    V_scB
            =
            π R_S^3 (1/24χ_S^6 - 1/2χ_S^4 + 2/3χ_S^3)
.

Inserting Eqs. (<ref>), (<ref>), and (<ref>) into Eqs. (<ref>) delivers a simple parametrization for the liquid volume


    V_L = 1/2πχ_S^4 R_S^3
.
The same way as for the cylindrical bridge, we can calculate the liquid volume fraction with


    
    c = V_L/V_L + V_S

    
    = 3χ_S^4/3χ_S^4 + 16
.
    

By rearranging the terms of Eq. (<ref>), we obtain a fourth order polynomial equation for χ_S


    χ_S^4 = 16c/3(1-c)
,
which has a real and positive solution


    χ_S
        = 2/√(3)√((1-c)√(c(1-c)))/1-c
.
Inserting Eq. (<ref>) into Eq. (<ref>) delivers the solution shown in Fig. <ref>.



  §.§.§ Large Spherical Capillary Bridge

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

The angles θ_1 and θ_2 are obtained in the same manner as above with


    
    Θ_S = θ_2 - θ_1

    
    θ_1 = π - arccos(h_S - R_S/χ_S R_S)

    
    θ_2 = π - arccos(h_S/R_S)
,
    

where θ_1 can be simplified to


    θ_1 = arccos(R_S - h_S/χ_S R_S)
.
The height h_S can be calculated in the same manner as in <ref>, and Eq. (<ref>) can be obtained again.



 §.§ Cylindrical bridge


Equation (<ref>) 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. <ref>.


To calculate the contact angle Θ_c, we consider the angle θ 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


    θ= arccos(R-h/R)
,
where h is half of the height of the cylindrical capillary bridge.
By considering the intercept theorem, we obtain the contact angle


    Θ_c = π/2 - arccos(1-χ_c)
,
where we have introduced the dimensionless scaling parameter χ_c∈[0,1] with h=χ_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,



    V_L
        = V_cylinder - 2 V_sphere-cap


    = 2 π R_S^2 h -
          2 π/6 h [ 3R_S^2 + h^2 ]


    = π R_S^2 h - π/3 h^3
.

By using the Pythagorean theorem

    R^2 = R_S^2 + (R-h)^2
,
the radius of the capillary bridge R_S in Eq. (<ref>) can be replaced by the radius of the solid sphere R_S.
With this replacement, we obtain for the fluid volume


    V_L
        = 2 π R h^2 - 4/3π h^3
.
Furthermore, we replace height h again with h= χ_cR so that the expression for the fluid volume


    V_L
        = (- 4/3χ_c^3 + 2 χ_c^2) π R^3

can be further simplified.
With Eq. (<ref>) and the volume of the solid,

    V_S = 8/3 πR^3
,
the liquid volume fraction can be calculated



    c
        = V_L/V_L + V_S


    = 2χ_c^3 - 3 χ_c^2/2χ_c^3 - 3 χ_c^2 - 4
.

Equation (<ref>) can be transformed in a third order polynomial


    χ_c^3 - 3/2χ_c^2 + 2c/1-c= 0
.
The polynomial equation (<ref>) has three real roots for c∈(0,0.2).
By substituting χ_c = t + 1/2, we can write Eq. (<ref>) in the form,


    t^3 + p t + q =0

with

    p = -1/4

and

    q = 2c/1-c -1/4
.
With the use of the cosine power rule


    cos^3(η)
        - 3/4cos(η)
        - 1/4cos(3η)
        = 0
,
we find three real solutions of Eq. (<ref>)


    t_n
        =
        √(-4p/3)cos(
            1/3arccos[
                -4q(-4p/3)^3/2]
            + 2nπ/3)

with n ∈{0,1,2}.
By back substitution, we also obtain three real solutions for the scaling factor


    χ_cn
        =
        -cos(
            1/3arccos(
                1-9c/1-c)
            + 2nπ/3)
        + 1/2
.
The physically meaningful solution of χ_c has to be limited to the interval [0,1], which is only fulfilled by χ_c2.
Inserting Eq. (<ref>) into Eq. (<ref>) delivers the result shown in Fig. <ref>.




