The receding contact line cools down during dynamic wetting

Wetting phenomena are ubiquitous in our daily life, in nature and in various scientific and engineering fields. In particular, the behavior of the contact line (CL), where a liquid-vapor interface meets a solid surface, has long been a topic of interest because it plays a key role in wetting properties.^1–4 For static wetting without CL motion, a common measure of wettability at the macroscopic scale is the contact angle (CA) described by Young’s equation,⁵ which was first proposed in 1805 based on the balance among solid–liquid, solid–vapor and liquid–vapor interfacial tensions. These interfacial tensions originate from the microscopic molecular interaction forces,⁶ and recent molecular simulation studies have provided significant advance in understanding of static wetting, e.g., mechanical^7–9 and thermodynamic^10–12 extraction of the interfacial tensions and analysis on line tension.^13–15 

The situation is more complex at the dynamic CL (DCL)—appearing typically during droplet spreading or moving on a substrate, where the advancing and receding CAs are different. To model the CA difference, a number of theoretical,^16–19 computational^20–25 and experimental^26–28 studies about the DCL have been carried out and have indicated that this dynamic effect is induced by the viscosity and friction in the vicinity of the DCL; however, the governing principle of the DCL motion still remains unclear mainly due to the lack of detailed information on the nanoscale thermal and flow fields around the DCL, and it is considered to be one of the long-standing unsolved problems of fluid dynamics.

In this article, we show a unique thermal phenomenon around the DCLs, cooling as well as heating at the DCLs. To elucidate its mechanism, we examined the heat flow field around the DCL using molecular dynamics (MD) simulations (indeed, as a method based on Newtonian mechanics, MD enables a clear connection between the microscopic molecular motion and the macroscopic conservation laws). To that aim, we have developed a heat transport analysis methodology applicable in multi-component MD systems based on the energy conservation law. We applied this methodology to a quasi-2D system with liquid–solid–vapor CLs, consisting of a Lennard-Jones (LJ) fluid between parallel solid walls moving in opposite directions as shown in the top panel of Fig. 1. Indeed, the LJ model captures the generic structural and dynamical properties of a large class of liquids, bound, e.g., by van der Waals or metallic interactions, and can also be used to describe more complex molecular liquids through the corresponding states approach. It is therefore useful to explore its behavior, as a reference to then explore the specificities of more complex liquids.^29–32 

The top panel of Fig. 1 shows the MD simulation system of a quasi-2D Couette-type flow, where the basic setups are the same as in our previous study.⁴⁰ The fluid–fluid and fluid–solid interactions were modeled by the 12-6 LJ potential $ΦLJ(rij)=4εijσijrij12−σijrij6$⁠, where r_ij is the distance between the particles i and j, while ɛ_ij and σ_ij denote the LJ energy and length parameters, respectively. Quadratic functions were added to this LJ potential so that the potential and interaction force smoothly vanished at a cut-off distance of r_c = 3.5σ.⁴³ We used the following parameters for fluid–fluid (ff) and fluid–solid (fs) interactions: σ_ff = 0.340 nm, ɛ_ff = 1.67 × 10⁻²¹ J, σ_fs = 0.345 nm, ɛ_fs = 0.646 × 10⁻²¹ J. The atomic masses of fluid and solid particles were m_f = 39.95 u and m_s = 195.1 u, respectively. Finally, the equations of motion were integrated using the velocity-Verlet algorithm, with a time step Δt of 5 fs.

Periodic boundary conditions were set in the x- and y-directions, and 20 000 LJ particles were confined between two parallel solid walls (dimension of x × y = 39.2 × 3.92 nm²) at a distance of $∼10.4$ nm, so that the LJ fluid formed two quasi-2D menisci with CLs on the walls upon the preliminary equilibration at a control temperature T_w = 85 K, by using a Langevin thermostat applied to the second outmost layers of the solid walls. The static CA on both top and bottom walls was ∼57°. After the equilibration, further relaxation runs to achieve a steady shear flow with asymmetric menisci were carried out for 10 ns by moving the particles in the outmost layers of both walls with opposite velocities of ±10 m/s in the x-direction. The above-mentioned temperature control only on the wall away from the fluid was applied so that the dissipation energy due to the flow in the liquid phase was removed, and a steady state was achieved as a result.

