Droplet impact onto a spring-supported plate: analysis and simulations
Michael J. Negus, Matthew R. Moore, James M. Oliver, Radu Cimpeanu
JJournal of Engineering Mathematics manuscript No. (will be inserted by the editor)
Droplet impact onto a spring-supported plate:analysis and simulations
Michael J. Negus · Matthew R. Moore · James M. Oliver · Radu Cimpeanu
Received: date / Accepted: date
Abstract
The high-speed impact of a droplet onto a flexible substrate is a highlynonlinear process of practical importance which poses formidable modelling chal-lenges in the context of fluid-structure interaction. We present two approachesaimed at investigating the canonical system of a droplet impacting onto a rigidplate supported by a spring and a dashpot: matched asymptotic expansions anddirect numerical simulation (DNS). In the former, we derive a generalisation ofinviscid Wagner theory to approximate the flow behaviour during the early stagesof the impact. In the latter, we perform detailed DNS designed to validate theanalytical framework, as well as provide insight into later times beyond the reachof the proposed mathematical model. Drawing from both methods, we observe thestrong influence that the mass of the plate, resistance of the dashpot and stiffnessof the spring have on the motion of the solid, which undergoes forced dampedoscillations. Furthermore, we examine how the plate motion affects the dynamicsof the droplet, predominantly through altering its internal hydrodynamic pressuredistribution. We build on the interplay between these techniques, demonstratingthat a hybrid approach leads to improved model and computational development,as well as result interpretation, across multiple length- and time-scales.
Keywords
Impact · Droplets · Interfacial flows · Asymptotic analysis · Directnumerical simulation · Fluid-structure interaction
Droplet impacts are a rich and ubiquitous phenomenon in both nature and indus-try, from inkjet printing [10] to pesticide spray deposition [31] and estimating the
Michael J. Negus · Matthew R. Moore · James M. Oliver · Radu CimpeanuMathematical Institute, University of Oxford, Radcliffe Observatory Quarter, Oxford OX26GG, UKE-mail: [email protected] CimpeanuMathematics Institute, Zeeman Building, University of Warwick, Coventry CV4 7AL, UKE-mail: [email protected] a r X i v : . [ phy s i c s . f l u - dyn ] S e p Michael J. Negus et al. early stages of oily aerosol dispersal in the atmosphere after large-scale spills [21].The dynamics of these processes are governed by the complex flow physics of thedroplet and the surrounding gas, as well as the properties of the substrate, suchas wettability [1] and roughness [6]. An additional layer of complexity is added tothe system if the substrate is deformable, meaning that the force of the impactorcauses the substrate to move or change shape. A common example of dropletimpact onto deformable substrates is that of rainfall onto leaves [8]. Previous the-oretical and experimental studies of this class of systems include droplet impactonto cantilever beams [7], silicone gels [13] and elastic membranes [25].Recent developments in experimental imaging techniques, improving in bothframe rate and spatial resolution, have reinvigorated investigative efforts in high-speed impact [14], revealing previously inaccessible features due to the small,rapidly developing regions upon impact, such as the impingement of micro-drops[38] and early azimuthal instabilities of the ejecta [17]. Similarly, the increase inhigh performance computing resources have allowed increasingly efficient directnumerical simulations (DNS) of these systems to be performed [4,19,34], whichare intensive due to the rapidly evolving interfaces, multi-scale flow features andhigh density and viscosity ratios present. These experimental and computationaldifficulties mean that rigorous mathematical modelling of such flows is of key im-portance to any comprehensive investigation, as analytical approximations to theflow can provide insight into the underlying physical processes, greatly enhancepredictive capabilities in regimes otherwise difficult to examine, and save computa-tional resources required for numerical simulations. In all cases, the deformabilityof the substrate makes studying these systems even more complex, from the diffi-culty in observing small deformations experimentally to the additional degrees offreedom required for numerical study.Studies of droplet impact usually focus on a specific timescale, as examiningthe rapidly evolving interfaces makes universal investigations challenging. Shortlybefore impact, the cushioning effect of the gas layer between the droplet and thesubstrate leads to high pressures which cause the bottom of the droplet to dimple.This interfacial deformation results in a gas bubble being entrapped inside thedroplet upon impact, which has been widely observed experimentally [35], as wellas reproduced in numerical studies [34], and modelled analytically [9]. At earlystages of the impact, close to the points of contact, the free surface rapidly turnsover and begins to spread across the substrate. During this timescale, instabilitiesin the free surface can cause micro-droplets to be rapidly expelled, which has beenobserved experimentally using high-speed photography [36]. Numerical schemeswith adaptive mesh refinement [28,29] allow for the computational study of thisearly timescale by concentrating resources in the small region close to the substrate[4,27,34]. Analytical approaches adapt Wagner theory [11,12,30,39], an inviscidfluid model for solid-liquid impact inspired by the study of aircraft landing onwater and ship slamming. Later in the impact process, once the droplet beginsto fully spread across the substrate, viscosity and surface tension, as well as the chemistry of the substrate, tend to play a more significant role. Various physicalparameters determine whether the droplet retracts, rebounds or splashes; see [14]for an extensive review on the experimental studies on late impact. Quantitiesof interest at this timescale include the minimum thickness and the maximumdiameter of the droplet as it spreads, and in particular, how these depend on the roplet impact onto a spring-supported plate: analysis and simulations 3 physical properties of the liquid and the substrate have been the focus of previoustheoretical studies [5,41].Experimental work focused on droplet impact onto elastic substrates is muchless common. Recent investigations have examined the impact of droplets onto theend of cantilever beams [7,40], drawing close parallels to the impact onto leaves.In both cases, the impact of the droplet excited oscillations at the end of thebeam, with characteristic time periods of the order of the late impact timescale.Different length beams were considered in order to show how stiffer beams resultedin oscillations with higher frequencies. Other studies concerning droplet impactonto elastic membranes [25] and silicone gels [13] show that the compliance of thesubstrate strongly affects the splashing threshold (regarded here as the minimumvelocity necessary to observe splashing), which is significantly increased when thestiffness of the substrate is reduced.Fluid-structure interaction problems are notoriously difficult to study numeri-cally and allowing the substrate to deform only exacerbates this. If the substrateonly exhibits translational motion (without accounting for bending), then theproblem can be simplified by considering a moving frame of reference, centredon the substrate [18], whereas substrates which exhibit bending need to be consid-ered using more complex techniques such as the immersed boundary method [26].Two-phase flow systems with immersed boundaries have not been studied exten-sively, with only a few noticeable exceptions [23,33]. Despite considering complexmoving boundaries (such as a twin screw kneader), the motion of the boundarieswere prescribed rather than resulting from fluid-structure interaction, although theproposed methods could in principle be extended to consider this. More recently,the impact of droplets in capillary-dominated regimes onto a flexible substrate hasbeen modelled using the Lattice-Boltzmann method [44], focusing on the spreadingand rebound of the droplets and how the contact time is affected by the bendingstiffness. The late-time spreading and rebound dynamics of an undamped plate-spring system have recently been studied numerically [45], inspired by the feathersof kingfishers. It was found that springs with certain stiffness values can shortenthe length of time the droplet is in contact with the substrate, as well as increasethe speed the droplet rebounds after impact.Analytical models for a liquid impacting a deformable substrate (or vice-versa),on the other hand, have been proposed for over half a century. One of the earliestmodels investigated a droplet impacting onto an elastic half-space by imposing aconstant uniform pressure over a circle whose radius increased in proportion tothe square root of time [2]. The full hydroelastic problem of the impact of a two-dimensional wave onto an Euler-Bernoulli beam has previously been studied usingWagner theory [16], and more recently this analytical model has been extended tostudy the axisymmetric impact of a droplet onto an elastic plate, where the elasticplate has a radius much smaller than the droplet and its deflection governed bythin-plate theory [24]. A thorough parameter study was conducted for differenttypes of plate, and regimes where the elasticity of the plate could cause splashingof the droplet at early times, defined as the detachment of the splash sheet from the surface of the plate, have been identified.To the best of our knowledge, a comprehensive study considering the fluid-structure interaction between an impacting droplet and a compliant substrate thatsystematically compares accurate numerical results to analytical or experimentalcounterparts has yet to materialise. One of the simplest types of deformable sub-
Michael J. Negus et al. strate systems which exhibits both elastic and damping effects is a rigid plate sus-pended by a Hookean spring and a linear dashpot, where the force of the impactingdroplet causes the spring to compress, with the dashpot damping the motion. Herewe present both a new analytical model extended from Wagner theory, as well asa direct numerical simulation platform for a droplet impacting onto a rigid platesupported by a spring and a dashpot. The system is used as a validation testbedfor the two approaches, as well as a framework for an extensive parametric studyfocused on the effect of surface compliance on the ensuing drop dynamics in thischallenging regime. We observe systematically how the influence of the substrateproperties (mass, spring stiffness and damping factor) affects the fluid-structureinteraction, emphasising how the resulting motion of the plate substantially altersthe pressure field of the droplet, and, in turn, the hydrodynamic force exerted ontothe plate, revealing a rich coupling of the different forces at play in the system.The rest of this paper is structured as follows. We outline the system geometryand general mathematical framework, as well as discuss assumptions at the level ofthe analytical and numerical models in §
2. In § §
4. The solutions from the analytical and numericalmodels are presented in §
5, where they are compared for a range of parametersof the plate-spring-dashpot system. We conclude by summarising our results anddiscussing the implications of this study and possible future extensions in § We consider the vertical impact of a droplet of incompressible, Newtonian liquidonto a rigid, planar, circular plate supported by a Hookean spring and lineardashpot. The droplet is initially spherical with radius R ∗ d and travelling uniformlydownwards with speed V ∗ , surrounded by an incompressible gas. Throughout thisstudy, dimensional quantities are denoted with a superscript ∗ . The plate hasradius R ∗ p and the plate-spring-dashpot system is initially in equilibrium. Thebottom of the droplet is introduced at a height δ ∗ > x ∗ , y ∗ , z ∗ ) is defined such that the surface of the plate lies inthe z ∗ = 0 plane, the droplet falls along the z ∗ > x ∗ , y ∗ , z ∗ ) = (0 , , δ ∗ ) at the onset of the dynamics, asillustrated in Fig. 1.The system is initialised at a time t ∗ = t ∗ = − δ ∗ /V ∗ . If the gas were absent, theplate would experience a zero net force until the droplet makes contact, and wouldtherefore remain in equilibrium. In this case, the droplet would make contact withthe plate at t ∗ = 0 with the plate stationary for t ∗ <
0. However the presence ofthe gas means there will be a pressure build-up prior to the impact [43,9], resultingin a net force that causes the plate to accelerate downwards so the droplet will make contact at a time t ∗ c > ρ ∗ l , ρ ∗ g and viscosities µ ∗ l , µ ∗ g respectively. The surface tension coefficient between theliquid and the gas is denoted by σ ∗ (taken to be constant) and the accelerationdue to gravity is g ∗ = g ∗ ˆn z , where ˆn z is the unit vector in the z ∗ direction. The roplet impact onto a spring-supported plate: analysis and simulations 5 Fig. 1: Schematic diagram of a spherical droplet of radius R ∗ d impacting with initialvelocity V ∗ onto a circular plate of radius R ∗ p and mass M ∗ . The plate is supportedby a spring with spring constant k ∗ and a dashpot with damping factor c ∗ .vertical position of the plate at time t ∗ is z ∗ = − s ∗ ( t ∗ ), where s ∗ ( t ∗ ) is referredto as the plate displacement.Denoting the variables in the liquid and gas with a subscript l and g respec-tively, the Navier-Stokes equations are assumed to hold in each fluid, ρ ∗ i (cid:18) ∂ u ∗ i ∂t ∗ + ( u ∗ i · ∇ ) u ∗ i (cid:19) = − ∇ p ∗ i + µ ∗ i ∇ u ∗ i − ρ ∗ i g ∗ , (1) ∇ · u ∗ i = 0 , (2)where i = l, g , u ∗ i is the velocity vector and p ∗ i represents the pressure in eachfluid.The impermeability condition on the plate states that the fluid velocity mustmatch the velocity of the plate along its surface, u ∗ i · ˆn z = − ˙ s ∗ ( t ∗ ) for z ∗ = − s ∗ ( t ∗ ) , x ∗ + y ∗ < R ∗ p , (3)where the overdot denotes differentiation with respect to time.The kinematic condition at the multivalued free surface z ∗ = h ∗ ( x ∗ , y ∗ , t ∗ ) is v ∗ n = u ∗ l · ˆn on z ∗ = h ∗ ( x ∗ , y ∗ , t ∗ ) , (4)where ˆn is the unit outward-pointing normal vector to the free surface and v ∗ n is the normal speed of the free surface. Continuity of velocity and normal stress across the free surface are given by u ∗ g = u ∗ l , ˆn · [ T ∗ g − T ∗ l ] = − σ ∗ κ ∗ ˆn on z ∗ = h ∗ ( x ∗ , y ∗ , t ∗ ) , (5)where T ∗ i is the Cauchy stress tensor in each fluid and κ ∗ = − ∇· ˆn is the curvatureof the free surface. Michael J. Negus et al.
Initially, at t ∗ = t ∗ , the liquid has a uniform downwards velocity V ∗ , u ∗ l ≡ − V ∗ ˆn z , (6)while the centre of the droplet is initially at z ∗ = δ ∗ + R ∗ d , meaning the free surface h ∗ ( x ∗ , y ∗ , t ∗ ) satisfies x ∗ + y ∗ + ( h ∗ ( x ∗ , y ∗ , t ∗ ) − δ ∗ − R ∗ d ) = R ∗ d . (7)The pressure p ∗ i is initially hydrostatic in each fluid, equal to p ∗ i = p ∗ atm − ρ ∗ i g ∗ z ∗ ,where p ∗ atm denotes atmospheric pressure. The gas far from the impact remainsundisturbed for all times, so that u ∗ g ∼ , p ∗ g ∼ p ∗ atm − ρ ∗ g g ∗ z ∗ as x ∗ + y ∗ + z ∗ → ∞ . (8)The plate, which has total mass M ∗ , is supported by a Hookean spring withspring constant k ∗ , and a dashpot with damping factor c ∗ . At t ∗ = t ∗ , the plate isin equilibrium. Hence by Newton’s third law, the force due to the compression ofthe spring balances the weight of the plate. Denoting the net hydrodynamic forceapplied to the plate in the downwards direction − ˆn z as F ∗ ( t ∗ ), the displacementof the plate from this equilibrium is governed by M ∗ ¨ s ∗ ( t ∗ ) = F ∗ ( t ∗ ) − c ∗ ˙ s ∗ ( t ∗ ) − k ∗ s ∗ ( t ∗ ) . (9)The net hydrodynamic force F ∗ ( t ∗ ) is equal to the sum of the contributionsfrom the hydrodynamic pressure and viscous stress above and below the plate.Assuming that the gas below the plate is at a constant pressure p ∗ atm and exerts anegligible amount of force due to viscous stress, then F ∗ ( t ∗ ) = (cid:90)(cid:90) √ x ∗ + y ∗ ≤ R ∗ p z ∗ =0 + ( p ∗ − p ∗ atm ) − µ ∗ ∂u ∗ z ∂z ∗ d x ∗ d y ∗ , (10)where p ∗ = p ∗ l , µ ∗ = µ ∗ l , u ∗ z = u ∗ l,z where the plate is wetted and p ∗ = p ∗ g , µ ∗ = µ ∗ g , u ∗ z = u ∗ g,z where the plate is unwetted.2.1 Non-dimensionalisationWe take the initial droplet radius, R ∗ d , and speed, V ∗ , as the characteristic lengthand velocity scales respectively. Then, choosing the advective and inertial timeand pressure scales, we non-dimensionalise by setting t ∗ = R ∗ d V ∗ t, ( x ∗ , y ∗ , z ∗ , h ∗ , s ∗ , R ∗ p ) = R ∗ d ( x, y, z, h, s, R p ) , u ∗ = V ∗ u , p ∗ = p ∗ atm + ρ ∗ l V ∗ p, F ∗ ( t ∗ ) = ρ ∗ l V ∗ R ∗ d F ( t ) , (11) where R p is referred to as the dimensionless plate radius.Under these scalings, the plate dispacement equation (9) becomes α ¨ s ( t ) + β ˙ s ( t ) + γs ( t ) = F ( t ) , (12) roplet impact onto a spring-supported plate: analysis and simulations 7 where α = M ∗ ρ ∗ l R ∗ d , β = c ∗ ρ ∗ l V ∗ R ∗ d , γ = k ∗ ρ ∗ l V ∗ R ∗ d . (13)The ratio between the mass of the plate and the mass of the droplet is equal to3 α/ π , hence α is referred to as the mass ratio. The damping factor, β , measuresthe strength of the resistance to motion due to the damping from the dashpot, andthe stiffness factor, γ , measures the strength of the restoring force due to elasticcompression of the spring.The relevant dimensionless parameters describing the flow dynamics are theReynolds, Weber and Froude numbers, defined respectively byRe = ρ ∗ l R ∗ d V ∗ µ ∗ l , We = ρ ∗ l R ∗ d V ∗ σ ∗ , Fr = V ∗ g ∗ R ∗ d . (14)Finally, the ratios between the densities and the viscosities of the gas and theliquid are given by ρ R = ρ ∗ g ρ ∗ l , µ R = µ ∗ g µ ∗ l . (15)2.2 Modelling assumptionsIn both § §
4, scenarios where the inertial effects of the impact are moresignificant than the effects of viscosity, surface tension and gravity are considered.Hence we assume throughout that the values of Re, We and Fr are large. As anillustrative example, consider the impact of a droplet of water with radius R ∗ d = 1mm and velocity V ∗ = 5 m/s, surrounded by air under atmospheric conditions.This gives Re ≈ , We ≈ , Fr ≈ , (16)and we make the assumption that they remain large throughout the early stagesof impact. Furthermore, the density and viscosity ratios for the air-water scenarioare ρ R ≈ . × − , µ R ≈ . × − , (17)which provides support for neglecting the effects of the gas phase in the analyticalmodel (but not in the DNS).As the system is initially radially symmetric about the z − axis, we assume thissymmetry remains in place throughout the impact and consider an axisymmetriccoordinate system ( r, z ), where r = x + y . This assumption restricts the applica-bility of the model away from systems that involve fully three-dimensional effects,such as prompt splashing [36] and azimuthal instabilities of the ejecta sheet [17]. The analytical model focuses on the dynamics shortly after the time of impact.The model follows the structure of previous works on axisymmetric Wagner theory, and the reader is directed to [11,12,20,22,27] for more details on the generalmethodology for impact involving stationary substrates. In this context the dropletimpacts the plate at t = 0, with the displacement and velocity of the plate equalto zero at t = 0. As the gas phase is ignored, all expressed quantities are in theliquid, and the subscript l is dropped for brevity. Michael J. Negus et al. t < t >
0. Therefore a velocity potential φ can be introduced, such that u = ∇ φ . The dimensional continuity equation (2)transforms to Laplace’s equation for φ , ∇ φ = ∂ φ∂r + 1 r ∂φ∂r + ∂ φ∂z = 0 , (18)and the dimensional momentum equation (1) results in the unsteady Bernoulliequation for φ and p , p + ∂φ∂t + 12 | ∇ φ | = C ( t ) , (19)for some C ( t ). The absence of viscosity and surface tension means the continuityof normal stress boundary condition (5) reduces to specifying that p = 0 at thefree surface. Finally, by neglecting the gas phase and liquid viscosity, the nethydrodynamic force (10) is just equal to the integral of the pressure p across thewetted part of the plate. Hence the governing equations are a set of non-linearequations for the velocity potential φ , pressure p , free surface location h and platedisplacement s .3.2 Asymptotic structureFollowing the structure of previous analytical models for droplet impact [11], weidentify that for t (cid:28)
1, the radial extent of the penetration region (where thedroplet would be below the r axis were the plate not present) is O ( √ t ). Giventhis, we introduce an arbitrarily small parameter 0 < (cid:15) (cid:28) t = (cid:15) ˆ t, (20)where ˆ t = O (1) as (cid:15) →
0. With the plate present, the free surface is violentlydisplaced and the liquid is ejected along the plate in a splash sheet. The curve atwhich the free surface is vertical is called the turnover curve , and for small timeswe assume its radial extent is close to that of the penetration region, meaning that r = (cid:15) ˆ d (ˆ t ) , ˆ d (0) = 0 , (21)where the turnover curve ˆ d (ˆ t ) = O (1) as (cid:15) →
0. The bottom of the penetrationregion is at z = − (cid:15) ˆ t , so if we assume that the plate acts to decelerate the verticalmotion of the droplet at early times, then the plate position z = − s ( t ) > − (cid:15) ˆ t ,which motivates the rescaling s ( t ) = (cid:15) ˆ s (ˆ t ) , (22) where ˆ s (ˆ t ) = O (1) as (cid:15) →
0. For brevity, the hat notation is dropped for the restof § (cid:15) →
0, the problem breaks down into four distinct spatial regions, as de-picted in Fig. 2. In the outer-outer region, for which ( r, z ) = O (1) × O (1), the bulkof the droplet is unaffected by the plate to leading-order, and is hence spherical roplet impact onto a spring-supported plate: analysis and simulations 9 Fig. 2: Schematic of the asymptotic structure of the system. The displacement ofthe plate − (cid:15) s ( t ) has been exaggerated for visibility.and moving downwards with unit speed. This means that C ( t ) = 1 / O ( (cid:15) ) × O ( (cid:15) ) region close tothe centre of the plate is referred to as the outer region, and here the splash sheetis neglected and the boundary conditions are linearised onto the undisturbed platelocation. This solution breaks down close to the turnover curve where the velocityand pressure become singular, which is corrected by introducing an inner regionof size O ( (cid:15) ) × O ( (cid:15) ) in which the free surface turns over, ejecting fluid into thesplash sheet. The thin, fast-moving splash sheet region is of size O ( (cid:15) ) × O ( (cid:15) ),emanating across the surface of the plate from the inner region. In this analysis,we assume the splash jet does not detach from the plate.In the present analysis we shall consider the outer, inner and splash sheetregions in detail, however we forgo an analysis for the outer-outer region as it doesnot contribute to the leading-order hydrodynamic force on the plate.3.3 Outer region Guided by a well-known scaling argument [11], in the outer region we set r = (cid:15) ˆ r, z = (cid:15) ˆ z, φ = (cid:15) ˆ φ, h = (cid:15) ˆ h, p = 1 (cid:15) ˆ p, and expand ( ˆ φ, ˆ h, ˆ p, d, s ) = ( ˆ φ , ˆ h , ˆ p , d , s ) + o (1) as (cid:15) →
0. The resultinggoverning equations in the outer region are ∂ ˆ φ ∂ ˆ r + 1ˆ r ∂ ˆ φ ∂ ˆ r + ∂ ˆ φ ∂ ˆ z = 0 for ˆ z > , (23) ∂ ˆ φ ∂ ˆ z = − ˙ s ( t ) on ˆ z = 0 , ˆ r < d ( t ) (24) ∂ ˆ φ ∂ ˆ z = ∂ ˆ h ∂t , on ˆ z = 0 , ˆ r > d ( t ) , (25)ˆ φ = 0 on ˆ z = 0 , ˆ r > d ( t ) , (26)where (23) is Laplace’s equation for ˆ φ , (24) is the kinematic boundary conditionon the plate, (25) is the kinematic boundary condition at the free surface and(26) is the dynamic boundary condition at the free surface, where (26) is foundby integrating the leading-order Bernoulli equation (19) with respect to time andapplying the initial condition for ˆ φ .The far-field conditions as we tend towards the outer-outer region state that toleading-order, the flow is travelling downwards with speed 1 and the free surfaceis parabolic in ˆ r , such that ˆ φ (ˆ r, ˆ z, t ) ∼ − ˆ z, as (cid:112) ˆ r + ˆ z → ∞ and ˆ h (ˆ r, t ) ∼ r − t as ˆ r → ∞ , (27)with subsequent initial conditionsˆ h (ˆ r,
0) = 12 ˆ r , s (0) = ˙ s (0) = 0 , d (0) = 0 . (28)In addition, the so-called Wagner condition is needed in order to match to theinner solution, ˆ h ( d ( t ) , t ) = − s ( t ) , (29)which states that, at leading-order, the free surface meets the plate at the turnoverpoint [12]. Finally, matching with the inner region reveals that ˆ φ = O ( (cid:112) d ( t ) − ˆ r )as ˆ r → d ( t ) − , as in the classical Wagner regime [12].Following [15], it is useful to consider the variational formulation of the ax-isymmetric problem by introducing the leading-order displacement potential, Υ ,as Υ (ˆ r, ˆ z, t ) = ˆ zt + (cid:90) t ˆ φ (ˆ r, ˆ z, τ ) d τ , (30)which is governed by the equations ∂ Υ ∂ ˆ r + 1 ˆ r ∂Υ ∂ ˆ r + ∂ Υ ∂ ˆ z = 0 for ˆ z > , (31) ∂Υ ∂ ˆ z = ( t − s ( t )) −
12 ˆ r on ˆ z = 0 , ˆ r < d ( t ) , (32) ∂Υ ∂ ˆ z = t + ˆ h (ˆ r, t ) −
12 ˆ r on ˆ z = 0 , ˆ r > d ( t ) , (33) Υ = 0 on ˆ z = 0 , ˆ r > d ( t ) , (34) roplet impact onto a spring-supported plate: analysis and simulations 11 such that the displacement potential is Υ = O (( d ( t ) − ˆ r )) / ) as ˆ r → d ( t ) − .The far-field condition for ˆ φ (27) implies that Υ is bounded as √ ˆ r + ˆ z → ∞ ,which means a separable solution for Υ can be found via Υ = (cid:90) ∞ ν ( λ, t ) e − λ ˆ z J ( λ ˆ r ) d λ , (35)where J is a Bessel function of the first kind and ν ( λ, t ) is a coefficient functionfound by solving the dual integral equations necessary to satisfy (32) and (34),namely ν ( λ, t ) = 2 π (cid:90) d ( t )0 σ (cid:20) σ − ( t − s ( t )) (cid:21) sin( λσ ) d σ . (36)We refer to previous studies [20,32] for further details.Evaluating (35) along the plate and expanding as ˆ r → d ( t ) − , we find Υ (ˆ r, , t ) = 2 π (cid:112) d ( t ) (cid:18) d ( t ) − ( t − s ( t )) (cid:19) (cid:112) d ( t ) − ˆ r − √ d ( t ) / π ( d ( t ) − ˆ r ) / + O (( d ( t ) − ˆ r ) / ) . (37) Hence, to enforce Υ = O (( d ( t ) − ˆ r ) / ) as ˆ r → d ( t ) − means that we must have d ( t ) = (cid:112) t − s ( t )) . (38)Note that if the plate is stationary, s ( t ) ≡ d ( t ) = √ t [42]. Clearly, when the plate is compliant, the dis-placement of the plate is expected to slow down the spreading of the droplet, atleast for times early enough that s ( t ) >
0. It is well known that the Wagner prob-lem is unstable under time-reversal [12], which means the solution breaks downif ˙ d ( t ) ≤
0. We therefore assume that ˙ s ( t ) < d ( t ) = √ − ˙ s ( t )) / (2 (cid:112) t − s ( t )), to remain positive.Finally, since it is necessary to calculate the hydrodynamic force on the plate F ( t ), we use the leading-order form of the Bernoulli equation (19) to find thepressure on the plateˆ p (ˆ r, , t ) = − ∂ Υ ∂t = 49 π d d t (cid:104) ( d ( t ) − ˆ r ) / (cid:105) , (39)where d ( t ) is given in terms of s ( t ) in (38).3.4 Inner regionSince the pressure is locally singular, there is an inner region moving with the turnover point at r = (cid:15)d ( t ) and the surface of the plate at z = − (cid:15) s ( t ). Theappropriate scalings are given by [12], r = (cid:15)d ( t ) + (cid:15) ˜ r, z = − (cid:15) s ( t ) + (cid:15) ˜ z,φ = (cid:15) (cid:104) ˙ d ( t )˜ r − (cid:15) ˙ s ( t )˜ z + ˜ φ (cid:105) , h = − (cid:15) s ( t ) + (cid:15) ˜ h, p = 1 (cid:15) ˜ p. (40) Under these scalings, it is straightforward to show that to leading-order, the plateis stationary in the inner region. Hence, as described in detail in [20], the leading-order inner problem is quasi-two-dimensional in each plane perpendicular to it,and is given by ∂ ˜ φ ∂ ˜ r + ∂ ˜ φ ∂ ˜ z = 0 for ˜ z > , (41) ∂ ˜ φ ∂ ˜ z = 0 on ˜ z = 0 , (42) ∂ ˜ φ ∂ ˜ z = ∂ ˜ φ ∂ ˜ r ∂ ˜ h ∂ ˜ r on ˜ z = ˜ h (˜ r, t ) , (43) (cid:18) ∂ ˜ φ ∂ ˜ r (cid:19) + (cid:18) ∂ ˜ φ ∂ ˜ z (cid:19) = ˙ d ( t ) on ˜ z = ˜ h (˜ r, t ) , (44)subject to appropriate matching conditions into the outer region,˜ φ ∼ − ˙ d ( t )˜ r + O ( (cid:112) ˜ r + ˜ z ) as (cid:112) ˆ r + ˆ z → ∞ , (45)and towards the splash sheet˜ h (˜ r, t ) → J ( t ) as ˆ r → ∞ , (46)where J ( t ) is referred to as the asymptotic sheet thickness.The solution to this problem is well known [39], and given parametrically by˜ φ = a ( t ) − ˙ d ( t ) J ( t ) π (1 + R [ ζ + log( ζ )]) , ˜ r + i ˜ z = J ( t ) π (cid:104) ζ + 4 i (cid:112) ζ − log( ζ ) + iπ − (cid:105) , (47)where ζ ∈ C , I ( ζ ) > a ( t ) ∈ C is an arbitrary function of time, the branch cutsfor log( ζ ) and √ ζ are taken along R ( ζ ) < , I ( ζ ) = 0 + and R , I denote the realand imaginary parts of a complex number.To match with the leading-order outer solution, we take | ζ | → ∞ in (47), whichyields ˜ φ ∼ − ˙ d ( t )˜ r + 4 ˙ d ( t ) J ( t ) π R (cid:104) i √ ˜ r + i ˜ z (cid:105) . (48)Thus, comparing (48) with (37) gives the leading-order sheet thickness J ( t ) = 2 d ( t ) π = 2 √ π ( t − s ( t )) / . (49)Again note that as d ( t ) = (cid:112) t − s ( t )), the displacement of the plate slows thespreading of the droplet, which leads to a thinner splash sheet. This is consistentwith the findings of [13], who showed that soft substrates inhibit splashing. Notethat the derivative of the sheet thickness ˙ J ( t ) = √ − ˙ s ( t )) (cid:112) t − s ( t ) /π is positive for all t , so the sheet thickness will still increase for all time within theWagner model.The leading-order pressure in the inner region is˜ p = − (cid:34)(cid:18) ∂ ˜ φ ∂ ˜ r (cid:19) + (cid:18) ∂ ˜ φ ∂ ˜ z (cid:19) − ˙ d ( t ) (cid:35) . (50) roplet impact onto a spring-supported plate: analysis and simulations 13 Along the surface of the plate, where ˜ z = 0, this solution is given parametricallyby˜ p (˜ r, , t ) = 2 ˙ d ( t ) e η (1 + e η ) , ˜ r = − J ( t ) π (cid:104) e η + 4 e η + 2 η + 1 (cid:105) for − ∞ < η < ∞ , (51)where d ( t ) and J ( t ) are given in terms of s ( t ) in (38) and (49) respectively.3.5 Splash sheet regionUpon impact, the fluid is ejected from the inner region into a thin, fast-movingsheet of fluid attached to the plate. In this region, we rescale r = (cid:15) ¯ r, z = − (cid:15) s ( t ) + (cid:15) ¯ z, h = − (cid:15) s ( t ) + (cid:15) ¯ h, φ = − (cid:15) ˙ s ( t )¯ z + ¯ φ, p = (cid:15) ¯ p. (52)As described in detail by [22], the leading-order splash sheet problem for the ra-dial velocity ¯ u = ∂ ¯ φ /∂ ¯ r and free surface height ¯ h reduces to the zero-gravityshallow-water equations. These equations can be solved using the method of char-acteristics and the solution is, as derived in [20,22],¯ r = 2 ˙ d ( τ )( t − τ ) + d ( τ ) , ¯ u = 2 ˙ d ( τ ) , ¯ h = ˙ d ( τ ) J ( τ )˙ d ( τ ) − d ( τ )( t − τ ) , (53)where 0 < τ < t .The subsequent solution for the pressure in the splash sheet region is foundby differentiating the Bernoulli equation (19) with respect to z , expressing in thesplash sheet variables, and expanding to leading-order, such that ∂ ¯ p ∂ ¯ z = ¨ s ( t ) − ˙ s ( t ) ∂ ¯ u ∂ ¯ r . (54)Integrating with respect to ¯ z and noting that ¯ p = 0 at ¯ z = ¯ h means the leading-order pressure along the plate is¯ p (¯ r, , t ) = (cid:18) ˙ s ( t ) ∂ ¯ u ∂ ¯ r − ¨ s ( t ) (cid:19) ¯ h (¯ r, t ) . (55)It is worth noting that in the classical case of a stationary plate, where s ( t ) ≡ p = O ( (cid:15) ). Thereforethe velocity ˙ s ( t ) and acceleration ¨ s ( t ) of the plate increase the magnitude of thepressure in the splash sheet region. In particular, if ˙ s ( t ) ∂ ¯ u /∂ ¯ r − ¨ s ( t ) <
0, the leading-order pressure would be below atmospheric pressure, which could providea possible mechanism for splash sheet detachment. However, the contribution thepressure in the splash sheet region makes to the leading-order force is still O ( (cid:15) ),which is lower in magnitude than the contributions from the outer and innerregions, so we shall neglect it henceforth in this analysis. Fig. 3: Analytical solutions for the pressure for (cid:15) = 0 . t = 5, with the left-handside black lines showing the stationary plate case, s ( t ) = 0 (reflected in r ), andthe right-hand side grey lines the moving plate case, s ( t ) = 0 . t . The outersolution for the pressure (39) is shown by the dashed lines, the inner solution (51)by the dotted lines and the composite solution (58) by the solid lines. The thinvertical dashed lines show the location of the turnover point r = (cid:15)d ( t ).3.6 Composite pressureClassically, the hydrodynamic force is determined by integrating only the outerpressure (39) for 0 ≤ ˆ r < d ( t ) [22]. However, as will be discussed in §
6, wefind better agreement with the numerical simulations at later times by using thecomposite expansion between the outer and inner regions. Van Dyke’s matchingprinciple [37] is used to find the overlap function between the outer and innersolutions for r < (cid:15)d ( t ). By evaluating the time derivatives in (39), the one-term-outer pressure is(1.t.o) p ( r, , t ) = 1 (cid:15) (cid:34) − r − d ( t ) ) ˙ d ( t ) π (cid:112) d ( t ) − ˆ r + 4 d ( t ) ¨ d ( t )3 π (cid:112) d ( t ) − ˆ r (cid:35) . (56)Expressing this in the inner variables and expanding to leading order gives theoverlap pressure p overlap ( r, , t ) = (1.t.i)(1.t.o) p ( r, , t ) = 2 √ d ( t ) / ˙ d ( t ) π(cid:15) √− ˜ r = 2 √ d ( t ) / ˙ d ( t ) π √ (cid:15) (cid:112) (cid:15)d ( t ) − r . (57)Therefore the composite expansion for p ( r, , t ) is p comp ( r, t ) = H ( (cid:15)d ( t ) − r ) (cid:20) (cid:15) ˆ p ( r/(cid:15), , t ) − p overlap ( r, , t ) (cid:21) + 1 (cid:15) ˜ p ( r/(cid:15) − d ( t ) /(cid:15) , , t ) , (58)where H is the Heaviside step function, ˆ p is given by (39) and ˜ p by (51).The composite pressure profile (58) depends on the plate displacement s ( t ),which is solved for in § roplet impact onto a spring-supported plate: analysis and simulations 15 order to illustrate the effects of a moving substrate on the pressure, we compare itsvalue in the case where the plate is stationary ( s ( t ) ≡
0) to that of a prescribedmoving plate case where s ( t ) = 0 . t . Note that this value for s ( t ) is chosen forillustrative purposes and is not a solution to (12), and only satisfies the assumptionthat ˙ s ( t ) < t < t = 5 in Fig. 3, where the left-hand side black linesshows the outer, inner and composite solutions to the pressure for the stationaryplate case (reflected in r ), and the right-hand side grey lines shown the correspond-ing values for the moving plate case. It can be seen that the turnover point in themoving plate case has advanced less than in the stationary plate case, as accordingto (38). The pressure in the moving plate case is significantly lower overall than inthe stationary plate case. This shows how the downwards motion of the plate doesnot just slow the spreading, but also decreases the hydrodynamic pressure insidethe droplet. Also noteworthy is that the inner solution under-estimates the solu-tion away from the turnover point in the stationary plate case, but over-estimatesit in the moving plate case.3.7 Hydrodynamic forceIn order to solve (12) for the displacement of the plate s ( t ), the value of the hy-drodynamic force, F ( t ), needs to be determined to leading order. We approximatethe force by integrating the composite expansion to the pressure (58) across theouter and inner regions. The composite expansion for the force is hence F comp ( t ) = 89 (cid:15)d ( t ) ((4 − √
2) ˙ d ( t ) + ¨ d ( t ) d ( t ))+ 8 (cid:15) ˙ d ( t ) J ( t ) π e η ( t ) (cid:20) πd ( t ) (cid:15) J ( t ) + 1 − e η ( t ) − e η ( t ) − η ( t ) (cid:21) , (59)where η ( t ) is defined implicitly by e η ( t ) + 4 e η ( t ) + 2 η ( t ) + 1 = πd ( t ) (cid:15) J ( t ) , (60)where the detailed derivation can be found in Appendix A.It is worth noting that the composite force F comp ( t ) differs from the force ona stationary plate if and only if the turnover point d ( t ) or sheet thickness J ( t )differ from their corresponding stationary plate values.3.8 Plate displacement solutionThe remaining unknown from § s ( t ), where s ( t ) appears in the solution (38) for the turnover point d ( t ) and in the solution (49) for the jet-thickness J ( t ). The plate displacement is foundby solving the second order ordinary differential equation (12), approximating theforce term F ( t ) by the composite force F comp ( t ) (59). The resulting equation is non-linear and implicit, and is solved using MATLAB’s ode15i solver in conjunctionwith the fsolve solver to find η ( t ) via (60) at each timestep. As the value of ˙ d ( t ) diverges at t = 0, the numerical scheme is initialised at a time t = t i = 10 − ,with zero initial guesses for s ( t i ) = ˙ s ( t i ) = ¨ s ( t i ) = 0. A small-time asymptoticanalysis of the plate displacement reveals that s ( t ) = O ( t / ) as t →
0, so theproblem is regular and we are hence justified in taking zero initial guesses. Theresults will be discussed in comparison to the full DNS in § (cid:15) = 0 .
1, the numerical solution for s ( t ) is found for 0 ≤ t ≤
100 on 1CPU in approximately 10 seconds. In comparison, the DNS results in § We build on the open-source, volume-of-fluid package
Basilisk [29] to implementthis complex multi-phase system, retaining effects due to viscosity and density inboth the liquid and the gas, as well as surface tension and gravity.
Basilisk , and itspredecessor
Gerris [28], have been used extensively to study interfacial flows, andin particular droplet impact, with great success over the past two decades, cross-fertilising investigative efforts within experimental, analytical and computationalcommunities alike [4,19,27,34,41].4.1 Moving frame coordinatesIn order to avoid including an embedded boundary in the quadtree-structured com-putational domain, we transfer the flow into a frame moving with the plate, fixingthe plate along the bottom computational boundary. The dimensionless movingframe coordinates are defined by x (cid:48) = ( x (cid:48) , y (cid:48) , z (cid:48) ) = x + s ( t ), u (cid:48) = u + ˙ s ( t ), where s ( t ) = (0 , , s ( t )), and the prime (cid:48) decorates all quantities in the moving frame. In-troducing the dimensionless variable density ρ (cid:48) ( x (cid:48) , t ) and viscosity µ (cid:48) ( x (cid:48) , t ), suchthat, following notation from § ρ (cid:48) = 1, µ (cid:48) = 1 in the liquid and ρ (cid:48) = ρ R , µ (cid:48) = µ R in the gas, the dimensionless governing equations in the moving frameare given by ρ (cid:48) (cid:18) ∂ u (cid:48) ∂t + ( u (cid:48) · ∇ (cid:48) ) u (cid:48) (cid:19) = − ∇ (cid:48) p (cid:48) + µ (cid:48) Re ( ∇ (cid:48) ) u (cid:48) + κ (cid:48) δ (cid:48) s We ˆn (cid:48) + ρ (cid:48) ¨ s ( t ) − ρ (cid:48) Fr ˆn (cid:48) z , (61) ∇ (cid:48) · u (cid:48) = 0 , (62) u (cid:48) z (cid:48) = 0 for z (cid:48) = 0 , x (cid:48) + y (cid:48) < R p , (63)where δ (cid:48) s is a Dirac distribution centred on the liquid-gas free surface and ˆn (cid:48) is theunit normal to the interface. Note that the kinematic condition (63) is that of a stationary plate, and the problem in the moving frame is equivalent to a dropletimpacting onto a stationary plate, with the liquid and gas under additional forcingequal to ρ (cid:48) ¨ s ( t ). In the far-field, we assume the pressure tends to 0 (neglectingvariations due to gravity) and the vertical velocity tends to u (cid:48) → ˙ s ( t ). The primenotation is dropped for the remainder of this section for brevity. roplet impact onto a spring-supported plate: analysis and simulations 17 Fig. 4: Direct numerical simulation setup at t = t = − . .
125 andtravelling with a vertical velocity −
1. The inset shows a snapshot of a simulationat t = 0 .
045 close to the surface, around the area indicated with the grey rectangle.The colour map illustrates the adaptive mesh refinement strategy, while the blackline depicts the location of the interface.4.2 Computational setupIn all our simulations, we consider a droplet with dimensional radius R ∗ d = 1 mminitially travelling vertically downwards at speed V ∗ = 5 m/s, where the values ofRe, We and Fr , ρ R , µ R , are given in § R p = 2, and the initial separationbetween the bottom of the droplet and top of the plate is δ ∗ = 0 . R ∗ d = 0 . α , damping factor β and stiffness factor γ arevaried across different parametric studies.The dimensionless governing equations (61)-(62) are solved using Basilisk within an axisymmetric domain, where the droplet initially has unit dimensionlessradius, travelling with uniform vertical velocity −
1. The axisymmetric computa-tional domain for the simulations is shown in Fig. 4. The domain is given by asquare box, with the r axis along the bottom boundary and the z axis along the left boundary. The side length of the domain is set to L = 6, which is sufficientlylarge so that far-field conditions do not artificially alter the target dynamics. Neu-mann conditions ∂ n u = 0 are specified along the top and right boundaries, where ∂ n is the partial derivative in the normal direction to the boundary. The verticalvelocity along the right boundary is specified as u z = ˙ s ( t ), to reflect the far-field velocity condition. The appropriate symmetry conditions are applied along the leftboundary. A mixed boundary condition is specified along z = 0, with u z = 0 for r ≤ R p and u z = ˙ s ( t ) for r > R p , the former representing the kinematic conditionalong the plate (63) and the latter representing the far-field condition. The no-slipcondition u r = 0 and a 90 ◦ contact angle are applied along the bottom bound-ary. Finally, the pressure p in the gas is set to to be zero along the top and rightboundaries in line with the far-field condition (8).The quadtree grid construction features in Basilisk allow for a high grid res-olution in areas of interest, varying from levels 5 to 13, where level n correspondsto 2 n square cells per dimension, if the grid were uniform. In this case, the largestcell has side length 6 / , corresponding to 0 .
188 mm in dimensional terms. Thesmallest cell has side length 6 / , corresponding to 0 . µ m. In order to accu-rately calculate the force on the plate, a region along the bottom boundary for0 ≤ r ≤ R p of height 24 / ( ≈ . µ m in dimensional terms) is held at level13, such that the bottom four grid cells in this region are at maximum level. Thisavoids numerical errors induced by multi-grid interpolation in a region which re-quires particular care due to the delicate fluid-structure interaction calculationsoutlined below. Adaptive mesh refinement is also used to refine the domain in re-gions where the velocities and interface location are rapidly changing. An exampleof the typical grid structure is shown in the inset of Fig. 4.At regular timesteps ∆t = 10 − (corresponding to 20 ns in dimensional terms),the hydrodynamic force applied to the surface of the plate, F ( t ), is determinedby numerically integrating the pressure, p , and the viscous stress, − µ∂u z /∂z ,along the bottom boundary for 0 ≤ r ≤ R p . This value of the force is then usedto solve the dimensionless plate displacement equation (12) using a second-orderfinite difference scheme, giving s ( t ) , ˙ s ( t ) and ¨ s ( t ). The boundary conditions arethen updated with the new value of ˙ s ( t ), and the vertical acceleration in all of thecells is incremented by ρ ¨ s ( t ).Several computational details are noteworthy in terms of ensuring a robustfluid-structure interaction calculation procedure. As observed in other studies ondroplet impact [27], numerical instabilities in the projection solver used for the re-sulting Poisson equation within Basilisk may cause the calculated pressure valuesto fluctuate between timesteps, thus causing the resulting force values to vary ar-tificially. These pressure spikes lead to artifacts in the finite difference scheme,which can ultimately result in the simulation breaking down due to diverging ac-celeration terms. To prevent this, we use a peak-detection algorithm [3] to identifynumerical spikes and smooth out the resulting force. Spatial filtering is also usedto manage the contrast in density and viscosity between the liquid and gas phases.Furthermore, any small gas bubbles or liquid drops that have a diameter smallerthan sixteen level 13 cells (corresponding to ≈ µ m in dimensional terms) aredeemed unphysical and dynamically removed, with the exception of the entrappedgas bubble centred at r = z = 0.The simulations span 0 . . reaches r ≈ .
9, close to the edge of the plate, and the turnover point reaches r ≈ .
3. The early impact stage can be considered over long before the turnoverpoint surpasses the initial droplet radius, hence we also capture timescales beyondwhen we expect the analytical results to be valid. Each individual simulationconsisted in approximately 60 ,
000 (dynamically adapted) degrees of freedom and roplet impact onto a spring-supported plate: analysis and simulations 19 was executed in parallel on 4 − § α = 100, with the damping and stiffness factor β = γ = 0. Formesh-independence, we ran multiple simulations varying the maximum refinementlevel from 9 to 13, finding that the calculated force did not vary noticeably beyondlevel 12. In dimensional terms, this means the smallest cells must have a side lengthof at most ≈ µ m. Similarly, the width of the computational domain was variedfor L = 3, 6, 12 and 24. We found that after L = 6, the calculated force value nolonger changed to a significant degree. We are thus confident in proceeding with acomprehensive parametric study, exploring the solution space with both analyticaland computational approaches. The aim of this section is two-fold: firstly, to systematically compare the predic-tions of the analytical model from § §
4, identifying timescales during which good agreement is observed and, sec-ondly, to provide insight into the physical mechanisms introduced once substratemotion is allowed, systematically showing how the mass ratio α , damping factor β and stiffness factor γ affect the dynamics of the system. To facilitate the com-parison of the analytical and numerical results, we re-express all quantities intothe original non-dimensional variables from § t = (cid:15) ˆ t in § § − . ≤ t ≤ . F ( t ) differsfrom the corresponding value for the stationary plate case, allowing us to identifywhere the plate motion has a strong influence on the dynamics of the droplet. The analytical and numerical predictions for the hydrodynamic force F ( t ) andthe plate displacement s ( t ) for α = 2, β = 0, γ = 500 are shown in Fig. 5, alongsidethe corresponding analytical and numerical predictions for the stationary platecase. Under no damping or hydrodynamic forcing, the plate displacement s ( t ) in(12) would oscillate about s ( t ) = 0 with a natural time period T = 2 π (cid:112) α/γ ≈ t = 0 .
015 (b) t = 0 . t = 0 .
535 (d) t = 0 . Fig. 5: (Top) Comparison of the hydrodynamic force F ( t ) between the stationaryplate case (black) and a moving plate (grey) with mass ratio α = 2, damping factor β = 0, stiffness factor γ = 500, with a dashed line for the analytical solution from(59) and a solid line for the corresponding numerical value. (Middle) Displacementof the moving plate case s ( t ), with the dashed line showing the analytical solutionto (12) and the solid line depicting the corresponding numerical results. (Lowerpanels) Comparison between the pressure p of the DNS between the stationary plate case (left) and the moving plate case (right) at the times labelled 1-4 in theplots above. roplet impact onto a spring-supported plate: analysis and simulations 21 . α, γ are chosen so T is of the same order of magnitude asthe timescale of the simulation − . ≤ t ≤ . ≈ k ∗ ≈ < .
1% of the pressure, hence the dominant contribution to the numerical resultsfor the force F ( t ) was due to the pressure itself.Point 1 ( t = 0 . F ( t ) for the moving plate is close to that of the stationary plate, as the platehas only deformed to within a distance of O (10 − ), and it can be seen in thepanels in Fig. 5 that the pressure distribution for both is similar. The plate dis-places downwards until around t = 0 .
25, at which point the strength of the elasticrestoring force, γs ( t ), causes the plate to move back upwards. Shortly after this,at point 2 ( t = 0 . t = 0 .
5, when the elastic force balanceswith the hydrodynamic force and the plate begins to accelerate downwards again.The hydrodynamic force reaches a local minimum at point 3 ( t = 0 . t = 0 . t = 0 . t (cid:28)
1, and that the radius of the turnovercurve remains small compared to the droplet radius, however the panels show theturnover curve at point 2 is close to r = 0 . s ( t ) and the difference in the locationof the fluid interfaces in the graphs and snapshots shown in Fig. 5 are small incomparison to the size of the droplet. Hence, upon experimental observation, thephysical system may not appear different to the stationary plate case. However theoscillations in the hydrodynamic force and the pressure differences in the snapshotsshow that flow inside the droplet is being significantly affected by the motion of theplate. This shows that just introducing substrate motion due to linear elasticityresults in a substantial change in the dynamics of the droplet.5.2 Plate parameter comparisonsIn § α , damping factor β and stiffness factor γ . In the following we aim to study physical mechanisms represented by these parameters individually.In order to systematically observe the effects of these physical mechanisms, weconducted a series of simulations for − . ≤ t ≤ .
675 with varying values of α , β , γ . The hydrodynamic force F ( t ) was calculated regularly and is shown bythe solid greyscale lines in Fig. 6, where darker lines correspond to higher values β = 0, γ = 0 (b) α = 2, γ = 100 (c) α = 2, β = 0 Fig. 6: Numerical values from the DNS for the hydrodynamic force on the plate F ( t ), with solid grey lines showing moving plate cases with varying values of themass ratio α , damping factor β and stiffness factor γ , and solid black lines showingthe stationary plate case. The black dashed lines indicate the analytical solutionto F ( t ) for a stationary plate given by (59). (a): β = γ = 0, and α ranges from1 to 100. (b): α = 2, γ = 100 and β = 0, 0 . β c , β c and 5 β c , where the criticaldamping value β c = 2 √ αγ ≈ .
28. (c): α = 2, β = 0 and γ ranges from 0 to 1000.of α , β or γ . For comparison, the solid black line and dashed line correspond tothe numerical and analytical hydrodynamic force for the stationary plate case.Analytical solutions for the rest of the cases are not shown for visual clarity on theplots, however an analysis similar to that presented in § α , β and γ .The mass ratio α = M ∗ /ρ ∗ l R ∗ d represents (up to a constant) the ratio betweenthe mass of the plate M ∗ and the mass of the droplet 4 πρ ∗ l R ∗ d /
3. Upon impact,the pressure of the droplet exerts a hydrodynamic force onto the plate, causing itto accelerate downwards. The downwards motion of the plate causes the pressureat the surface of the plate to decrease, in turn decreasing the hydrodynamic force.For lighter plates (smaller α ), this downwards motion will be faster, and hence weexpect the hydrodynamic force will be lower in lighter plates than for heavier ones.The mass ratio α is varied from 1 to 100 in Fig. 6a, with β = γ = 0. In these cases,the only force acting on the plate is the hydrodynamic force of the droplet fromabove, hence the plate accelerates downwards at a rate depending on the mass ratio α . It can be seen from Fig. 6a that increasing α causes the overall force to increase,tending towards the stationary plate value for large values of α . In addition, thetime at which the hydrodynamic force reaches a maximum increases as α increases, which happens once the plate has accelerated to its maximum velocity, resultingin a lower hydrodynamic force. It takes longer for this to happen the heavier theplate is, hence the time at which the maximum is reached increases as α increases.The damping factor β determines the amount of resistance to motion the dash-pot exerts. We note that the ODE for the plate displacement (12) under no external roplet impact onto a spring-supported plate: analysis and simulations 23 forcing has a critical damping value of β = β c = 2 √ αγ . If this unforced systemwere displaced from its equilibrium position and released, the undamped system( β = 0) would experience oscillations of a fixed amplitude about the equilibriumpoint. For 0 < β < β c the system would be underdamped, and the amplitude of theoscillations would decay at a rate increasing with β . If β > β c , the system wouldbe overdamped and would exponentially return to equilibrium, returning moreslowly with increasing β . Finally, if β = β c , the system would be critically dampedand would return to the equilibrium in the fastest time. However inclusion of thehydrodynamic forcing will alter these dynamics, and, in particular, we expect thathigher values of β would lead to smaller displacements from equilibrium due tothe resistance to motion. In Fig. 6b, the mass ratio and stiffness factor are fixedat α = 2, γ = 100 such that the critical damping value is β c = 20 √ ≈ .
28. Thegreyscale lines show the values of force for β = 0, 0 . β c , β c and 5 β c . For β = 0,we can clearly see oscillations in the force, and the amplitude of these oscillationsdecreases when the system is underdamped for β = 0 . β c . These oscillations aresuppressed in the case of critical damping β = β c , where the force follows a trendthat is initially lower than the stationary plate value, whereas the force approachesthe stationary plate value for the overdamped case β = 5 β c . Since the force de-pends predominantly on the hydrodynamic pressure in the droplet, the fact thatthe force follows the same behaviour of under-, over- and critical damping showsthe strong influence the dashpot has on the behaviour of the droplet.The strength of the elastic force from the compression of the spring is rep-resented by the stiffness factor γ . In the absence of damping and external forc-ing, the solution of (12) for s ( t ) would oscillate with a dimensionless time period T = 2 π (cid:112) α/γ , hence an increase in the stiffness factor results in the oscillationshaving a shorter time period. The spring does work on the droplet via a verti-cal force equal to γs ( t ), so the droplet loses kinetic energy depending on how farbelow the z axis the plate is displaced. Unlike damping, the elastic force is con-servative, meaning the loss of kinetic energy from the droplet due to the elasticityis converted into potential energy in the spring, which is then in turn convertedinto kinetic energy via oscillations. Fig. 6c shows the hydrodynamic forces in sys-tems for mass ratio and damping factor α = 2 and β = 0, with stiffness factor γ varying from 0 to 1000. As expected, we observe that the time periods of theoscillations decrease as γ increases in Fig. 6c. For the two highest γ values (500and 1000), it can be seen that the values of F ( t ) oscillate centred on the forcevalue for the stationary plate case. This suggests that as γ increases, although thefrequency of oscillations increases, the amplitude of the oscillations would decreaseand eventually the force would tend to the stationary plate value.Although the analytical solution was only shown for the stationary plate case inFig. 6, it is worth noting the good agreement this solution has with the numericalvalues at early times for the majority of the moving plate cases shown. At earlytimes, the velocity of the plate is still small, hence it does not significantly alterthe hydrodynamic force. It is only once the plate has been accelerated that itsmotion affects the hydrodynamic force by doing work on the liquid. In this section, we have shown the rich variety in behaviours that the sys-tem exhibits as a result of the individual physical contributions due to the massof the plate, strength of the dashpot and stiffness of the spring. These physicalmechanisms result in previously unreported changes in the droplet dynamics, suchas pressure oscillations which can be suppressed by energy losses due to damp- ing. Although the magnitude of the plate displacement is small in comparison tothe length scale of the droplet in all cases, the strong coupling observed betweenthe plate displacement s ( t ) and the hydrodynamic force F ( t ) justify making useof models where the fluid-structure interaction is retained in order to accuratelypredict the dynamics of the droplet over these timescales. In this paper we have presented two models for the vertical impact of a droplet ontoa plate supported by a spring and a dashpot: an analytical model using matchedasymptotic expansions and a full computational framework based on DNS. Al-though droplet impact onto elastic beams has been considered previously [24], theanalytical model we present is the first to consider the Wagner theory formulationwhere the substrate experiences both elastic forcing and damping. As opposed toprevious axisymmetric models [22,24], we approximate the hydrodynamic force onthe substrate using the leading-order composite expansion of the pressure betweenthe outer and inner region (rather than just the outer region). Significantly, wefound that the composite force shown in Fig. 5 is within 10% of the numericalsolution up to t ≈ .
2, in contrast to the force contribution due to the outer regiononly remains within 10% of the numerical solution up to t ≈ .
04, which justifiesconsidering the contributions from the inner region in order to extend the timescalein which such analytical models are valid more generally. Previous numerical in-vestigations involving a plate-spring system [45] do not take into account forcingdue to damping, and focus on the late-time dynamics of spreading and rebound,whereas we focus on the influence the plate motion has on the delicate early stagesof impact in a high-speed context. Finally, the response of an elastic substrate onan impacting droplet has very recently been modelled using an effective boundarycondition on the pressure in order to consider a stationary computational domain[13]. By contrast, in our model the substrate motion is resolved by using a movingframe of reference centred on the surface of the substrate, fully representing thefluid-structure interaction.The two proposed methodologies are distinct in their approach to understand-ing the system, and yet they are stronger in combination. In a problem with suchviolent topological changes over short scales, it is vital to have analytical resultsto both validate and inform our DNS platform. The analytical model providesguidance into key quantities, such as the location of the turnover point, whichcan be used as a prediction for the simulation duration and refinement strategy.In addition, the analytical model can be used to rapidly search for parameterswhere interesting coupling between the droplet and plate can be observed, ratherthan spending considerable computational resources searching for these regimesnumerically. However, the analytical model relies on a series of assumptions, suchas neglecting viscosity, surface tension and gravity, and is limited to early impacttimes. On its own, it is impossible to assess the consequence of these assump- tions, and where they break down. By systematically comparing the analyticalpredictions to the numerical model, we can identify the regimes where these as-sumptions are valid and support the use of the analytical model as opposed tothe costly DNS. If desired, the DNS can then be used to go beyond those regimesand study timescales inaccessible to the analytical model. It is only when used in roplet impact onto a spring-supported plate: analysis and simulations 25 conjunction that the predictive power and robustness of these models reach theirfull potential.Not only have the methods presented in this paper extended existing analyticaland numerical models, they have also allowed us to provide physical insight intothe dynamics of a novel, complex multi-phase system. We recognise that the dis-placement of the plate and perturbation of the free surface of the droplet are small,hence these models provide insight into a physical regime that would otherwise bevery difficult to study experimentally. Through an extensive parameter study, weidentified the influence that the mass ratio α , damping factor β and stiffness factor γ have on the hydrodynamic force exerted by the droplet. In particular, we foundthat lighter plates (smaller α ) result in a lower value of the force; stiffer springs(higher γ ) result in oscillations of higher frequency and that resistive dashpots(higher β ) suppress the oscillations due to the elasticity from the spring.The plate-spring-dashpot system is one of the simplest models for a flexiblesubstrate, as it only allows for vertical motion. Droplet impact onto the end ofa cantilever, such as a leaf, is more complex as the bending of the beam breaksthe axisymmetry. However, as considered in [7], the deflection of the end of thebeam can be modelled using the second order differential equation (12) when thedeflection is small (hence negligible bending). Therefore the model for the substratemotion considered in this paper could provide insight into the early time dynamicsof droplet impact onto the end of cantilever beams. In addition, there is much scopeto extend these models to more complex substrates, such as elastic membranesunder tension, as previously studied experimentally [25], further guided by recentanalytical [24] and computational [23,33] progress. In conclusion, we believe thatthe proposed mathematical framework embodies productive co-development andinvestigative interplay between rigorous state-of-the-art methodologies, providinga general and highly efficient route to studying complex systems involving fluid-structure interaction in the future. Conflict of interest
The authors declare that they have no conflict of interest.
Appendix A Composite hydrodynamic force
For the analytical model, the hydrodynamic force is determined by integratingthe pressure p across the surface of the plate. As the leading-order solution in theouter region (39) diverges at the turnover point, the force contribution close tothe turnover point must be an over-estimate. Hence the leading-order compositepressure between the outer and the inner region was found in § determine the resulting composite force, where the force contribution from thesplash sheet region was determined to be negligible. When viscosity is neglected,the dimensionless hydrodynamic force on the plate (10) is given by F ( t ) = 2 π (cid:90) R p rp ( r, − (cid:15) s ( t ) , t ) d r . (64) The composite force is found by integrating the composite pressure (58) alongthe surface of the plate, which is determined by splitting the range of the integralacross the respective asymptotic regions, F comp ( t ) = F outer ( t ) + F inner ( t ) − F overlap ( t ) , (65)where F outer ( t ) and F inner ( t ) are the result of integrating the leading-order outerpressure (39) and inner pressure (51) respectively, and F overlap ( t ) is the result ofintegrating the overlap pressure (57). The resulting force due to the outer regionis hence, F outer ( t ) = 2 π (cid:90) (cid:15)d ( t )0 r (cid:15) ˆ p ( r/(cid:15), , t ) d r = 89 (cid:15) (cid:90) d ( t )0 d d t (cid:104) ˆ r ( d ( t ) − ˆ r ) / (cid:105) dˆ r = 89 (cid:15) d d t (cid:90) d ( t )0 ˆ r ( d ( t ) − ˆ r ) / dˆ r = 845 (cid:15) d d t (cid:104) d ( t ) (cid:105) = 89 (cid:15)d ( t ) (4 ˙ d ( t ) + ¨ d ( t ) d ( t )) . (66)For F inner ( t ), the integration variable is changed to the parameter η from (51).Then, η ( t ) is defined such that r = 0 for η = η ( t ), i.e. e η ( t ) + 4 e η ( t ) + 2 η ( t ) + 1 = πd ( t ) (cid:15) J ( t ) , (67)where η ( t ) must be solved for numerically for each t . Note that η → −∞ as r → ∞ , where ˜ p decays exponentially. Hence we take the upper limit of theintegral to be r = ∞ , which introduces exponentially small errors. The resultingforce due to the inner region is hence F inner ( t ) = 2 π (cid:90) ∞ r (cid:15) ˜ p ( r/(cid:15) − d ( t ) /(cid:15) , , t ) d r = 2 (cid:15)π (cid:90) ∞− d ( t ) /(cid:15) ( (cid:15)d ( t ) + (cid:15) ˜ r )˜ p (˜ r, , t ) d˜ r = 8 (cid:15) ˙ d ( t ) J ( t ) (cid:90) η −∞ (cid:18) (cid:15)d ( t ) − (cid:15) J ( t ) π ( e η + 4 e η + 2 η + 1) (cid:19) e η d η = 8 (cid:15) ˙ d ( t ) J ( t ) π e η ( t ) (cid:20) πd ( t ) (cid:15) J ( t ) + 1 − e η ( t ) − e η ( t ) − η ( t ) (cid:21) . (68)Finally the overlap force is F overlap ( t ) = 2 π √ d ( t ) / ˙ d ( t ) π √ (cid:15) (cid:90) (cid:15)d ( t )0 r (cid:112) (cid:15)d ( t ) − r d r = 16 √ (cid:15)d ( t ) ˙ d ( t ) . (69)Combining (66), (68) and (69), the composite force is F comp ( t ) = 89 (cid:15)d ( t ) ((4 − √
2) ˙ d ( t ) + ¨ d ( t ) d ( t ))+ 8 (cid:15) ˙ d ( t ) J ( t ) π e η ( t ) (cid:20) πd ( t ) (cid:15) J ( t ) + 1 − e η ( t ) − e η ( t ) − η ( t ) (cid:21) , (70)where d ( t ) and J ( t ) are given by (38) and (49) respectively. roplet impact onto a spring-supported plate: analysis and simulations 27 References
1. Antonini, C., Amirfazli, A., Marengo, M.: Drop impact and wettability: From hydrophilicto superhydrophobic surfaces. Physics of Fluids (10), 102104 (2012)2. Blowers, R.M.: On the response of an elastic solid to droplet impact. IMA Journal ofApplied Mathematics (Institute of Mathematics and Its Applications) (2), 167–193 (1969)3. van Brakel, J.P.: Robust peak detection algorithm (using z-scores). URL https://stackoverflow.com/questions/22583391/peak-signal-detection-in-realtime-timeseries-data
4. Cimpeanu, R., Moore, M.R.: Early-time jet formation in liquid-liquid impact problems:theory and simulations. Journal of Fluid Mechanics , 764–796 (2018)5. Eggers, J., Fontelos, M.A., Josserand, C., Zaleski, S.: Drop dynamics after impact on asolid wall: Theory and simulations. Physics of Fluids (6) (2010)6. Ellis, A.S., Smith, F.T., White, A.H.: Droplet Impact on to a Rough Surface. The Quar-terly Journal of Mechanics and Applied Mathematics (2), 107–139 (2011)7. Gart, S., Mates, J.E., Megaridis, C.M., Jung, S.: Droplet Impacting a Cantilever: A Leaf-Raindrop System. Physical Review Applied (4), 044019 (2015)8. Gilet, T., Bourouiba, L.: Rain-induced ejection of pathogens from leaves: Revisiting thehypothesis of splash-on-film using high-speed visualization. Integrative and ComparativeBiology (6), 974–984 (2014)9. Hicks, P.D., Purvis, R.: Air cushioning and bubble entrapment in three-dimensional dropletimpacts. Journal of Fluid Mechanics , 135–163 (2010)10. Hoath, S.D. (ed.): Fundamentals of Inkjet Printing. John Wiley & Sons, Weinheim, Ger-many (2016)11. Howison, S.D., Ockendon, J.R., Oliver, J.M., Purvis, R., Smith, F.T.: Droplet impact ona thin fluid layer. Journal of Fluid Mechanics , 1–23 (2005)12. Howison, S.D., Ockendon, J.R., Wilson, S.K.: Incompressible water-entry problems atsmall deadrise angles. Journal of Fluid Mechanics (-1), 215–230 (1991)13. Howland, C.J., Antkowiak, A., Castrej´on-Pita, J.R., Howison, S.D., Oliver, J.M., Style,R.W., Castrej´on-Pita, A.A.: It’s Harder to Splash on Soft Solids. Physical Review Letters (18), 184502 (2016)14. Josserand, C., Thoroddsen, S.: Drop Impact on a Solid Surface. Annual Review of FluidMechanics (1), 365–391 (2016)15. Korobkin, A.A.: Formulation of penetration problem as a variational inequality. DinamikaSploshnoi Sredy , 73–79 (1982)16. Korobkin, A.A.: Wave impact on the center of an Euler beam. Journal of Applied Me-chanics and Technical Physics (5), 770–781 (1998)17. Li, E.Q., Thoraval, M.J., Marston, J.O., Thoroddsen, S.T.: Early azimuthal instabilityduring drop impact. Journal of Fluid Mechanics , 821–835 (2018)18. Li, L., Sherwin, S.J., Bearman, P.W.: A moving frame of reference algorithm forfluid/structure interaction of rotating and translating bodies. International Journal forNumerical Methods in Fluids (2), 187–206 (2002)19. L´opez-Herrera, J.M., Popinet, S., Castrej´on-Pita, A.A.: An adaptive solver for viscoelas-tic incompressible two-phase problems applied to the study of the splashing of weaklyviscoelastic droplets. Journal of Non-Newtonian Fluid Mechanics , 144–158 (2019)20. Moore, M.R.: New mathematical models for splash dynamics. Ph.D. thesis, University ofOxford (2014)21. Murphy, D.W., Li, C., D’Albignac, V., Morra, D., Katz, J.: Splash behaviour and oilymarine aerosol production by raindrops impacting oil slicks. Journal of Fluid Mechanics , 536–577 (2015)22. Oliver, J.M.: Water entry and related problems. Thesis, Oxford University p. 150 (2002)23. Patel, H.V., Das, S., Kuipers, J.A., Padding, J.T., Peters, E.A.: A coupled Volume ofFluid and Immersed Boundary Method for simulating 3D multiphase flows with contactline dynamics in complex geometries. Chemical Engineering Science , 28–41 (2017)24. Pegg, M., Purvis, R., Korobkin, A.: Droplet impact onto an elastic plate: a new mechanismfor splashing. Journal of Fluid Mechanics , 561–593 (2018)25. Pepper, R.E., Courbin, L., Stone, H.A.: Splashing on elastic membranes: The importanceof early-time dynamics. Physics of Fluids (8), 082103 (2008)26. Peskin, C.S.: The immersed boundary method. Acta Numerica (January 2002), 479–517(2002)8 Michael J. Negus et al.27. Philippi, J., Lagr´ee, P.Y., Antkowiak, A.: Drop impact on a solid surface: short-time self-similarity. Journal of Fluid Mechanics , 96–135 (2016)28. Popinet, S.: Gerris: A tree-based adaptive solver for the incompressible Euler equations incomplex geometries. Journal of Computational Physics (2), 572–600 (2003)29. Popinet, S.: A quadtree-adaptive multigrid solver for the Serre-Green-Naghdi equations.Journal of Computational Physics , 336–358 (2015)30. Scolan, Y.M., Korobkin, A.A.: Three-dimensional theory of water impact. Part 1. InverseWagner problem. Journal of Fluid Mechanics , 293–326 (2001)31. Smith, D.B., Askew, S.D., Morris, W.H., Shaw, D.R., Boyette, M.: Droplet size and leafmorphology effects on pesticide spray deposition. Transactions of the American Societyof Agricultural Engineers (2), 255–259 (2000)32. Sneddon, I.N.: Mixed boundary value problems in potential theory. North-Holland (1966)33. Sun, X., Sakai, M.: Numerical simulation of two-phase flows in complex geometries by usingthe volume-of-fluid/immersed-boundary method. Chemical Engineering Science , 221–240 (2016)34. Thoraval, M.J., Takehara, K., Etoh, T.G., Popinet, S., Ray, P., Josserand, C., Zaleski, S.,Thoroddsen, S.T.: Von K´arm´an vortex street within an impacting drop. Physical ReviewLetters (26), 264506 (2012)35. Thoroddsen, S.T., Etoh, T.G., Takehara, K., Ootsuka, N., Hatsuki, Y.: The air bubbleentrapped under a drop impacting on a solid surface. Journal of Fluid Mechanics (-1),203–212 (2005)36. Thoroddsen, S.T., Takehara, K., Etoh, T.G.: Micro-splashing by drop impacts. Journalof Fluid Mechanics , 560–570 (2012)37. Van Dyke, M.: Perturbation methods in fluid mechanics. Parabolic (1975)38. Visser, C.W., Frommhold, P.E., Wildeman, S., Mettin, R., Lohse, D., Sun, C.: Dynamics ofhigh-speed micro-drop impact: Numerical simulations and experiments at frame-to-frametimes below 100 ns. Soft Matter (9), 1708–1722 (2015)39. Wagner, H.: ¨Uber Stoß- und Gleitvorg¨ange an der Oberfl¨ache von Fl¨ussigkeiten. ZAMM- Zeitschrift f¨ur Angewandte Mathematik und Mechanik (4), 193–215 (1932)40. Weisensee, P.B., Tian, J., Miljkovic, N., King, W.P.: Water droplet impact on elasticsuperhydrophobic surfaces. Scientific Reports (1), 30328 (2016)41. Wildeman, S., Visser, C.W., Sun, C., Lohse, D.: On the spreading of impacting drops.Journal of Fluid Mechanics , 636–655 (2016)42. Wilson, S.K.: The mathematics of ship slamming. Ph.D. thesis, University of Oxford(1989)43. Wilson, S.K.: A mathematical model for the initial stages of fluid impact in the presenceof a cushioning fluid layer. Journal of Engineering Mathematics (3), 265–285 (1991)44. Xiong, Y., Huang, H., Lu, X.Y.: Numerical study of droplet impact on a flexible substrate.Physical Review E (5), 1–15 (2020)45. Zhang, C., Wu, Z., Shen, C., Zheng, Y., Yang, L., Liu, Y., Ren, L.: Effects of eigenand actual frequencies of soft elastic surfaces on droplet rebound from stationary flexiblefeather vanes. Soft Matter16