Modeling liquid droplet impact on a micropillar-arrayed viscoelastic surface via mechanically averaged responses

Droplet impact on a substrate is an intriguing phenomenon that widely exists in our daily life and a broad range of industrial processes. However, droplet impact dynamics on soft textured surfaces are less explored and the underlying mechanisms remain elusive. Here, we report numerical simulation of droplet impact dynamics on a micropillar-arrayed soft surface using BASILISK, which involves a multiscale geometric domain containing the micropillars and droplet that are in the order of $ \rmu {\rm m} $ μm and $ {\rm mm} $ mm, respectively. As such, the volume of fluid (VOF) method is coupled with the finite volume method (FVM) to build the fluid fields and track their interface. From a conceptual point of view, the micropillared substrate is formed by imposing interstitial gaps into the otherwise intact soft material, whose viscoelastic properties can be quantified by gap density ϵ. Via a five-parameter generalized Maxwell model, the viscoelastic properties of the micropillared substrate can be approximated by its equivalent elastic response in the Laplace–Carson (LC) space, and the averaged bulk strain of the micropillared substrate in the real space is obtained by the inverse LC transform. Moreover, through parametric studies of splash extent, it turns out that for a specific ϵ, the splash is dramatically intensified with increasing impact velocity $ U_i $ Ui. The splash also turns more violent with increasing ambient pressure $ P_a $ Pa, which is evidenced by a larger splash angle of $ 114.44^\circ $ 114.44∘ between the ejected sheet and the horizontal substrate at 5 atm. Conversely, the splash becomes more depressed with increasing surface tension σ. Overall, the splash magnitudes of our simulations agree well with those predicted by the Kelvin-Helmholtz instability theory. By leveraging the LC transform in the fluid-viscoelastic solid interactions, our simulation methodology captures the main features of droplet impact dynamics on microstructured viscoelastic surfaces by means of the mechanically averaged responses while avoiding the predicament of domain scale inconsistency.


Introduction
Droplet impact dynamics and the related splash phenomena on various substrates have been extensively studied over the past several decades.Contingent on droplet properties, ambient conditions, and surface topography, six distinct procedures such as deposition, prompt splash, corona splash, receding break-up, partial rebound, and complete rebound may take place in sequence during the entire droplet impinging process on a solid surface (Rioboo et al., 2001).Within the spreading stage, the droplet behaviour is mainly determined by the conversion of kinetic energy to surface energy.Wildeman et al. (2016) found that for a droplet with high impact velocities, nearly half of its kinetic energy would be converted to the surface energy of the spreading liquid film, regardless of the impact parameters or the exact energy budgets.Furthermore, capillary force impact velocity, and surface conditions (Lee et al., 2015;Thoroddsen & Sakakibara, 1998;Xu et al., 2005).Mandre and Brenner (2012) probed the formation and deflection of the ejected splash sheet on a smooth surface, and quantitatively delineated the splash threshold as a function of impact velocity and droplet size.Rioboo et al. (2002) found that there exists a dimensional similarity of the droplet impact shape during the initial impinging phase on a dry solid surface.
Since most of the droplet impact studies were focused on rigid substrates, many outstanding problems of droplet impingement on soft surfaces involving energy transport and force imposition have been less explored.Also, the droplet impact on soft surfaces presents some unique characteristics, such as an extended contact time and different evolutions of spreading and splashing.Langley et al. (2020) investigated the effects of substrate stiffness on droplet spreading and found that the transition from gliding to ring contact is delayed at lower stiffnesses, characterized by a smaller Weber number We (We = ρ l U 2 i D/σ , where ρ l and D are the density and diameter of the droplet, respectively, U i is impact velocity, and σ is surface tension).As surface wettability also affects contact time, Weisensee et al. (2016) studied the droplet impact on elastic nanostructured superhydrophobic substrates and found that the springboard effect could conversely lead to an impact time reduction.Moreover, compared to a rigid substrate, the spreading diameter develops distinctively on a soft surface and the initiation of splash generally requires more kinetic energy.Howland et al. (2016) studied droplet impact and splash on soft surfaces with varying stiffness, droplet size and viscosity, and air pressure.They concluded that droplet splash could be substantially mitigated by reducing the stiffness of the substrate, and it needs approximately 70% more kinetic energy to initiate the droplet splash on a soft surface.
As a unique soft material owning elasticity and viscosity simultaneously, the viscoelastic materials such as polymer-matrix composites, physiological tissues, and polymeric compounds have received increasing interest in a wide spectrum of researches (Christensen, 2012).Specifically, Ghezelbash et al. (2022) examined the timedependent viscoelastic behaviours of blood clots through the compression and shear stress relaxation tests.Rambausek et al. (2022) developed a numerical framework to model viscoelastic behaviours of finite strain magnetorheological elastomers which consist of magnetically hard or soft phases.Furthermore, a variety of mathematical models of viscoelasticity have been developed to describe the deformation, relaxation, and creep responses of viscoelastic materials.Using atomic force microscopy (AFM), Chyasnavichyus et al. (2014) measured the probe tip force-displacement curve on poly(n-butyl methacrylate) and analysed its properties via Sneddon's model for the elastic region and Johnson's model for the viscoelastic zone, which was carried out across the whole glass-transition process.Boisly et al. (2016) experimentally investigated the rate dependence and the characteristic force-deformation tendency of the cutting process of viscoelastic materials.Then, the generalized Maxwell model and the Mooney-Rivlin model were utilized in their finite element method (FEM) simulations to describe the linear viscoelasticity and the force-displacement relation, respectively.
Facing an external impetus, the responses of a viscoelastic bulk material are influenced not only by the viscoelastic properties of itself but also by its geometric configuration.Nguyen et al. (2017) developed a methodology that analogizes the bulk and shear moduli of the intact viscoelastic solids obtained by the generalized Maxwell model to the corresponding moduli of microcracked materials in the Laplace-Carson (LC) space.Then, they validated this analogical deduction in the real space with a simple loading case through the inverse LC transform.Regarding the contact properties of a viscoelastic micropillar, Gong et al. (2021) analysed the contact properties between a spherical asperity and a soft polymer substrate, which was validated by the FEM simulations, and subsequently studied the viscoelastic contact between an asperity and a composite micropillar.On a single viscoelastic polydimethylsiloxane (PDMS) micropillar, Du et al. (2013) calculated the cell contraction force in the frequency domain by FEM, in which the complex modulus was obtained by fitting the generalized Maxwell model with the least square method, and the reaction force was deduced by the cellular contraction data transformed to the Fourier series.
Regarding droplet impingement on a structured soft surface such as a micropillar-arrayed viscoelastic substrate, the aforementioned viscoelastic models need to be revisited to predict the elastic moduli and time dependent stress-strain relationship.According to the soft material's viscoelastic responses, linear viscoelasticity and nonlinear viscoelasticity models have been put forward respectively.In the case of the linear responses, stress could be scaled or superposed in the same way as the corresponding composite strain is processed mathematically (Wineman, 2009).Regarding the nonlinear viscoelasticity, it could result from the material's nonlinear characters or its geometric structure (Drapaca et al., 2007).
With sufficient accuracy and modest complexity, linear viscoelasticity is adopted by us to investigate the viscoelastic properties of the micropillar-arrayed substrates.In this respect, several typical models have been developed to describe the stress-strain relationship, including the Kelvin-Voigt model where the viscoelastic system consists of a spring and a dashpot in parallel, and the Maxwell model which is composed of the same components rather connected in series (Findley et al., 1977).As opposed to the Kelvin-Voigt model that has limitations in elucidating the stress relaxation and the Maxwell model that is deficient to handle the creep behaviour (Findley et al., 1977), the generalized Maxwell model, in which one additional spring is connected in parallel with other Maxwell elements (i.e. each Maxwell element comprises an elastic spring and a viscous damper in series), can account for both the stress relaxation and the creep behaviour simultaneously (Kohlrausch, 1863).In general, the more the Maxwell elements are integrated, the more accurate the generalized Maxwell model would be.However, this strategy might lead to a prohibitive computational cost along with implementation difficulty.Thus, a three-parameter generalized Maxwell model (i.e. the standard linear solid model (Lazan, 1968)) has been developed as a trade-off.Likewise, the five-parameter model has also been put forward to render a more accurate prediction of viscoelastic behaviours with moderate computational cost.Mishra and Patra (2018) utilized both linear and non-linear fiveparameter models to investigate the creep behaviour of pile groups in soil study, finding that compared with other numerical and experimental results the maximum errors induced by these two models are only 18% and 3%, respectively.
Because of the complexity in determining the viscoelastic properties of the microstructured soft substrate and the difficulty in describing the associated dropletsubstrate interactions, in which a multiscale geometric domain is involved including micropillars in the order of µm and the droplet in the order of mm, studying droplet impact on such a micropillar-arrayed viscoelastic substrate becomes rather challenging and arduous.Consequently, only few studies have been conducted on this topic.For instance, Parizi et al. (2007) conducted a parametric study of the molten nickel and zirconia droplets impacting on rigid micro-cube patterned surfaces, in which the finite difference method and the VOF approach were coupled to track the fluid-fluid interface.Choi et al. (2017) performed numerical simulations of the droplet impact and evaporation on a porous surface by solving vapour fraction-related equations, in which the effects of impact velocity, porosity, and particle size of the porous medium on the impingement and evaporation were carefully appraised.Guo et al. (2014) leveraged the level set method and the VOF approach to probe the two-dimensional (2D) impact dynamics on a liquid film and analysed the effects of impact velocity and film thickness on spreading diameter.Because of the complexity in force implementation and response determination in the three-dimensional (3D) space, there have been only few studies of droplet impact dynamics for the 3D scenario.Bussmann et al. (1999) proposed a volume tracking algorithm delineating the free surface of the droplet and used the software RIPPLE to determine the surface tension force as a volumetric force exerted on the free surface.Therefore, it is more feasible to study droplet impact on a micropillar-arrayed viscoelastic substrate in the 2D domain under reasonable assumptions.
In this work, we adopted the five-parameter viscoelastic model for the micropillar-arrayed substrate and then simulated droplet impact dynamics thereon by means of the mechanically averaged responses in the LC space.This work can be divided into several parts: (i) construct the discretized model for the fluid fields including the liquid droplet and the ambient gas as well as the interface between them by the coupled FVM and VOF method; build the mathematical model of the micropillar-arrayed viscoelastic substrate in the LC space, and realize the real-space coupling between the fluid field evolution and the substrate response via the inverse LC transform; (ii) based on the elastic moduli and viscosities of micropillared substrates characterized by gap density , conduct parametric simulations of droplet impact dynamics regarding the effects of diverse factors such as impact velocity U i , ambient pressure P a , and surface tension coefficient σ on the substrate deformation magnitude and the splashing intensity; and (iii) study the impact and splash mechanisms by comparing the simulation results with the Kelvin-Helmholtz (K-H) instability theory of the interface growth.Based on the LC transform, our simulation methodology can reconcile the geometric scale inconsistency between the impinging droplet and the micropillar-arrayed viscoelastic substrate by way of the mechanically averaged responses.