After the relaxation run, the main calculation was conducted for an average time of 400 ns. We calculated the external body force and volume averaged velocity, i.e., streaming velocity, on the RHS of Eq. (3) using cuboid bins of size Δx ×Δy ×Δz = 0.150 × 3.92 × 0.149 nm³, while we calculated the energy flux, velocity, stress and the specific energy on the RHS of Eq. (3) using the MoP with the faces of each local bin. Regarding the calculation of the energy flux and the specific energy, see details in supplementary material. We also obtained the temperature distribution by calculating the thermal energy in the local bin, i.e., the kinetic energy by using the relative velocity of molecules to the streaming velocity.

The middle panel of Fig. 1 shows the density distribution and velocity field obtained by the volume average. Due to the shear applied by the wall, a caterpillar-like flow was induced, and DCLs with different CAs, i.e., advancing and receding CLs, appeared. In addition to the CA difference, we showed the stress inhomogeneity in the bulk liquid induced by this flow in our previous study.⁴⁰ In the present study, we report a distinct thermal difference in the DCLs as shown in the temperature distribution in the bottom panel of Fig. 1: temperature rises around the advancing CLs (bottom right and top left), and temperature drops around the receding CLs (bottom left and top right). Quantitatively, in the bulk liquid away from the interfaces, the temperature is around 86.5 K, which is slightly higher than the control temperature of the wall due to viscous dissipation, whereas that around the advancing CLs is about 2 K higher and that around the receding CLs is about 2 K lower than in the bulk, as shown in Fig. 1. The cooling at the receding CLs is especially intriguing, because viscous dissipation can only induce temperature rise through heat production.

To elucidate the mechanism of the heat production/absorption around the DCLs, we conducted a heat flux analysis. The top panel of Fig. 2 shows the heat flux field superimposed on the temperature distribution. Note that the heat flux field was depicted only for the fluid sufficiently away from the wall, where the effect of potential field from the wall on the fluid was negligibly small; the MoP methodology of the heat flux calculation is shown in supplementary material. Heat flow from the high temperature area to the low temperature area can be observed, meaning that the heat produced around an advancing CL induces temperature rise there, and flows to the cold neighboring receding CL due to the heat absorption. To quantitatively evaluate the heat production/absorption around the DCLs and in the liquid area sufficiently away from the CLs, we set three control volumes (CVs) as shown with magenta rectangles in Fig. 2; a CV surrounding the receding CL, a CV surrounding the advancing CL, and a CV between them. We calculated the integral of the divergence of the heat flux in each CV by using Eq. (3) and the corresponding values are shown inside the CVs with a unit of mW/m, calculated as the heat production/absorption rate divided by the system depth, in the top panel of Fig. 2. According to Fig. 2, heat is produced and absorbed at the CVs surrounding the advancing (97 mW/m) and receding (−89 mW/m) CLs, respectively, and the absolute values are approximately twice as large as the heat production 42 mW/m in the middle CV even though its volume is twice as large as the others. From this, we see that viscous dissipation is not the main cause of heat production/absorption at the DCLs.

We also show the internal energy distribution and the velocity field as the background of bottom panel of Fig. 2, and one can observe that the internal energy changes along the streamline specifically near the DCLs, where heat is produced/absorbed. At the advancing CL, heat is produced when the fluid flows from the solid–vapor and liquid-vapor regions to the solid–liquid region, whereas at the receding CL heat is absorbed when the fluid flows from the solid–liquid region to the solid–vapor and liquid–vapor regions. During these processes, the internal energy of the fluid changes due to the surrounding density change as well as due to the potential field induced by the solid surface, and it leads to the cooling and heating at the DCLs. This heating/cooling phenomenon can be interpreted due to the change of the potential field along the pathline, meaning that along the pathline into/out of the solid–liquid interface. In that sense, this phenomenon is reminiscent of latent heat, which induces the heat production/absorption upon the phase change.

