The effect of contact angle hysteresis on a droplet in a viscoelastic two-phase system

We investigate the dynamic behaviour of a two-dimensional (2D) droplet adhering to a wall in Poiseuille flow at low Reynolds numbers, in a system where either the droplet is viscoelastic (V/N) or the surrounding medium (N/V). The results reveal that the deformation of the viscoelastic drop over time is changed due to the presence of polymeric molecules. In the first stage, the viscoelastic droplet speeds up and deforms faster, while in the second stage, the Newtonian counterpart accelerates and its deformation outpaces the viscoelastic droplet. The deformation of viscoelastic drop is retarded significantly in the second stage with increasing Deborah number $De$. In the V/N case, the viscous bending is enhanced on the receding side for small $De$, but it is weakened by further increase in $De$, and this non-monotonic behavior brings about an increase in the receding contact line velocity at small $De$ and a decrease at large $De$. On the advancing side, the viscous bending is decreased monotonically for $Ca<0.25$, and hence the advancing contact line velocity is decreased with increasing $De$; however, the viscous bending presents a non-monotonic behavior for $Ca=0.25$ at the advancing side. The non-monotonic behavior on the receding side is attributed to the emergence of outward pulling stresses in the vicinity of the receding contact line and the inception of strain-hardening at higher $De$, while the reduction in the viscous bending at the advancing side is the result of just strain-hardening due to the presence of dominant extensional flow on the advancing side. Finally, in the N/V system, the viscoelasticity of the medium suppresses the droplet deformation on both receding and advancing sides, and this effect is more pronounced with increasing De; the weakening effect of viscous bending is enhanced significantly at the advancing side by increasing the Giesekus mobility parameter.


