Closed-Form Expressions for Contact Angle Hysteresis: Capillary Bridges between Parallel Platens

: A closed form expression capable of predicting the evolution of the shape of liquid capillary bridges and the resultant force between parallel platens is derived. Such a scenario occurs within many micro-mechanical structures and devices, for example, in micro-squeeze flow rheometers used to ascertain the rheological properties of pico- to nano-litre volumes of complex fluids, which is an important task for the analysis of biological liquids and during the combinatorial polymer synthesis of healthcare and personal products. These liquid bridges exhibit capillary forces that can perturb the desired rheological forces, and perhaps more significantly, determine the geometry of the experiment. The liquid bridge has a curved profile characterised by a contact angle at the three-phase interface, as compared to the simple cylindrical geometry assumed during the rheological analysis. During rheometry, the geometry of the bridge will change in a complex nonlinear fashion, an issue compounded by the contact angle undergoing hysteresis. Owing to the small volumes involved, ascertaining the bridge geometry visually during experiment is very difficult. Similarly, the governing equations for the bridge geometry are highly nonlinear, precluding an exact analytical solution, hence requiring a substantial numerical solution. Here, an expression for the bridge geometry and capillary forces based on the toroidal approximation has been developed that allows the solution to be determined several orders of magnitude faster using simpler techniques than numerical or experimental methods. This expression has been applied to squeeze-flow rheometry to show how the theory proposed here is consistent with the assumptions used within rheometry. The validity of the theory has been shown through comparison with the exact numerical solution of the governing equations. The numerical solution for the shape of liquid bridges between parallel platens is provided here for the first time and is based on existing work of liquid bridges between spheres.


Introduction
Capillary forces arising from surface tension is a frequently encountered phenomenon. This tension is due to the imbalance of forces that act on individual molecules that reside at the interface between two media. These forces may include London dispersion forces, hydrogen bonds, dipoledipole interactions, dipole-induced dipole interactions, π-bonds, donor-acceptor bonds, and electrostatic interactions [1] and may include others. How these forces interact with the molecules and each other is complex and is not completely understood. This is because they depend very much on the substances under investigation and so many conclusions are based on empirical relationships [2].
The study of these forces is of great practical importance especially in the powder technology and related industries [3][4][5][6]. This is because moisture leads to the formation of liquid bridges between particles and so significantly affects the strength of agglomerates. Because of this, studies in this field have largely concentrated on liquid bridges between spheres [3][4][5][6][7][8][9] or between a sphere and a plane [10][11][12]. A few solutions have been proposed each using a different approximation. However, in general these solutions are not easily adaptable for modelling capillary bridges between two planes. The study of this geometry is of increasing importance as it can explain stiction in MEMS devices [13,14] and head/disk systems [15], be used for particle transfer in microfluidics channels [16], and is the geometry increasingly being used for squeeze flow rheology [17][18][19][20][21][22]. To this end, several investigators have studied capillary bridges between flat parallel surfaces [23][24][25][26][27].
It is the effect that the shape of the capillary bridge has on squeeze flow rheology that is of interest here. Generally, analysis of squeeze flow assumes that the height of the fluid layer is much smaller than its radius. This is known as the lubrication approximation [28] and allows for the assumption that the fluid flows radially. In using this approximation, it is then assumed that edge effects can be neglected. This is because dealing with the boundary conditions caused by the free surface is not a trivial matter [2,29]. However, in rheometry, the radius of the fluid is finite and so surface tension effects need to be taken into account as they influence the geometry of the fluid flow as well as imparting an additional force to the rheometer [30]. The change in the geometry of the liquid bridge caused by the movement of the platens of the rheometer manifests itself as a change in curvature of the bridge and the movement of the three-phase boundary between the liquid, the surrounding vapour phase, and the solid platen. This movement is further complicated by the effect of contact angle hysteresis which prevents the three-phase boundary from moving smoothly [26,27,31].
In order to predict the fluid flow and the forces acting on the rheometer, it is necessary to know how the shape of the liquid bridge evolves whilst it is being acted upon by a squeeze flow rheometer. To this end, this paper focuses on the derivation of a closed form expression capable of predicting the evolution of the shape of the liquid bridge during squeeze flow rheometry, including the effects of contact angle hysteresis and hence the capillary force. The validity of the theory will then be shown through comparison with an exact numerical solution.

