Stress distribution and surface shock wave of drop impact

Drop impact causes severe surface erosion, dictating many important natural, environmental and engineering processes and calling for substantial prevention and preservation efforts. Nevertheless, despite extensive studies on the kinematic features of impacting drops over the last two decades, the dynamic process that leads to the drop-impact erosion is still far from clear. Here, we develop a method of high-speed stress microscopy, which measures the key dynamic properties of drop impact responsible for erosion, i.e., the shear stress and pressure distributions of impacting drops, with unprecedented spatiotemporal resolutions. Our experiments reveal the fast propagation of self-similar noncentral stress maxima underneath impacting drops and quantify the shear force on impacted substrates. Moreover, we examine the deformation of elastic substrates under impact and uncover impact-induced surface shock waves. Our study opens the door for quantitative measurements of the impact stress of liquid drops and sheds light on the origin of low-speed drop-impact erosion.


Continuous displacement field
We implemented the moving least squares (MLS) interpolation method to obtain a continuously differentiable displacement field from a discrete DIC displacement field of experiments [1]: where u(x) = [u r , u z ] is the target continuous displacement field, P T (x) is a polynomial basis and a(x) is the corresponding coefficients of the basis. x is the continuous position vector. Because of the cylindrical symmetry of the dropimpact geometry, x = (r, z) with r = 0 and z = 0 at the initial impact point.
The negative z direction points along the direction of the impact velocity of drops. We adopted a cubic basis P T = [1, r, z, r 2 , z 2 , rz, r 3 , z 3 , r 2 z, rz 2 ]. The coefficients of the target function at x was obtained by minimizing the weighted least-square error, where w is the discrete DIC displacement field from experiments. b i indicates the discrete DIC coordinates with a total n = 540 points in our experiments.
is the weighted function defined as [2]: . m determines the fitting range, which was fixed at m = 43 in this study. The minimization of Eq. (2) gives where 2 Impact pressure PDMS gels are nearly incompressible with Poisson's ratio close to 0.5, which leads to a large λ. On the other hand, the bulk strain ε b is close to 0. Therefore, the impact pressure cannot be accurately determined from the product of λε b in the linear constitutive equation To overcome the difficulty, we adopted the quasi-steady-state assumption [1].
Specifically, λε b is replaced by an undetermined term −µ in Eq. (5), Under the quasi-steady-state assumption, the stress tensor obeys: As a result, µ can be calculated as: where x 0 sets a reference point for integration. We chose the reference point at the top layer z = 0 with r as large as possible away from the impact point at short times. Taking p(x 0 ) = σ zz (x 0 ) = 0, we have µ(x 0 ) = 2Gε zz (x 0 ). Finally, the pressure can be obtained Here, the integration path goes from the reference point x 0 to x, as shown in Supplementary Fig. 1.
To verify the quasi-steady-state assumption, we compared the ratio of the inertial force to the elastic force per unit volume for solid-sphere impact: where ρ s is the density of the PDMS gels. R is less than about 2% near the impact axis where the impact pressure is above the noise level. Consistent with the estimate, the impact pressure of solid spheres calculated based on the quasi-steady-state assumption quantitatively matches that of the finite-element simulations without the assumption (Fig. 2b of the main text). As the time scale of solid-sphere impact is much shorter than that of drop impact under comparable impact conditions, the quasi-steady-state assumption should work even better for drop impact. Lastly, it is worth of noting that the calculation of shear stress-the quantity most relevant to surface erosion-does not rely upon the quasi-steady-state assumption and therefore is immune to the errors associated with the assumption.

Theory and numerical solutions
The deformation of an elastic medium at small strains is described by the Navier-Lamé equation: the shear modulus. ρ s , E and ν are the density, Young's modulus and Poisson's ratio of the medium, respectively. u = (u r , u θ , u z ) is the displacement of the medium in a cylindrical coordinate. F is the body force per unit volume. Since F is small compared with the impact stress, we ignore the term. In the cylindrical coordinate, an axisymmetric Navier-Lamé equation reduces to  12) and (13): where ρ, U and D are the density, the impact velocity and the diameter of the liquid drop, respectively. r t = √ 6U Dt/2 is the position of the turning point, where the pressure and shear stress diverge and exhibit the finite-time singularity [3]. Note that the shear stress σ rz | z=0 is a factor 1/ √ Re ∼ 7 × 10 −3 smaller than the impact pressure σ zz | z=0 . Thus, the impact pressure is the leading factor controlling the deformation of the impacted substrate ( Fig. 4d-f and the associated discussion in the main text).
To make the equations dimensionless, we propose the following scaling, t ∼ The proposed scaling yields the following boundary value problem of partial differential equations: and where all the quantities are now dimensionless with the turning point at r t = √ 6t/2. The resulting equations and boundary conditions are independent of M in terms of the scaled variables. The scaling suggests that the shear force should scale as F d ∼ σD 2 ∼ E 1/2 (ρU 2 ) 1/2 D 2 , matching our experimental measurements at different E (Fig. 5b of the main text). Note that the density ratio between the liquid and the gel ρ/ρ s ∼ O(1) in our experiments.
We numerically solve Eqs. (16) and (17) with the boundary conditions Eqs. (18) and (19) using the finite element method (Methods). To avoid the singularity at r t , we impose a small cut-off δ: and We choose δ = 0.1 for our numerical simulations, as a good convergence of solutions is achieved for the chosen spatial resolution. To compare with the surface wave induced by the impact of a solid sphere, we apply the Hertzian contact force in the dimensionless form [4] F = 2 √ 2 3 on the surface of the elastic medium over a small region of r 0 = 0.1 around the impact axis. Here, we assume that the deformation of the solid sphere is negligible compared with that of the elastic medium and the vertical displacement x ≈ U t at short times based on the experimental observation [5]. In comparison with the impact stress of liquid drops (Eqs. (18) and (19)), Eq. (22) is stationary in space and does not exhibit either the finite-time singularity or the shock front associated with the turning point at small t. With a stationary