{'fraction_non_alphanumeric': 0.025629273137007753, 'fraction_numerical': 0.009396126962300905, 'mean_word_length': 4.759348878134624, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Contact angle measurement on curved wetting surface in multiphase lattice Boltzmann method

Contact angle is an essential physical quantity that characterizes the wettability of a substrate. Although it is widely used in the studies of surface wetting, capillary phenomena and moving contact lines, measuring contact angles in experiments and simulations is still complicated and time-consuming. In this paper, we present an efficient scheme for the real-time and on-the-spot measurement of contact angles on curved wetting surfaces in lattice Boltzmann simulations. The measuring results are in excellent agreement with the theoretical predictions by the spherical cap method without considering the gravity effect. A series of the simulations with various drop sizes and surface curvatures confirm that the present scheme is grid-independent. Then, it is verified in gravitational environments by simulating the deformations of sessile and pendent droplets on the curved wetting surface. The numerical results are highly consistent with experimental observations and support the theoretical analysis that the microscopic contact angle is independent of gravity. Furthermore, the scheme is applied to capture the dynamic contact angle hysteresis on homogeneous or chemically heterogeneous curved surfaces. Importantly, the accurate contact angle measurement enables the mechanical analysis at moving contact lines. The present measurement is simple and efficient, and can be extended to implement in various multiphase lattice Boltzmann models.


I. INTRODUCTION
Contact angle is an important characteristic quantity used to express surface wettability and has a wide range of applications in nature and industrial production, such as wetting, microfluidics, capillary phenomena, coating technology and moving contact line [1-3].
Experimenters have researched and invented various schemes in order to measure contact angles. One of the earliest and widely used methods is the technique of measuring the contact angle of sessile droplets using a telescopic goniometer, which is based on the principle of measuring the line tangency of the three-phase contact points of a droplet profile on a smooth surface [4]. The angle measured in this way is usually very close to advancing contact angles. By applying explicit vibration [5], equilibrium contact angles can be obtained. McDougall et al. modified the sessile drop method and obtained the advancing and receding contact angles by tilting the solid surface until the droplet just started to move [6].
Subsequently, Extrand and Kumagai used this method to study the contact of liquids on various polymer surfaces [7]. Axisymmetric drop shape analysis-profile is a technique used to measure liquid-fluid interfacial tension and contact angles, and has high precision [8].
Kwok et al. control droplet injection rate or extraction rate and then use the technique to measure low-velocity dynamic contact angle [9]. Langmuir and Schaeffer used the specular reflection of the droplet surface to measure the contact angle [10]. Later Fort and Patterson improved the method and used it for static drops and meniscus on flat plates or inner tubes [11]. Different from the methods where a sessile droplet is formed above a solid sample, the captive bubble method provides a direct measurement of the contact angle of bubble formation in a liquid by forming a bubble below the solid sample and then immersing it in the test liquid [12]. The method of tilting the plate is to immerse one end of the solid plate in the liquid, and rotate the other end toward the liquid surface until it is immersed in the liquid, forming the meniscus on both sides of the plate [13]. The plate is slowly tilted until the meniscus on one side becomes horizontal. The angle between the plate and the horizontal plane is the contact angle. In addition, researchers are not only able to calculate the contact angle by direct measurement, but also by indirect methods. In the Washburn capillary rise method, for example, the contact angle is derived from the rate at which the liquid rises through the powder-filled bed by capillary action [14]. The Wilhelmy balance method is also one of the common methods for measuring contact angles [15], which is an indirect force method that reduces the measurement of angle to a measurement of weight and length. The result of this method is highly accurate and not subjective, and it is also suitable for studying the advancing or receding contact angle and contact angle hysteresis.
In numerical simulations involving surface wettability, the scheme that used a goniometer to measure the contact angle from images generated by the simulated data is subjective and rough. More precisely, image analysis can be employed to obtain the contact angle from images [16]. Sakugawa et al. obtained contact angles by using image processing and polynomial fitting [17]. In order to improve the accuracy of contact angle, Scanziani et al.
and Klise et al. both used X-ray microtomography images to calculate the contact angle [18,19]. Measuring contact angles in low-resolution images is cumbersome, especially as the need to derive fluid images prior to measurement is tedious and time-consuming, and may introduce subjective bias, so is not optimal for in-situ measurements. Without the influence of gravity, surface tension makes the droplet appear as a spherical cap shape on a plain surface. The contact angle can be calculated by measuring the height and bottom width of the droplet [20,21] , and this theoretical method is known as the spherical cap method. When the droplet is reduced to the nanometer scale, since there is no stable interface between gas and liquid, the descending contour needs to be fitted by the least square method [22,23]. The spherical cap method is simple and achievable, but it cannot be used in gravity or nonequilibrium environments. Subsequently, researchers have gradually expanded the study of contact angles from droplets to fibres and porous media, and have proposed several methods to measure their contact angles from an energy perspective. Amrei et al. studied the variation of rough fibre contact angle with fibre roughness by means of an energy minimization method [24]. Blunt et al. determined the contact angle of three-phase flow in porous media by using energy balance [25]. Jasper proposed a general variational method for predicting contact angles considering the Laplace pressure case [26]. However, these are domain-specific methods and are not generic.
There are also many studies in numerical simulation on the modification of boundary conditions to improve the accuracy of contact angle measurements. For the simulation of diffusion interfaces, prescribed contact angles can be obtained by using geometric formulas for the wetting conditions [27]. To improve the accuracy and stability of the contact angle boundary conditions, Lee et al. use characteristic interpolation to obtain contact angles [28].
Dong further extended the boundary conditions of the contact angle after considering the relaxation of the dynamic contact angle to simulate dynamic wall-confined gas/liquid flows with large density ratios [29]. Leclair et al. used Dirichlet boundary conditions to study incompatible two-phase pore-scale suction and discharge forces in porous media using the desired contact angle imposed at the boundary [30,31]. These methods aim to modify the wetting boundary condition to impose an accurate contact angle, rather than improve the algorithm of the contact angle measurement. These imposing procedures of contact-angle boundary conditions are computationally complex and nonlocal. Especially, they involve the intervention to the evolution of flow field. Essentially, a contact angle is a geometrical concept. The Youngs equation can only theoretically explain some special cases, such as a sessile drop on a flat substrate at zero-gravity mechanical equilibrium. In dynamic or nonequilibrium environments, the contact angle should be measured through a geometrical method. Recently, Wen et al. proposed a geometry-based contact angle measurement on a plain substrate. The simulation results showed that the method was accurate and efficient [32]. Nevertheless, experiments, natural phenomena and industrial applications often involve complex boundary shapes and even soft substrates, which appeal to an effective and realtime scheme to measure the contact angle on curved surfaces.
In this paper, we design a in situ method for contact angle measurement on curved wetting surface. In Section II, we introduce the lattice Boltzmann method and the chemical-potential multiphase model. Section III describes in detail the measuring method contact angle on curved wetting surface and the chemical-potential boundary condition. In Section IV, we verify the measurement accuracy and the grid-independence. A series of simulations of sessile and pendent droplets under the effect of gravity indicate that the contact angle is microscopic and independent of gravity. Subsequently, the dynamic hysteresis phenomena of contact angle of droplets on homogeneous surfaces as well as chemically patterned surfaces are computed and analyzed. Based on the accurate contact angle, we can perform the mechanical analyses at the contact line regions of droplets on chemically patterned surfaces.
Finally, Section V briefly summarizes the work.
The lattice Boltzmann equation (LBE) is fully discretized in space, time, and velocity.
The multiple-relaxation-time (MRT) version of LBE improves the numerical stability and computational accuracy, and it can be expressed as [44] where M is a transformation matrix that linearly transforms the distribution functions to the where ω i is the weighting coefficient and u is the fluid velocity.
The lattice Boltzmann Eq. (1) is decomposed into two basic steps of collision and advection, revealing the phenomenon of fluid movement at the meso level.
In the MRT model, its biggest feature is that multiple relaxation times are used in the collision process, and different moments can use different relaxation times. S is a diagonal matrix of non-negative relaxation times: S = diag (0, s e , s ε , 0, s q , 0, s q , s v , s v ) . In this paper, the relaxation times are given by s e = 1.64, s ε = 1.54, s q = 1.9, s v = 1/τ for the simulations with the MRT LBE.

B. Chemical-Potential multiphase model
The chemical potential is the partial differential of the Gibbs free energy to the composition [45]. For a nonideal fluid system, following the classical capillarity theory of van der Waals, the free energy functional within a gradient-squared approximation is written as [46,47] where the first term represents the bulk free-energy density and the second term describes the contribution from density gradients in an inhomogeneous system, and κ is the surface tension coefficient. The general chemical potential can be derived from the free energy density functional [45,48], Gradients in the chemical potential act as a thermodynamic force on the fluid. With respect to the ideal gas pressure c 2 s ρ, the nonideal force can be evaluated by a chemical potential The general equation of state can also be defined by the free energy density, Solving the linear ordinary differential Eq. (8) gives the general solution of the free-energy where C is a constant. When the general expression of equation of state (EOS) is selected, substituting Eq. (9) into Eq. (8) will solve the relevant chemical potential, and the constant is eliminated. For example, the famous Peng-Robinson (PR) EOS and its chemical potential are, and where R is the gas constant, a is the attraction parameter, b is the volume correction parameter, and the temperature function is In our simulations, the parameters are given by a = 2/49, b = 2/21, and R = 1. The acentric factor ω is 0.344 for water. To make the numerical results closer to the actual physical properties, we define the reduced variables T r = T /T c and ρ r = ρ/ρ c , in which T c is the critical temperature and ρ c is the critical density.
A proportional coefficient k is introduced to decouple the dimension unit of the length between the momentum space and the mesh space, namely δx = kδx . Here the quantities in the mesh space are marked by a superscript. Following dimensional analysis, the chemical potential in the mesh space can be evaluated by [49,50] We further apply the central difference method with fourth order accuracy to calculate the gradients. These approaches greatly improve the stability of the chemical-potential multiphase model, and the transformation has no loss of accuracy holding the mathematical equivalence.
In addition, we chose the exact difference method proposed by Kupershtokh et al. to incorporate the nonideal force F into LBE [51]: where δu = δtF /ρ . The body force term F i is simply equal to the difference of the equilibrium distribution functions before and after the nonideal force acting on the fluid during a time step. Correspondingly, the macroscopic fluid velocity is redefined as the average momentum before and after the collision: v = u + δtF /(2ρ) .

III. MESOSCOPIC MEASUREMENT OF CONTACT ANGLE ON CURVED SUR-FACE
A. Chemical-potential boundary condition The chemical potential plays a fundamental role in driving a phase transition or indicating the wettability of a solid surface. In the previous work, the chemical-potential boundary condition is implemented on a plane surface [32,50,52]. Here, we improve the scheme to a curved wetting surface. Fig. 1 presents the three-phase contact region of a drop on a curved wetting substrate. Since the central difference method with fourth order accuracy involves the neighboring nodes whose distances are less than or equal to two lattice units, the chemical-potential boundary condition must treat two layers of solid nodes at least. A specific chemical potential is assigned to these solid nodes in order to regulate the wettability of the substrate. It influences on the gradient calculation of chemical potential on the fluid nodes adjacent to the substrate, and this chemical-potential gradient reflects the interaction between the fluid and the wetting substrate. On the other hand, the boundary condition needs to estimate the densities of the solid nodes in the two layers in order to calculate the density gradient. Fig. 1 marks the first and second layers of solid lattice nodes by cyan and orange, respectively. The densities on the two layers of solid nodes can be calculated based on the nearest neighbor nodes, where x s + e i δ t indicates the adjoining nodes, and s w is a switching function. For the first layer of nodes (in cyan), s w = 1 when x s + e i δ t is a fluid node; for the second layer of nodes (in orange), s w = 1 when x s + e i δ t is in the first layer; otherwise, s w = 0 . The miss distribution functions on the boundary, which stream, in concept, from a solid node to a fluid node, are calculated by the multiphase curved boundary condition with mass conversation [53]. measurement are highly desirable for the local mechanical analyses. As shown in Fig. 1, the liquid-gas interface is defined by the contour line where the density is equal to the mean density of the gas and liquid phase, which is very close to the theoretical interface defined by the equimolar division and has a much simpler calculation [40]. This surface distinguishes the liquid and gas nodes in the transition region of the drop. The intersections between the drop surface and the liquid-gas links can then be obtained by the linear interpolation, where x l and x g represent the liquid and gas nodes of a liquid-gas link respectively, and ρ m is the mean density of the liquid and gas. The two intersections, whose distances are  analysis of droplets on surfaces, the calculation field Dx is extended to 3000 lattice units and the droplet radius is r 0 = 100 lattice units. In this paper, droplets on a curved solid surface with a specific chemical potential are simulated using the PR EOS. The density of the flow field is initialized as follows [54]: where ρ g and ρ l are the gas and liquid coexistence densities obtained using Maxwells equal area method of construction, the initial interface width is W = 10 , r 0 is the initial radius of the droplet, and r = (x − x 0 ) 2 + (y − y 0 ) 2 .
A. Accuracy of contact angle measurement The droplet will have a perfect spherical cap under gravity-free conditions. If the length and height of the bottom of the droplet are L and H, then the radius of the droplet can be calculated R 0 = (4H 2 + L 2 )/8H and then the horizontal angle tan β = L/2(R 0 − H), and then the contact angle can be calculated by the spherical cap method. The spherical cap method is generally used as a benchmark to validate the proposed measurement method for calculating contact angles on a mesoscopic scale, as the base length and height of the droplet can be readily calculated. Fig. 3 shows the contact angles measured by the present method and the spherical cap method at two temperatures from hydrophilic to superhydrophobic surfaces. It can be seen from the Fig. 3 that the results of the present method (black star) agree with those of the spherical cap method (green line). The overall trend for both the present method and the spherical cap method is linear, with the contact angle increasing as the chemical potential increases, but the spherical cap method is intuitively seen to bend in the superhydrophilic and superhydrophobic cases. We therefore fitted the linear type by least squares for the accuracy of this scheme.
To investigate the accuracy of the present scheme at temperature and for different wettability, the contact angles are measured on the hydrophilic and hydrophobic surfaces, and are drawn in Fig. 3  angle and the surface chemical potential is very useful in multiphase simulations, because the surface wettability can be readily adjusted according to actual requirements.

B. Verifications of grid independence
Verification of the grid-independence is crucial in the results of numerical simulations.
We first measure the contact angles of a series of droplets from r 0 = 30 to r 0 = 100 lattice units on five different wettability substrates and find that the contact angles are the same for different droplet sizes. On the other hand, we measure the contact angles of droplets on substrates with the curvature radii from 150 to 400 lattice units. As can be seen in Fig. 4, the contact angles measured for different droplet sizes are the same and the contact angles for the same wettability on substrates with different radii of curvature are highly consistent.

C. Deformations of sessile and pendent drops
Both theory and experiment have verified that gravity does not affect the equilibrium contact angle of a droplet on a uniformly smooth surface [55] . If the effect of gravity is considered, the droplets on the solid surface will undergo deformation gradually deviating from the spherical cap shape. As the radius of the droplet increases, the effect of gravity becomes significant gradually, and the drop deformation is larger and larger. Because the deformation extends the drop footprint radius and lower the drop height, the contact angle computed by the spherical method, which is based on the height and footprint radius, decreases inevitably under gravity. A series of droplets are simulated and their diameters vary from 0 mm to 3 mm. These droplets were located on curved wetting substrates with the contact angles 70 • , 100 • and 140 • , and the corresponding chemical potentials took -0.01, 0.04, and 0.1, respectively. The macroscopic diameter 0 is equivalent to the case of zero gravity. Fig. 5 presents that both of the present scheme and the spherical cap method obtain almost the same contact angles when the droplet diameter is less than 1 mm. This confirms the theoretical prediction by Picknett and Bexon that a droplet resting on a smooth homogeneous surface takes the shape of a spherical cap and the gravity effect is negligible provided that its mass is less than about 1 mg [56]. When the macroscopic diameter of the droplet exceeds 1 mm, the contact angle calculated using the spherical cap method decreases significantly due to a decrease in height and an increase in width, and gradually deviates from the gravity-free value. However, the contact angles measured by the present method is remain the same all the time. This confirms that the present scheme obtains the microscopic contact angle, which is independent of gravity [55]. The deformations of droplets under gravity are shown in Fig. 6. The macroscopic diameters of these two sets of droplets are 2 and 3 mm. Since the surface tensor of the water/vapor system is constant, the larger droplet suffers a larger gravity force and displays more apparent deformation. It can be clearly seen in Fig. 6 that the droplets with the initial droplet diameter of 3 mm are much flatter and shorter than those of 2 mm. The simulation results are consistent with the those in the literature by Xie et al. [57].
A pendent droplet adsorbed on the undersurface of a homogeneous curved substrate is stretched and its footprint radius is contacted by the gravity effect, thus the contact angle calculated by the spherical cap method inevitably increases with the growth of the droplet.  are -0.01, 0.04, and 0.1, respectively. Fig. 7 shows that on the same homogeneous surface, the contact angles of the droplets with different macroscopic diameters remain the same all the time, whereas those calculated by the spherical cap method gradually increase as the growth of the drops. Similar to the sessile drops, the independence between the microscopic contact angle and the gravity effect is verified again.
The obvious difference between a pendent droplet and a sessile droplet is that the pendent droplet will fall off, when its size is large enough so that the gravity force is greater than the adhesion force. Therefore, as shown in Fig. 8, droplets cannot be stretched unceasingly.
When the substrate is more hydrophobic, the pendent droplet is stretched and the footprint radius is narrower. The droplet falls off when the diameter exceeds 2.4, 1.95, 1.65 mm for the substrates with the contact angle 70 • , 100 • and 140 • , respectively.

D. Dynamic contact angle hysteresis
In general, the roughness and chemical heterogeneity of a solid surface can lead to contact angle hysteresis [3]. In this section, we investigate the dynamic contact angle hysteresis on chemically homogeneous and heterogeneous curved surfaces, which is too difficult to capture in experiments and theoretical calculations. As shown in Fig. 9, the contact angles of the where γ represents the liquid-gas surface tension, and θ is the dynamic contact angle and θ eq is the equilibrium contact angle in relation to the surface property. When the liquid/gas transition region is on a homogeneous region, θ eq is equal to the contact angle of the region.
If it locates on a heterogeneous region that is composed of two components, the equilibrium contact angle is evaluated by the modified Cassie-Baxter equation [58] : where θ s1 and θ s2 are the intrinsic equilibrium contact angles for the two components, r represents the roughness of the wetting surface (r = 1 for a smooth surface). ϕ d and 1 − ϕ d indicate the area ratios of the liquid/solid and liquid/gas interfaces. In the context of twodimensional diffuse interface model, ϕ d and 1 − ϕ d refer to the ratio of the length of the isodensity line on the surfaces with contact angles of θ s1 and θ s2 , respectively, to the total length of the isodensity lines. When the contact line is located at the border of hydrophilic and hydrophobic surfaces, we use Eq. (18) to calculate the corresponding equilibrium contact angle. anced Young's force also fluctuates around zero, as shown in Fig. 13. The force exhibits a large jump when the advancing contact line crosses the border of hydrophilic and hydrophobic regions. In the hydrophilic region, the advancing angle is greater than the equilibrium angle because the hydrophobic surface resists the liquid spreading; thus, the Youngs force on this region is much greater than zero. oscillates and is often greater than zero. When the receding contact line crosses the border of hydrophilic and hydrophobic regions, the force also exhibits a large jump. It should be note that since the right is the positive direction, θ A > θ eq leads to a positive unbalanced Youngs force, whereas θ R > θ eq leads to a negative one. The simulation results manifest that the present scheme can accurately evaluate the dynamic contact angle and conduct the in-situ mechanical analysis at moving contact line. We expect to gain further insight into capillary phenomenon and dynamic hysteresis through microscopic contact angle and real-time mechanical analysis.

V. CONCLUSIONS
Wetting and capillarity are ubiquitous in nature. As the most important physical quantity in this field, contact angle illustrates the competitive interactions between liquid, gas and solid surface. It is significant to obtain the accurate contact angle in scientific researches and industrial applications. This paper presents a geometry-based scheme to measure the real-time contact angle on curved wetting substrates. The accuracy and gird independence of the scheme are carefully verified. The theoretical prediction that the microscopic contact angle does not depend on gravity is confirm, even though drops are deformed under gravity.
With the accurate measurement of microscopic contact angle, the dynamic contact angle hysteresis can be captured readily, and the mechanical analyses at moving contact line are implemented in situ. These simulations manifest that the present scheme can be a powerful tool to investigate the issues involving surface wetting, capillary phenomena and moving contact lines. Since it is based on the geometry of contact angle, the present scheme can be used in other multiphase lattice Boltzmann models, such as pseudopotential models [21,43,59], field phase models [48,60], etc.