Lattice-Boltzmann simulations of the dynamics of liquid barrels

We study the relaxation towards equilibrium of a liquid barrel—a partially wetting droplet in a wedge geometry—using a diffuse-interface approach. We formulate a hydrodynamic model of the motion of the barrel in the framework of the Navier–Stokes and Cahn–Hilliard equations of motion. We present a lattice-Boltzmann method to integrate the diffuse-interface equations, where we introduce an algorithm to model the dynamic wetting of the liquid on smooth solid boundaries. We present simulation results of the over-damped dynamics of the liquid barrel. We find that the relaxation of the droplets is driven by capillary forces and damped by friction forces. We show that the friction is determined by the contribution of the bulk flow, the corner flow near the contact lines and the motion of the contact lines by comparing simulation results for the relaxation time of the barrel. Our results are in broad agreement with previous analytical predictions based on a sharp interface model.


Introduction
Droplets in wedge geometries appear in many natural environments. For example, many shorebird species have wedge-shaped beaks that allow them to feed on water-bound organisms [1]; water striders have arrays of tapered bristles that help them brush off droplets from their legs [2]; and the material properties of wet granular media depend on the adhesion and lubrication provided by capillary bridges wedged between solid grains [3][4][5].
Understanding the motion of droplets in wedges is important to improve technologies that use the geometry of the confinement for purposes of transport, positioning or actuation of small volumes of liquid. Wetting droplets inside tapered capillary tubes [6] or wedge-shaped channels [7] self-propel towards regions of stronger confinement, while non-wetting droplets trapped in non-parallel channels migrate to regions of weaker confinement [7,8]. Such principles have been used to transport capillary bridges spontaneously [9], using mechanical [10] or photo-induced [11] actuation, and even to separate droplets formed by two immiscible liquids [12].
For the ideal situation of a droplet trapped between perfectly flat and smooth walls that form a wedge, the shape of the droplet depends on two parameters: the contact angle that the liquid makes with the solid, θ e , and the opening angle of the wedge, β (see figure 1) [13]. In this paper, we are interested in a specific interface configuration which we call a liquid barrel [14]. Liquid barrels form when the wetting angle of the droplet satisfies the relation θ e > 90 • + β [15]. They differ from edge spreads (θ e 90 • − β), edge blobs (90 • − β < θ e 90 • + β), and free drops (θ e = 180 • ) in that they equilibrate into a truncatedsphere shape away from apex of the wedge [15][16][17]. When displaced from their equilibrium position, they relax back driven by capillary forces (mediated by the geometry of the wedge) and damped by frictional forces [18].
Despite their ubiquity, the dynamics of liquid barrels has not been explored in detail. Experimentally, Ruiz-Gutiérrez et al studied the relaxation of water droplets in smooth wedges formed by slippery liquid infused surfaces, where a contact line is absent [18]. The authors found an exponential relaxation towards equilibrium. By comparing their measurements of the relaxation time to a sharp-interface model, they showed that the dominant source of dissipation is the flow within the bulk of the barrel. Then, assuming superposition of sources of dissipation, they expanded the sharp-interface model to include the effect of the flow near the contact line and of the motion of the contact line itself [14]. These assumptions, however, have not been verified. Therefore, a more detailed study of the motion of the liquid barrels, which includes the evo lution of the hydrodynamic fields and the contact line, is needed.
Diffuse-interface numerical simulations are good candidates to study the dynamics of liquid barrels, as they have the capacity of modelling capillary phenomena, including the dynamic wetting of smooth solid surfaces [19][20][21][22][23][24][25]. They do this by coupling the Navier-Stokes equations [26] with a mesoscopic equation of motion, such as the Cahn-Hilliard equation [27,28]. The main advantage of this approach is that the interface dynamics occur naturally through convection and diffusion-the latter driven by chemical potential gradients [19,29,30]. This contrasts with sharp-interface models, where one needs to track the evolution of the interface [31] and to specify a boundary condition for the contact line in an ad hoc manner [32].
In this paper we study the relaxation of a liquid barrel using numerical simulations of a diffuse-interface model; we use a lattice-Boltzmann algorithm, which is a suitable method of integrating the diffuse-interface equations [33][34][35][36]. We begin, in section 2, by presenting the diffuse-interface model and formulate the relaxation dynamics of the liquid barrel in this framework. Then, in section 3, we specify the lattice-Boltzmann method. We introduce a new algorithm to model a smooth wedge geometry, but which can be applied to implement boundaries of arbitrary shape. In section 4 we present our simulation results. We validate our simulation method by comparing the equilibrium configuration of the barrels to previous analytical results. We study the relaxation of the barrels to their equilibrium position. We investigate the role that hydrodynamic fields play during the translational motion of the barrel. We analyse the dissipation of energy and the relaxation time, comparing our simulation results with those of the sharp-interface model proposed by [14]. Finally, in section 5 we present the conclusions of this work.