Theory and discretization method for liquid droplet and ambient gas phases
The simulation of droplet impact on a structured viscoelastic substrate necessitates modelling multiphase fluid flow and fluid-solid interactions in a coupled manner.In this work, we used the open-source codes BASILISK to numerically simulate the droplet impacting, spreading, and splashing processes on micropillararrayed viscoelastic surfaces by solving the incompressible Navier-Stokes (N-S) equations (Equations ( 1)-( 3)) (Temam, 2001) and the related VOF equations (Equations (4)-( 10)).The N-S equations are first discretized by FVM and then solved through the obtained algebraic equations.As represented in Equation (3), the velocity divergence is equal to zero within each computation step.At the impacted surface, the deformation velocity of the viscoelastic substrate is treated as a Dirichlet boundary condition to two distinct phases, i.e. the liquid droplet and the surrounding gas.And the volume of fluid (VOF) method is used to obtain the concrete volume fraction of each phase at each computation step by solving the advection equation (Equation ( 4)), which can be mathematically derived from the continuity equations based on the individual phase density (Equation ( 5)) and the weighted density (Equation ( 6)).According to the VOF theory, both phases share the same value of either the velocity field or the pressure field in each discretized cell of the computation domain.By substituting the weighted density and the weighted viscosity (Equations ( 7) and ( 8), respectively) into the N-S equations (Equations ( 1)-( 2)), the velocity and pressure fields can be obtained.Regarding the two-phase case, if a specific grid cell is completely filled with one phase, volume fraction f is assigned as unity for that phase and zero for the other phase.Therefore, for two-phase flow (i = 2 in Equations ( 7)-( 8)), Equations ( 7) and ( 8) are simplified to be Equation ( 9) and (10), respectively.
where u is the velocity vector shared by all phases, and the deformation tensor D is defined as D mn ≡ (∂ m u n + ∂ n u m )/2, ρ and μ are the weighted density and the weighted dynamic viscosity of various phases, respectively.ρ i and μ i (i = 1, 2 for the two-phase case) are the density and the dynamic viscosity of phase i, respectively.Additionally, p stands for the pressure, a is the acceleration vector considering both the gravity and the external forces except for surface tension, and F s represents the capillary force caused by surface tension.f i is the volume fraction of phase i (only f is adequate for the two-phase flow).
During the computation, volume fraction f and pressure p are centre staggered in each grid cell.By contrast, velocity u, viscosity μ, and acceleration term are face staggered.The detailed temporal discretization of the N-S equations is realized by BASILISK, in which the momentum equation in vector form is discretized with a second-order time scheme using an intermediate velocity u (Equation ( 11)) (Popinet, 2009).The velocity at the next time step can be expressed as Equation ( 12), where the pressure gradient needs to be calculated in advance by the Poisson equation as shown in Equation ( 13) (Popinet, 2009).Here, Equation ( 13) is solved by the quad/octree-based multilevel solver (Popinet, 2003).
where σ is the surface tension coefficient; κ and n are the curvature and the normal vector of the interface, respectively; δ s is the Dirac distribution function related to the interface.Comparing Equation (1) with Equation (11), the capillary force could be written as: This surface tension force can be transformed into the specific control volume within the discretized domain as a volumetric force (Dupont & Legendre, 2010): To get the exact value of u , Equation ( 11) is reorganized and can be expressed as: As can be seen, the right-hand side of Equation ( 16) only relies on the values at time steps n and n + 1/2.This Helmholtz-type equation (Equation ( 16)) could be solved by an alternative form of the multilevel Poisson solver.The viscous term is discretized by the Crank-Nicholson method which is second-order convergence and unconditionally stable.Also, the Bell-Colella-Glaz second-order unsplit upwind scheme (Bell et al., 1989;Popinet, 2003) is utilized to discretize the advection term . When the CFL number (a scale to assure that the time step size is within a reasonable range) is less than one, the Bell-Colella-Glaz scheme is relatively stable.

Initial construction of the droplet-gas interface by level set method
In addition, the level set method (Osher & Sethian, 1988) is utilized to initially construct the interface between the liquid droplet and the ambient gas so that the periphery of the impinging droplet can be accurately built at the beginning of the simulation.Rather than by either explicit parameterization or direct specification, the level set method could build the profile of an object in a fast and reliable fashion and obtain the topology representation with excellent adaptation (Wang et al., 2003).According to the level set method, a scalar function of an iso-surface in implicit form is used to get the object boundary: where scalar function is the level set function, i.e. (x) : R n → R, y is the iso-value to define the specific iso-surface, and x is the set of points that comprise the iso-surface associated with the iso-value y.In this way, the set of vertices of the initial grid cells are applied to derive the values of the level set function so that the volume fraction of each cell is obtained to delineate the primary interface between the two phases.