Introduction
The wetting of a solid substrate by a polymeric drop is an important process in many applications such as inkjet printing, coating, and self-cleaning surfaces.Among these applications, superhydrophobic surfaces have recently drawn the attention of many researchers due to their application in drag reduction and self-cleaning (Cheng and Rodak, 2005;Neinhuis and Barthlott, 1997).These surfaces are inspired by natural surfaces, such as gecko skin (Watson et al., 2015) and lotus leaves that give rise to the so-called lotus effect (Shirtcliffe et al., 2010).Superhydrophobic surfaces are usually characterized (on a macroscopic length scale) by the advancing contact angle θ A and the receding contact angle θ R , and the difference between these two contact angles called Contact Angle Hysteresis (CAH).The properties of a surface such as microscopic roughness and inhomogeneities may lead to CAH (Joanny and Robbins, 1990), and the movement of a droplet on the surface can be significantly affected by CAH.Furthermore, the non-Newtonian properties of a droplet affect its motion on a solid surface.For example, the elasticity of a droplet moderately speeds up the contact-line motion on a hydrophilic flat and chemically homogeneous surface at late stages of wetting and for a high-viscosity liquid (Yue and Feng, 2012;Wang et al., 2015Wang et al., , 2017)), but a recent study (Yada et al., 2023) revealed that increasing the polymer concentration does not modify the contact-line friction for aqueous polymer (shear-thinning) solutions on the smooth surface in the initial wetting regime.Thus, the elasticity of a drop has a negligible effect on the early stages of the wetting process, and the contact line motion at early times can be controlled by changing the viscosity of the solvent or adding roughness/inhomogeneities to the surface.
In contrast to the smooth hydrophilic and hydrophobic surfaces, the elasticity of a droplet plays an important role on the superhydrophobic surfaces.Experiment results suggest that the superhydrophobic surface may lose its selfcleaning property when they interact with highly elastic fluids (Xu et al., 2018).The motion of a Boger fluid drop, which shows viscoelasticity without shear thinning, on the superhydrophobic surface is strikingly different from a Newtonian drop.The velocity of the viscoelastic drop is notably reduced in comparison with the Newtonian drop, and complex branch-like patterns are left on the tilted superhydrophobic surface.Previous studies also showed that the polyacrylamide drops at higher polymer concentrations were noticeably more elongated when sliding down a hydrophilic surface with small CAH (Varagnolo et al., 2017).The wetting area of a viscoelastic drop in a Newtonian surrounding fluid has been studied by numerical simulations using the Oldroyd-B model (Wang et al., 2023).They found that the wetting area decreased when the elasticity of the fluid increased, and the same was observed when the polymeric viscosity ratio decreased.Also droplet motion on lubricant surfaces has been analysed.A lubricant surface is constructed by trapping a suitable low surface tension lubricating liquid inside a texture (Lafuma and Quéré, 2011), and these surfaces hardly pin drops.An oscillatory motion has been reported for the viscoelastic drop having both shear-thinning and elasticity on lubricant surfaces (Sartori et al., 2022).
Our aim is to model the effect of CAH on viscoelastic fluids on a continuum level, neglecting the microscopic details of the roughness, using a model that can handle a dynamic contact angle, but still allows receding and advancing contact angles to have different values.The most popular approach has been developed by Spelt (2005) for a level-set method.An intermediate contact angle is computed, and the contact line pins if it is within the hysteresis window [θ R , θ A ], otherwise it moves.This approach has then been extended to the phase-field method (Ding and Spelt, 2008), volume-of-fluid method (Dupont and Legendre, 2010), lattice Boltzmann method (Ba et al., 2013), and fronttracking method (Shin et al., 2018).It should be noted that the above methods need ghost cells outside the boundary to impose the contact angle condition and a slip model to circumvent the well-known singularity in the contact line (Batchelor et al., 2003); however, the Cahn-Hilliard equation automatically regularizes the singularity at the contact line (Jacqmin, 2000).A CAH model capable of capturing pinning, advancing and receding automatically without needing any explicit knowledge of contact line velocity was proposed by Yue (2020), and it has been further extended to level-set method (Zhang and Yue, 2020).Recently, an angle-dependent line friction model for the Cahn-Hilliard equation has been proposed by Amberg (2022), and the numerical results of the proposed model were in a very good agreement with the experimental data.
In the present work, we investigate the effect of a fluid's viscoelastic properties on the deformation and movement of a droplet with contact line hysteresis in a 2D Poiseuille flow at low Reynolds number (Re) by direct numerical simulation.The deformation and yield condition (depinning of contact lines) of both 2D and three-dimensional (3D) Newtonian drops were the subject of many previous studies at low and high Re, which are summarized below.
The yield condition of a 2D droplet in a shear flow and stokes regime was examined for a broad range of capillary numbers (Ca), viscosity ratios (λ µ ), and CAH (∆θ = θ A − θ R ) (Dimitrakopoulos and Higdon, 1997).It is found that increasing ∆θ results in an increase in the critical capillary number (Ca cr ), above which the droplet starts to move.Their study also demonstrated that a droplet with a higher viscosity (higher λ µ ) on a surface with a specific CAH starts to move at a smaller Ca cr than a droplet with a lower viscosity.
The effect of Ca, λ µ , slip velocity, and surfactants on the 2D droplet deformation and yield condition in Stokes flow regime was studied by Schleizer and Bonnecaze (1999) in both shear-driven and pressure-driven flows.They demonstrated that surfactants reduce the deformation of a drop due to the generation of marangoni stresses along the droplet surface; in addition, a decrease in the deformation of the drop when the surface is hydrophilic and an increase in the deformation when the surface is hydrophobic have been observed in the presence of a slip velocity without considering CAH.
Later, 2D droplet deformation in shear flow with pinned and moving contact lines was studied at moderate Re using a level-set method by Spelt (2006).The author showed that the critical Weber number (W e cr ) approaches a constant value when Re increases.The droplet becomes unstable for W e > W e cr , and the transient behaviour is very different at different Reynolds numbers.It was demonstrated that the slip at the wall increases the W e cr , and the slip length has a significant effect on the contact line speed.It should be mentioned that the behavior of the droplet with slip velocity was reported to be dependent on the slip length.
Afterward, the critical conditions for the motion of a 3D droplet at both low and moderate Re in the presence of CAH were studied (Ding and Spelt, 2008).The inertial effects led to a complicated structure around the deformed drop with a strong 3D effect for the moderate Re, whereas their results were in agreement with previous 2D studies at low Re.
In the present study, we conduct direct numerical simulations to study the effect of fluid viscoelasticity on the deformation and depinning of a 2D drop either in V/N or N/V systems over a surface in the presence of CAH at low Re.
In particular, we focus on the physical effects of polymers on the deformation and depinning processes.We apply the Cahn-Hilliard equation with a dynamic contact angle boundary condition for capturing the interface between the two phases, and the Giesekus constitutive equation to model the viscoelasticity of the drop.First, the qualitative effect of the elasticity is investigated over a wide range of the Ca at a low constant De = 1.We compare the time evolution of both dynamic contact angle and velocity of the receding and advancing contact lines for a viscoelastic droplet (V/N) and a Newtonian droplet in a Newtonian medium over surfaces with a varying but significant degree of contact line hysteresis (∆θ = [100 • , 130 • , 140 • ]).Next, the effect of De is studied in the V/N system.Finally, the effect of De and α has been investigated on the drop's deformation and the depping of contact lines in the N/V system.

Governing equations
Consider an incompressible two-phase system composed of a Newtonian fluid with viscosity µ n = µ s and a viscoelastic (Giesekus) fluid with solvent viscosity µ s , and polymeric viscosity µ p .A phase-field variable ϕ varies smoothly from phase one (ϕ = 1) to another phase (ϕ = −1) with ϕ = 0 denoting the fluid/fluid interface.This two-phase system is modeled by the following coupled system of equations (Yue et al., 2004;Abels et al., 2012): the linear momentum and continuity equations the Giesekus constitutive model for the viscoelastic fluid and the Cahn-Hilliard equation for the movement of the interface In the above equations, u(x, t) is the velocity vector, p(x, t) is the pressure, and τ p is the extra stress due to the polymers.In the Giesekus equation, λ H is the polymer relaxation time, and α is the Giesekus mobility parameter.
In the Cahn-Hilliard equation, G is the chemical potential, M is the mobility parameter, λ is the mixing energy density which is related to the surface tension in the sharp-interface limit (Yue et al., 2004) as The function h(ϕ) in Eq. ( 5) is defined as with F (ϕ) = λ 4η 2 (ϕ 2 − 1) 2 being the double-well potential, and where η is the capillary width of the interface.The density ρ and the dynamic viscosity µ of the mixture are a linear function of the phase-field variable as: The total viscosity of the non-Newtonian phase is µ = µ s + µ p .The flux of density due to the diffusion of components J (Abels et al., 2012) in Eq. ( 1) is given by The equations (1), (3), and (4) should be supplemented with appropriate boundary and initial conditions.The following boundary conditions are imposed: Equations ( 11) and ( 12), respectively, impose the no-slip and zero diffusive flux at the wall, where n is the outward pointing normal vector to the boundary.Equations ( 13) and ( 14) correspond to the wall energy relaxation (Jacqmin, 2000;Carlson et al., 2009) and neutral wettability condition, respectively.∂Ω\∂Ω w denotes the boundaries other than the solid walls intersected with contact line, and µ f is the contact line friction.The surface chemical potential L(ϕ, ∇ϕ) is defined as where cosθ is the wall energy.