Diffuse-interface model
We describe the system using a diffuse-interface model where the liquid barrel, referred to as the inner phase, and the surrounding fluid, or outer phase, are identified using a phase field φ(x). We define the phase field such that the inner phase corresponds to the φ > 0 phase, while the outer phase corresponds to the φ < 0. The fluids are contained in a 3D domain, Ω ⊆ R 3 , and are enclosed by a 2D surface boundary, ∂Ω, which represents the solid wedge geometry.
We introduce the Helmholtz free energy, as the relevant thermodynamic potential for constant volume and temperature situations. The first term is a volume contribution, where the free-energy density, ψ, is expressed as [37,38] The first two terms in equation (2) allow the coexistence of the inner and outer phases, while the last term gives rise to a diffuse interfacial region of typical thickness ξ and surface tension γ. The second term in equation (1) is the surface contrib ution to the free energy. The constant h controls the energy cost incurred when the fluids come into contact with a solid boundary, and is related to the contact angle of the interface with the solid, θ e , by where α = arccos sin 2 θ e [20,27]. From equation (1), it is possible to calculate the chemical potential field, and the pressure tensor field [38], where I is the identity matrix. From equation (5), it is possible to show that the jump of the normal stress across a gently curved interface obeys the Young-Laplace condition (see, e.g. [28]), where κ = −∇ ·n i /2 is the mean curvature of the diffuse interface, and n i is the local unit normal vector to it. Out of equilibrium, a gradient in the chemical potential will cause a diffusive current −M∇µ, where the constant M is called the mobility. In addition, the phase field will be advected by the velocity field, u. Therefore, the local conservation of φ is given by the Cahn-Hilliard equation [37]: The local conservation of momentum is governed by the incompressible Navier-Stokes equation [26], i.e.
where ρ and η are the local density and dynamic viscosity of the fluid [28]. Following [39], we define the local fluid density and kinematic viscosity in terms of the phase field via the profiles and where ρ in , ρ out , η in and η out are the saturation mass densities and dynamic viscosities of the inner and outer phases. Note that ρ and η are required to be positive definite quantities regardless of the sign in φ.
The system of partial differential equations (7) and (8) is subject to the boundary conditions at the solid walls: where n is the unit normal vector to the solid surface. Equation (11) determines the normal component of ∇φ at the boundary according to the surface energy parameter, h, and gives rise to the wetting behaviour of the fluid. Equation (12) imposes a vanishing flux of φ in the normal direction to the solid boundary. Equation (13) combines the impenetrability and no-slip conditions. Equation (14) imposes the normalstress balance at the solid boundary.
Within the diffuse-interface model, the motion of the contact line occurs by virtue of diffusive currents caused by a local imbalance in the chemical potential field [21,24]. This is because, while the velocity field vanishes at the solid-fluid interface by virtue of equation (13), the diffusive term in equation (7) does not. This regularises the singularity that stems from the no-slip boundary condition [23]. As shown in [20] by Briant et al, the combination of both features leads to a region of characteristic size m where the contact line 'slips' past the solid surface. The connection between the diffuse-interface contact-line dynamics and the corresponding sharp-interface description can be established by identifying m with the cutoff length-scale below which the hydrodynamic description breaks down (see, e.g. [24,40]); this is the contact-line slip length within the Cox-Voinov relation [41] θ where and v cl is the projection of the velocity of the interface in the direction perpendicular to the contact line and λ ≡ η out /η in is the viscosity ratio. The Cox-Voinov relation links the hydrodynamic distortion of the fluid-fluid interface at a typical macroscopic length scale M , characterised by the apparent contact angle θ, to its velocity, v cl . Apart from the hydrodynamic effect, the contact line will give rise to a local friction force, −ζ 0 v cl , where ζ 0 is a contact-line friction coefficient. Several authors have studied the dependence of m and ζ 0 on the diffuse-interface model parameters [19,22,23,29]. Recently, Kusumaatmaja et al, [24] showed the existence of two scaling regimes of m with the fluid mobility and viscosity: a diffuse-interface regime occurs if ξ (Mη) 1/2 , and leads to m ∼ ξ 1/2 (Mη) 1/4 ; a sharp-interface regime occurs if ξ (Mη) 1/2 , and leads to m ∼ (Mη) 1/2 . In both regimes, however, the slip length obeys where we define the dimensionless mobility as M * ≡ Mη in /ξ 2 . Jacqmin studied the motion of the contact line using a diffuseinterface model, and found that the speed of the contact line increases with the strength of diffusion [19]. Therefore, we expect that the friction coefficient obeys