Viscoelastic model for the micropillar-arrayed substrate
As mentioned in the introduction, there are three main models to describe the stress-strain relationships for viscoelastic materials, namely the Kelvin-Voigt model (Figure 1(a)), the Maxwell model (Figure 1(b)), and the generalized Maxwell model (Figure 1(c)).For the first two models, the dashpot (representing viscosity) and the spring (standing for elasticity) are connected in parallel and in series, respectively.Nonetheless, the generalized Maxwell model is composed of one spring and a load of Maxwell elements (a dashpot and a spring connected in series) in the parallel configuration.
Compared to the other two models, the generalized Maxwell model gives rise to better description of the stress-strain relationship by taking both relaxation and creep into account.That said, the Kelvin-Voigt model fails in interpreting stress relaxation and the Maxwell model falls short of accounting for creep response of the viscoelastic materials.However, the higher accuracy of the generalized Maxwell model usually incurs a prohibitive computational cost, which can be optimized by simplifying the generalized Maxwell model to a fiveparameter model, i.e. two Maxwell elements in parallel with an additional spring which stands for the long-term elastic effect.
Regarding the viscoelastic micropillar-arrayed substrate whose elastic and viscous properties are different from those of the intact bulk of the viscoelastic material, we need to build a viscoelastic model for the micropillared configuration.According to Nguyen et al. (2017), the temporal response of a non-aging viscoelastic material to external stress could be approximated by the elastic temporal response of an equivalent elastic material in the LC space.After this LC transform, the viscoelastic material can be modelled by the comparable elastic material whose stress-strain relationship is expressed in a linear tensor format: σ s * = C(q) : ε * , where C(q) is the stiffness tensor in the LC space, * represents the LC transform, and q denotes the LC variable.Consequently, for the generalized Maxwell model in the LC space, the bulk and shear moduli of the intact viscoelastic material, k * and μ * , could be given by: where k sub and μ sub (sub = i or ∞) indicate the bulk and shear moduli of either the specific spring in one certain Maxwell element or the additional spring.η s i and η d i are the bulk and shear viscosities of the corresponding Maxwell element, respectively.
In accordance with Nguyen et al.'s elastic response study of micro-cracked materials in the LC space (Nguyen et al., 2017), all the parameters in Equations ( 18) and ( 19) should be modified for cracked (i.e.micropillararrayed in this study) materials considering the effects of crack density = Na, which is also deemed as gap density for the micropillar-arrayed substrate in this paper.Here, a is the interstitial gap (or gap length for shapes other than circle) for the micropillared surface and N is the number of gaps in a unit volume of the micropillared substrate.Thus, the above two equations (Equations ( 18) and ( 19)) are modified correspondingly to obtain the bulk and shear moduli for the micropillared configuration, i.e. k hom * and μ hom * : Also, it should be noted that the modified moduli in Equations ( 20) and ( 21) are identical to those of the media with randomly distributed cracks (i.e.gaps for the micropillar-arrayed substrate) as predicted by the classical Mori-Tanaka approach (Mori & Tanaka, 1973).Therefore, these two governing equations can be effectively applied to the micropillared case regardless of the shape and distribution of the micropillar gaps when employing the five-parameter generalized Maxwell model in the LC space.As for the gas entrapped in the gaps between the micropillars, its cushion effect between the droplet and substrate is neglected due to the relatively low gap density in this study (all the cases are listed in Table 1).In our simulated cases, the equivalent shear modulus is not dominant because the droplet impacts Table 1.Effectiveness j of the bulk moduli and viscosities of the micropillar-arrayed substrates compared to the corresponding parameters of the intact material with varying gap density (Nguyen et al., 2017).along the normal direction where the shear deflection of the micropillars plays a relatively minute role in the response.Therefore, only the bulk strain caused by the bulk modulus (Equation ( 20)) is evaluated in our analysis, resulting in the bulk deformation velocity of the substrate.This deformation velocity is then applied as a Dirichlet boundary condition to the fluid fields of the liquid droplet and the ambient gas.Furthermore, the bulk strain for 2D cases in the LC space could be calculated by: where σ s0 and ε 0 (q) represent the stress and strain in the normal direction under a specific loading, respectively.However, instead of directly using the obtained strain in the LC space to calculate the boundary velocity, the inverse LC transform should be performed in advance so that the equivalent strain in the real space could be implemented.Due to the analytical complexity of this inverse transform, the direct inversion method proposed by Schapery (1962) is adopted by us to get the approximate value of the strain as: According to Schapery (1962), the direct inversion method has a great adaptivity when the logarithmic time doesn't change quickly as in our case.Therefore, the strain in the real space can be approximately deduced after this direct inversion operation.

Schematic representation of the micropillar-arrayed geometry
The schematic of the droplet impacting process on a micropillar-arrayed substrate is shown in Figure 2. The 2D simulation domain is a 5 unit lengths × 5 unit lengths (unit length is the basic length scale in BASILISK) region meshed with 512 × 512 discretized points.The height h of the micropillars is set at 100 × 10 −6 unit length, and the diameter D of the droplet is 1 unit length.The computational process of droplet impact by BASILISK starts at the droplet's falling moment from its initial position, during which the progress is characterized by time t.As stated above, the micropillar matrix configuration is characterized by gap density , and only the strain in the normal direction is considered in this study.Due to the axial-symmetry of the spherical liquid droplet, only half of the droplet is modelled to get the representative scenario of the impact process.

A simplified five-parameter model for the generalized Maxwell model
In Equations ( 20) and ( 21), the bulk and shear moduli for the micropillared material in the LC space have been established.Nonetheless, in order to reduce the computational cost while achieving a sufficient accuracy, a five-parameter model is adopted by us to represent the elastic response in the LC space as shown in Figure 3.In this model, two parallel Maxwell elements are included and the shear moduli and shear viscosities are omitted as stated above.Moreover, according to Nguyen et al.'s validation study of the creep case (Nguyen et al., 2017), the five-parameter model has an excellent agreement with the analytical solution of the governing equations in the LC space in terms of the temporal strain under a constant stress load.
Meanwhile, instead of getting the exact values of k sub ( ) (sub = i or ∞) and η s i ( ) in Equation ( 20) by the analytical polynomial suggested by Nguyen et al. (2017), the approximate effectiveness j of the viscoelasticity   20) with i = 2.The values of effectiveness j for the bulk moduli and viscosities of the micropillar-arrayed substrate with varying gap density are given in Table 1, which are used in the subsequent simulations.Also, the bulk moduli and viscosities of the adopted intact material are listed in Table 2, which are modified to become relatively softer compared to the primary intact material (Nguyen et al., 2017).

Algorithm chart for the computation process
The algorithm of the whole computational procedure is listed in Table 3, including the N-S solver and the VOF solver adopted by BASILISK.From the ad hoc point of view, the deformation velocity of the micropillared substrate, which is acquired by averaging the strain rates at all the discretized points on the boundary, functions as a Dirichlet boundary condition to the N-S equations of the fluid fields.Additionally, the corrected pressure, calculated by the pressure at the current step and the pressure at the previous step, is applied to replace the analogous stress, i.e. σ s0 in Equation ( 22), in the LC space.This correction approach could avoid the race condition (i.e. a disorder caused by the sequence of computation steps).For this case, the race condition is that the substrate deformation velocity needs to be computed while being used as the initial boundary condition at the same time.The flow chart of the computation process is showcased in Figure 4.
Table 3. Algorithm of the N-S equation and VOF solvers.