Numerical method
The governing equations are discretized by a second-order finite difference method on a dual-resolution grid.The numerical method used in the present study has been described in detail and validated in (Bazesefidpar et al., 2022), except for the Cahn-Hilliard equation which is described in the following section.

Cahn-Hilliard equation
Different approaches have been proposed to deal with the non-linear term in the Cahn-Hilliard equation (Eq.4) embedded in the chemical potential G (Gao and Wang, 2014;Shen and Yang, 2010;Guillén-González and Tierra, 2013;Shen et al., 2018), but scalar auxiliary variable (SAV) proposed by (Shen et al., 2018) has the advantages that the resulting discretized equations are linear, unconditionally energy stable, and relatively straight-forward to implement.Challenges still remain when using this method, and the accuracy of SAV scheme declines if too large time steps are chosen so that the results are generally no longer accurate.One strategy to overcome this issue is to use adaptive time stepping strategy (Shen et al., 2019) or the relaxed-SAV (RSAV) method (Jiang et al., 2022;Zhang and Shen, 2022), and the RSAV method proposed by Jiang et al. (2022) has been used in the present study.The Equations (4-5) with the boundary conditions (12-14) can be reformulated and discretized in time by using RSAV-BDF2 scheme (Jiang et al., 2022) as where E(t) = Ω F (ϕ), g * ,n+1 = 2g n − g n−1 represents the second-order explicit approximation of g n+1 , ĝ = 2g n − 1 2 g n−1 , and γ 0 = 3 2 .The L n+1 A,R on the righthand side of Eq. 20 is defined as (Yue, 2020) where the explicit second order approximation ϕ * ,n+1 of the phase field function has been used in computing L A and L R with A and R corresponding to advancing and receding contact angles respectively.If the contact angle is between the hysteresis window [θ R , θ A ], the contact line is pinned; otherwise it moves.The speed of the receding contact line (when it is not pinned) depends on the θ d,R , and the advancing contact line speed on the θ d,A , where θ d is the dynamic contact angle defined on the wall (y = 0).The above system of equations (16-21) can be transformed into two coupled equations (Yang et al., 2019) by introducing new field functions as resulting in the following final set of equations.
The first coupled system of equations for ϕ 1 and G 1 is with the following boundary conditions and the second coupled system of equations for ϕ 2 and G 2 is with the following boundary conditions First, the coupled equations (26-27) with the boundary conditions (28-29) are solved, and then the coupled equations (30-31) with the boundary conditions (32-33) are solved.The resultant coupled equations are non-symmetric linear system so that GMRES iterative solver with geometric multigrid preconditioner with the Message-Passing Interface (MPI) for distributed-memory parallelization provided by PETSc library (Balay et al., 2019) has been used to solve them.Then, ϕ n+1 and G n+1 is computed by using Eq.(24-25).At the end, the relaxation technique Eq. ( 22) is used to improve the accuracy and consistency of the SAV scheme.The interested readers are referred to (Jiang et al., 2022) for the definition of ζ 0 .) is imposed at the inlet and outlet of the channel.The droplet initially is semicircle with contact angle θ 0 = 60 • , radius R = 0.9023H, and height h = 0.4511H.This test case is defined by the length scale L ref = H, which is the half channel width; u ref = ( 3U H )h, which is the velocity scale.This results in the following dimensionless numbers (Schleizer and Bonnecaze, 1999;Dupont and Legendre, 2010;Liu et al., 2015;Yue, 2020).From the physical parameters of the problem, we obtain the Reynolds number, Re = , the ratio between inertial and viscous forces; the Capillary number, Ca = µ1u ref σ , the ratio between viscous and surface tension forces, density ratio, λ ρ = ρ2 ρ1 , the ratio between ambient density to droplet density, and the viscosity ratio, λ µ = µ2 µ1 , which is the ratio between ambient viscosity to droplet viscosity.Furthermore, from the phase field model, we obtain the Peclet number P e = 2   (Schleizer and Bonnecaze, 1999) and (Yue, 2020) Following Yue (2020), we consider the following values of the relevant dimensionless numbers: The comparison of the droplet interface ϕ = 0 with the results by Schleizer and Bonnecaze (1999); Yue (2020) at different Capillary numbers when they reach steady state represents a good agreement except at Ca = 0.15 at the receding contact angle, see Fig. 2; the same agreement and error at Ca = 0.15 is indeed reported by (Yue, 2020), and it is due to a very small receding contact angle.

