Nanosecond shock wave-induced surface acoustic waves and dynamic fracture at fluid-solid boundaries

We investigate the generation and propagation characteristics of leaky Rayleigh waves (LRWs) caused by a spherical shock wave incident on a water-glass boundary both experimentally and numerically. The maximum tensile stress produced on the solid boundary is attributed to the dynamic interaction between the LRWs and an evanescent wave generated concomitantly along the boundary. The resultant tensile stress field drives the initiation of crack formation from pre-existing surface flaws and their subsequent extension along a circular trajectory, confirmative with the direction of the principal stress on the boundary. We further demonstrate that this unique ringlike fracture pattern, prevalent in damage produced by high-speed impact, can be best described by the Tuler-Butcher criterion for dynamic failure in brittle materials. The orientation of the ring fracture extension into the solid also follows closely with the trajectory of the local maximum tensile stress distribution.


I. INTRODUCTION
Surface acoustic waves (SAWs), such as Rayleigh wave and leaky Rayleigh wave (LRW) [1], have been widely used in nondestructive testing, telecommunication, and microfluidic applications [2][3][4][5][6]. SAWs are also prevalent in many destructive nature phenomena, such as earthquakes and tsunami [7]. Further, SAWs have long been speculated to contribute to surface damage produced by cavitation on ship propellers [8], high-speed impact on wind turbine blades [9], supersonic flights through rain forests [10], and fragmentation of rocks by jetting streams [8,10,11]. Despite this, the fundamental mechanism whereby SAWs are generated under various shock wave or impact conditions is largely unknown, and the correlation between SAW, resultant stress field, and damage produced in the solid has not been well established [12][13][14]. More recently, in medical applications such as shock wave lithotripsy (SWL), SAWs have been postulated to play a vital role in ensuring the success of noninvasive disintegration of kidney stones in patients [15]. Altogether, a fundamental knowledge of the mechanism underpinning SAW generation, propagation characteristics, and dynamic interaction with solid materials will be critical for either preventing or harnessing SAW-induced fracture under impulsive shock loadings [16].
In this work, we present new insights into the physical processes responsible for the generation and propagation of LRW and resultant damage produced by a nanosecond sparkgenerated shock wave in water near a glass boundary. First, we combine experimental and computational techniques to dissect the relationship between various stress waves excited by the incident shock wave at the fluid-solid boundary and the resultant transient stress field generated in the solid. Second, we establish an experimental system to create reproducible ringlike fracture patterns on glass surfaces under different shock wave conditions. Third, we characterize the surface crack extension into the solid by confocal microscopy and compare the results with the stress field distribution in the solid. The uniqueness of this experimental system to initiate, under well-controlled experimental conditions, and drive the progression of ringlike fractures, prevalent in damage produced by high-speed impact, overcomes many of the primary limitations in previous studies [9][10][11][12]. The results provide us for the first time a mechanistic understanding of the characteristic ringlike fractures that is consistent with the prediction by the Tuler-Butcher criterion for dynamic failure in brittle materials.

Dynamic photoelastic/shadowgraph imaging
A dynamic photoelastic/shadowgraph imaging system [17] was used to capture the generation and propagation of the shock waves in water produced by an 3.6 Fr. Nano Pulse lithotripsy (NPL) probe (Lithotech Medical, Israel) operated at 10 kV, and the subsequent shock wave interaction with a glass sample (50 × 50 × 3.3mm in LxWxH, BOROFLOAT® Borosilicate Windows, Edmund Optics Inc.) at different standoff distances. The transient shock wave-glass sample interaction was recorded by using an ultrahigh-speed camera (Kirana-M5, Specialised Imaging) operated at 5 million frames per second (fps) with backlighting provided by a 10-ns pulsed laser illumination system (SI-LUX-640, Specialised Imaging). The imaging system was triggered by the output signal from a photo-detector (DET110, THORLABS), which detected the spark discharge produced by the NPL probe.