Algorithm 1 Computation for cell-centred velocity u and cell-centred pressure p and volume fraction f
The time step size is controlled by CFL number, which is 0.8 by default.Part A -solve for u and p: 1.For Equation ( 16): in order to get u , the face velocity u n+1/2 f at t + t/2 needs to be predicted first; 2. Divergence-free face velocity u n+1/2 f at t + t/2 can be obtained after successive iterations, using the Poisson-Helmholtz equation solver similar to Equation (13), i.e.: (ρ is constant for incompressible flow, the default tolerance for this solver is 10 −3 which must be reached within the specified maximum iterations); 3. Utilize this divergence-free face velocity u n+1/2 f to calculate the advection term u n+ 1 2 • ∇u n+ 1 2 via the standard Bell-Collela-Glaz advection scheme, which is a second-order unsplit upwind scheme; 4. Before leveraging the implicit viscosity solver, which is the same multilevel Poisson solver in step 2, to gain the viscous term ∇ the acceleration and pressure gradient effects can be added to u n by a correction function where a combined term c consisting of cell-centred acceleration a and pressure gradient ∇p at the previous time step is adopted; 5. Remove c dt from the previously modified u n after step 3 and step 4.Then, the provisional face velocity u n+1 f at t + t is obtained by adding the face velocity u n f interpolated from the centred velocity u n and the product of the acceleration term a and t; 6. Make one more projection to obtain the final face velocity u n+1 f at t + t along with the cell-centred pressure p, which could be used to calculate c at the current time step; 7. Use correction function mentioned in step 4 to add the effect of the combined acceleration and pressure gradient term, c dt, at current step to u n , then u n+1 is derived; 8. Face velocity u n+1 f obtained in step 6 could be applied to predict the next half step's face velocity as step 1; Part B -solve for f : 9.For the following equation similar to Equation (4): from step 6 could be employed in the advection solver of volume fraction; 10. f n+1 i is acquired through time integration of the discretized equation resulting from the equation in step 9, mainly considering the advection term.

Validation through modeling droplet impact on rigid and soft surfaces
Prior to investigating droplet impact dynamics on a viscoelastic surface, we need to validate the simulation protocols and algorithms developed by us (Figure 4) as well as the mathematical models built for both the fluid fields and the soft substrate.Since droplet impact on viscoelastic structures has been less explored, we first simulated droplet impact on a rigid surface and compared the BASILISK simulation results with the available experimental data for the validation purpose.As such, the stress-strain relationship σ s = Eε, where σ s and ε are the stress and the strain in the bulk direction, respectively, E is the bulk elastic modulus, was used in the case of the   1 of Liu et al. (2010).See Movie S1 in the supplemental materials.relatively rigid substrate to directly acquire the substrate deformation velocity.In specific, we compared the experimentally captured images provided by Liu et al. (2010) regarding the impact process of an ethanol droplet on a Plexiglas surface with our simulation contours of volume fraction f under the same conditions.The Weber number We was set as 1870, and the physical properties of ethanol such as density, viscosity, and surface tension are given in Table 1 of Liu et al. (2010).In this validation simulation, the bulk elastic modulus of the Plexiglas surface E = 2.95 GPa and its thickness is 2 mm.As shown in Figure 5, our simulation results can capture the main features of the impact process and are consistent with the experimental images at each frame.Moreover, the nondimensionalized radial trajectory r t /a r as a function of the dimensionless time t e /τ for the ethanol droplet impact process is showcased in Figure 6, where r t and a r are the travelling distance of the outmost rim and the droplet radius, respectively; t e is the elapsed time after the initial impingement, τ = a r /U i is a characteristic falling time, and U i is the impact velocity.As can be seen in Figure 6, the evolution of the radial trajectory obtained by our modelling has a good agreement with both the simulation results of Wu and Cao (2017) and the experimental data provided by Stow and Hadfield (1981).
To make further validation, the impact dynamics of ethanol and water droplets over soft PDMS silicone gels with various stiffnesses have been simulated by us, which are then compared with the experimental results obtained by Basso and Bostwick (2020).The spreading factor α = D max /D (where D max is the maximum spreading diameter) is selected as the parameter for Figure 6.Nondimensionalized radial trajectory r t /a r of the ejected sheet against dimensionless time t e /τ for the ethanol droplet impact simulation shown in Figure 5.The experiment data are from Stow and Hadfield (1981), and the referred simulation results of ethanol droplet impact were obtained by Wu and Cao (2017).The lines are guides to the eye.
comparison, which demarcates the maximum spreading capability of the generated liquid film.As can be seen in Figure 7(a), the simulation results are generally consistent with the experimental data.For the case of water droplet, the maximum difference is only 1.72% at E = 47.4 kPa compared to the corresponding experimental result.Furthermore, the contour of the maximum spreading diameter for ethanol droplet impact on different soft substrates with We = 270 is illustrated in Figure 7(b).
Hence, our simulation results are consistent with the experimental data on both the rigid substrates and the soft surfaces, indicating that the mathematical model of the entire droplet impact process implemented in BASILISK is reliable and robust.Consequently, the general protocols developed by us can be applied in the subsequent simulations of the micropillared cases.

Effect of viscoelasticity of the micropillar-arrayed substrate on droplet impact
By directly leveraging the strain-stress relationship as discussed in Section 4, droplet impingement on either the stiff substrate or the soft surface has been simulated by us for validation purpose.However, whether there exists a critical point for viscoelasticity of the viscoelastic micropillar-arrayed substrate, at which its deformation behaviour is analogous to that of the stiff solid substrate, deserves further investigation.In this respect, the effects of viscoelasticity of the micropillar-arrayed substrate on the surface deformation need to be investigated.
Because the Plexiglas surface is treated as a stiff substrate in the droplet impingement modelling of Section 4, we performed an asymptotic analysis to find out at what viscoelasticity the deformation of the micropillar-arrayed substrate can be approximated to that of the Plexiglas surface.We chose the primary viscoelasticity of the bulk material from Nguyen et al. (2017) with only the viscosity terms modified, as listed in Table 4.The effectiveness   1 are applied to form the basic micropillar-arrayed substrate.Moreover, to represent various micropillar-arrayed substrates, each of the elasticity terms (i.e.k 1 ( ), k 2 ( ) and k ∞ ( ), with = 0.20) of the basic micropillared substrate is multiplied by a pre-factor m, which ranges from 0.005 to 0.05.For the Plexiglas surface, its elastic modulus is 2.95 GPa, and its thickness is equal to the height of micropillars, i.e. 100 × 10 −6 unit length.
In addition, the same conditions used in the Plexiglas surface case of Section 4 are employed to simulate ethanol droplet impact on the micropillar-arrayed substrates, i.e. droplet diameter D = 1.7 mm, We = 1870 and ambient pressure P a = 1 atm.By comparing the dimensionless deformation magnitude of different substrates as shown in Figure 8, it turns out that the larger the pre-factor is, the smaller the deformation tends to be.When taking the deformation of the Plexiglas surface as a benchmark, with a critical pre-factor m = 0.00875, the viscoelastic response of the micropillar-arrayed substrate Figure 8. Dimensionless deformation displacement d avg /h of various substrates against time t with We = 1870 and ambient pressure P a = 1 atm for ethanol droplet impact simulations.The pre-factor m is applied to represent diverse micropillar-arrayed substrates.The droplet diameter D = 1.7 mm, the gap density of the micropillar-arrayed substrates = 0.20, and the physical properties of ethanol droplet are from Table 1 of Liu et al. (2010).During the impingement process, stages I and II indicate the going-down stage and the lift-up stage of the substrates, respectively.
to the impinging ethanol droplet can be deemed almost the same as the elastic behaviour of the stiff Plexiglas surface.Regarding the maximum dimensionless indentation depth, the difference is merely 0.76% between the micropillar-arrayed substrate with m = 0.00875 and the Plexiglas surface.