A simple model of the dynamics of a liquid barrel close to equilibrium
The equilibrium configuration of a liquid barrel corresponds to a spherical interfacial shape truncated by the walls of the wedge [15]. For a given droplet volume, equilibrium contact angle, and wedge angle, the radius and centre of the truncated sphere (relative to the apex of the wedge) read and The equilibrium surface energy of the liquid barrel reads [14] E 0 = 3γV 2/3 π 6 (cos 3θ e − 9 cos θ e ) Therefore, equation (21) can be used as the zero-point energy of the system. The truncated-sphere equilibrium states are stable [14,15,17], i.e. E 0 F for any deformation of the interface that keeps the volume constant. Consequently, a deformed droplet will relax back to the equilibrium barrel configuration. In general, this relaxation process will include the reshaping of the interface to a truncated sphere and the translation of the droplet along the wedge. We are interested in the long-time translational motion of the barrel where short-wavelength perturbations have already been damped [26,43]. In such a situation, the free energy and the energy dissipation can be expressed as F = F(X) and Ė =Ė(X,Ẋ), where X and Ẋ are the position and translational velocity of the barrel. As the droplet relaxes, the free energy is consumed by dissipation [42], thus, during its motion we obtain, which leads to the equation of motion. The free energy can be expanded in a Taylor series around equilibrium where k ≡ d 2 F/dX 2 | X=Xe > 0 is the restitution constant. In equilibrium the linear term vanishes, and thus, does not appear in equation (23). The dissipation function can be expressed aṡ where ν ≡ − 1 2 ∂ 2Ė /∂Ẋ 2 |Ẋ =0 > 0 is the translational friction coefficient. In the expansion, we have dropped the constant and linear terms because dissipation of energy cannot take place in equilibrium in a closed system (constant term), and because spontaneous creation of energy cannot occur (linear term). In the overdamped regime, the kinetic energy is assumed negligible, implying that the interfacial energy is dissipated by friction. Therefore, we can construct the equation of motion of the barrel by equating (23) with (24), and using the chain rule, dF/dt =Ẋ dF/dX. This yields which has the solution is the relaxation time.
Ruiz-Gutiérrez et al [14] formulated a sharp-interface analytical model of the exponential relaxation of a liquid barrel in a wedge geometry where η in η out , i.e., in the limit λ → 0. They assumed a quasi-spherical barrel shape out of equilibrium, obtaining that the restitution constant can be expressed as with the constant of proportionality being a function of the contact angle.
To compute the translational friction coefficient ν they considered three main sources of dissipation: the bulk flow within the barrel (ν bulk ), the corner flow near the contact line (ν corner ), and the motion of the contact line itself (ν cl ), We shall follow a similar approach and include some generalisations that consider the effect of a viscous outer phase. The bulk coefficient is composed of the dissipation from the inner and outer phases ν bulk = ν in + ν out where we take, from [14]. Analogously, we expect that the friction coefficient of bulk from the outer phase is The corner flow dissipation coefficient is described by equation (15), for small deviations of the contact angle, where we have used equation (17) to estimate m and M ∝ V 1/3 . The dissipation coefficient related to the contact line is given by where we have used the scaling for ζ 0 as in equation (18). We thus expect, by combining equation (27) through equation (33), that the relaxation time of the liquid barrel immersed in a viscous flow is generalised to the scaling form: .., 4, are functions of the contact angle yet to be determined.

