Dynamic capillary wetting studied with dissipative particle dynamics

We present a study on dynamic capillary wetting in the framework of dissipative particle dynamics (DPD) based on a novel wall model for wetting on solid boundaries. We consider capillary impregnation of a slit pore in two situations: (i) forced (piston-driven) steady state flow and (ii) capillarity driven imbibition out of a finite reservoir. The dynamic contact angle behavior under condition (i) is consistent with the hydrodynamic theories of Cox under partial wetting conditions and Eggers for complete wetting. The flow field near the contact line shows a region of apparent slip flow which provides a natural way of avoiding a stress singularity at the triple line. The dynamics of the capillary imbibition, i.e. condition (ii), is consistently described by the Lucas–Washburn equation augmented by expressions that account for inertia and the influence of the dynamic contact angle.

3 field itself might influence the wetting dynamics rendering the dynamic contact angle problem non-local.
Simulation methods such as molecular dynamics (MD) are in principle capable of representing dynamic wetting problems as a whole, and have extensively been used in the past to shed more light on the molecular processes close to the CL region [20]- [22]. Their application to realistic problems is numerically, however, very costly and therefore limits the size and timescale of the problems to be studied. On the other hand pure continuum methods such as computational fluid dynamics (CFD) methods based on finite elements can conveniently be used to solve the global flow problems in general. However, in the specific problem of free surface flows there is no unique way of coupling the relevant molecular microscopic processes to continuum fluid dynamics [6], [23]- [25].
In recent years novel simulation schemes such as the grid-based Lattice-Boltzmann method (LBM) (see [26] and references therein), kinetic Monte Carlo approaches such as direct simulation Monte Carlo (DSMC) [27,28] or dissipative particle dynamics (DPD) [29,30] have emerged that offer alternatives to address the above mentioned problems. The DPD method was originally invented to model the dynamics of solvents on mesoscopic time and length scales, hydrodynamics is represented by the collision dynamics of an ensemble of pseudo-atoms. Fluid properties but also flow boundary conditions arise as a consequence of the collision dynamics or particle interactions. One may think of these methods as coarse-grained versions of an MD, that still retain some molecular characteristics. Thus, when applied to dynamic wetting problems, they are likely to provide some reasonable model for CL motion, such as Navier-slip [31]. But equally important, they are much more efficient than, e.g. MD in solving the global flow problem. In the present work, we apply a multi-body-DPD (MDPD) approach to a nontrivial transient problem of capillary dynamics, namely, the impregnation of a slit channel out of a finite reservoir (section 5). For representing the sl-system, we employ a recently developed model for solid boundaries, which is reviewed in section 3, and some important results on capillary flow induced by a piston at constant capillary number, referred to as forced wetting are explained in section 4 in greater detail. In the following section, we outline the MDPD method used in this work originally proposed by Warren [32].