Results and discussion
In general, impinging droplet over the solid surface starts spreading upon the initial contact, which is followed by the potential splashing afterward.By considering various factors such as impact velocity U i , surface tension coefficient σ , and ambient gas pressure P a , we have conducted parametric simulations of droplet impact dynamics involving spreading and possible splashing processes on micropillar-arrayed viscoelastic substrates with gap density varying from 0.05 to 0.20 as listed in Table 1.Correspondingly, we compared the splashing magnitudes obtained from our simulations with the splashing extents predicted by the Kelvin-Helmholtz instability theory, which is a type of interfacial instability perturbed by small-amplitude variation originating from the substantially large velocity difference across the interface (Lee & Kim, 2015).
In essence, the Kelvin-Helmholtz instability is caused by the incompetence between the stabilizing effect of gravity and surface tension and the destabilizing effect of shear stress across the interface.Jepsen et al. (2006) found that the Kelvin-Helmholtz instability might serve as the main mechanism for the occurrence of splashing and finger formation during the spreading process of the impinging droplet.The dispersion relation of the Kelvin-Helmholtz instability could be expressed as: where ω is the interface growth rate, k is the wave number, and k max corresponds to the maximum point of ω (Equation ( 24)) which is achieved by setting the derivative of ω function as zero; U rel is the relative velocity between the gas phase and the liquid droplet phase; ρ g and ρ l are the densities of the ambient gas and liquid droplet, respectively.A series of droplet impact processes in the ambient environment have been simulated by us under a variety of impact velocities and gap densities for parametric study.In addition, while = 0.15 and U i = 2 m/s or 3 m/s, droplet impact dynamics with varying ambient gas pressure P a or surface tension coefficient σ have also been simulated in this study.

Grid demonstration and mesh independence test
Regarding numerical simulations of droplet impact on the micropillar-arrayed substrate, a schematic of the meshed domain is shown in Figure 9(a), which is a 512 × 512 grid network as mentioned in Section 3.1.Along with the computation by BASILISK, a self-mesh adaptation process is implemented to adaptively adjust the mesh density especially in the vicinity of the fluid-fluid interface.Due to the wide gap of the geometric dimensions between the droplet and the micropillars, which span several orders of magnitude from µm to mm, there is no physical mesh built for the micropillars in our modelling.Instead, via the five-parameter model for the micropillararrayed substrate in the LC space, each discretized node at the bottom boundary serves as a coupling position between the substrate and the main fluid fields by transmitting the corresponding deformation velocity.In this way, the inconsistency of geometric scale between the droplet and micropillars has been reconciled by means of the mechanically averaged responses.
Mesh independence test needs to be conducted before the subsequent parametric studies.Hence, we added two more meshes to observe and compare the dimensionless displacement of the micropillar-arrayed substrate, i.e. one is 256 × 256 mesh and the other is 1024 × 1024 mesh, in addition to the adopted 512 × 512 mesh.The averaged dimensionless deformation d avg /h of the micropillararrayed substrate against time t with impact velocity U i = 2 m/s, gap density = 0.15, ambient pressure P a = 1 atm, and surface tension σ = 1/50 N/m for various mesh densities is illustrated in Figure 9(b).As can be seen in 9(b), all the deformation results are quite uniform with negligible deviations during the entire impact process.The difference of the maximum indentation depth between 512 × 512 mesh and 1024 × 1024 mesh is only 0.55%, indicating that the 512 × 512 mesh is accurate enough for subsequent analyses.This independence test strongly confirms the high efficiency of the self-mesh adaptation method implemented in BASILISK.
In the spirit of mechanically averaged responses, the equivalent elastic properties of the entire micropillararrayed viscoelastic surface are governed by the fiveparameter generalized Maxwell model in the LC space, as stated in Section 2.3.We elucidate this averaging procedure by illustrating the dimensionless deformation d l /h of each local discretized point at the base boundary at a certain moment (for instance, t = 1.15 s) as shown in Figure 10 with impact velocity U i = 2 m/s, gap density = 0.15, ambient pressure P a = 1 atm, and surface tension σ = 1/50 N/m.As such, the averaged deformation at one moment is the arithmetic mean of all the local momentary deformations.It is noteworthy that point 1 in Figure 10 coincides with the droplet centre.And the discretized points ranging from node 60 to node 100 represent the rim of the spreading liquid laminate by droplet impingement.

Contours of spreading and splashing processes and droplet impact behaviours
Based on the parametric simulations of droplet impact processes, we observed the evolutions of both the velocity field and the volume fraction field as shown in Figures 11-15.Also, the movies for the contours of volume fraction f are showcased in the supplemental materials.For all the cases in Sections 6.1-6.5, as stated in Section 3.1, the height of micropillars h = 100 × 10 −6 unit length and the droplet diameter D = 1 unit length, corresponding to 0.1 µm and 1 mm, respectively, if the unit length is set to be 1 mm.While = 0.15, P a = 1 atm, and σ = 1/50 N/m, the contours of volume fraction f for the droplet impact processes with U i = 1 m/s, 2 m/s, 3 m/s, and 4 m/s, respectively, are shown in Figure 11.With increasing impact velocity, the splashing magnitude is enhanced and more satellite droplets are generated.Additionally, with = 0.15, σ = 1/50 N/m, and U i = 3 m/s, the contours of volume fraction f and vertical velocity U v are showcased in Figures 12 and 13, respectively, when splashing occurs with varying ambient pressure P a = 2 atm, 3 atm, and 5 atm.As can be seen in Figure 12, the splash angle changes from 92.65 • to 114.44 • as P a increases from 2 atm to 5 atm.In general, more satellite droplets are ejected with increasing ambient pressure.Meanwhile, with a higher ambient pressure, each of the generated satellite droplets tends to become smaller.Moreover, while = 0.15, P a = 1 atm, and U i = 3 m/s, the splash processes with diverse surface tension coefficients, i.e. σ = 1/25 N/m, 1/100 N/m, and 1/200 N/m, are displayed in Figures 14 and 15, respectively.As displayed in Figure 14, with σ decreased from 1/25 N/m to 1/200 N/m, the splash angle increases from 58.76 • to 84.96 • , and more satellite droplets are ejected at the rear end of the splash tip with a tendency to break into even smaller ones.Thus, as opposed to the declining trend of splash intensity with increasing surface tension coefficient σ , the splash extent gets enhanced with either increasing impact velocity U i or rising ambient pressure P a .