Theory
In squeeze flow rheometry, when a fluid sample is placed between two parallel platens it forms a liquid bridge with a specific geometry, for instance, see Figure 1. This is because surface tension acts to minimise the surface area of an interface. This minimisation of the surface area is described by the Young-Laplace equation [31]. This equation equates the pressure difference, ∆ , across an interface between two fluids to the resulting curvature due to surface tension. For the geometry shown in Figure 1, this is most simply expressed in terms of the principal radii of curvature at the plane of symmetry of the liquid bridge, i.e., at the 'neck' of the bridge: where is the surface tension between liquid and its vapour, and the other terms were defined in Figure 1. This pressure difference and surface tension determine the capillary force that will act on the rheometer: It can be seen that this force is simply the sum of the tension acting on the surface, denoted by surface tension multiplied by the circumference of the liquid bridge, and the pressure difference acting on the cross-sectional area, see Figure 2 for a schematic.

The Toroidal Approximation
The Young-Laplace equation, as given by eq. 1, cannot be solved analytically except in a few special circumstances, such as planar and cylindrical geometries [10,32]. As such, it is common practice to use numerical methods to find the solution [6]. However, a few closed form approximate solutions exist. These include the gorge approximation [6], the parabolic approximation [8], and the toroidal approximation [33]. As the parabolic approximation requires direct physical measurement of the liquid bridge (which is impractical in this case) and the gorge method leads to high inaccuracies [6], the toroidal approximation is considered the option most worth investigating. The toroidal approximation assumes the profile of the liquid bridge to be a circular arc when the Bond number is small. It has been noted that the error of the toroidal approximation can approach 10% when the bridge gap is large compared to the neck radius [34] and so a further assumption ⁄ > is sensible. It can be seen in Figure 3 that this approximation is a reasonable one. Here, as an example, the profile of a liquid bridge between two flat parallel platens with a given contact angle (45° in this case), volume ( × nL) and Bond number (0.002) for a range of curvatures (9500 to 52,500 m −1 ) is shown. This profile was calculated through the numerical integration of eq. 2 using the method provided in the appendix. Also shown are circular arcs fitted to the 'true' profile representing the toroidal approximation. The circular arcs of the toroidal approximation are a good representation of the bridge profile geometry, which provides justification for the approach used here. It is necessary to remember that the profile of the liquid bridge is not actually a circular arc. For the geometry to be valid, the mean curvature must be constant at all points on the surface of the bridge. This is not true of circular arc geometry where the curvature varies along the axis of the bridge, suggesting that a pressure difference exists within the bridge. This, of course, cannot happen or else the fluid will flow to negate the pressure difference. This means the bridge profile is not circular, but it is still a close approximation. The toroidal approximation has been used to describe the geometry of liquid bridges between spheres [35] but has not been solved completely to deal with the flat on flat geometry as is of interest here.
To use the toroidal approximation, consider the geometry for a liquid bridge between two flat platens with a gap, , shown in Figure 4: From simple geometry, the bridge profile, assumed to be a circular arc with radius , can be expressed by the equation: where the radius of curvature of the bridge profile is defined as: The neck radius, which is the narrowest point of the liquid bridge in the vertical axis, can therefore be defined as: By using the disc method for calculating the volume of a revolution, the volume of the liquid bridge can be shown to be: We can solve eq. 6, and so the neck and contact radii for a given gap, contact angle, and volume for the bridge can be shown to be: This means we can calculate the capillary force as a function of gap for a given volume and contact angle using the definition given by eq. 9: where the neck radius and the radius of curvature are those approximately given by eq. 7 and 4.