Numerical simulations
The numerical model was constructed in COMSOL Multiphysics 5.3 (Burlington, MA) using the Acoustic-Solid Interaction, Transient Multiphysics Interface provided in the Acoustic module. An axisymmetric model was considered with the z axis of the cylindrical coordinate system aligned with the NPL probe axis [18]. Both the fluid and solid domains were discretized by the second order triangular Lagrangian elements. The mechanical properties of the borosilicate glass were taken from the product (BOROFLOAT®, SCHOTT North America, Inc) website. The detailed description of the model and calibration against a known analytical solution of the acoustic waves generated by an impulsive point source near a flat fluid-solid boundary [19] are shown in [18]. Figure 1(a) shows schematically the experimental setup (left) and a finite element multiphysics model simulation in COMSOL (right), in which a nanosecond spark discharge was used to create a spherically divergent shock wave in water. The expanding shock front in water was captured by high-speed shadowgraph imaging and the pressure waveforms along the critical incident angle for LRW (θ LRW = ~28° from the probe axis) were measured at different radial distances (r) using a fiber-optic probe hydrophone (FOPH-500, RP acoustics). As shown in Fig. 1(b), the leading shock front has a rise time (t r ) about 40 ns, followed by a compressive wave of ∼300 ns in pulse duration (t p ), and a peak pressure (p+) approximately scaled inversely with r, which is characteristic of a point source. The transient shock wave interaction with a borosilicate glass plate (50 × 50 × 3.3mm in L×W×T) and various resultant stress waves were visualized by photoelastic imaging (movie 1) and recorded by using an ultrahigh-speed camera (Kirana-M5, Specialised Imaging). To vary the incident pressure of the shock waves, the standoff distance (S d ) between the spark source and the glass surface was adjusted using a 3D translational stage. In this work, p+ varies in the range of 168 MPa at S d = 0.5mm to 28 MPa at S d = 3.0mm.