Substrate deformation during droplet impact with varying impact velocity
Since the viscoelastic properties of the micropillared substrate are substantially determined by gap density in the five-parameter viscoelastic model, we investigated its equivalent elastic response under a certain impact velocity U i with varying .To study the effects of impact velocity U i on the droplet impinging process, we conducted simulations by setting U i equal to 1 m/s, 2 m/s, 3 m/s, and 4 m/s (corresponding to We = 50, 200, 450, and 800), respectively, for each gap density , while P a = 1 atm and σ = 1/50 N/m.
Instead of using each individual micropillar's indentation depth, we adopted the average indentation d avg across all the discretized mesh points in our analysis.Then, the average dimensionless indentation depth d avg /h of the micropillared substrate against time t with changing gap density for each impact velocity U i is showcased in Figure 16.As stated in Section 3.1, the initial height of micropillars h = 100 × 10 −6 unit length.As can be seen in Figure 16, the deformation response of the micropillared substrate to impact can be divided into two distinct stages, i.e. the going-down stage I (including the pre-impact stage) and the lift-up stage II (containing the post-impact stage).In addition, the maximum indentation depth is enlarged with increasing impact velocity U i for each gap density .For instance, for = 0.15 and U i = 2 m/s, the maximum dimensionless indentation depth reaches 2.565 × 10 −8 at t = 1.193 s.In comparison, the maximum indentation depth rises to 6.743 × 10 −8 at t = 0.88 s in the case of = 0.15 and U i = 3 m/s.Also, the duration of the lift-up stage II is shortened with increasing U i for each .Therefore, there exists a larger displacement variation within a shorter period for the entire lift-up stage under the higher impact velocity.Consequently, with increasing U i , a more intensified splash would be incurred, and more satellite droplets get ejected from the mother droplet, as shown in Figure 11.
Furthermore, the averaged dimensionless deformation velocity d avg /|U i | (here d is the deformation velocity for a discretized cell, i.e. the derivative of d with respect to the time step size t) of the micropillared substrate against time t with varying gap density for each impact velocity U i is displayed in Figure 17.It is noteworthy that the negative sign of d avg /|U i | just indicates the going-down phase (indentation) of the micropillars.As can be seen in Figure 17, with increasing impact velocity U i , both the maximum indentation velocity and the maximum lift-up velocity are enhanced.In specific, for = 0.15 and U i = 3 m/s, the maximum indentation velocity reaches 1.500 × 10 −10 at t = 0.867 s and the maximum lift-up velocity gets to 2.982 × 10 −11 at t = 0.902 s.In comparison, they increase to 2.794 × 10 −10 at t = 0.738 s and 4.841 × 10 −11 at t = 0.765 s, respectively, for the case of = 0.15 and U i = 4 m/s.Meanwhile, the duration of either the going-down stage or the lift-up stage (marked as regions II and III in Figure 17, respectively) of the micropillared array is dramatically shortened as the impact velocity rises for each gap density.This trend indicates that in the case of higher impact velocity U i , the kinetic energy transferred from the droplet into the micropillars during the impingement is at least partially converted back to the droplet within a much shorter period throughout the lift-up process.In other words, an      even more drastic variation of deformation velocity of micropillars ensues as U i increases, which would facilitate the splash occurrence with an enhanced intensity, as reflected in Figure 11.

Substrate deformation during droplet impact with distinct ambient pressure
Since the ambient gas density ρ g also serves as a dominant factor in determining the maximum wave number k max (Equation ( 25)), we simulated splash circumstances under gap density = 0.15 and impact velocity U i = 2 m/s and 3 m/s, respectively, with varying ambient pressure P a .The dimensionless indentation depth d avg /h of the micropillared substrate against time t under two distinct impact velocities with increasing ambient pressure P a is displayed in Figure 18(a,b), respectively.Furthermore, the dimensionless deformation velocity d avg /|U i | of the micropillared substrate against time t under the same conditions is shown in Figure 18(c,d), respectively.
With increasing ambient pressure, a similar decreasing trend of indentation depth exhibits for both U i = 2 m/s and 3 m/s, as displayed in Figure 18(a,b), respectively.In particular, for U i = 3 m/s, the maximum indentation depth reduces from 7.058 × 10 −8 to 4.628 × 10 −8 as the ambient pressure P a is increased from 3 atm to 5 atm .Moreover, the displacement regions labelled with the red lines in Figure 18(a,b) correspond to the velocity ranges marked with the red solid lines in Figure 18(c,d), respectively.A close look at Figure 18(a,b) reveals that the slope of the displacement range demarcated by red line is not consistent with the slope of its prior continuous range of the lift-up phase.The slope inconsistency of the red line may be due to an intense velocity instability, which would facilitate a more intensified splashing.In (d), the ejection and rejoining behaviours of satellite droplets, which induce the drastic oscillations during the lift-up stage as P a = 2 atm and 3 atm, can be seen in Movies S6 and S7 of the supplemental materials, respectively.
As shown in Figure 18(c,d), both the maximum indentation velocity and the maximum lift-up velocity decrease with increasing ambient pressure P a .With P a increased from 2 atm to 3 atm while the impact velocity U i = 2 m/s, ignoring some minor oscillations, the maximum indentation velocity reduces from 4.824 × 10 −11 to 3.651 × 10 −11 and the maximum lift-up velocity declines from 1.783 × 10 −11 to 1.604 × 10 −11 .Consequently, a lower P a results in a velocity profile of larger values around the main splash tip as illustrated in Figure 13(a).However, this observation doesn't mean that the splashing under a lower P a occurs in a more intensified fashion, i.e. multiple satellite droplets ejected and a bigger splash angle formed (the angle gauged between the main splash tip and the horizontal substrate surface), because the largest splash angle of 114.44 • appears at P a = 5 atm as shown in Figure 12(c).On the contrary, the splash magnitude is mainly determined by the instability of the dimensionless deformation velocity d avg /|U i |.In the cases of U i = 2 m/s and 3 m/s while P a = 5 atm, there exists a consecutive range whose slope is not consistent with its prior range within the lift-up phase, which is marked by the red solid lines as shown in Figure 18(c,d) (corresponding to the red lines in Figure 18(a,b)), respectively.This slope inconsistency is indicative of a substantial instability at the interface between the liquid droplet and the ambient gas, which would induce a more intense splashing with satellite droplets generated and a larger splash angle formed, as displayed in Figure 12(c).

Substrate deformation during droplet impact with different surface tension coefficients
As indicated by Equation ( 25), surface tension coefficient σ has a significant influence on the maximum wave number k max as well.Therefore, we simulated and investigated the splashing behaviours with a variety of surface tension coefficients while gap density = 0.15, ambient pressure P a = 1 atm, and droplet impact velocity U i = 2 m/s and 3 m/s, respectively.With a varying surface tension coefficient σ , the dimensionless deformation displacement d avg /h and the dimensionless deformation velocity d avg /|U i | against time t for impact velocity U i = 2 m/s and 3 m/s are depicted in Figure 19(a-d), respectively.
As can be observed in Figure 19(a,b), only a minor increment of the maximum indentation depth exhibits with decreasing surface tension coefficient for both impact velocities, i.e.U i = 2 m/s and 3 m/s.Specifically, with U i = 3 m/s, the maximum indentation depths with surface tension coefficient σ = 1/100 N/m and 1/200 N/m are 6.758 × 10 −8 and 6.768 × 10 −8 , respectively.Although there is a rather small increment in the maximum lift-up height with decreasing surface tension σ for each impact velocity, the maximum lift-up heights under the three different surface tension coefficients can be deemed almost the same for both U i = 2 m/s and 3 m/s.Therefore, the variations of both the maximum indentation depth and the maximum lift-up height under distinct values of σ are relatively minute, indicating that the variation of surface tension σ is not the dominant factor determining the splash extent.
Not only the maximum indentation velocity but also the maximum lift-up velocity does not alter too much when surface tension coefficient σ is decreased from 1/25 N/m to 1/200 N/m for impact velocity U i = 2 m/s and 3 m/s, as shown in Figure 19(c,d), respectively.Therefore, surface tension has a minor influence on the deformation velocity of the micropillared substrate.For instance, under U i = 2 m/s, the maximum indentation velocity has a minor variation between 5.782 × 10 −11 and 5.740 × 10 −11 and the maximum lift-up velocity slightly varies between 1.385 × 10 −11 and 1.370 × 10 −11 as σ changes from 1/100 N/m to 1/200 N/m.Nonetheless, the velocity instability, which is reflected by a consecutive range whose slope of its dominant trend, i.e. without considering intermittent peaks pointing up and down, is not consistent with its prior range during the lift-up phase as demarcated by the red solid lines in Figure 19(c,d), is strengthened with decreasing surface tension coefficient.This strengthening trend is evidenced by the fact that the inconsistent range of slope generally deviates further away from its prior range under the lower surface tension condition.Because the strengthened velocity instability in the case of a lower surface tension may lead to intensified splash, a larger extent of splash consequently ensues in this situation with a bigger splash angle formed and more accompanying satellite droplets would be generated, as displayed in Figure 14.However, as opposed to the aforementioned slope inconsistency of the dominant and consecutive trend, there exist some intermittent oscillations or peaks of the dimensionless velocity d avg /|U i | during the lift-up stage for σ = 1/100 N/m under U i = 2 m/s and for σ = 1/25 N/m under U i = 3 m/s, as shown in Figure 19(c,d), respectively.It is noteworthy that these intense intermittent oscillations themselves do not stand for the entire velocity instability, which is signified by the slope inconsistency as discussed above.Instead of resulting in a more intense splash, these intermittent oscillations would just induce a large velocity divergence at the main splash tip, as displayed in Figure 15(a) for the case of U i = 3 m/s.