A brief account of the scheme
DPD is a mesh-free particle based simulation method that was first devised by Hoogerbrugge and Koelman [29] in 1992 as a generic method for modeling fluid dynamics. In the DPDscheme point particles interact via short range forces [29,33] such that momentum and mass are conserved, which is responsible for the hydrodynamic behavior of a liquid at large scales [34]. A common interpretation is that each particle is considered as a lump of molecules [34]. In this picture the DPD method can be regarded as a coarse-graining of MD and as such is capable of modeling larger length and timescales.
The dynamics of the DPD particles is governed by Newton's second law: F i = m iri , where m i is the mass of the DPD particle. Particles interact via pairwise central forces consisting of a random, dissipative and conservative part: If r i denotes the particle position, the conservative force in standard DPD reads: F C i j = A w C (r i j )e i j , with the repulsion parameter A > 0, r i j = r i − r j , r i j = |r i j | and e i j = r i j /r i j . The weight function w C (r ) = (1 − r/r c ) vanishes for an inter-particle distance r i j larger than a cut-off radius r c . The random force reads: F R i j = qw R (r i j )ξ i j r i j , while the dissipative force damps the relative 4 velocity normal to the connecting line between the particles, v i j · e i j = (v i − v j ) · e i j , which leads to: The random and dissipative force act as a thermostat if the amplitudes q of the random white noise ξ i j and the dissipation constant γ satisfy a fluctuation-dissipation theorem: q 2 = 2γ k B T and w R (r ) 2 = w D (r ). The usual choice for the weight functions is w R = w C . In this context k B is the Boltzmann constant and T is the temperature determined by the DPD thermostat. The resulting force F i on the particle i is the sum of all interactions with the neighbors within the cut-off range r c : . Over the past years several authors have investigated and refined this numerical scheme. Español showed that the Navier-Stokes equations are correctly represented by the DPDmethod [35] and also showed how energy conservation can be incorporated if needed [36]. Marsh [37] used kinetic theory in order to relate macroscopic properties of the DPD-fluid to the model parameters. As a recent development several groups [30,32,38] devised an extension of the simple conservative force law in order to generate more complex equations of state (EOS) showing a van der Waals loop. The latter is commonly realized by augmenting the conservative force F C i j by a density-dependent contribution. This class of schemes is termed MDPD. The density-dependent part can, e.g. be derived from a mean-field approximation to the excess free energy [38]. In an analogous manner though technically different (eventually density gradient terms must be evaluated), also LBM or DSMC can be extended to represent two-phase fluids [27,39].
In the approach of Warren [32] adapted in this work the density-dependent contribution is rather introduced empirically, and has a different cut-off range r d : with an additional weight function w d (r ) = 15/(2πr 3 d )(1 − r/r d ) 2 . In contrast to standard DPD, the density-independent force (with an interaction range r c ) is an attractive force, i.e. A i j < 0, while B > 0 applies for the density-dependent part, rendering it repulsive. The local densitȳ ρ i at the location of particle i is estimated by the instantaneously weighted averageρ i = i = j w d (r i j ) [32]. Warren showed that this approach produces stable free capillary surfaces. An appropriate choice of parameters results in sharply defined liquid-vapor (lv) interfaces that can easily be kept track of in simulations. There are virtually no particles in the vapor phase, which is thus effectively neglected. In the simulation, the numerical cost is only on the dense phase. In this respect, the MDPD approach is easier to handle than, e.g. LBM as there extra effort is required to deal with two-phase flows when the phases have significantly different densities [26]; for a solution to represent an air/water-like interface with a density ratio of 1 : 1000, see [40]. In addition to the original scheme, we introduce a species-dependent force constant A sl , defining the interaction between liquid (l) particles and those of the solid (s) wall in order to establish adjustable adhesive properties of the sl-interaction. All relevant quantities of the MDPD liquid (similar to those in [32]) are summarized in table 1. As in many MD schemes, in DPD dimensionless units are used, i.e. the thermal energy T = k B T , the length scale r c and the mass of the DPD particles m i are all set to unity. Reference to dimensional units can be made by defining the length scale r c , the mass of a DPD particle and the mean thermal energy T of a DPD particle. Different procedures for matching length and timescale have been discussed [43]. One way of matching the length scale is to match the dimensionless isothermal compressibility κ = κ T /κ (0) , which assures that relative density fluctuations within the DPD-fluid are matched to the real liquid. Here, κ T = −1/V (∂ V /∂ p) T is the isothermal compressibility of the real fluid computed from the change of the volume V under variation Table 1. Parameters used in the simulations. The surface tension has been obtained in a planar geometry, by integrating the lateral part of the pressure tensor in normal direction [41]. The viscosity η was obtained by a separate simulation following [42]. We use a velocity Verlet algorithm for the numerical integration of the MDPD equations using a time step t = 0.01.

Parameter
Symbol mu Fluid particle density ρ 6.00 Interaction range (attr.) of the pressure p using isothermal conditions and κ (0) is the isothermal compressibility of an ideal gas. An alternative way of matching the DPD-simulation to a real system is to calculate the system specific dimensionless numbers, e.g. the Reynolds number Re = ρvl/η and capillary number Ca = ηv/σ and to generate the corresponding situation in the DPD-simulation. In this work all results are given in model specific, i.e. DPD units, and are used for calculating the dimensionless numbers wherever needed.