Contact Angle Hysteresis
The above discussion is only valid when the contact angle remains constant. When the contact line is in motion, such as when the liquid in the bulk phase is flowing, as in during squeeze flow rheology, different mechanics apply. Many aspects of this dynamic wetting are poorly understood, and the subject is of great scientific interest [36]. What is known is that when the contact line moves due to flow within the bulk of the liquid, hysteresis is observed in the contact angle. As the contact line moves towards the vapour, the contact angle will become larger than the equilibrium angle; this is known as the advancing contact angle. Similarly, when the contact line contracts, the contact angle is smaller than the equilibrium value, and this is known as the receding contact angle.
It has been found that the contact line does not move smoothly, but instead pins and slips [27,37]. This is because there is a transition period as the contact angle varies between its advancing and receding state. In the context of the rheometer, this contact line will undergo four stages of motion:

1)
When the platen moves downwards and the gap decreases, initially, the contact line remains stationary until the contact angle equals the advancing contact angle (see Figure 5a). This is known as pinning.
2) If the gap decreases further, the contact angle remains at the advancing angle and the contact line expands (see Figure 5b). This is known as slipping.
3) When the gap increases, the contact angle first decreases while the contact line remains constant (pins-see Figure 5c). 4) When the contact angle reaches the receding angle, provided the gap is still increasing, it will remain at that angle while the contact line declines (slips-see Figure 5d). And so on.
When the contact line slips, the contact angle remains constant and so the discussion above is valid. However, during the pining regime, the contact angle varies and so eq. 6 would need to be solved for as a function of the gap, . This is not trivial, and an analytical solution cannot be derived. To find a solution, a numerically derived expression needs to be found. It is easy to see that the variation of the contact angle as a function of the gap will depend on the volume of the bridge and the contact line radius. Therefore, in order to obtain a single expression that will cover all situations, the equations will be non-dimensionalised by dividing all lengths and radii by the cube root of the bridge volume: It is also necessary to normalise the gap between the theoretical extreme gaps for a liquid bridge with constant contact line radius and volume. These extremes correspond to the situation where the contact angle is at a maximum (180°) or a minimum (0°). These cases relate respectively to the (i) minimum gap, where any further decrease would cause the contact line to slip, and (ii) the maximum gap, where any further increase would cause the contact line to slip, as illustrated in Figure 6. A real liquid will generally oscillate between limits in-between the theoretical limits. where: The normalised dimensionless gap between the platens, * , is thus defined as: * = * − * * − * There is also another case of special interest, which corresponds to a contact angle of 90°. When this is the case, the bridge assumes the shape of a cylinder, and is the point where the profile changes from being concave to convex and visa-versa. It is useful to note that this case represents a singularity in the solution of eq. 6, as it means that eq. 4 equals infinity (i.e., the bridge profile is straight).
There are multiple solutions to eq. 6 for a given contact line radius, but the solution which coincides with the smallest gap corresponds to that of the most stable bridge, where the effects of drainage is minimal [38]. However, when the dimensionless contact line radius, * , is less than 1, the only solutions that exist correspond to potentially unstable bridges with gaps between the platens much larger than for bridges with larger contact line radii. Due to the intractability of the solution of eq. 6, the contact angle as a function of the normalised dimensionless gap was found numerically by using a bisection method. The result can be seen in Figure 7: The curves in Figure 7 have the shape of a generalised logistic curve, which has the form of eq. 14 (in units of degrees). * = + + * ⁄ where is the lower asymptote, is the distance between the asymptotes, affects near which asymptote maximum growth occurs, is the growth rate and is the point of maximum growth. These variables were all found to be as functions of * . is given in terms of a polynomial, and , , , and can be given by Gompertz curves, the coefficients of which were found, using a standard curve fitting algorithm, to be: * = . .
The quality of the fit of the coefficients can be seen in Figure 8, where the coefficients of determination for each of the fits is 0.9987 or better.
The geometry of a liquid bridge undergoing a change in height can now be calculated for the two circumstances, namely when the contact line slips and when it pins. The procedure for any subsequent calculations would be as follows: • State initial conditions • We need to know four things: gap, initial contact angle, contact line radius, and volume.
• Gap is a design/experimental parameter and is known, equilibrium (and dynamic for later) contact angle is a fluid property and can be found easily beforehand. • The volume or the contact line radius needs to be measured. The other can be calculated using eqs 6 or 8, respectively. However, they only need to be measured once. The values can then be used throughout the calculation.
• Change the gap • It is assumed that subsequent motion will be continuous and no relaxation towards equilibrium conditions is observed. • The change in geometry depends on the direction the upper platen is moved.
• If it is moving upwards, initially, the receding contact angle will be lower than the equilibrium contact angle. Here, the contact line radius (given) will stay constant until the instantaneous contact angle (given by Equation 16) equals the receding contact angle. If the motion continues to increase, the contact angle stays constant and the contact line radius changes according to Equation 8.
• Change the direction • At this point, the contact angle equals the receding contact angle and the contact line radius is known. Decreasing the gap, the contact line radius will remain constant and the contact angle will increase, as in Equation 16, until it reaches the advancing contact angle. Decreasing the gap further will mean the contact angle stays at the advancing angle and the contact line radius changes.
Using this procedure, the capillary forces can be found, as demonstrated in Figure 10 for two scenarios covering a wide range of contact angle hysteresis. The maximum relative error between the force as calculated using the toroidal approximation and that calculated numerically from the exact expression using the algorithm given in the Appendix, is 0.21% for the data shown in Figure 10(left) and 1.3% in Figure 10(right).
Significantly, the data in Figure 10(left) (that is, the 120 data points used to describe one cycle here) took c.a. 880 s to generate using MATLAB R2018a on a HP Omen Laptop with Intel ® Core™ i7-9750H CPU @ 2.60GHz, with 6 Cores, 12 Logical Processors and 16 GB of RAM. The equivalent data generated on the same system using Equation 14 took just c.a. 0.12 s. This demonstrates the significant speed-up the use of these closed-form expressions can give the user as compared to having to solve the exact equations numerically.
It can be clearly seen that the form of these curves is very similar to those measured experimentally by De Souza et al. [26] and by Shi et al. [27]. Direct comparison is not possible between these studies and the theory proposed here as those studies used asymmetric bridges with large bond numbers and large extensions/velocities in which viscous forces and bridge collapse is significant and not of interest here. For the case shown in Figure 10, the contact line radius and contact angle, as calculated using the toroidal approximation, vary as a function of gap and are as shown in Figures  11 and 12: This relationship may be better appreciated as a function of time:  Figure 10 (left). In this case, a frequency of 10 Hz is arbitrarily assumed.