Considering this feature, it is expected that this cooling and heating effect should be increased with the flow rate around the DCLs, i.e., the faster wall speed. Here, we additionally conducted the heat analysis for the CVs with various wall speeds: u^w = 1.0, 2.5, 5.0, 7.5, and 12.5 m/s (the density and velocity fields and temperature distribution with each condition are shown in supplementary material). Top, middle and bottom panels of Fig. 3 show the volume integral values of the heat flux divergence, stress work term in Eq. (6) and internal energy change along the streamline in Eq. (7), respectively. Blue and red solid lines denote the values of CV including the receding CL (RCL) and the advancing CL (ACL), and green one denotes the CV between them. Note that the CV arrangement for all wall speeds is the same as that in Fig. 2. Similar to Fig. 2, the internal energy change is the main part of heat production/absorption in CVs including DCLs while the stress work term is dominant in the middle CV (referred to as “Bulk” in Fig. 3). The internal energy change appears to be proportional to the wall speed, implying that the spatial distributions of the density and the specific energy do not largely change depending on the wall speed. On the other hand, the stress work term appears to be proportional to the square of the wall velocity in the middle CV since the shear stress, i.e., viscous stress, is proportional to the shear rate in the bulk where the shear rate can be roughly proportional to the wall velocity. Also in the CVs including the DCLs, the work done by the solid–fluid interaction force in Eq. (6) should largely depend on the wall speed because that frictional force is supposed to be proportional to the slip velocity. Under the present wall velocity range where the steady-state caterpillar-like flows with DCLs are achieved, the internal energy change is always dominant over the stress work term in the CVs including DCLs: and thus the temperature rise/drop near the DCLs should always exist. In addition, we observed these cooling/heating at the DCLs induced by the same mechanism also on less wettable walls as shown in supplementary material. Note that this quasi-latent heat around the DCLs is not a dissipation energy, meaning that it cannot be included in the dissipation terms of existing macroscale DCL models.^1,17,19 However, it indeed induces the temperature change in the vicinity of the DCLs, which should be included in the DCL models.

In this article, we have presented a heat transport analysis methodology applicable in multi-component MD systems, which we have used to investigate the heat transport features of the DCL. The heat analysis revealed that heat is not only produced but also absorbed around DCLs, mainly due to the quasi-latent heat induced by the internal energy change of fluid along the pathline, when the fluid moves among the interfaces, which is accompanied by a change in fluid–fluid and fluid–solid interaction energy. In addition, this latent heat is not a dissipation energy, thus almost the same heat is absorbed and produced at receding and advancing CLs, respectively, while heat is only produced in bulk liquid due to viscous dissipation. Overall, these results provide new insights into the molecular mechanisms controlling the dynamics of the CL. Moreover, the framework for analyzing the heat transport at the molecular scale should be useful for investigating various nanoscale systems such as the flow in carbon nanotubes or in nanoporous media.^44–51 

The supplementary material contains the calculation methods of the energy density, energy flux and heat flux by the Method of Planes and the derivation of the Lagrangian derivative of the internal energy in Eq. (4). We also show therein the density, velocity and temperature distributions with various wall speeds on the lyophilic walls corresponding to Fig. 3, and the density and temperature distributions around the DCL on the lyophobic walls.

H.K., T.O. and Y.Y. are supported by JSPS KAKENHI Grant Nos. JP23KJ0090, JP23H01346, and JP22H01400, respectively. Y.Y. is also supported by JST CREST Grant No. JPMJCR18I1, Japan. Numerical simulations were performed on the Supercomputer system “AFI-NITY” at the Advanced Fluid Information Research Center, Institute of Fluid Science, Tohoku University.

The authors have no conflicts to disclose.

Hiroki Kusudo: Conceptualization (equal); Data curation (lead); Funding acquisition (equal); Investigation (lead); Methodology (lead); Resources (equal); Software (lead); Visualization (lead); Writing – original draft (lead). Takeshi Omori: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Supervision (supporting); Writing – original draft (equal). Laurent Joly: Methodology (supporting); Supervision (supporting); Writing – original draft (supporting). Yasutaka Yamaguchi: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Project administration (lead); Resources (equal); Supervision (lead); Visualization (supporting); Writing – original draft (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.