A model for solid boundaries with adjustable adhesive properties
Modeling sl-interfaces in mesoscopic simulation approaches is particularly difficult: even on a molecular scale, the transition from a liquid to the solid is very sharp and abrupt. Naturally, the true behavior of a sl-interface cannot be represented, but it is possible to capture selected aspects phenomenologically. For instance, in LBM stick boundary conditions can be imposed by the so-called 'back reflection' prescription [31]. The solid boundary is represented as a mathematical plane, and particles crossing it are reflected by reversing their momentum. This principle has also been adopted in DPD approaches [44], but care must be exercised in order not to introduce numerical artefacts [45,46]. Complications also arise if, in order to represent wetting phenomena [41], also the boundaries are constructed from particle assemblies, analogous to MD [47]. A dense impenetrable particle configuration may lead to spurious density oscillations [48] or artificial slip [49]. An elegant and simple solution for flat interfaces has been presented by Merabia and Pagonabarraga only recently [50]. They replace the solid by yet another liquid phase with similar properties. For both phases, bounce-back conditions apply at the boundary plane so that they remain separate. Otherwise, the inter-particle interaction is left unchanged even across the boundary. This way, distortions in the fluid phase directly at the sl-interface are avoided.
With this system they investigated the spreading of droplets on solid boundaries, whereas they did not investigate the behavior under forced wetting as considered in the next section of this work. Similar in spirit, we have recently suggested a wall model that provides a smooth transition from the liquid to the solid with respect to the forces specified in equation (1) [51]. The main difference in our approach is that the solid is not defined by a mathematical description (e.g a predefined plane), but instead consists of an ensemble of particles, which are pinned to the initial position via a harmonic spring which largely enhances the possibility to create complicated channel geometries or textured surfaces. For this to be achieved, we allow for a thin diffuse interfaces between both phases to develop, and we shall now briefly summarize the most important aspects. The wall regions consist of an amorphous configuration of particles at the same density as in the liquid. The particles are attached to fixed sites by harmonic forces, while the interparticle interaction is a 'normal' liquid-liquid interaction given by equation (1) (A ss = A ll = −40.00). Thus, in our approach solid particles oscillate around the initial position over time leading to thermal roughness. In order to prevent particles of the liquid phase from diffusing into the wall indefinitely, an additional soft repulsive force normal to the wall surface is applied. Since the solid phase is made of DPD particles of the same density as the liquid phase a smooth transition from the liquid to the solid, see figure 1, is established with constant density of the dissipative particles in the liquid region. In contrast to a rigid wall, spurious density oscillations as well as temperature oscillations near the wall can essentially be reduced. Figure 1 illustrates the difference between a rigid wall (with a (100)-fcc arrangement of the wall particles) and the 'thermalized' wall presented here. In summary, the wall model used here provides a numerically simple scheme avoiding local distortions of liquid properties at the sl-interface.
Note that the only property imposed on the wall model is the affinity of the solid surface, governed by the parameter A sl , which is larger than A ll and A ss (A sl A ll , A sl A ss ). Qualitatively this is similar to the methodology in LBM, where an attractive external potential is used to tune the interfacial free energy in a van der Waals/mean field framework [52,53]. In a pure liquid, a value A sl > −40.0 would lead to a separation into immiscible phases (see also the work of Pagonabarraga and Frenkel on mixtures of non-ideal fluids [30] or Diaz-Herrera et al [54], who also achieve a broad miscibility gap in a binary mixture of Lennard-Jones fluids by varying the attractive part of the interaction). In the boundary model described above the choice of A sl thus assists in shaping the sl-interface accordingly. Interestingly, even at large A sl , when the repulsive (i.e. hydrophobic) interaction dominates, still no particle layering develops. This seems similar to the behavior of water close to hydrophobic surfaces, as studied by Grigera [55], where also a very smooth transition from the bulk density to zero directly at the interface is found. In that situation, it is believed that the depletion of liquid originates on account of entropic reasons, due to a reorganization of hydrogen bonds [56]. It might be of some interest to consider generic interface models as developed here in order to explore under which other conditions layering/depletion may arise (see also the discussion in [32] on layering at a free capillary surface).
There is a rather simple dependence of the static contact angle on the parameter A sl . The set-up for extracting the static contact angle is illustrated in figure 2. The meniscus was segmented into 20 slices, and the center of mass was determined for each slab separately in order to extract the shape of the meniscus and thus the static contact angle θ 0 . After achieving stationary conditions the center of mass was averaged 10 000 time steps for each slice separately in order to determine the static contact angle. Figure 3 shows the static contact angle as a function of the attraction parameter A sl . The wetting behavior of the wall can be tuned arbitrarily between hydrophilic and hydrophobic by changing the mutual attraction A sl . Furthermore, the resulting static contact angle is consistent with the expected value given by Young's law: θ 0 = arccos((σ sv − σ sl )/σ lv ), with σ sv , σ sl , σ lv being the interfacial tensions of the solid-vapor (sv), sl-and lv-interface. The surface tensions were measured independently in a planar geometry, by integrating the lateral part of the pressure tensor p in normal, i.e. z-direction, to the interface, σ i j = ( p zz − 1 2 ( p x x + p yy ) − φ) dz, (σ i j ∈ {σ sl , σ sv , σ lv }) as described in detail in [33]. φ as defined by Nijmeijer et al [47] accounts for the contribution of an external potential. In our case this is the soft repulsive force preventing particles of the liquid phase from diffusing into the wall over time. The contribution to the surface tension is less than 5% for all cases of A sl . Within the numerical accuracy the two data sets agree implying that the 'measured' contact angle θ 0 is in agreement with the thermodynamic prediction.