Linearization and Applicability to Squeeze Flow Rheometry
For this theory to be useful for squeeze flow rheometry, it is important that it is consistent with the fluid flow theory. In the fluid flow theory, it is assumed that the fluid has cylindrical geometry, the flow is purely radial, and that the no-slip condition applies [30]. These assumptions may seem to be at odds with the theory presented here. The reason for this is that, here, it is stated that the boundary is curved due to surface tension, precluding cylindrical geometry, and that the contact line can move, which will cause infinite shear stresses at the interface if the no-slip condition is applied and the flow is purely radial. This is because dealing with the exact boundary conditions in squeeze flow rheometry is not a trivial matter. This discrepancy is negated by ensuring that the lubrication approximation is applicable. The lubrication approximation states that if the fluid gap is much smaller than its radius, perturbations in the radial flow caused by edge effects can be neglected [17].
This means the individual force contributions from the surface tension and viscous flow can be treated separately and summed together [30].
One way to maintain consistency is to ensure the volume of liquid under consideration in the two theories is the same. To do this, the radius of fluid used in the fluid flow theory must be the rootmean-square of the radius of the bridge given in Equation 3, or: To see why this is consistent, recall that the volume of the bridge is defined as: And the volume of a cylinder with diameter and height is defined as: If the definition of the root-mean-square radius, , from Equation 17 is substituted into Equation 19 as the radius of the cylinder so that = , it can be seen that Equation 18 and Equation 19 are in fact the same (see Equation 10).
Therefore, matching the value for the bridge radius does ensure that the volume of liquid used in both theories are consistent.
It was shown in [30] that if the capillary force can be described as a linear function then the equations of motion for a squeeze flow rheometer can be found. To linearise Equation 9, the change in the gap, , must be sufficiently small as to not cause the contact line radius to slip. While the capillary force is still non-linear, if the change in gap is small enough, this nonlinearity will be small. This means the capillary force can be approximated by the linear relationship given in Equation 21, as demonstrated in Figure 13.
The coefficients used in Equation 21 are: = And: where is the mean or nominal bridge gap. While these expressions are difficult to solve analytically, they are easily solved numerically. For instance, for the data shown in Figure 13, = 179.4 N/m and = -0.0017 N. The maximum relative error between the linear function given by Equation 21 and the data shown in Figure 13 is 0.11%. Therefore, Equation 21 can be used to determine the correction to the squeeze force measured during rheometry due to capillary forces and to check that these forces are linear, facilitating analysis of the rheometry data. Figure 13. Linearised capillary force. Note that, as an example, the amplitude has been here reduced to 25 nm and, as a result, the contact radius is not able to slip, hence the capillary force is essentially linear.