Droplet impacting with varying droplet diameter and micropillar height
In addition, the impact processes of a droplet of 2 unit lengths in diameter with impact velocity U i = 1 m/s and 2 m/s, respectively, for each gap density listed in Table 1 were simulated, as shown in Figure 20.Here, the ambient pressure P a and surface tension σ are maintained at the same levels as in Section 6.3, i.e.P a = 1 atm and σ = 1/50 N/m.Also, the height of the micropillars is increased to 3000 × 10 −6 unit length.In the realistic scale, with the unit length being millimeter, the droplet diameter D = 2 mm and the initial micropillar height h = 3 µm.
For the cases with an enlarged droplet diameter (D = 2 mm) and an increased micropillar height (h = 3 µm), the averaged dimensionless displacement d avg /h of the micropillared substrate against time t with changing gap density for impact velocity U i = 1 m/s and 2 m/s is showcased in Figure 20(a,b), respectively.Compared to the displacement variation in the case of the smaller droplet diameter D = 1 mm and shorter micropillar height h = 0.1 µm (as shown in Figure 16(a,b)), the trend of the dimensionless displacement for the corresponding case of D = 2 mm and h = 3 µm acts very similarly under the same impact velocity U i and gap density , as shown in Figure 20(a,b).However, an increment of the maximum indentation depth under the same U i and becomes prominent in the cases of the larger droplet diameter and taller micropillar height.And the averaged dimensionless deformation velocity d avg /|U i | of the micropillared substrate against time t with varying gap density for impact velocity U i = 1 m/s and 2 m/s is displayed in Figure 20(c,d), respectively, while D = 2 mm and h = 3 µm.Regarding each gap density, there is much similarity between Figure 20(c,d) and Figure 17(a,b) (D = 1 mm and h = 0.1 µm) for the same impact velocity U i .Nonetheless, due to the larger impact pressure caused by an enlarged droplet size, the dimensionless deformation velocity in the case of D = 2 mm and h = 3 µm generally increases under the same U i and .In comparison with prior cases, more violent oscillations of d avg /h and d avg /|U i | have emerged in the cases of a larger droplet size and taller micropillars, resulting from the stronger impact pressure associated with the bigger impinging droplet.

Splash magnitude compared with the Kelvin-Helmholtz instability
After simulating the splashing process and analysing the tendency of splash magnitude under the effects of impact velocity U i , ambient pressure P a , and surface tension coefficient σ , we conducted a quantitative estimation of flow instability, which may result in potential splash, via the Kelvin-Helmholtz instability theory (Equations ( 24)-( 25)) to determine how intense the splash would be.By substituting the values of impact velocity U i (approximately treated as U rel in Equations ( 24)-( 25)), gas density ρ g , ambient pressure P a , and surface tension coefficient σ used in one certain simulation into Equations ( 24) and ( 25), the splash magnitude of the corresponding case can be quantified by the interface growth rate ω, which owns a maximum point at k max .As a result, for the various cases with = 0.15, the maximum interface growth rate ω k=k max against wave number k max is shown in Figure 21.It can be seen that U i plays an even more dominant role in determining ω k=k max than either P a or σ , which is due to the higher order (the second order) of the velocity term in Equation ( 25).
The trend of the interface growth rate ω k=k max shown in Figure 21 is consistent with the splash circumstance displayed in Figures 12 and 14, which are the contours of volume fraction f with varying P a and σ , respectively, as = 0.15 and U i = 3 m/s.With P a in the range of 2 atm to 5 atm, the magnitudes of ω k=k max for points b-d shown in Figure 21 conform to the splash intensity results shown in Figure 12.In specific, a larger splash angle of 114.44 • and more satellite droplets are present at P a = 5 atm in Figure 12(c), which is consistent with the higher ω k=k max of point d as shown in Figure 21.In addition, with σ decreased from 1/25 N/m to 1/200 N/m, the values of ω k=k max for points e-g displayed in Figure 21 match the splash extent contours represented in Figure 14.For instance, a larger splash angle of 84.96 • and more satellite droplets of smaller volume become prominent at σ = 1/200 N/m in Figure 14(c), which agrees with the higher ω k=k max of point g as represented in Figure 21.Consequently, the contour results are in an excellent agreement with the K-H instability theory, indicating the conducted simulations of droplet impact on the micropillared substrate still abide by the interface growth mechanism.

Conclusions
In this work, we performed the numerical simulations of liquid droplet impact dynamics on micropillar-arrayed viscoelastic surfaces using BASILISK.As such, a level set method is used to initially construct the profile of the droplet.As to the mathematical model of the micropillararrayed viscoelastic substrate, a five-parameter generalized Maxwell model in the Laplace-Carson (LC) space is adopted by us to obtain the equivalent elastic response related to its viscoelastic properties.Subsequently, the bulk strain in the real space, which is used for calculating the substrate deformation velocity, is acquired by the inverse LC transform.The deformation of the micropillar array in the horizontal shear direction is omitted in this study due to its secondary influence on the droplet's normal impingement.Eventually, this bulk deformation velocity at each time step is treated as a Dirichlet boundary condition to the flow fields of both the liquid droplet and the ambient gas.
By tuning a variety of factors such as impact velocity U i , ambient pressure P a , and surface tension coefficient σ in the parametric study, the mechanisms of how these factors affect the splash extent on the micropillar-arrayed surface are revealed.For instance, with increasing impact velocity U i , both the maximum indentation velocity and the maximum lift-up velocity increase under each gap density .In addition, both the going-down period and the lift-up period of micropillars are dramatically shortened with impact velocity U i enhanced.Also, a comprehensive view of the fluid fields has been obtained by virtue of contours of various parameters.These parametric studies shed light on how the viscoelasticity of the micropillar-arrayed substrate influences the droplet's spreading and splashing behaviours.
Furthermore, the splash magnitude has been quantitatively analysed by the Kelvin-Helmholtz instability theory for different cases, which reaches a good agreement with our simulation results regarding the splash extent under a wide range of impact velocity U i , ambient pressure P a , and surface tension coefficient σ .With increasing U i , the splash is intensified and there are more satellite droplets generated.This increasing trend of splash also holds for the cases with rising P a , which is verified by the experimental results (Xu et al., 2005) that the splashing would be strengthened by the increased ambient pressure.On the contrary, the declining σ can enhance splashing possibility and magnitude.And these splash tendencies are consistent with the experimental observations of droplet impact on a rigid solid surface conducted by Liu et al. (2010).Hence, our numerical simulation methodology captures the main features of the droplet impact dynamics on microstructured viscoelastic surfaces and underscores the importance of the LC transform in studying fluid-viscoelastic solid interactions by virtue of the mechanically averaged responses to reconcile the geometric inconsistency of multiscale domain.
Regarding the air entrapped within the gaps between the micropillars, its cushion effect on the droplet impact dynamics is neglected due to the relatively low gap density in this study.However, with larger gap density, the cushion effect deserves further investigation in future work.

Figure 1 .
Figure 1.Schematic configurations of mathematical models for viscoelastic materials: (a) the Kelvin-Voigt model; (b) the Maxwell model; and (c) the generalized Maxwell model.

Figure 2 .
Figure 2. Droplet impact process on a micropillar-arrayed substrate (not in scale), the unit of the denoted lengths is one unit length of the grid in BASILISK.