Behaviour of dynamic contact angles under forced flow in a thin slit
In the following, we demonstrate that the DPD scheme provides a global solution for the moving CL-problem, which can be interpreted in terms of effective continuum theories [4]. For investigating the dynamic contact angle behavior, we consider a situation of driven capillary flow in a thin slit, similar to the experiments of Hoffman [57]. This is termed here as forced wetting, as the liquid front is made to move with constant velocity across the solid without changing shape, the velocity is determined by the piston. This facilitates experimental but also the theoretical or numerical analysis, since one is essentially dealing with a stationary situation, in contrast to many natural wetting phenomena such as droplet spreading, where the time varying shape of the capillary meniscus must be taken into account. For the channel system, we generate a stationary situation as illustrated in figure 4(a). The liquid is placed between two plates that are moved with a relative velocity v wall , while the plug is kept fixed by a piston. In the reference frame of the liquid the walls move with velocity −v wall relative to a piston that keeps the liquid in place. v wall determines the capillary number Ca = ηv wall /σ . We may thus extract the shape of the meniscus by the procedure described above for the static case, see figure 2, and obtain a relation θ d (Ca). Note that this procedure will in general yield an apparent dynamic contact angle θ d , as we do not measure the local inclination of the meniscus with the boundary and the data points at the interface are not considered for the   [4]. The continuum model predicts pronounced convection (rolling) next to the CL in agreement with experimental observations [8] and the numerical result obtained from the MDPD-simulation. spherical fit. This set-up is also suitable for extracting the flow field after sufficient averaging. Figure 4(b) illustrates the flow field resulting from the MDPD-simulation in the fixed reference frame, i.e. the frame where the wall has a velocity −v wall . Along the solid surface we observe significant slip over a distance of ∼1.5 r c .
The slip occurs in combination with pronounced caterpillar motion of the fluid, which has been demonstrated in an experiment by Dussan [8]. Figure 4(c) shows the streamlines as obtained from solving the Navier-Stokes equations for an infinite fluid wedge, with constant velocity at the free surface and a no-slip condition at the sl-interface: with the origin of the reference frame being the position of the CL. The coefficients C A , D A , E A and F A are given in [4]. Equation (2) leads to a logarithmic singularity of the viscous stress as r → 0. In contrast figure 4(b) shows that the DPD-method provides slip as a mechanism for removing the stress singularity. According to Cox [4] the singularity can be avoided by identifying at least two distinct regions: (i) a region close to the CL, where slip predominates the flow and (ii) a no-slip region Figure 5. Comparison of the constitutive law for dynamic contact angles according to Cox [4] (equation (4)) and the corresponding characteristics obtained by the MDPD-simulations. Vertical axis in rad 3 .
(and (iii) sometimes an intermediate region is also needed). The asymptotic matching of the different regions provides a solution for the constitutive law of the dynamic contact angle to the leading order of Ca [4]: with the ratio = L micro /L macro independent of Ca and g(θ) = θ 0 dx. L micro is a microscopic length scale over which a slip boundary condition is assumed to dominate the flow field. L macro signifies a macroscopic length. By studying forced wetting of silicone oils in thin glass tubes, Hoffman [57] found empirically a relationship that almost quantitatively agrees with equation (3) and that complies with the Taylor expansion of equation (3) for small angles (θ d < 100 • ) [9]: In the case of partial wetting, see figure 5, the data sets can well be fitted to equation (4). If we take 9 · ln( −1 ) as a compound fit parameter we find 9 · ln( −1 ) = 20.8 ± 0.3.
In contrast, the case of complete wetting has a more complex characteristics, see figure 6. Although for Ca > 0.04 the slope attained is comparable to the case of partial wetting, a significant deviation from the linear behavior occurs for small Ca. Eggers [9] pointed out that in the case of θ 0 = 0 • the microscopic length L micro actually depends on Ca. Eggers showed that the solution of the lubrication equation yields an additional dependency such that L micro = 0.54λCa −1/3 , if slip is assumed as the mechanism responsible for removing the stress singularity. In this context λ is the Navier slip length such that, v − v wall = λ · (∂v/∂ x), if v is the fluid velocity at the wall and v wall is the velocity of the moving boundary. This leads to the constitutive law for θ 0 = 0 • : Figure 6. Comparison of the constitutive law for dynamic contact angles according to Eggers [9] (equation (5)) and the corresponding characteristics obtained by the MDPD-simulations. For reference the model of Cox [4] is depicted. Vertical axis in rad 3 .
where α = L macro /(0.54λ). The additional weak logarithmic dependency only becomes relevant for small Capillary numbers. We compare the constitutive law, i.e. equation (5), with the result obtained by the MDPD-method. Figure 6 shows that within the numerical accuracy the MDPDmethod reproduces the behavior given by equation (5) for θ 0 = 0 • . The numerical value obtained by the fit is α = (16.7 ± 0.3).
We should stress that all that really is required for the exponent of 1/3 to apply is the absence of an extended precursor foot, present only for strongly wetting systems with a large excess interfacial free energy. A numerical value of 1/3 is not strictly related to the appearance of slip, but instead may originate from some other mechanism, that removes the stress singularity. In this context the diffuse nature of the sl-interface (see also for use of diffuse Cahn-Hillard interfaces [58]) has been pointed out as one possible mechanism to avoid a stress singularity. Furthermore, it must be expected that along the slip region, especially at the foremost liquid front, the interface needs some time to form, which yet would be a refinement of a diffuse-interface mechanism. As a matter of fact, a possible interface relaxation plays a major role in an elaborate continuum theory put forward by Shikhmurzaev [10,59]. Exploring the present boundary model along these lines further could be of considerable interest.