Conclusions
In this work, the geometry of a pendular liquid capillary bridge between two flat parallel platens has been determined, taking into account the contact angle hysteresis exhibited by the three-phase contact line of the bridge as the height of the bridge varies. The solution was based on the toroidal approximation, which is the assumption that the profile of the bridge is well approximated by a circular arc. This approximation has been extended for all contact angles from 0-180° for liquid bridges characterised by low Bond numbers. A closed form expression has been developed that allows for the fast approximation of the contact angle and contact radius at any bridge height without the need for significant experimentation or numerical analysis. All that is needed to be known before analysis is the bridge volume, the bridge height and the advancing and receding contact angles. This means that the effects of contact angle hysteresis, such as the resultant nonlinear capillary forces, can be calculated very conveniently whilst still maintaining accuracy. The inconsistency between the complex geometry of the liquid bridge and the simple geometry assumed during micro-squeeze flow rheology has been resolved, and the techniques for minimising the nonlinearities in the capillary forces during rheometry have been determined. The validity of the theory has been shown through comparison with the numerical solution of the exact governing equations.

Appendix A
The numerical scheme used in this paper to check the validity of the toroidal approximation, as seen in Figure 3, is similar to that used in [6]. As the boundary conditions were modified to be consistent with the flat-on-flat geometry, the theory is re-encapsulated here. First, the axisymmetric Young-Laplace equation is non-dimensionalised (capitalised letters denote nondimensionalised terms) by dividing all the relevant geometric distances by the characteristic length given by Equation where * is twice the dimensionless mean curvature. The solution of this equation gives the profile of the liquid bridge, which can be approximated by a truncated Taylor series: Eq. A.2 can be integrated to give: where: where * is the dimensionless contact radius. Therefore, the first and second differential, as used in Equation A.3, can be given as: However, in the solution of Equation A.3, the curvature and contact line radius/angle are unknown. Here, a bisection method was used to find the correct geometry. For a maximum and minimum chosen curvature, the volume of the bridge is calculated using the geometry defined by Equation A.3 for a given contact line radius or contact angle. The curvature range is then bisected, and, in each iteration, the calculated volume is compared to the known correct volume until a satisfactory convergence is reached. This in turn defines the gap between the platens. Therefore, if the procedure is repeated for each contact line radius, we obtain the geometry of the liquid bridge for a range of gaps (see Figure A1). Given the geometry for a range of bridge gaps, the curvature and neck radius are known, allowing for the determination of the capillary force through the solution of Equation 2 (see Figure  A2).
(a) (b) Figure A2. How the capillary force varies with the gap as given by Equation 2 for the geometry shown in Figure A1. (a) The contact angle is constant; (b) the contact angle varies.