Figure 3 .
Figure 3. Schematic configuration of the five-parameter model.

Figure 4 .
Figure 4. Flow chart of the simulation procedures realized in BASILISK.

Figure 5 .
Figure 5.Comparison between the experimentally captured images (top) provided byLiu et al. (2010) for the impact process of an ethanol droplet on a Plexiglas surface and the contours of volume fraction f (bottom) simulated by BASILISK under the same conditions, i.e. droplet diameter D = 1.7 mm, We = 1870 and ambient pressure P a = 1 atm.The physical properties of ethanol droplet used in the simulations are the same as Table1ofLiu et al. (2010).See Movie S1 in the supplemental materials.

Figure 7 .
Figure 7. Validation of simulations of droplet impact on soft silicone substrates with various elastic moduli, i.e. 15.1 kPa, 22.2 kPa and 47.4 kPa: (a) spreading factor α against Young's modulus E, where the experimental data of ethanol and water droplets with different We are from Basso and Bostwick (2020); (b) contour of the maximum spreading diameter when ethanol droplet impinges on silicone substrates with diverse stiffnesses while We = 270.The lines are guides to the eye.

Figure 9 .
Figure 9. Representative meshed domain and independence test of mesh density: (a) schematic of the 512 × 512 node mesh used in this study; (b) dimensionless deformation displacement d avg /h of the micropillar-arrayed substrate against time t with impact velocity U i = 2 m/s, gap density = 0.15, ambient pressure P a = 1 atm, and surface tension coefficient σ = 1/50 N/m for different meshes.

Figure 10 .
Figure 10.Dimensionless deformation d l /h at t = 1.15 s of each local discretized point at the bottom boundary with impact velocity U i = 2 m/s, gap density = 0.15, ambient pressure P a = 1 atm, and surface tension σ = 1/50 N/m.The droplet contact radius on the substrate approximately ranges from node 1, which is aligned with the right boundary shown in Figure 2, to node 100 at the corresponding moment.

Figure 11 .
Figure 11.Contours of volume fraction f during droplet impact processes with impact velocity U i increased from 1 m/s to 4 m/s while = 0.15, P a = 1 atm, and σ = 1/50 N/m: (a) U i = 1 m/s; (b) U i = 2 m/s; (c) U i = 3 m/s; and (d) U i = 4 m/s.See Movies S2-5 for the cases of (a)-(d) in the supplemental materials, respectively.

Figure 12 .
Figure 12.Contours of volume fraction f during droplet impact processes with ambient pressure P a increased from 2 atm to 5 atm while = 0.15, σ = 1/50 N/m, and U i = 3 m/s: (a) P a = 2 atm; (b) P a = 3 atm; and (c) P a = 5 atm.The inset in each figure represents the zoom-in view of the splash process for the corresponding case.See Movies S6-8 for the cases of (a)-(c) in the supplemental materials, respectively.

Figure 13 .
Figure 13.Contours of vertical velocity U v during droplet impact processes with ambient pressure P a increased from 2 atm to 5 atm while = 0.15, σ = 1/50 N/m, and U i = 3 m/s: (a) P a = 2 atm; (b) P a = 3 atm; and (c) P a = 5 atm.

Figure 14 .
Figure 14.Contours of volume fraction f during droplet impact processes with surface tension coefficient σ decreased from 1/25 N/m to 1/200 N/m while = 0.15, P a = 1 atm, and U i = 3 m/s: (a) σ = 1/25 N/m; (b) σ = 1/100 N/m; and (c) σ = 1/200 N/m.The inset in each figure represents the zoom-in view of the splash behaviour for the corresponding case.See Movies S9-11 for the cases of (a)-(c) in the supplemental materials, respectively.

Figure 16 .
Figure 16.Dimensionless deformation displacement d avg /h of the micropillared substrate against time t with changing gap density while P a = 1 atm and σ = 1/50 N/m for varying impact velocity U i : (a) U i = 1 m/s; (b) U i = 2 m/s; (c) U i = 3 m/s; and (d) U i = 4 m/s.The initial height of micropillars h = 100 × 10 −6 unit length.Two stages, i.e. stage I (going-down including pre-impact process) and stage II (lift-up containing post-impact process), denote the entire droplet impact process.

Figure 17 .
Figure 17.Dimensionless deformation velocity d avg /|U i | of the micropillared substrate against time t with changing gap density while P a = 1 atm and σ = 1/50 N/m for varying impact velocity U i : (a) U i = 1 m/s; (b) U i = 2 m/s; (c) U i = 3 m/s; and (d) U i = 4 m/s.The initial height of micropillars h = 100 × 10 −6 unit length.The entire impact process is composed of four stages, i.e.I (prior), II (goingdown), III (lift-up), and IV (post).

Figure 18 .
Figure 18.Dimensionless deformation displacement d avg /h and dimensionless deformation velocity d avg /|U i | of the micropillared substrate against time t with ambient pressure P a increased from 2 atm to 5 atm while σ = 1/50 N/m, gap density = 0.15, and U i = 2 m/s, 3 m/s, respectively: (a) d avg /h vs t, U i = 2 m/s; (b) d avg /h vs t, U i = 3 m/s; (c) d avg /|U i | vs t, U i = 2 m/s; and (d) d avg /|U i | vs t, U i = 3 m/s.The initial height of micropillars h = 100 × 10 −6 unit length.In (d), the ejection and rejoining behaviours of satellite droplets, which induce the drastic oscillations during the lift-up stage as P a = 2 atm and 3 atm, can be seen in Movies S6 and S7 of the supplemental materials, respectively.

Figure 19 .
Figure 19.Dimensionless deformation displacement d avg /h and dimensionless deformation velocity d avg /|U i | of the micropillared substrate against time t with surface tension coefficient σ decreased from 1/25 N/m to 1/200 N/m while P a = 1 atm, gap density = 0.15, and U i = 2 m/s, 3 m/s, respectively: (a) d avg /h vs t, U i = 2 m/s; (b) d avg /h vs t, U i = 3 m/s; (c) d avg /|U i | vs t, U i = 2 m/s; and (d) d avg /|U i |, U i = 3 m/s.The initial height of micropillars h = 100 × 10 −6 unit length.

Figure 20 .
Figure 20.Dimensionless deformation displacement d avg /h and dimensionless deformation velocity d avg /|U i | of the micropillared substrate against time t while P a = 1 atm, σ = 1/50 N/m, and U i = 1 m/s, 2 m/s, respectively, with varying gap density : (a) d avg /h vs t, U i = 1 m/s; (b) d avg /h vs t, U i = 2 m/s; (c) d avg /|U i | vs t, U i = 1 m/s; and (d) d avg /|U i | vs t, U i = 2 m/s.The initial height of micropillars h = 3000 × 10 −6 unit length.In (a) and (b), stages I and II are consistent with Figure 16.In (c) and (d), stages I-IV have the same meanings as in Figure 17.

Figure 21 .
Figure 21.Maximum interface growth rate ω k=k max against wave number k max for cases with gap density = 0.15.
η s

Table 2 .
Nguyen et al. (2017)osities of the intact material used in this study (based on the data byNguyen et al. (2017)and modified).Effective Viscoelastic Properties Figures provided by Nguyen et al.With a specific gap density , the bulk modulus or viscosity of the micropillared substrate is obtained by multiplying the corresponding modulus or viscosity of the intact substrate material by the correlated effectiveness value.Then, k hom * for the micropillared substrate is acquired by substituting thus-obtained moduli and viscosities of all the components (dashpots and springs) in the five-parameter viscoelastic model into Equation (

Table 4 .
Nguyen et al. (2017)osities of the intact material used for the asymptotic analysis (based on the data byNguyen et al. (2017)and with only viscosities modified by us).