1 Introduction

A liquid droplet in contact with a soft substrate, such as elastomers, gels and rubbers, can deform the substrate due to the localized force acted upon by the liquid on the solid at the contact line.^1,2 This morphological change in soft solids has a significant impact on the shape of the droplet,³ as well as the dynamical process when the droplet is in motion.^4–8 Apart from the scenario of droplets interacting with a soft solid, studying how soft solids respond to external localized stresses is important to the understanding of industrial and biological systems such as soft adhesion,^9,10 and the functions of cells and tissues.^11,12

The deformation of soft solids below the contact line appears as a ridge-like shape, which has been measured, for example, using X-ray microscopy.¹ The equilibrium shape of the ridge is critically determined by the elastic properties of the soft solid,^2,13–15 as well as the forces acting on it. An apparent normal force pulling the solid upwards by the liquid at the contact line is the tension with a magnitude equal to γ sin θ,¹⁶ here γ is the liquid–air surface tension and θ is the liquid equilibrium contact angle. A length scale named the elastocapillary length γ/G can be constructed, with the shear modulus G serving as a resistance to deformation.¹⁵ Earlier computations of the equilibrium shape of the wetting ridge date back to the studies by Lester¹⁷ and Rusanov,¹⁸ in which a droplet is in contact with a soft semi-infinitely thick substrate. Later, more investigations were done for a soft solid layer of finite thickness H attached to a rigid substrate.^19,20 All these studies point to a conclusion that both length scales γ/G and H play a crucial role in solid deformation.

Remarkably, recent studies^15,21–24 have demonstrated that when including the contribution of solid surface tension γ_s to the elastic model, the agreement with the experimentally measured ridge profiles turns out to be promising. The solid surface tension suppresses the deformation and regularizes the singularity at the contact line.²¹ The solid angle formed at the contact line follows the balance between the surface forces, or named Neumann's relation.^3,25 Neumann's relation has also been validated by studies using nonlinear elastic constitutive models²⁶ and experimentally by measurements using confocal microscopy.⁹ When γ/G is very small, i.e. rigid solids, one expects the liquid contact angle at the contact line to follow the Young's relation. By including a microscopic length scale l_m characterizing the width of the liquid–air interface,²⁷ it is shown that the transition from Neumann's relation to Young's relation occurs when γ/G ∼ l_m.^28,29 Young's relation is recovered for γ/G ≪ l_m. The crossover from the Neumann's relation to Young's law has also been demonstrated in studies using the molecular dynamics simulations.³⁰

Apart from the elastic properties, many soft solids like tissues and gels also demonstrate viscoelastic features. When a droplet is placed in contact with such smooth soft substrates, the flat surface of the solid around the contact line starts to deform and gradually form a wetting ridge. Both experimental^1,6 and theoretical⁶ studies find that a sharp tip has developed even at the early stage of ridge formation. Another study³¹ investigated the relaxation of a deformed viscoelastic half-space subsequent to the removal of the point load. It was shown that solid surface tension reduces the local elastic deformation and stress field near the load point.

As we mentioned above, the dynamic response of a viscoelastic solid to a time-dependent localized force has been demonstrated by some previous studies. Nevertheless, how different regions of a viscoelastic solid respond in time during the growth and the decay of the ridge has not been fully investigated and described in the literature, in particular the role played by the solid surface tension, which is crucial for the understanding of phenomena such as stick-slip motion of contact lines;^7,32 droplet evaporation^5,33,34 and soft electrowetting.³⁵ In this study, we examine both the growth and the decay of a ridge and show how the solid interface and the characteristic features of the ridge change with time. Importantly, we focus on the tip region of the ridge and show that the solid material in that region relaxes in a different dynamical way from the remaining part of the ridge. To describe the viscoelastic properties of the soft material, the standard solid model will be implemented. The solutions are computed using the approach of the boundary element method (BEM), which has been used widely to study droplet^36–38 and solid deformation.^39,40 The BEM can be extended to study problems involving viscoelastic deformation coupled with fluid flow such as in wet adhesion.