III. RESULTS AND DISCUSSION
Theoretical analysis pertinent to this problem may be carried out using the conventional approach based on decomposition of a spherical wave into harmonic plane waves by applying a Fourier transform with respect to time, followed by a Fourier-Bessel transform in cylindrical coordinate system [20]. Alternatively, the problem can be solved using the modified Cagniard's technique, widely adapted in geophysics analyses [21][22][23]. For example, the reflected pressure waves in the fluid (p r ) can be determined by a convolution of the input signal of the source and the space-time Green's function of the boundary configuration [19]: where ρ f is the fluid density, ϕ V (t − τ) is the source strength, and G f r ( ⋅ , τ) is the time Laplace transform of the space-time Green's function. In both methods, however, the solution has to be obtained numerically for each point in the field, although the latter approach is computationally more efficient [19]. Therefore neither approach is practical for stress field characterization due to the extensive computational times required to map out the entire field near the NPL probe tip. To circumvent this obstacle, we model the spark-induced shock wave as a monopole source [24], based on the pressure waveforms measured in water, and numerically calculate the resultant shock wave-glass interaction using a finite element multiphysics model constructed in COMSOL. We have validated the COMSOL model by comparison with the de Hoop's analytical solution [23]. Since the spark-generated shock wave is spherically divergent with the pressure amplitude decreasing inversely with the propagation distance, and furthermore, part of the LRW energy will be reradiated constantly into the fluid during propagation, nonlinear effects during the shock wave-solid interaction in NPL is presumed to be insignificant. In fact, it is generally believed that linear elastic model is sufficient to analyze the stress field produced in the target stones during SWL [25].
The generation and propagation of LRW by the incident shock wave on the water-glass surface was confirmed both experimentally and numerically [ Fig. 1(c) and also see movie S1]. Excellent agreements were observed regarding the wave front positions of the incident (P i ) and reflected (P r ) shock (or pressure) waves in water, the transverse (T) wave in the glass, and the LRW, which is characterized by a leading transverse branch (LRW T ) jointed asymptotically by a trailing longitudinal branch (LRW L ) on the boundary [26]. The experimentally estimated LRW speed (c LRW = 3028m/s) matches well with the theoretical value (3163 m/s) at the water-glass boundary [1]. While propagating on the boundary, the LRW reradiated part of its energy as a Schmidt head wave [20] back into the water, connecting tangentially to the P r wave front (i.e., the leaky pressure wave front was aligned at about θ LRW with respect to the boundary normal).
Next, we quantified the transient stress field excited by the incident shock wave inside the glass using the COMSOL model. In particular, the variations in the first principal stress (σ T ) profile along five different refraction rays of the T wave at progressively increased radial distances from the boundary were plotted in Fig. 2. As the refraction angle (θ r ) increased from 0° (ray 1) to 5° (ray 2), we observed a transition from compressive stress produced by the longitudinal (L) wave (c L = 5498m/s) to tensile stress generated by the T wave (c T = 3479m/s). Note that at θ r = 0° the refracted L and T rays overlap with each other. As θ r increased further to 42° (ray 3), the tail ends of the two LRW branches started to emerge.
The leading LRW T had a shorter pulse duration than the trailing LRW L ; yet both were comparable in σ T magnitude to that produced by the T-wave sandwiched between them. At θ r = 85° (ray 4), the dominance of LRW-induced σ T with significantly increased magnitude became clear when the two branches gradually merged and obscured the T-wave-induced σ T . Finally, the LRW T and LRW L jointed together at the water-glass boundary (ray 5), producing the maximum σ T (i.e., σ T,max ) an order of magnitude stronger than its counterpart induced by the T wave inside the glass at small θ r . Near the boundary (rays 4 and 5), the σ T profile consists of a leading tensile component followed by a compressive tail. It is also worth noting that during convergence, the direction of the σ T produced by LRW T and LRW L branches rotate in opposite directions and eventually align with σ rr on the surface along the tangential direction of the boundary.
Figure 3(a) shows the depth-dependent variations of the six stress components produced by the LRW near the glass surface at r = 0.5mm under S d = 0.5mm where the global peak in σ T,max is produced. On the boundary σ rr and σ ϕϕ are the only nonzero components with σ rr (= σ T on the boundary) significantly greater than σ ϕϕ , indicating that the solid surface is under biaxial tension with the σ T aligned exclusively with σ rr . As the depth (in the −z direction) increases, σ rr decays exponentially and converts from tension to compression at about 0.08 λ 0 where λ 0 = c LRW t p . Within the same range of shallow depth, σ ϕϕ decreases monotonically while σ zz and σ rz increase from zero to a maximum, which is still yet an order of magnitude lower than the peak value of σ rr , before a gradual decay. These characteristic features of the stress field are consistent with those produced by a surface acoustic wave [1]. Figure 3(b) shows the changes in the temporal profiles of σ T and the corresponding pressure in the fluid near the boundary at two different radial distances. At S d = 0.5mm, the leading tensile peak in σ T first appears immediately after r = r i = 0.16mm, which is smaller than r LRW (= 0.27mm) -the corresponding radial distance when the shock wave incident angle (θ i ) equals to θ LRW . In addition, the spatial peak of the maximum tensile stress σ T,max is achieved at r = r m = 0.5mm, far exceeding r LRW . In comparison, at S d = 3.0mm, r i , r LRW and r m shift to 0.89, 1.62, and 2.2 mm, respectively. Therefore the σ T,max produced by the incident shock wave on the solid boundary varies significantly with S d both in terms of amplitude and radial position. Overall, however, the radial distribution of σ T,max at different S d [ Fig. 3(c)] shows an initial rapid buildup to a peak value especially at small S d , followed by a gradual decay along the propagation direction.
Several critical features can be more clearly observed in Fig. 4(a) when σ T,max is normalized by P i *, the peak acoustic pressure incident on the water-glass boundary at θ LRW , and plotted against the normalized radial distance r S d (or equivalently θ i shown in the top horizontal axis). First, the initiation of a leading tensile pulse in σ T appears at θ i = ∼17°, slightly above the critical incident angle for the L wave θ L * = 15°, independent of S d . Second, as S d increases, the peak of σ T, max P i * (appeared beyond θ LRW ) will increase exponentially in amplitude while shifting asymptotically toward θ LRW . Third, beyond the peak, σ T, max P i * decays initially at different rates depending on S d (e.g., r −1.63 at S d = 3mm). However, as θ i becomes greater than ∼60° (i.e., in the far field of the source region for the generation of the LRW), the decay of σ T, max P i * converges to the same rate of r −1.12 , independent of S d . It is worth noting that all these decay rates are greater than the corresponding decay rate of r −0.5 associated with Rayleigh waves produced by a point source at an air-solid boundary [27]. This result is consistent with the gradual reradiation (or loss) of LRW energy by the Schmidt head wave into the surrounding fluid [20,28]. Also importantly, the characteristic dual branch feature of the LRW first appears at about θ i = 26° [ Fig. 4(b)], slightly beyond the critical incident angle for the T wave θ T * = 25°. Therefore it is clear that the first tensile peak in σ T is initiated by the T wave at θ i = 17°, and boosted subsequently by the ensuing LRW. Altogether, these results suggest that θ LRW represents a resonant condition [1], around which LRW of varying launching efficiency will be generated by the sweeping incident spherical shock wave to gradually build up the intensity. As S d becomes larger, the projection area within the resonant θ LRW bandwidth will increase [29], leading to a stronger LRW buildup through more efficient superposition of individual wavelets and thus higher peak of σ T, max P i * . This interpretation is consistent with the theory that for a plane wave incidence (equivalent to S d → ∞ in our experiemnt) σ T, max P i * will approach to infinity at θ LRW [20], as shown by the asymptotic dash line in Fig. 4(a).
Equally important, the modeling results also reveal that a compressive pressure wave is produced by the advancing evanescent wave (EW) along the boundary on the fluid side when θ i exceeds θ T * [see Fig. 3(b)]. The speed of the advancing EW on the boundary ( = c 0 sinθ i ) [29] decreases progressively with θ i , reaching eventually to the sound speed in water (c 0 ) at glancing incidence. Therefore, within the resonant bandwidth of θ LRW two simultaneously generated waves (i.e., LRW and EW) will initially superimpose destructively with each other until the faster LRW has advanced sufficiently to move out of the compressive region covered by the EW [see Fig. 3

Author Manuscript
Author Manuscript

Author Manuscript
Author Manuscript interaction is manifested by a rapid decay of σ T, max P i * immediately after the peak, compared to a gradual decay downstream when the characteristics of the LRW propagation on the boundary become clear without the influence of the EW. It is worth noting that since the propagation distance of the newly generated LRW on the boundary increases with S d , more high frequency components of the LRW will be reradiated into the fluid [20,28], leading to a longer rise time in the resultant σ T pulse at large S d , as shown in Fig. 3(b).
Exposures of the borosilicate glass plates to NPL-generated nanosecond shock waves produce concentric, ringlike fracture patterns on the glass surface, not directly underneath the probe, but at progressively increasing radial distances with S d . The ring cracks produced by NPL bear some qualitative similarities to the cone cracks formed by the Hertzian contact [30,31]. However, they differ significantly in the mechanism and the loading rate under which the cracks are formed. Figure 5(a) shows a representative example at S d = 1.5mm, in which a circumferential crack was first observed after the 18th shock at a radial distance R 1 = 1.6mm from the projected center of the probe axis through the glass surface. With sufficient tensile stress buildup in this region [see Fig. 3(c)], the fracture was likely initiated from a pre-existing surface flaw of critical size, which grew stepwise by the coalescence of microcracks to form a macrocrack as described by the cohesive zone model of brittle failure under impulsive shock loading [32,33].The initiated crack extended in both directions by a short arc length on the order of λ 0 [10]. It has been shown that by removing pre-existing surface flaws using an etching process with 10% hydrofluoric acid, the ring crack formation produced by high-speed liquid droplet impact on glass surfaces can be effectively eliminated [10,11]. After the initiation, the crack propagated circumferentially along a trajectory orthogonal to σ T , and progressively developed into a ringlike fracture through dynamic fatigue after the 39th shock. During this period, a second circumferential crack was also observed after the 23th shock at a smaller radius (R 2 = 1.4mm), which then propagated along a different circular trajectory (movie 2). Altogether, these observations suggest that while the pre-existing surface flaws dictated the location of fracture initiation, the extension of each crack was driven primarily by the local σ T (= σ rr ) via the first mode of brittle fracture [34]. The circumferential crack propagation pattern is consistent with the principle of maximum hoop stress (or tensile stress in this case) or local symmetry [33,35], which favors crack extension along arbitrary paths with a symmetric stress distribution near their tips. In addition, at small S d (= 0.5mm) where σ T was sufficiently high over a large area of the glass surface [see Fig. 3(c)], multiple fracture rings were observed at different radii following the shock wave treatment, presumably due to the existence of multiple randomly distributed critical flaws on the glass surface. In contrast, at large S d (= 2.5 mm) where σ T became sufficiently low, only a single ring fracture was observed within 30 shocks. At further increased S d (⩾ 3.0mm), no fractures could be produced even after 100 shocks, suggesting that the tensile stress field was below a critical threshold to initiate fracture. Overall, as S d became larger, the radii of the ringlike fractures increased [ Fig. 5(b)], while the width of the fractured area and number of fracture rings would decrease [ Fig. 6(a)]. Moreover, while the mean radius of the ring cracks normalized by S d remains almost unchanged, the number of shocks (NS) required to initiate a ring crack decreases significantly with the incident pressure (P i ) on the solid boundary [ Fig. 6(b)].
We have further calculated a stress integral (SI) value in the fractured area of the glass sample based on the Tuler-Butcher criterion for brittle material failure under dynamics loading [36] for different S d (Fig. 7). SI = ∫ 0 T σ T − σ 0 2 dt for σ T > σ 0 (2) in which σ 0 = 44 MPa, is taken as the practical tensile strength of glass samples, 1 and T = 6 μs is the simulation running time. Physically, SI scales with the strain energy density impulse in the deformed material [37]. The result indicates the existence of a nominal SI threshold, corresponding to the peak value produced at S d = 3.0mm (i.e., 2.4 × 10 −4 MPa 2 s). All the observed ring-fractures were produced in regions with SI values exceeding this threshold. This finding is consistent with the Tuler-Butcher criterion for dynamic fracture of brittle materials under high strain-rate loading, which implies that both the amplitude and duration of the tensile stress pulse play a critical role in initiating the fracture [36,38,39]. In contrast, no damage threshold can be identified based solely on the magnitude of σ T,max produced at S d = 3.0mm [see Fig. 3(c)].
Moreover, the increased magnitude with extended range of radial distance above the SI threshold at small S d is consistent with the higher number of ring fractures and broader damage zone observed experimentally [Figs. 5(b) and 6(a)]. Furthermore, above the damage threshold the number of shocks (NS) required to initiate the ring-type fracture were found to be inversely correlated with SI [ Fig. 7(b)], which is characteristic of dynamic fatigue processes through coalescence of adjacent microcracks driven by impulsive tensile stress pulses [32,34]. Finally, confocal microscopy examinations of the ring fractures revealed that the circumferential (or ring) cracks were initiated on the surface along the boundary normal direction, and then extended into the glass at an inclination angle θ in the range of 41° to 63° at various radial distances (Fig. 8). The trajectory of the ring fracture extension was found to first decrease and then increase with radial distance-a feature that is consistent with the zaxis variation of LRW L , which is behind the LRW front with a broader pulse duration than LRW T . This observation suggests that LRW L is likely to play a dominant role in determining the direction of the crack propagation. It is worth noting that a large specimen size (50mm × 50mm × 3.3mm) was chosen in this work to eliminate the influence of wave reflection and interaction near the boundaries. More complex stress wave interaction and resultant fracture patterns in small-sized samples with different geometry and material properties, and their implications to stone fragmentation in NPL will be investigated in the future.

IV. CONCLUSIONS
In conclusion, we have analyzed the initiation and propagation characteristics of LRW on a water-glass boundary produced by a nanosecond spherical shock wave. The maximum tensile stress generated on the glass surface was attributed to the LRW interaction with a concurrently generated evanescent wave in a frequency and standard distance dependent (a) Schematic diagram (left) and COMSOL simulation results (right), illustrating nano pulse lithotripsy (NPL)-induced shock wave interaction at a water-glass boundary and various pressure and elastic stress waves generated in water and glass, respectively.P i : incident shock wave; P r : reflected shock wave; L: longitudinal wave; T: transverse wave; LRW: leaky Rayleigh wave. (b) Pressure waveform of the NPL-generated shock wave measured by a fiber-optic probe hydrophone (FOPH) at about 30° (which is near the critical incident angle for the excitation of a LRW at water-glass boundary) from the NPL probe axis and 9 mm radial distance from the source, as well as the shadowgraph images of the spherically expanding shock waves in water at 6 and 8 μs. The FOPH probe was further tilted by about 13° from the shock front normal direction to avoid picking up the flash from the spark charge produced at the tip of the NPL probe. (c) High-speed photoelastic/shadowgraph images (left) and COMSOL simulation results (right) showing the propagation of P i , P r and Schmidt head wave in water, T in the glass, and LRW at the water-glass boundary, respectively.

FIG. 2.
The COMSOL model-calculated pulse profiles for the first principal stress σ T (positive value indicates tension) detected at different radial distances (r) along five T wave rays (rays 1-5) of different refraction angles. Ray 1 overlaps with the L wave ray at 0° refraction angle,    A confocal microscopy examination of ring fractures induced at S d = 0.5mm and the inclination angles of ring fractures at various radius observed in experiments (blue dots) and predicted from COMSOL simulations (red dots and curve).