Results
We scrutinize the effect of the viscoelasticity on both deformation and depinning of a drop over a surface with different degrees of the hysteresis window, ∆θ = 140 • with the choice of θ R = 30 • and θ A = 170 • , ∆θ = 130 • with the choice of θ R = 30 • and θ A = 160 • , and ∆θ = 100 • with the choice of θ R = 40 • and θ A = 140 • .These hysteresis windows are larger than encountered in physical experiments, and are chosen to prolong the pinning of the contact lines, following previous numerical works of (Dimitrakopoulos and Higdon, 1997;Schleizer and Bonnecaze, 1999;Dimitrakopoulos, 2007;Ding et al., 2010;Dupont and Legendre, 2010).A rectangular domain, Ω = [0, 2H]×[0, 8H], is used for the V/N system with the droplet placed on the middle of the bottom wall, and a larger domain, Ω = [0, 2H] × [0, 16H], for the N/V system in order to allow a fully developed polymeric stresses in the inlet and outlet boundaries.It should be noted that the velocity inside the droplet has been set to zero in the beginning of the simulation, but a parabolic velocity profile is imposed elsewhere.
In addition to the Newtonian dimensionless numbers introduced in section 4.1, the following dimensionless numbers have been used to study the effect of the viscoelasticity : the dimensionless cross-sectional area A = (4A d /L 2 y ) representing the ratio between the droplet area (A d ) and the square of channel half height, the Deborah number De = (λ H u ref /H) representing the ratio between elastic forces to the viscous forces, and the polymeric viscosity ratio β = (µ s /µ) which is the ratio between the solvent viscosity and the total viscosity.The dimensionless cross-sectional area is kept constant at A = 0.5 in the following sections (section 5.1 and section 5.2), and the semi-circular droplet's initial radius r 0 = A d H/(θ 0 − sin θ 0 cos θ 0 ) and the droplet's initial height h = r 0 (1 − cos θ 0 ) are computed for a specific initial contact angle (θ 0 = 90 • ).Furthermore, the following dimensionless numbers are kept constant in all simulations: Re = 0.5 (so that the inertial effects are negligible), λ µ = 1, λ ρ = 1, and β = 0.1.It should be noted that λ µ is the ratio between total viscosities of the two fluids; Since polymeric viscosity increases with polymer concentration, the viscosity of the Newtonian fluid also needs to increase if the polymer concentration changes (for example by adding more viscous fluids).
Regarding numerical parameters related to the phase field model, we choose an affordable but small capillary length (Cn = η/H = 0.00375), and the Peclet number is set to P e = (2 √ 2/3)(Ca/Cn) according to the guidelines (Yue et al., 2010;Yue, 2020) in order to reach the so-called sharp-interface limit; the contact line friction is set to µ f = (Λ/Cn)µ 1 (Carlson et al., 2011;Yue, 2020) with Λ = 0.1.The contact line is pinned when ∂ϕ ∂t | t n+1 = 0, so the droplet's deformation and the contact line's pinning/depinning are independent of the choice of µ f prior to depinning.
To analyse the results further, we have extracted different quantities in the simulation, e.g., the locations of the receding contact line X cl,R , and the advancing contact line X cl,A , and the values of the dynamic contact angle at the receding side θ D,R , and at the advancing side θ D,A .Here, the dynamic contact angle has been computed using the geometric expression cos θ D = (n • ∇ϕ)/(|∇ϕ|) at the interface (ϕ = 0) at the height y = 2 Cn above the wall (Wang et al., 2015).The velocity of the contact line (U cl = σU * cl µ f ) is computed by using second-order central difference with extracted X cl in the simulation, and the average contact line velocity (U cl ) is computed by taking average over a nominal interval [0, t end ].
In the following sections, we present results for non-Newtonian cases, where either the droplet is laden with polymers and therefore viscoelastic (V/N), or the doplet is Newtonian but the ambient fluid is viscoelastic (N/V).In both cases, polymers affect the droplet deformation and dynamics in two ways: the polymeric stresses affect the interfacial deformation, and they also induce changes in the flow field which can affect the droplet dynamics (Ramaswamy and Leal, 1999;Yue et al., 2005b;Aggarwal and Sarkar, 2007).The stress tensor excluding the term due to mixing energy can be written as (Yue et al., 2005a) Components of the stress tensor normal (n s ) and tangential (t) to the interface (ϕ = 0), see Fig. 1, for the N/N, V/N, and N/V systems are shown along the interface in order to perceive the dominant forces there.Furthermore, our numerical experiments showed that the polymeric stresses inside the droplet were growing without bound for De ≥ 5 when using Oldroyd-B model (α = 0).It is known that results from Oldroyd-B model may not be accurate in extensionally dominated flows, and it will be shown later that there are such regions in these flows.Therefore, we employ the Giesekus model by setting α = 0.05 for the V/N and α = [0.05,0.1] for N/V systems in order to prevent this issue while producing sufficiently large polymeric stresses to be able to observe viscoelastic effects.