2 Formulation

We consider that a droplet is placed in contact with a layer of isotopic viscoelastic solid with thickness H which is attached to a rigid solid. We denote the displacement at any position x = x + yŷ + zẑ of the undeformed viscoelastic solid and time t as u(x,t) = u_x(x,t) + u_y(x,t)ŷ + u_z(x,t)ẑ. In index notations, we write the component as u_i where i = x, y or z. In the following, the subscripts i, j and k take either one of the components in rectangular coordinates. The repeated indices conform to the Einstein summation convention, i.e., they are summed over the three components. For viscoelastic materials made of polymers such as elastomers and gels, we assume that the bulk modulus K can be considered as a constant, independent of time. We further assume that the strain is small, namely the strain tensor ε_ij(x,t) = (∂u_i/∂x_j + ∂u_j/∂x_i)/2 ≪ 1. Hence, the stress tensor σ_ij(x,t) and the strain tensor follow the linear constitutive relation given as (1) where the deviatoric part of the strain tensor , ė_ij ≡ ∂e_ij/∂ and is the shear relaxation function which depends on the viscoelastic model characterized by the parameters G, G₁,… and η, η₁,…. We assume the dynamics of the solid deformation is slow enough so that the net force on a soft solid element vanishes which gives the governing equation (2) when body forces such as gravity are neglected.

The range of solid deformation depends on the elastocapillary length γ/G and the thickness of the soft layer H. In this study, we assume that both the elastocapillary length and the thickness of the soft layer are much smaller than the contact radius R of the droplet, i.e., γ/G ≪ R and H ≪ R, meaning that the solid deformation is localized around the edge of the droplet. In this limit, it is justified to assume that the deformation of a solid behaves as the case of plane strain, for which u_x and u_y are a function of x and y only, and u_z is a constant, and so the governing eqn (1) can be simplified due to the conditions of the vanishing of the following quantities, namely (3)

In the plane strain situation, the domain of the soft solid layer that we consider is a rectangle with thickness H and length L ≫ H. Fig. 1 shows the schematic of the setup. The soft layer is attached to a rigid solid, hence at the bottom boundary (y = −H) of the domain where the soft and the rigid solids are in contact, the displacement u_i(x,y = −H,t) = 0. (4) At the two other ends of the domain, we also assume the displacements are zero, i.e. u_i(x = −L/2,y,t) = u_i(x = L/2,y,t) = 0. (5) The boundary conditions at the solid/fluid interfaces require more discussions. As it has been shown that the solid surface tension plays an important role in the deformation of the wetting ridge, we include this specific property in our model. For soft solids, we have to distinguish between surface energy and surface stress due to the Shuttleworth effect.^41–43 In this study, for the simplicity of modelling, we neglect this effect and hence the tension force at the solid interface is independent of the stretching of the solid. We also assume that the surface tensions are the same for the solid/liquid and the solid/air interfaces, hence the tension is described by a constant γ_s along the whole solid interface, and the liquid contact angle θ = 90°. We now consider the forces acting on an element of the interface. We first define the traction in the soft solid as f_i ≡ σ_ijn_j, where n_j is the normal vector of a surface pointing outward from an element. We assume that the external force acting on the soft solid only has the out-of-plane component (i.e., the y-component). Suppose the external force per unit area acting on the interface is given by f ^ext_i(x,t), for small interface slope due to small deformation ∂u_y/∂x ≪ 1, we have the force balance on the interface as f_x = −f ^ext_x = 0 (6) and (7) where the normal force balance eqn (7) includes the contribution of the solid surface tension γ_s.^15,22,23 Here we have adopted the approximation for the curvature of the interface, i.e., κ = ∂²u_y/∂x², when the slope is small.