Dynamic capillary impregnation
Finally, we address the wetting behavior of a transient problem, the capillary impregnation of a slit pore. The asymptotic limit (t → ∞) of this problem can effectively be described by a onedimensional equation, the so-called Washburn equation [60], which describes the penetration depth z over time t, i.e. z(t) ∝ √ t. This behavior has been confirmed experimentally over a wide range of time and length scales [60]- [62] down to the nano-scale. Also MD simulations have shown that the characteristics can be described in the framework of a one-dimensional continuum description [63,64] even for short time and length scales. We pursue this path by studying a system as shown in figure 7. We consider the case where θ 0 = 0 • , i.e. A sl = −40.0. The asymptotic limit for t → ∞ cannot directly be applied since inertial effects and the presence of the dynamic contact angle, i.e. θ d = θ 0 , alter the dynamics of the wetting process.
Before investigating the dynamics of the impregnation, the analytical solution of this problem is derived by extending the Washburn law [60] in order to account for the dynamic contact angle contribution and inertia. The dynamics of the capillary impregnation is governed by balancing capillary forces F cap , viscous forces F η and inertial forces due to the momentum change d p(t)/dt generated by the mass M(t) which is moved in the z-direction during the process of impregnation: Assuming that the energy is dissipated in the capillary equation (6) gives: 1 bh where b denotes the slit width, see figure 7, h denotes the width in the y-direction (12r c ) and z is the penetration depth into the pore. Equation (7) corresponds to the Washburn law [60] (adapted to a slit) if inertia is set to zero, i.e. d p(t)/dt = 0, and the dynamic contact angle is replaced by the static contact angle, i.e. θ d = θ 0 . This corresponds to the special case in the asymptotic limit t → ∞. The momentum change in z-direction can directly be expressed in terms of the momentum of each particle, d p(t)/dt = d( m iżi )/dt. While the momentum change of the mass M slit in the slit can be calculated analytically by d p slit /dt = d(M slitż (t))/dt = d(ρbhz(t)ż(t))/dt, the momentum change in the reservoir is more difficult to determine (due to the complex flow field) and therefore was measured during the course of the simulation and used as input to equation (7). Also the dynamic contact angle θ d was monitored during the simulation. Figure 8 shows the observed dynamic contact angle over time during the MDPD-simulation. During the timescale studied here the resulting dynamic contact angle clearly exceeds the t inertial t Figure 8. Dynamic contact angle behavior θ d obtained by the MDPD-simulation versus time during the capillary filling process. Vertical axis in degrees. Figure 9. Dynamics of the penetration into the pore over time z(t) obtained by the MDPD-simulation compared to equation (7). The dynamic contact angle is measured online during the simulation and inserted into equation (7) for the numerical integration.
static contact angle (θ 0 = 0 • ). Over time the dynamic contact angle slowly decreases as the penetration velocity into the slit decreases, see figure 9. Interestingly, the inertial timescale [64] τ = ρb 3 /σ ≈ 80 roughly corresponds to the timescale needed for the meniscus to relax from its initial configuration (starting conditions corresponds to the static contact angle θ 0 = 0 • ). Since inertia has a considerable influence on the dynamics of the impregnation, the dynamic contact angle is directly extracted from the simulation and inserted into equation (7) in order to check if the MDPD-simulation can be consistently described by equation (7).
Integrating equation (7) numerically with the data set of figure 8 results in excellent quantitative agreement with the MDPD-simulation, see figure 9. As a reference in figure 9 the Washburn law [60] is shown, which clearly overestimates the penetration dynamics as inertia and the reduction of the capillary pressure due to the dynamic contact angle are missing.

Summary and conclusions
We presented a study on capillary wetting under static as well as dynamic conditions in the framework of DPD. The measured static contact angle is in agreement with Young's law. The characteristics of the dynamic contact angle was investigated under stationary conditions and compared to the hydrodynamic models of Cox [4] and Eggers [9], respectively. The dynamic contact angle is in agreement with the theoretical prediction of Eggers [9] for a static contact angle θ 0 = 0 • with an additional logarithmic dependency on Ca compared to the model of Cox [4]. We find that for a static contact angle θ 0 > 0 • the microscopic length L micro is independent of the capillary number in agreement with theory [1,4,9]. As a last step a transient problem, the wetting into a slit was studied. The dynamics can be described by a simple approach balancing capillary and viscous forces under consideration of inertia. This analytical solution was compared to the result of the MDPD-simulation resulting in good quantitative agreement. In summary, we have shown that the MDPD-method can successfully be applied for studying dynamic stationary as well as transient capillary wetting phenomena. The MDPDapproach therefore has the potential to become a valuable explorative tool as an alternative or supplement to CFD or analytical approaches.