Governing equations
The lattice-Boltzmann method numerically integrates the discretised Boltzmann equation [44] where f q (x, t) is a particle distribution function that represents the average number of fluid particles with position x and velocity c q at time t. Space and time are discretised, and the velocity space is restricted to a set {c q } Q−1 q=0 where Q is the number of directions in which particles can move. We use the D3Q15 model [34], which consists of a 3D lattice with 15 velocity vectors.
The time evolution of the distribution function in equation (35) consists of a collision step where f q relaxes to an equilibrium value f e q over a timescale determined by the collision parameter τ f (right-hand side of the equation), followed by a streaming step where f q is propagated along the direction of c q over a unitary time increment (left-hand side).
The Navier-Stokes equation, equation (8), is recovered by means of a Chapman-Enskog expansion of equation (35) [33]. The local momentum density is related to the first moments of the distribution function, i.e.
The equilibrium distribution function, f e q , is constructed to convey the thermodynamic behaviour of the fluid, and to ensure local mass and momentum conservation. This is done by requiring that the moments of f e q obey q f e q = ρ, q f e q c q = ρu, and q f e q c q c q = P + ρuu. A suitable choice for the equilibrium distribution is, which, in combination with equation (10), specifies τ f . To integrate the Cahn-Hilliard equation, equation (7), we introduce a second discretised Boltzmann equation, where g q is a distribution function with a collision parameter τ g = 1, whose zeroth moment is related to the phase field, The corresponding equilibrium distribution function, g e q , is defined as thus satisfying q g e q = φ, q g e q c q = φu and q g e q c q c q = Mµ(τ g − 1/2) −1 I + φuu. For a set of values of the distribution functions f q and g q , equations (36) and (40) determine the fluid velocity, and phase field. To model fluids of different densities, we calculate the local mass density using equation (9). We then evaluate f e q and g e q using equations (37) and (41), and the collision terms in equations (35) and (39). To compute the pressure tensor and the chemical potential in equations (37) and (41) we use equations (5) and (4). These require computing the gradient and Laplacian of the phase field, which we approximate using the finite-differences stencils where the w q are used as weighting factors to optimise the accuracy of the approximation [46]. The collision parameter τ f is determined by combining equations (10) and (38). Once the collision terms are determined, we iterate equations (35) and (39). The lattice-Boltzmann method is known to break Galilean invariance in situations where the fluid has density inhomogeneities [33,47]. Following [48], we added a correction term in the pressure tensor when calculating the equilibrium distribution function f e q to reduce this effect.