To model the localized external traction, we use the Gaussian function to describe the spatial distribution of the normal traction, i.e., , where _m is the microscopic length such that _m ≪ H. In the limit that _m → 0, the traction f ^ext_y(x,t) = Γ(t)δ(x), where δ(x) is the Dirac delta function, which is also used in some other studies.⁶ The strength of the external force Γ(t) is a function of time. In this study, we will use two different functional forms, respectively, to investigate the growth and the decay of the ridge, namely case 1: Γ(t) = γ(t) and case 2: Γ(t) = γ[1 − (t)] = γ − γ(t). Here (t) is the Heaviside step function. As the governing equations are linear, the solution for case 2 (decay case) will be obtained by combining the steady solution (at t → ∞) of case 1 and the solution obtained when using Γ(t) = −γ(t).

There are different types of viscoelastic models which have been used to describe the behavior of different viscoelastic materials. The simplest models are the Kelvin–Voigt model or the Maxwell model which consists of a spring and a viscous damper connected in parallel or in series respectively. In this study, we assume the response of the soft material follows the standard solid (Zener) model which consists of a spring (characterized by G) connected in parallel with a Maxwell element, i.e., a spring (characterized by G₁) connected in series with a viscous damper (characterized by η), which is relatively simple but can describe both stress relaxation and creep. The standard solid model has been used to describe the viscoelastic behavior of elastomers. We note that our approach can be used for other linear viscoelastic models. Please refer to the Appendix for the details of the description of the viscoelastic model. Here we point out that for our viscoelastic model, the initial instantaneous shear modulus is G + G₁ and the static (long-term) shear modulus is G. Since many soft solids can be considered as incompressible, we will take the incompressible limit, i.e. K ≫ (G + G₁).

2.1 The boundary element method (BEM)

To solve the governing eqn (1) and (2) with the plane strain conditions (eqn (3)) and the boundary conditions given by eqn (4)–(7), we implement the approach of the boundary element method (BEM) in the time domain.^39,44–46 We denote the closed contour of the boundary of the domain as D, the position vector of a point on the contour as and as the differential length on the contour. The governing eqn (1)–(3) are formulated as boundary integral equations (see the Appendix for the details of derivation) which read (8) where (9)

The expression of the tensor components P_ij(x,,t) and U_ij(x,,t) depends on the specific viscoelastic model as for the shear relaxation function . In the Appendix, we derive the expressions for the standard solid (Zener) model. We note that P_ij and U_ij are singular at x = ; the integral over the singular point is thus computed analytically by expanding the tensor components in series around x = .

To solve the boundary integral equations numerically, we implement a time marching scheme^39,44–46 and discretize the boundary of the domain D as follows: we have evenly distributed marker points for the two side boundaries at x = −L/2 and x = L/2, and the bottom boundary at y = −H. The number of marker points, respectively, equals to 50, 50 and 200. For the top boundary at y = 0, we use 400 marker points and distribute them with a separation increasing exponentially from the center. The distance d_m between two marker points at the center is as small as d_m/H = 10⁻⁶ so that we can resolve the fine structures around the contact line where the localized force is imposed. We take a constant time step δt. At any time t = nδt, the boundary integral eqn (8) is discretized in the following way (using trapezoidal rule for the time integral): where _ij(x,,t) = ∂U_ij(x,,t)/∂t, Ṗ_ij(x,,t) = ∂P_ij(x,,t)/∂t and S_i(x) is the term corresponding to the history of responses, which can be written as (10)

Regarding the rescalings, unless specified differently, we will divide all the lengths by the thickness of the soft layer H, the stresses by γ_s/H and the time by η/G; we end up with the following independent dimensionless control parameters (with a tide): (11)

3 Results

We consider that the width of the soft layer L is much larger than its thickness H, and the microscopic length scale _m is much smaller than the thickness of the soft layer. We take = 50 and (e.g. consider a liquid/air interface of nanometric thickness and a 10 micron thick soft solid layer). We present the results in three main parts. In Section 3.1, we focus on the evolution of the global features of the ridge. In Section 3.2.1, we discuss at how the static solid interface profile depends on the parameters. In Sections 3.2.2 and 3.2.3, the role of the solid surface tension and the relaxation at the tip region is explored in detail.