Elastic effects in the V/N system
In this section, we investigate the case where a viscoelastic droplet is surrounded by a Newtonian fluid (V/N system), while the capillary number, contact line hysteresis window and Deborah number are varied.

Effect of the capillary number
In this section, we investigate the effect of the capillary number on the deformation of the viscoelastic droplet, at a fixed Deborah number (De = 1) and at different values of receding and advancing contact angles (θ R and θ A ). Figure 3 presents the time evolution of the instantaneous contact angle θ D at both the receding and advancing sides over two different surfaces with ∆θ = 100 • and ∆θ = 140 • , and we observe that viscoelasticity of the droplet affects them differently.Figures 3(a counterpart for small Ca when the contact line is pinned, and θ D,R gains negligibly larger value when contact line is moving for (De ≤ 1), see Fig. 5(a).The advancing contact angle θ D,A reaches a smaller equilibrium value for Ca ≤ 0.2 and a larger value for Ca > 0.2, see Figs. 3(b) and 3(c); It can be observed that θ V E D,A is approaching θ N D,A and surpassing it as the advancing contact line moves with considerable velocity which depends on both Ca and ∆θ, see Fig. 5(b).It is worth mentioning that for larger values of De, the behaviour of θ D,R is different, as will be discussed in section 5.1.2.
The transient behavior of the θ D reveals that the difference between dynamic contact angles (θ V E D − θ N D ) for a specific Ca is not monotonic in time, and this trend is more pronounced for the (θ D,R ) in the early times, see Figs. 3(a) and 3(c) for t ≤ 2 s.The reason is that the Newtonian and viscoelastic droplets deform at different rates.The viscoelastic drop speeds up and deforms faster than the Newtonian counterpart at early times, and this behavior has been  Thus, neither the viscous stress nor the polymeric stress is influential in the faster deformation of the viscoelastic droplet at the early times.The reason is instead found from the pressure distribution at the interface (Fig. 4 (e)).The pressure pulls the droplet inward for all values of the φ, and its magnitude is larger in comparison to other components due to the presence of the surface tension in the interface (Young-Laplace equation).The magnitude of the pressure is larger for φ ≲ 20 • and smaller for φ ≳ 158 • for the viscoelastic droplet in comparison to the Newtonian counterpart; thus, the pressure is responsible for the faster decrease of the dynamic contact angle in the receding side, and faster increase on the advancing side, both of them increasing the deformation of the viscoelastic droplet at t = 0.1 s.The change in the pressure is the consequence of the modification of the overall flow field by polymeric stresses.
Next, we will look at the effect of contact line hysteresis window ∆θ = θ A −θ R on pinning and velocity of the contact lines.Firstly, we observe that as hysteresis increases, both Newtonian and viscoelastic droplet are pinned at higher capillary numbers.Secondly, for droplets that are not pinned, the steady state contact line speed is higher for viscoelastic droplets when Ca is larger (Ca ≥ 0.25).Figs.5(a,b) display the time evolution of the contact line velocity for different Ca over surfaces with ∆θ = 100 • , and Figs.5(c,d) present the variation of averaged dimensionless contact line velocity (U * cl ) vs. Ca over surfaces with different ∆θ.To characterise the depinning process, the time for the inception of motion (U cl (t in ) ≈ 10 −4 m s ) in both contact lines have been marked on their corresponding velocity curve for the hysteresis window ∆θ = 100 • , Figs. 5(a,b).Both of the contact lines are pinned when Ca < 0.2.When Ca ≥ 0.2, the advancing contact line starts to move first (when the θ D,A | y=0 ≤ θ A at t ≤ t in ) for the Newtonian drop; the receding contact line of the Giesekus drop moves first in this cases.At this moderate De, the both advancing and receding contact lines of the viscoelastic droplet move faster over all surfaces (different ∆θ) with Ca = 0.25, but the advancing contact line of the Giesekus drop moves faster when Ca = 0.25 and slower when Ca = 0.2 in comparison with its Newtonian counterpart, see Figs. 5(c,d).
The contact line speed can be approximated by the difference between the dynamic contact angle θ and the equilibrium contact angle θ e (θ A and θ R in the present of CAH) for minor variation of θ from θ e (Yue and Feng, 2011;Amberg, 2022;Yada et al., 2023): where the contact line friction µ f is a constant determined by surface and fluid properties.Hence, the contact line speed is enhanced if the droplet experiences more bending causing the dynamic contact angle to differ from the equilibrium contact angle.At the receding contact line, more bending implies a decrease in θ D,R .The elasticity (De = 1) of the droplet enhances the viscous bending at the receding contact line (θ V E D,R ≤ θ N D,R ) after its depinning over a large span of time, see Figs. 3(a,c), for the simulations in this section and therefore θ However, the viscous bending at the advancing contact line is initially weakened by the elasticity of the droplet (θ V E D,A ≤ θ N D,A ).Thus, the elasticity of the droplet causes In this subsection, we have studied the effects of moderate viscoelasticity on the deformation of the droplet and movement of both contact lines by comparing a Newtonian droplet and a viscoelastic droplet with De = 1.Some of the conclusions will change when viscoelasticity is further increased in the next subsection.

Effect of increasing Deborah number
The effect of increasing the droplet's elasticity at a fixed capillary number has been investigated by changing De over the surfaces with ∆θ = 100 • and ∆θ = 140 • .The main effect observed is that bending of both contact lines is decreased significantly with increasing Deborah number to De = 5 − 10.Thus, we have a non-monotonic behavior with respect to De in the receding contact line; the elasticity of the droplet enhances the viscous bending for droplets with small elasticity, but the viscous bending is decreased for droplets with larger De.This non-monotonic behaviour in the receding contact angle also changes the receding contact line's velocity trend, and the velocity of the receding contact line decreases for De = [5,10] in comparison to the Newtonian counterpart in all cases, see Fig. 7. Further increasing the elasticity of the droplet depending on the hysteresis window of the surface and capillary number brings the contact lines to a halt, see Figs.The advancing dynamic contact angle of the viscoelastic droplet generally attains a smaller steady-state value with further increasing the elasticity to De = 10 (Fig. 6), and the velocity of the advancing contact line decreases with increasing the elasticity of the droplet, see Figs.To gain a better understanding of the elastic effects, we have visualized the trace of the conformation tensor c defined as c = (λ H /µ p )τ p + I, which is a measure of stretching of the polymer molecules, in figure 8 for the highest Deborah number De = 10 at Ca = 0.2 and ∆θ = 100 • (corresponding to drops in Figs.6(c) and 7(a,b)).The results for De = 5 are very similar to De = 10.The interface (defined by ϕ = 0) of the viscoelastic drop is shown with a black contour, and the corresponding Newtonian drop by a white contour, with the trace of the conformation tensor in the background, see Fig. 8. Furthermore, we show by arrows the principal direction of the conformation tensor, representing the orientation of the dumbbells.The pull of the polymers affects the interface position most in regions where the arrows are perpendicular to the interface ( (Bird et al., 1987b), pp. 64-69).
At early times (Figs.8a), the viscoelastic drop speeds up and deforms faster than the Newtonian counterpart.The polymer molecules are most stretched in the vicinity of the receding and advancing contact lines at t ≈ 0.15 s, but the dumbbells have acute angle with the interface there, and hence do not influence the interface.At t ≈ 0.6 s (Fig. 8(b)), the receding and advancing contact lines of the Newtonian drop start moving already.However, the polymeric stresses have now built up (which takes a finite time proportional to λ H ). Hence, deformation and time evolution of the viscoelastic droplet have been retarded,  and both contact lines of the viscoelastic droplets remain pinned for De = 5−10.Furthermore, the dumbbells are becoming perpendicular to the interface in φ ≈ [140 • , 160 • ] and because polymers are also very elongated here, the stresses should have a large influence on the interface.
To compare the relative importance of polymeric and flow-induced stresses, we again show them (Eq.34) along the interface ϕ = 0, separated into normal and tangential components, in Fig. 9. Increasing viscoelasticity to De = 5 − 10 suppresses the tangential component (Fig. 9(b)) of the viscous stress along the interface ( 3), and this decrease is mainly responsible for the delay in the deformation of the viscoelastic droplets.The change in viscous stresses due to viscoelasticity has a larger magnitude than the previously analysed polymeric stresses (Fig. 9(c,d)) which are only active in the segment φ ≈ [110 • , 170 • ].Further, the polymeric stresses are largest for De = 5 and decrease at De = 10, while the contact line stays pinned longest for De = 10.Hence, the retardation and pinning of the viscoelastic droplet is mainly caused by velocity modification in viscoelastic flows around the interface, rather than by the polymeric stresses themselves.
Moreover, another important rheological property of the viscoelastic fluid is the strain-hardening which usually occurs in extensional flow, and can be To further contrast the effects at small and large elasticity, we compare the flow topology of the viscoelastic droplets at De = 0 (Newtonian case, Fig. 10 (a)), De = 1 (Fig. 10 (b)) and De = 5 (Fig. 10(c)).The flow type inside the viscoelastic droplet with De = 5 is more extensional (red color) at the front where viscoelastic stresses pull the droplet interface, in comparison to the other droplets.This is the extensional flow region, where extensional viscosity may be higher.The extensional viscosity as a function of the extensional strain rate at steady-state for the Giesekus model (α ̸ = 0) can be expressed as: ((Bird et al., 1987a), p.p 368) Fig. 10(d) presents µ + e as a function of strain rate1 ε for three Deborah numbers.The magnitude of µ + e for De ≥ 5 increases steeply already at small elongational strain rates.The strain rate for De = 5 shown in figure 10 (f) is highest in sheardominated regions at the top, but also moderate in extensional flow areas at the front, and therefore we should expect a large extensional viscosity resisting deformation at the front.As observed, the shear strain rate is always higher in the viscoelastic fluid than in the Newtonian fluid, and this is a consequence of keeping the total viscosity ratio constant (λ µ = 1); the solvent viscosity of the viscoelastic droplet is smaller.The magnitude of N 1 is largest in the top side where shear flow dominates and increases with De,but those regions of the interface do not directly influence the contact line movement.
The magnitude and tendency (direction) of the total stress vector in Eq. 34 along the interface ϕ = 0 is depicted by arrows on top of figures 10 (a,b,c):(n s •τ , where n s is the unit vector normal to the interface.The stresses are mostly pulling the interface inwards for both Newtonian (10a) and viscoelastic (10b,c) droplets.However, in the vicinity of the receding contact line marked with circle, the stresses are pulling the interface outward for the viscoelastic droplets.The magnitude of these outward stresses increases with De (elasticity of the droplet) so that they are bringing about a larger θ D,R .The outward stresses around the receding contact line were not observed when capillary number was lowered to Ca = 0.15, at least until De = 10.
Concluding, the non-monotonic behavior in the receding dynamic contact angle with increasing De can be attributed to the combination of two effects: the outward stresses in the vicinity of the receding contact line, and the inception of the strain-hardening occurring for De ≥ 5. Also the smaller advancing dynamic contact angle at steady-state might be attributed to strain-hardening.

Effect of matrix viscoelasticity on a Newtonian droplet
In this section, we consider the effect of matrix viscoelasticity on the deformation and depinning of a Newtonian droplet (N/V system) in the presence of contact line hysteresis.The same boundary and initial conditions have been applied as previously, with the polymeric stresses are set to zero in the beginning of the simulation in the Giesekus medium.Figs.11(a,c,e) present the time evolution of θ D,R of a Newtonian drop in the Giesekus medium over a surface with ∆θ = 140 • , and it can be seen that the elasticity of the matrix has suppressed the drop's deformation in the receding side significantly for De ≥ 5, in other words, the viscous bending in the receding side has weakened drastically.The effect of changing capillary number on the steady-state value of θ D,R is limited and small in the range Ca = [0.15,0.25] used in our simulations when De ≥ 5.The same trend can be observed on the advancing side, where the viscoelasticity of the medium decreases the viscous bending in particular for De = 5 − 10; Furthermore, the shear-thinning also plays an important role at higher De by further reducing the viscous bending, see Figs. 11(b,d,f).
It is worth mentioning that the decreased viscous bending results in that both contact angles are insensitive to capillary number for the Newtonian drop in the viscoelastic matrix with high elasticity.For a Newtonian droplet in a Newtonian matrix, when increasing capillary number from Ca = 0.15 to Ca = 0.25, both receding and advancing contact angles vary by ∝ 30 • .For a Newtonian droplet in a viscoelastic matrix (at De = 5 and De = 10), the change with capillary number is less than 10 • .Both contact lines are pinned at Ca = 0.15 − 0.25 over the surface with ∆θ = 140 • in a viscoelastic matrix (Fig. 12(a)), while the receding contact line of the Newtonian drop in the Newtonian medium moves with a significant velocity at Ca = 0.25.If the hysteresis window is decreased to ∆θ = 100 • (Fig. 12(b)), the advancing contact line of the least elastic medium (De = 1) starts moving at t ≈ 0.85 s and the receding contact line remains pinned, while both contact lines start to move simultaneously in a Newtonian matrix with appreciable velocities.The trace of conformation tensor along with principal stress orientation is shown at De = 1 13(a,c) and De = 5 13(c,d), respectively, over a surface with ∆θ = 100 • .The stretching of polymer molecules is small in the vicinity of the contact lines, while large stretching occurs around the interface for φ ≈ [54 • , 138 • ] since dumbbells are no longer parallel to the interface in this portion (as they were in the V/N case, fig.8), this implies that large polymeric stresses is expected over a large portion of the interface.
The normal component of the viscous stress τ n vi has been lessened significantly due to the elasticity of the matrix at t = 0.6 s for φ ≈ [0, 120 • ] responsible , is the consequence of the small strain rate in these regions.Summarizing, also viscoelasticity in the matrix phase in the N/V system suppresses the deformation of the droplet and pins the contact lines, and does so more efficiently than the viscoelasticity of the droplet in a V/N system.In both cases, the effect of viscoelasticity mainly stems from changes in the flow field and weakening the flow gradients around the droplet such as tangential viscous shear, which results in decreasing viscous bending and hence decreased movement of the contact lines.

Conclusions
In the present study, two-dimensional numerical simulations have been conducted to investigate the effect of viscoelasticity either in the drop or in the matrix on the depinning and deformation of a drop adhering to the wall in Poiseuille flow at low Reynolds number, under the influence of contact angle hysteresis.The results indicate that the viscoelastic properties have a significant impact on the drop's deformation and depping of the contact lines.
We first compare Newtonian and viscoelastic droplets at a low constant Deborah number (De = 1) and different capillary numbers.The elasticity of the drop causes the viscoelastic drop to deform faster than the Newtonian counterpart at the early initial times, but the Newtonian drop then catches up.The viscoelastic drop continues to deform, and its receding dynamic contact angle reaches a smaller final value than the Newtonian counterpart prior to depinning; the advancing dynamic contact angle of a viscoelastic drop reaches a smaller value if the receding contact angle is pinned, otherwise it attains a larger final value.This final stage of the deformation depends on both Ca and ∆θ for a fixed De.
Secondly, we investigated the effect of increasing elasticity (De = 5 and De = 10) over the surfaces with ∆θ = 140 • and ∆θ = 100 • .In the first stage of the deformation, the effect of increasing De is negligible, but in the second stage, the time evolution of the dynamic contact angle at both receding and advancing sides is hindered significantly with non-monotonic behavior in the receding dynamic contact angle.This trend is attributed to the emergence of the pulling outward stresses changing to inward along the interface depending on both Ca and De and the inception of the strain-hardening inside the drop at both receding and advancing sides.
Finally, we scrutinized the effect of matrix viscoelastic properties in N/V system over surfaces with ∆θ = 140 • and ∆θ = 100 • .The viscoelasticity of the matrix lessens the deformation of the Newtonian drop, in other words, the viscous bending has been weakened in both receding and advancing sides due to the presence of polymer molecules in the medium, and these reductions are more pronounced for De ≥ 5.These reductions are bring the drops in the N/V system to halt in the matrix with large elasticity (De ≥ 5).Furthermore, the shear-thinning effect reduces the the viscous bending further in the advancing side by increasing α from 0.05 to 0.1.
Eq. A8 anticipates a large variation in the velocity with a change in θ d due to the dependency on cot function over superhydrophobic surfaces since these surfaces usually own large contact angles.

Figure 1 :
Figure 1: The Initial configuration for the droplet in the 2D with θ 0 .θ D has been extracted at the height y = 2 Cn above the wall.
the ratio between the advection and diffusion; the Cahn number, Cn = η H , which represents the ratio between the interface width and the characteristic length scale; Λ = µ f Cn µ1 , which characterizes the wall energy relaxation rate.

Figure 4 :
Figure 4: Different components of the stresses along ϕ = 0 at time t = 0.1 s for drops with Ca = 0.2 and De = 1 over the surface with ∆θ = 140 • .(a) τ n vi (b) τ t vi (c) τ n p (d) τ t p (e) τ n pre .

Figure 5 :
Figure 5: The time evolution of U cl,R (X t=0 cl,R = 3.436) and U cl,A (X t=0 cl,A = 4.564) for the Newtonian and Giesekus drops with De = 1.N and VE refer to Newtonian and Giesekus drops respectively.(a) U cl,R for the surface with ∆θ = 100 • (b) U cl,A for the surface with ∆θ = 100 • (c) U * cl,R vs Ca (d) U * cl,A 7(a) and 7(b) and for De = 10, Ca = 0.2, and ∆θ = 100 • .
7(b), 7(d) and 7(f).The decrease in the velocity of the advancing contact line is in line with the decreasing θ D,A .The same trend for U cl,R exists in the receding contact line for De = [5 − 10] but θ cl,R gains a larger value, see Figs. 7(a), 7(c) and 7(d).

Figure 8 :
Figure 8: Time evolution of viscoelastic drops with Ca = 0.2 over the surface with ∆θ = 100 • .The interface ϕ = 0 of the Newtonian and Giesekus drops are depicted with a white and black, respectively.The trace of the conformation tensor has been visualized on the background, and the principal eigenvector of the conformation tensor is represented by purple line segments.(a)De = 10 at t ≈ 0.15 s, (b) De = 10 at t ≈ 0.6 s, (c) De = 10 at t ≈ 1.8 s, (d) De = 10 at t ≈ 6 s.

Figure 9 :
Figure 9: Different components of the stresses along ϕ = 0 at time t = 0.6 s for drops with Ca = 0.2 over the surface with ∆θ = 100 • .(a) τ n vi (b) τ t vi (c) τ n p (d) τ t p (e) τ n pre

Figure 10 :Figure 11 :Figure 12 :
Figure 10: Different flow quantities of the Newtonian and Giesekus drops with Ca = 0.2 over ∆θ = 140 • at t ≈ 6 s.(a)The flow parameter (ξ) on the background and the total stress vector (n s • τ in Eq. 34) along the interface ϕ = 0 for the Newtonian droplet with Ca = 0.2 and ∆θ = 140 • at t ≈ 6 s, (b) the same for a Giesekus drop with Ca = 0.2, De = 1, (c) the same at Ca = 0.2, De = 5.(d) The steady-state dimensionless elongational viscosity as a function of the elongational strain rate for the Giesekus model with different De corresponding to Ca = 0.2 and ∆θ = 140 • .(e)The magnitude of strain rate (| γ|) for the droplet at De = 1, and (f) the same at De = 5.(g) The first normal stress difference at De = 1, and (h) the same at De = 5.

Figure 13 :
Figure 13: Time evolution of the Newtonian drop in viscoelastic matrix with Ca = 0.25 over the surface with ∆θ = 100 • .The interface ϕ = 0 of the Newtonian and Giesekus drops are depicted with a white and black lines, respectively.The trace of the conformation tensor has been visualized on the background, and the principal eigenvector of the conformation tensor is represented by purple line segments.(a) De = 1 at t ≈ 0.6 s, (b) De = 5 at t ≈ 0.6 s, (c) De = 1 at t ≈ 1.8 s, (d) De = 5 at t ≈ 1.8 s.

Figure 14 :Figure 15 :
Figure 14: Different components of the stresses along ϕ = 0 at time t = 0.6 s for drops with Ca = 0.25 over the surface with ∆θ = 100 • in N/V system.(a) τ n vi