Boundary conditions
Boundary conditions in the lattice-Boltzmann method arise in the streaming steps of equations (35) and (39), and in the spatial derivatives of the phase field needed to compute the equilibrium distribution functions in equations (37) and (41).
In the following, we will refer to these as kinetic boundary conditions and finite-differences boundary conditions. We define a near-boundary node, of position vector x b ∈ Ω, as a node that has at least one lattice vector that crosses the boundary of the simulation domain (see figure 2). These lattice vectors define a set of cut links, Γ c [49]. Each cut link is characterised by its length, |δ q c q |, where 0 < δ q 1, and by a local normal vector to the boundary, n q =n(x b + δ q c q ).
The kinetic boundary conditions consist of specifying the particle population fq streaming into the simulation domain opposite to the cut link, where q ∈ {q | c q = −c q ; q ∈ Γ c }. For no-slip walls at rest, fq is given by a bounce-back algorithm [49][50][51] with linear interpolation, is called the post-collision distribution function. The same algorithm is used to determine the kinetic boundary condition for gq. The bounceback algorithm satisfies the no-flux boundary conditions, equations (12)- (14).
The finite-differences boundary conditions consist of calculating the derivatives of the phase field at a boundary node in such a way that equation (11) is satisfied. Specifically, we write the Taylor expansion [52] if q ∈ Γ c , and otherwise. In equations (45) and (46) the gradient vector, ∇φ(x b ), and the Hessian matrix, ∇∇φ(x b ), are unknown. In 3D, the gradient vector and the Hessian matrix comprise 3 + 6 independent components, forming a set of 9 unknowns. Equations (45) and (46), however, give Q − 1 = 14 equations. Therefore, the system is over-specified.
To solve the system of equations we introduce a pseudoinverse algorithm. First, we express equations (45) and (46) in the same units by multiplying every instance of equation (45) by δ q |c q | 2 /n · c q . Then, we express the system of equations in the matrix form: where is a 9 × 1 vector containing the unknown entries of the gradient vector and the Hessian matrix, Φ is a 14 × 1 vector of known field values and boundary conditions whose entries read and G is a 14 × 9 matrix of coefficients that reflects the local structure of φ, including the boundaries. The pseudo-inverse algorithm consists of estimating the solution, Λ = G −1 Φ, computing G as In equation (50), E is a 9 × 14 matrix which projects G into a 9 × 9 matrix. This can be thought of as a weighting of the entries of G while preserving linear independence; in the spirit of equations (42) and (43), we define columns of the projection matrix as, where we have generalised δ q = 1 for q ∈ Γ c . Although the expression of E is not unique, we found that equation (51) produces the expected interface profile near the solid boundaries, which we quantified using the equilibrium contact angle.
Because the matrix G stores the structure of the lattice and of the solid boundary (which do not change over time), the pseudo-inverse algorithm, equation (50), is applied numerically at the initialisation of a simulation, and is therefore not more expensive than the usual application of a finite-differences stencil.

Simulation setup
We set up a simulation domain contained in a box of dimensions N x × N y × N z , as shown in figure 3. We fix the walls of the A near-boundary lattice node (large square) has at least one lattice link that intersects the solid boundary. Each of such cut links is defined by its direction, and by the fractional distance to the wall, δ q . At the intersection, the boundary is defined by a local normal vector n q . wedge at (x − b) ·n = 0, where n(β) = (− sin β, 0, cos β) is the normal vector of the upper (+β) and lower (−β) wall. We introduce the offset b to avoid that the walls intersect within the simulation domain, which we found gave rise to numerical errors. We use b ·n = 1.72 which allows at least two nodes of separation between the boundaries. We use periodic boundary conditions along the y direction, and impose two solid planes at x = 1/2 and x = N x − 3/2.
As an initial condition the fluid is at rest. The interface has a spherical configuration with radius R 0 and centre X 0 as shown in figure 3. This corresponds to setting the following initial velocity and phase field profiles We have used two initial values of the centre of sphere, X 0 = (0.75N x , 0, 0) and X 0 = (0.24N x , 0, 0), which, for the range of tapering and contact angles considered, correspond to the droplet initial positions, X 0 > X e and X 0 < X e . These initial conditions allow us to study the inwards or outwards motion of the barrel. The instantaneous volume, V(t), and position of the liquid barrel, X(t), are computed as follows. To calculate the volume of the droplet we use To calculate X(t), we take a slice of the phase field at the bisector plane, x = (x, y, 0). We then use a linear interpolation scheme to find the interface, defined as the level curve φ = 0.
We take X(t) as the x coordinate of the centre of the circle that best fits the interface profile using a least mean squares algorithm. The values of the simulation parameters are summarised in table 1. For these parameter values, a simulation time of ∼ 5 × 10 6 iterations was found to be a sufficient timescale for the system to relax to equilibrium.