3.1 Shape evolution of the wetting ridge

3.1.1 Profiles of the wetting ridge during the growth and the decay

We look at how a wetting ridge grows when and how the ridge decays when , with = 0.2. As our viscoelastic model is linear, which assumes small solid interfacial slopes (small strains), we thus pick the value of to be significantly smaller than 1. Note that the solid interfacial slope around the contact line increases with based on the Neumann's relation. The time series of is plotted in Fig. 2(a) and (b). We take = 10⁻¹ and ₁ = 10². Fig. 2(c) and (d) show the out-of-plane displacement ũ_yo ≡ ũ_y(,ỹ = 0,) at different times. Fig. 2(e) and (f) show the corresponding in-plane displacement ũ_xo ≡ ũ_x(,ỹ = 0,). At = 0⁺, there is an instantaneous elastic deformation with shear modulus equal to + ₁. As ₁ is large, the initial solid surface deflection is tiny. The ridge shape, described by ũ_yo, is symmetric along the central axis ( = 0) while the in-plane displacement ũ_xo is anti-symmetric. A significant difference between the growth and the decay of the ridge is that, during the growth, a sharp tip (around = 0) is maintained from the early stage of the growth. For the decay, when the localized external stress vanishes (at = 0), the sharp tip quickly relaxes to round shapes with a much smaller curvature as can be seen in Fig. 2(d). Note that the profiles at = 0⁻ in Fig. 2(d) and (f) are the steady profiles before the external stress vanishes.

3.1.2 Evolution of the characteristic quantities

There are several features of the ridge that we can describe here. First, the maximum out-of-plane displacement ũ_y,max ≡ ũ_yo( = 0,) of the solid surface or called the ridge height is at the center = 0. Second, the out-of-plane displacement decreases from the center and a minimum value ũ_y,min ≡ ũ_yo( = _min,) is achieved at = _min (and = −_min), where a clear dimple is formed with a negative out-of-plane displacement. Third, for the in-plane displacement ũ_xo, we see from Fig. 2(e) and (f) that around the center = 0, the solid material is compressed (with a negative slope shown in the plot, i.e., ∂ũ_xo/∂ < 0, and there is a maximum ũ_x,max ≡ ũ_xo( = _max,) at the negative side of . These characteristic quantities are also indicated in the schematic of the setup in Fig. 1.

Now we look at the evolution of the characteristic quantities. In Fig. 3(a) and (b), we plot ũ_y,max, |ũ_y,min| and ũ_x,max as a function of the rescaled time , respectively, during the growth and the decay of the ridge. All these quantities increase/decrease monotonically with time and saturate to a steady value during the growth/decay of the ridge. Next, we look at the time series of the magnitudes of the positions |_min| and |_max| as plotted in Fig. 3(c) and (d). Both |_min| and |_max| equal to unity at = 0⁺ (instantaneous elastic response). This means that for a large ₁, these two positions are equal to the thickness of the soft substrate. They increase with time during the growth. However, during the decay of the ridge, instead of decreasing, the positions still increase with time. This feature can also be observed from the profiles shown in Fig. 2(d) and (f).

To characterize the dynamics of the whole wetting ridge, we here focus on the evolution of the ridge height ũ_y,max. For both the growth and the decay, we find that, a rescaled ridge height relaxes exponentially to the steady states at late times. The rescaled ridge heights are defined as for the growth and for the decay. In Fig. 3(e) and (f), the rescaled ridge heights are respectively plotted as a function of for different values of . We find that both ũ_g() and ũ_d() relax the same way, i.e., ũ_g() or ũ_d(), scales as ∼exp(−1.15) at late times. The late time ridge evolution possesses the same time scale as the relaxation of the viscoelastic material in our standard solid model, i.e., η/G.

3.2 The role of the solid surface tension

3.2.1 Static profiles: from the soft regime to the stiff regime

Before we look at the dynamic response of the soft solid at the tip region of the ridge, it is insightful to see how the static shape (purely elastic) of solid deformation at the tip region varies when is changed. Note that the static shape is the same as the shape at → ∞ in the growth case. From the scaling with , we will also deduce the dependence on the shear modulus G and the soft layer thickness H, which are dimensional quantities of physical and practical importance. Before we focus on the tip region, let us first look at the ridge height ũ_y,max( → ∞), which is plotted in Fig. 4(a) as a function of . As one might expect, the ridge height increases with decreasing . Moreover, when ≫ 1 or in other words, when the soft layer thickness is much larger than the elastocapillary length γ_s/G, we find that ũ_y,max( → ∞) ∼ ⁻¹. This scaling implies that the dimensional ridge height u_y,max( → ∞) scales linearly with the elastocapillary length and is independent of the soft layer thickness. The confinement (thickness) of the soft layer plays no role in this regime. When becomes smaller (≲10), we see that ũ_y,max( → ∞) scales as ^α locally with −1 < α < 0. Hence u_y,max( → ∞) increases with decreasing G and increasing H.

In the following, we focus on the tip region. In Fig. 4(b), we plot the static out-of-plane displacement ũ_yo(, → ∞) in the region of −0.01 < < 0.01, for a wide range of = 10⁻²–10⁴. Note that the profiles are shifted by the ridge height ũ_y,max( → ∞) so that all the tips are at the same position. For ≤ 10, we see the solid interfaces form the same solid angle θ_s = 168.6°. The solid angle can be understood by the force balance between the solid surface tensions and the external liquid–air surface tension, namely γ = 2γ_s cos(θ_s/2). For > 10, the interface profiles have a solid angle getting closer to 180° as increases. This shows how the solid deformation at the tip region varies from a soft regime to a stiff regime which are consistent with Neumann's relation and Young's relation, respectively, in the two different limits.

To show clearly how the slope of the solid interface varies with the distance from the center at different length scales, we plot the interfacial slope as a function of (in log scale) in Fig. 4(c). For the stiffer case we have computed here, e.g. = 10⁴, the magnitude of the slope is small. The whole solid interface maintains close to a flat surface even though being exerted by a localized external force. As decreases, the magnitude of the slope increases. There appear overlapping regions of the curves for different values. The size of the overlapping region increases with decreasing . The overlapping of the curves means that the solid deformation in that region is independent of , and hence the solid surface tension dominates over the elastic stress. Moreover, we also see that there is a region of a constant slope −0.1 for < 10. Note that the slope −0.1 means that the solid angle θ_s = 168.6°. Indeed the constant-slope region is what we have plotted in Fig. 4(b). At distances from the center smaller than , the distribution of external stress plays a role, the slope of the interface is not −0.1, although the solid surface tension dominates over the elastic stress (for ≤ 10). It is important to point out that plays a role for the existence of a region where Neumann's relation is valid. When taking , there always exists a region close enough to = 0 where the Neumann's relation is valid for any values. For example, in the inset of Fig. 4(d), we plot as a function of for both and , and a fixed value of = 10³. We find that when approaching the contact line ( = 0) from a large distance, the curves of for and overlap until ≈ 3 × 10⁻⁴ where starts to play a role. In the case of , the magnitude of the interfacial slope keeps increasing until ≈ 3 × 10⁻⁷ for which and a constant-slope region is observed around ≈ 3 × 10⁻⁷. In reality, the thickness of the interface (i.e. _m) is in the nanometric scales, that means when is larger than a certain value, Neumann's relation cannot be observed at any distance from the contact line. Instead, the solid interface maintains as a flat surface as we have shown in Fig. 4(b) and (c).

How does the size of the constant-slope region depend on ? To address this issue, we introduce the boundary of the constant-slope region as follows. We define _γ as the largest value of || such that . Thus when , we have , which defines the range of the constant-slope region. Since the constant-slope-region disappears for large when is not small enough, we here use in our computations. In Fig. 4(d), we plot _γ as a function of . As we have demonstrated already in Fig. 4(c), we find that _γ increases with decreasing . Similar to the scaling of the ridge height, when ≫ 1, the size of the constant-slope region _γ ∼ ⁻¹, which implies that the dimensional quantity x_γ scales linearly with the elastocapillary length γ_s/G and is independent of the soft layer thickness. When ≲ 1, we see that _γ scales as ^β locally with −1 < β < 0. Hence the dimensional size x_γ increases with decreasing G and increasing H. We end this section with an example to estimate the dimensional size of the constant-slope region. Taking the shear modulus G = 1 kPa, the thickness of the soft layer H = 10 μm, the solid surface tension γ_s = 0.01 N m⁻¹ and γ = 0.002 N m⁻¹, we get = 1. From Fig. 4(d), we obtain _γ ≈ 10⁻³, i.e., x_γ ≈ 10 nm. No constant-slope-region can be observed when x_γ ≲ 3_m.

3.2.2 Dynamic response at the tip region during the growth

We first look at the growth case. We take = 10⁻¹ and ₁ = 10². In Fig. 5(a), the interfacial slope as a function of (in log scale) is plotted for different times. At = 0⁺, the solid responds by an instantaneous elastic deformation. As time goes on, we observe that there is a region around the tip where the interfacial slope is independent of the time after. The size of this region grows with time. For example, this region has a range 0 < ≲ 10⁻⁴ at = 0⁺ and 0 < ≲ 10⁻³ at = 0.02. To understand this result, we note that the shear modulus relaxes from + ₁ to for our viscoelastic model. We have also shown in Fig. 4(c) and (d) that there exists a surface-tension-dominated region with a size increasing with decreasing . Hence, once a solid element relaxes to within the surface-tension-dominated region, the solid element is in an equilibrium state. The time taken to relax to an equilibrium state thus depends on the distance from the center. The solid relaxes to an equilibrium state in a shorter time for materials closer to the center as we have shown in Fig. 5(a). In contrast, for stiffer solid materials for which no surface-tension-dominated region exists, all the solid elements relax in the same time scale η/G as shown in one example in Fig. 5(b) in which = 10³ and ₁ = 10⁴.

How does the relaxation at the tip region depend on material parameters? To address this question, we pick a position = _o = 10⁻² and study how the solid interface slope at this position varies with time for different values of , ₁ and . We find that the slope relaxation is independent of when plotting as a function of the rescaled time / as shown in Fig. 5(c). In other words, the relaxation at the tip region is much shorter than the long-term relaxation of the ridge. Note that / = tγ_s/ηH. On the other hand, the short-term shear modulus ₁ does play a role as shown in Fig. 5(d). It is interesting to see how the surface tensions determine the tip relaxation. In Fig. 5(e), we plot recalled by /2 as a function of / for different values of ; all curves collapse on each other. Hence we find that scales linearly with . From the results we have shown, the relaxation of the slope can be written as , with the functional form of F undetermined. The physical interpretation is as follows: the ratio γ/γ_s determines the equilibrium slope (solid angle). The relaxation time scale is determined by the solid surface tension γ_s and the short-term relaxation time of the viscoelastic material. The relaxation is characterized by a spreading velocity γ_s/η.

3.2.3 Dynamic response at the tip region during the decay

Once the external force disappears, the sharp tip relaxes immediately. How does the tip region relax during the decay of the ridge? In Fig. 6(a), we show the slope as a function of (in log scale) for different times during the decay. As there is an instantaneous elastic response at , the constant-slope-region may disappear as early as at = 0⁺ as shown in Fig. 6(a) for which ₁ = 10². Since the static shape at the tip region is determined by the surface tensions (in the soft regime), one may expect the tip region relaxes differently from the rest part of the ridge during the decay, as we have observed for the growth case. To address this question, we pick two positions, i.e. = _o = 10⁻² and = _o = 1, and study how the solid interface slope at these two positions varies with time for two different values of , namely = 10⁻¹ and = 10⁻². Note that from Fig. 4(d), _o = 10⁻² ≲ _γ and _o = 1 ≫ _γ. In Fig. 6(b), we see that at _o = 10⁻², the evolution of the interfacial slope is the same for = 10⁻¹ and = 10⁻². However, at _o = 1, the slope evolves differently as we can see in Fig. 6(c), and the time scale depends on η/G, which is much larger than the relaxation time at the tip region. On the other hand, the slope scales with for both = _o = 10⁻² and = _o = 1 as we have shown in Fig. 6(d) and (e), in which the y-axis is rescaled by /2. This is because the out-of-plane displacement depends linearly with γ for the whole ridge, and hence the ratio γ/γ_s determines the interfacial slope.

4 Discussion and conclusions

In summary, we study the dynamic response of an incompressible viscoelastic solid being acted upon by a localized surface force. The responses of the global features of the wetting ridge, e.g., the growth and the decay of the ridge height, depend on the rheological properties of the solid material. For the standard solid model we implement in this study, the ridge height relaxes exponentially at late time with a time scale given by η/G. Some other studies use the Chasset-Thirion model to describe the viscoelastic behavior of silicone gels;⁶ the relaxation of the ridge is found to follow a power law with time. The dynamic around the tip, on the other hand, is different from the remaining part of the soft solid during both the growth and the decay. For the growth, the equilibrium state of a solid element is achieved when it is within the surface-tension-dominated region, which has a size increasing with time due to viscoelastic properties of the solid. Hence an equilibrium can be achieved within a time scale much shorter than the relaxation time of the shear modulus of the viscoelastic material (η/G). Similarly for the decay, the tip region relaxes with time scales that depend significantly on the solid surface tension, namely the relaxation has a characteristic spreading velocity γ_s/η. Although the results presented in this study are obtained using the standard solid model, the existence of the surface-tension-dominated region at the ridge tip is expected to be universal. Hence the features of the dynamics at the tip region shown in this study can be compared with experimental measurements given that the viscoelastic properties of the solid material are determined.

How does the dynamic response of a soft solid play a role when a droplet is in motion? Current models^4,6,47 assume dissipation being dominated by the viscoelastic effects in the soft solid and neglect flow in the fluids, and consider a constant contact line velocity. On the other hand, experimental studies observe the stick-slip motion of the contact line.⁷ The interplay between the fluid flow and the time-dependent morphological changes of the wetting ridge is expected to be crucial for the dynamics of the contact line and the droplet motion, which remains a challenging topic to be explored. Our findings show the features of the local morphological changes at the ridge tip during both growth and relaxation, which provide an important fundamental knowledge for investigating a contact line moving over a soft surface.

Conflicts of interest

There are no conflicts to declare.

Appendix A1. Derivation of the boundary integral equation in time domain for linear viscoelastic models (plane strain cases)

We here decompose the stress σ_ij(x,t) and the strain ε_ij(x,t) tensors into the hydrostatic and deviatoric parts. Namely, we write (12) and (13) with the strain tensor defined as (14)

For isotropic linear viscoelastic materials, the constitutive relations can be written in the generalized differential forms which read (15) and (16) where p_l, q_m, a_l and b_m are coefficients that depend on the specific viscoelastic model. Applying Laplace transform to the constitutive relations (15) and (16), with Laplace transform of a function f(t) defined as , the constitutive relation can be written as⁴⁴ Σ_ij(x,s) = (s)E_kk(x,s)δ_ij + 2Ĝ(s)E_ij(x,s), (17) where , E_kk(x,s) ≡ [ε_kk(x,t)] and E_ij(x,s) ≡ [e_ij(x,t)]. The expressions in Laplace space of the bulk modulus and the shear modulus . In Appendix A2 we show how to obtain them for the standard solid (Zener) model. Eqn (17) has the same form as the linear elastic constitutive relation. Hence by following the same approach of deriving the boundary integral equation for a purely linear elastic material in the quasi-static state (∂σ_ij/∂x_j = 0), we can obtain the boundary integral equations for linear viscoelastic models in the Laplace space. Here we consider plane strain cases, for which the boundary integral equation in Laplace space is given as⁴⁵ (18) where (19) û_j(x,s) ≡ [u_j(x,t)], _j(x,s) ≡ [σ_ij(x,t)n_j(x)], with n being the normal vector pointing outward the soft solid domain, and the line integrals on the right hand side are carried out along the closed contour D of the boundary of the domain with dl() denoting the differential length of D and = + ȳŷ being the position vector of a point on the contour. The functions Û_ij(x,,s) and _ij(x,,s) are given by^45,46 (20) (21) where r = [(x − )² + (y − ȳ)²]^1/2, , and . The functions Ĵ₁(s), Ĵ₂(s), Ĵ₃(s) and Ĵ₄(s) are given as (22)

To obtain the boundary integral equation in the time domain,^39,44–46 we perform an inverse Laplace transform ⁻¹ to eqn (18) and obtain (23) where û_ij(x,,) ≡ ⁻¹[Û_ij(x,,s)] and P_ij(x,,) ≡ ⁻¹[_ij(x,,s)]. Here we have used the convolution and the derivative properties of the Laplace transform.^45,46

Appendix A2. The fundamental solutions in the time domain for the standard solid model

Here we show how to obtain the expressions of the fundamental solutions in the time domain, namely U_ij(x,,) and P_ij(x,,). As shown in Appendix A1, the stress σ_ij and the strain ε_ij tensors are decomposed into the hydrostatic and deviatoric parts. For viscoelastic materials made of polymers such as elastomers and gels, the response to the hydrostatic pressure acts like purely elastic materials. We thus write σ_kk = 3Kε_kk (24) for the hydrostatic part.

We consider the standard solid (Zener) model which consists of a spring connected in parallel with a Maxwell element (a spring connected in series with a viscous damper). That means we can construct the following relations between the different components of stresses and strains. The superscript ‘Ma’ is used to represent the Maxwell element; ‘e2’ and ‘η’ are respectively for the elastic and the viscous component of the Maxwell element, and ‘e1’ for the elastic part that is connected in series with the Maxwell element. Hence for the deviatoric part we have s_ij = s ^Ma_ij + s ^e1_ij (25) and e_ij = e ^Ma_ij = e ^e1_ij. (26) The Maxwell element gives the following relations as s ^Ma_ij = s ^e2_ij = s ^η_ij (27) e ^Ma_ij = e ^e2_ij + e ^η_ij. (28)

The constitutive relations for each component are given as s ^e1_ij = 2Ge ^e1_ij, s ^e2_ij = 2G₁e ^e2_ij, s ^η_ij = 2ηė ^η_ij. (29)

Using the above relations we end up with the constitutive relation for the whole deviatoric part of the Zener model given as (30) In terms of the generalized differential form of the constitutive relation in eqn (15), we get for the Zener model (31) and the other coefficients of p_l and q_m being vanishing. Next, we impose the Laplace transform to the differential form of the constitutive relation (30), we end up with the constitutive relation in Laplace space given as (32) where (33) and (34)

Including the hydrostatic part, we have the full constitutive relation in the Laplace space (35) where (36) Hence we can define the bulk modulus and the shear modulus in Laplace space as (37)

To obtain the fundamental solutions U_ij(x,,t) ≡ ⁻¹[Û_ij(x,,s)] and P_ij(x,,t) ≡ ⁻¹[_ij(x,,s)], what we need to do is to perform the inverse Laplace transform of Ĵ₁(s), Ĵ₂(s), Ĵ₃(s) and Ĵ₄(s) using the expressions of (s) and Ĝ(s) in (37), we then obtain (38) (39) (40) (41) At t = 0, the functions are given as (42) (43) (44) (45) where G₁₂ ≡ G + G₁.

When t → ∞, we have (46) (47) (48) (49)

Acknowledgements

TSC gratefully acknowledges financial support from the Research Council of Norway (Project No. 315110). TSC thanks Christian Pedersen for the fruitful discussions and the referees of the paper for the inspiring comments to help to improve the manuscript.

References

Footnote