Equilibrium
We first focus on the equilibrium configuration of the droplet. As expected from the analytical model, the interface adopts the shape of a sphere truncated by the solid planes of radius and position given by equations (19) and (20) (figure 4(a)). The simulations agree with the analytical prediction over the whole range of equilibrium and wedge angles considered (see figures 4(b) and (c)); this agreement holds regardless of the initial conditions (compare symbols to solid lines in the figure).
Beyond the configuration of the interface, the simulations provide details of the equilibrium hydrodynamic fields. The velocity field vanishes everywhere, except for small currents (∼ 10 −5 in lattice-Boltzmann units; see figure 5(a)) which arise due to spurious effects in the lattice-Boltzmann method a Unless otherwise specified, values are reported in lattice-Boltzmann units.  [36]. Despite being augmented by the lower density and viscosity of the outer phase, however, they have a negligible effect on the liquid barrel. The scalar pressure profile, p ≡ trP/3, along the bisector line, dips at two points (see figure 5(b)). These correspond to the interfacial regions (compare to the phase-field profile shown as a dashed line). As pointed out by [39], this is due to the free-energy density contribution to the pressure and gives rise to the surface tension effect. The pressure is higher in the bulk of inner phase than in the outer phase. This is the combined effect of surface tension and the curvature of the interface, as expected from the Young-Laplace relation, equation (6).

Relaxation towards equilibrium
We now focus on the translational motion of the droplet during equilibration. Figure 6(a) shows a sequence of simulation snapshots at intervals of 8 × 10 5 time steps for droplets moving inwards and outwards to the same equilibrium position. In both cases, the position of the droplet as a function of time follows an exponential relaxation (figure 6(b) and inset).
To better understand the mechanism driving the relaxation process, we consider the instantaneous pressure and chemical potential distributions within the droplet. The slope of the pressure profile, ∂p/∂x, is negative for droplets that move outwards and positive for droplets that move inwards (see figure 7). This indicates the action of a net capillary driving force resulting from Laplace pressure differences.
We now turn our attention to the flow pattern during the equilibration process. We observe two recirculating vortices at the sides of the droplet (see streamlines in figure 8(a)). We found that these structures persist over the whole set of simulations, regardless of the direction of motion or the velocity of translation. This is reasonable, as the dynamics in the simulations always fall in the exponential regime. Therefore, we expect that while the magnitude of these features decreases as the system approaches equilibrium, the overall structure of the flow remains the same [26].
The outer flow contrasts with the flow in the bulk of the barrel, which shows a remarkably laminar structure (see figure 8(a)). From the simulation results, we can identify two regions where the inner flow pattern is qualitatively distinct, corresponding to the bulk flow and the flow near the contact line. The laminar flow occurs in the bulk of the barrel, as shown in the cross-section shown in figure 8(b). The flow points in the direction of the apex, growing in magnitude    from the walls to the bisector plane. This is reminiscent of the pressure-driven flow within a wedge, also known as a Jeffery-Hamel flow [26]. The structure of the flow near the contact line, shown in the close-up of figure 8(c), is consistent with the generic corner flow of wetting dynamics predicted by Cox and Voinov [40,53], which results in a tread-milling motion of the interface as documented in experiments by [54]. Thus, the structure of the bulk and corner flows in the simulations is in good qualitative agreement with the theoretical model of the barrel dynamics introduced by [14].

Relaxation time
The exponential relaxation of the barrels implies that Ẋ = −(X − X e )/τ, where τ is the relaxation time. We found that this linear relation is satisfied for barrels over the whole range of contact angles considered, and that, for constant values of all other parameters (β, M and η in ), the data collapse onto a single line (see figure 9(a)). The data collapse implies that the relaxation time does not change significantly with the contact angle. While at lower contact angles we expect a stronger restitutive force from equation (28), we also expect a higher hydrodynamic dissipation near the contact lines due to the development of vortices in the inner (more viscous) phase. We can infer that the variation of the dissipative and conservative forces with the contact angle is therefore approximately equal, thus cancelling out a net dependence of τ on θ e . From equation (34), we expect τ ∝ 1/ sin 2 β , and thus − sin 2 βẊ ∝ (X − X e ), which is confirmed by the simulation results (see figure 9(b)). Physically, this functional dependence arises from the projection of surface tension and pressure forces along coordinate of translational motion, which, as first noted by [17], balance out in equilibrium. Because the restitution coefficient scales as 1/ sin 2 β, we can conclude that the friction coefficient-and therefore the dissipation itselfis independent of the wedge geometry (at least for the small wedge angles considered here).
We found that τ increases linearly with the viscosity of the inner phase, η in (see figure 10(a)). This is in agreement with equation (34), implying that all contributions to the drag coefficient arising from the inner phase scale with η in . The extrapolation of the relaxation time to a finite value as η in → 0 is due to the non-zero viscosity of the outer phase.
With the intention of giving generality to our results, we now turn to the dimensionless representation of the relaxation time, τ * . The dependence of τ * on the dimensionless mobility, M * , is presented in figure 10(b). We found a monotonic decrease, which can be reasoned in terms of both a larger microscopic length-scale m , and a smaller contactline friction coefficient, ζ 0 , at higher M. More quantitatively, we compared the simulation results to the theory by varying simultaneously M * and λ to fit the parameters b i , i = 1, ..., 4 and α to equation (34). The functional form is in good agreement with the simulation data for the set of parameter values b 1 = 4.1, b 2 = 11.4, b 3 = −1.26, b 4 = 3.46 × 10 −7 and α = −0.07. In fact, we found that none of the terms in equation (34) are negligible; a fit cancelling the bulk (fixing b 1 = 0 and b 2 = 0, shown as a dash-dotted line in figure 10(b)), the corner flow (b 3 = 0, dashed line) or the contact line dissipation (b 4 = 0, dotted line) gives a qualitatively poor representation of the simulation data. Therefore, we can conclude that all sources of dissipation, corresponding to the bulk, corner flow and contact line, contribute to determining the global friction coefficient of the barrel, and thus to the relaxation timescale.

Conclusions
We have studied the translational relaxation of liquid barrels in wedges towards equilibrium using a diffuse-interface model.
Within our lattice-Boltzmann simulation method, we have introduced a smoothing algorithm to model the interaction of a liquid barrel in contact with a solid wedge geometry. This algorithm can be used to model the wetting behaviour of two fluid phases in contact with solid boundaries of an arbitrary shape. We have introduced an alternative method to compute the energy landscape of a capillary system by integrating the dissipation function in time. This is advantageous because it avoids a direct computation of the surface area of the interfaces as is done in sharp-interface approximations, or the free energy at the solid walls in the diffuse-interface model.
We have validated our simulation method by analysing the equilibrium state of the liquid barrels. We find a quantitative agreement with previous theoretical [15] and experimental results [17]; the shape of the barrel in equilibrium is a truncated sphere, whose position and radius obey the expected dependence on the wetting angle and the opening angle of the wedge.
We have analysed the relaxation of the liquid barrels towards equilibrium. The motion of the liquid is driven by a distribution of the curvature of the interface which creates a pressure gradient. The resulting flow field in the bulk of the droplet is laminar. Near the contact lines, the flow field changes to the treadmill pattern described by [54]. The motion of the contact lines is due to differences in the chemical potential caused by the out-of-equilibrium interface curvature [21].
We have compared our simulations to the model of Ruiz-Gutiérrez et al [14], who derived an expression of the relaxation time of the liquid barrel including the effect of the hydrodynamic dissipation of the bulk flow and the corner flow near the contact line, and the dissipation arising from the motion of the contact line itself. We have identified the scaling of the relaxation time in terms of the diffuse-interface model parameters. Our results confirm the presence of the three contributions to the relaxation time. The relative contribution of the contactline and corner-flow dissipation (to the bulk dissipation) is governed by the size of the interface. In our simulations this length scale is two orders of magnitude smaller than the size of the droplet, which contrasts with millimetre-sized droplets of molecular liquids (such as water) where the interface thickness is several orders of magnitude smaller than the typical macroscopic length scale of the flow. Qualitatively, however, the scaling relation for the relaxation time is expected to hold.
We hope that our results will help guide experiments to identify the contribution of the different sources of dissipation in this system. We acknowledge that, in practice, two important phenomena are likely to occur: contact angle hysteresis and gravity. On the one hand, hysteresis emerges from either physical or chemical inhomogeneities that produce weak pinning points on the solid surface. We expect that these have the effect of truncating the relaxation to equilibrium in a stochastic relaxation. On the other hand, strong body forces, such as gravity, that distort the shape of larger droplets would require investigation beyond the liquid barrel approximation. In both cases, these phenomena deserve their own treatment and we hope to have inspired further investigation on these lines.