3.1. Governing equations

This study addresses the full FSI problem using a one-way coupling framework, which comprises the simulation of a droplet impacting a rigid wall and the resulting response of the solid material. The fluid impact is modelled as a compressible, multi-component, two-phase flow problem, involving complex physical phenomena such as phase transitions and strong shock waves. Here, the five-equation model (Allaire, Clerc & Kokh Reference Allaire, Clerc and Kokh2002; Saurel, Petitpas & Abgrall Reference Saurel, Petitpas and Abgrall2008; Wu et al. Reference Wu, Xiang and Wang2018) is employed to solve the fluid flow problem, and the governing equations are

(3.1a) \begin{align} \frac {\partial \alpha _l\rho _l}{\partial t}+\boldsymbol{\nabla }\boldsymbol{\cdot }(\alpha _l\rho _l \boldsymbol{u})=\dot {S}_{\rho ,l}, \\[-28pt] \nonumber \end{align}

(3.1b) \begin{align} \frac {\partial \alpha _v\rho _v}{\partial t}+\boldsymbol{\nabla }\boldsymbol{\cdot }(\alpha _v\rho _v \boldsymbol{u})=\dot {S}_{\rho ,v}, \\[-28pt] \nonumber \end{align}

(3.1c) \begin{align} \frac {\partial \alpha _g\rho _g}{\partial t}+\boldsymbol{\nabla }\boldsymbol{\cdot }(\alpha _g\rho _g \boldsymbol{u})=0, \\[-28pt] \nonumber \end{align}

(3.1d) \begin{align} \frac {\partial \rho \boldsymbol{u}}{\partial t}+\boldsymbol{\nabla }\boldsymbol{\cdot }(\rho \boldsymbol{uu}+p\boldsymbol{I})=0, \\[-28pt] \nonumber \end{align}

(3.1e) \begin{align} \frac {\partial E}{\partial t}+\boldsymbol{\nabla }\boldsymbol{\cdot }[(E+p)\boldsymbol{u}]=0, \\[-28pt] \nonumber \end{align}

(3.1f) \begin{align} \frac {\partial \alpha _l}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\alpha _l=\dot {S}_{\alpha ,l}, \\[-28pt] \nonumber \end{align}

(3.1g) \begin{align} \frac {\partial \alpha _v}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\alpha _v=\dot {S}_{\alpha ,v}, \\[-28pt] \nonumber \end{align}

where $\alpha _k$ and $\rho _k$ denote the volume fraction and the density of component $k$ , respectively. Three components, namely liquid water, water vapour and air, are considered and denoted by subscript $l$ , $v$ and $g$ , respectively. The volume fractions satisfy the relation $\alpha _g=1-\alpha _l-\alpha _v$ . The variables $\rho$ , $\boldsymbol{u}$ , $p$ , $E=\rho e+\rho \boldsymbol{u}^2/2$ and $e$ represent the density, velocity, pressure, total energy and specific internal energy of the mixture, respectively.

The thermodynamic state of the three components is described by the stiffened-gas equation of state (Saurel et al. Reference Saurel, Petitpas and Abgrall2008):

(3.2a) \begin{align} e_k(p,\rho _k)=\frac {p+\gamma _k p_{\infty ,k}}{\rho _k(\gamma _k-1)}+q_k, \\[-28pt] \nonumber \end{align}

(3.2b) \begin{align} \rho _k(p,T)=\frac {p+p_{\infty ,k}}{C_{v,k}T(\gamma _k-1)}, \\[-28pt] \nonumber \end{align}

(3.2c) \begin{align} h_k(T)=\gamma _k C_{v,k} T+q_k, \\[-28pt] \nonumber \end{align}

(3.2d) \begin{align} c_k(p,\rho _k)=\sqrt {\frac {\gamma _k(p+p_{\infty ,k})}{\rho _k}}, \\[-28pt] \nonumber \end{align}

(3.2e) \begin{align} g_k(p,T)= \big(\gamma _k C_{v,k}-q'_k \big)T-C_{v,k}T {\rm {log}} \frac {T^{\gamma _k}}{(p+p_{\infty ,k})^{\gamma _k-1}} +q_k, \\[-28pt] \nonumber \end{align}

(3.2f) \begin{align} \mu _v=g_v(p,T)+(\gamma _v-1)C_{v,v}T\mathrm{log} \frac {\alpha _v}{1-\alpha _l},\mu _l=g_l. \\[0pt] \nonumber \end{align}

Here $e_k$ , $h_k$ , $c_k$ , $g_k$ and $\mu _k$ denote the specific internal energy, the specific enthalpy, the speed of sound, the specific Gibbs free energy and the chemical potential of component $k$ , respectively. The equation parameters, including $\gamma _k$ , $p_{\infty ,k}$ , $C_{v,k}$ , $q_k$ and $q'_k$ , can be determined referring to previous studies (Han, Hantke & Müller Reference Han, Hantke and Müller2017; Wu et al. Reference Wu, Xiang and Wang2018). Then, the mixture quantities can be calculated as follows (Saurel et al. Reference Saurel, Petitpas and Abgrall2008):

(3.3a) \begin{align} \rho =\alpha _l\rho _l+\alpha _v\rho _v+\alpha _g\rho _g, \\[-28pt] \nonumber \end{align}

(3.3b) \begin{align} \rho e=\alpha _l \rho _l e_l+\alpha _v \rho _v e_v+\alpha _g \rho _g e_g, \\[-28pt] \nonumber \end{align}

(3.3c) \begin{align} \rho c^2=\alpha _l\rho _lc_l^2+\alpha _v\rho _vc_v^2+\alpha _g\rho _gc_g^2. \\[-12pt] \nonumber \end{align}

The source term caused by the phase transition phenomenon in the governing equations can be calculated by

(3.4a) \begin{align} \dot {S}_{\rho ,l}=-\dot {S}_{\rho ,v}=\nu (\mu _v-\mu _l), \\[-28pt] \nonumber \end{align}

(3.4b) \begin{align} \dot {S}_{\alpha ,l}=\frac {\nu }{\varrho _l}(\mu _v-\mu _l),\dot {S}_{\alpha ,v}=\frac {\nu }{\varrho _v}(\mu _l-\mu _v). \\[0pt] \nonumber \end{align}

Here, the relaxation parameter $\nu$ can be taken as positive infinite once the phase transition is triggered and as zero otherwise. A detailed expression of the parameter $\varrho$ is given in Zein, Hantke & Warnecke (Reference Zein, Hantke and Warnecke2010).

The response of the solid material is modelled using the classical elastic mechanics. Under this framework, the material is assumed as continuous, homogeneous and isotropic, with a small deformation during the entire simulation (Gould Reference Gould2013). These assumptions are well justified by the material properties of the W4600 coating in table 1 and the descriptions from the supplier (3M 2008). The governing equations for elastic mechanics consist of the kinematic equation for continuum (Slaughter Reference Slaughter2002),

(3.5) \begin{align} \frac {\partial \sigma _{\textit{ij}}}{\partial x_j} + f_i=\rho _{\kern-1pt s} \ddot {u}_i, \end{align}

where $\sigma _{\textit{ij}}$ , $f_i$ and $u_i$ denote the stresses, body forces and displacements, respectively, of component $i$ in local coordinates ( $x_i$ , $x_j$ ). Here $\ddot {u}_i$ denotes the second-order time derivative $\partial ^2 u_i/\partial t^2$ of the displacement variables. The stresses in (3.5) relate to the strain tensor, $\epsilon _{\textit{ij}}$ , according to the physical (constitutive) equation

(3.6) \begin{align} \sigma _{\textit{ij}}=C_{\textit{ijkl}}\epsilon _{kl}, \end{align}

where $C_{\textit{ijkl}}$ denotes the fourth-order stiffness tensor (Belen’kiĭ et al. Reference Belen’kiĭ, Salaev and Suleĭmanov1988). For isotropic solid materials, this simplifies to the Lamé constants of the volumetric deformation parameter $\lambda$ and the shear deformation parameter $\mu$ (Salencon Reference Salencon2001), i.e.

(3.7) \begin{align} C_{\textit{ijkl}}\epsilon _{kl}=\lambda \epsilon _{kk}\delta _{\textit{ij}} + 2\mu \epsilon _{\textit{ij}} \end{align}

for volumetric strain $\epsilon _{kk}$ and the Kronecker delta $\delta _{\textit{ij}}$ . The parameters of $\lambda$ and $\mu$ are functions of the material’s Young’s modulus $E_s$ and Poisson’s ratio $\nu _s$ , and can be calculated referring to the elastic theories (Timoshenko & Goodier Reference Timoshenko and Goodier1951):

(3.8a) \begin{align} \lambda =\frac {2\mu \nu _s}{1-2\nu _s}, \\[-28pt] \nonumber \end{align}

(3.8b) \begin{align} \mu =\frac {E_s}{2(1+\nu _s)}. \\[0pt] \nonumber \end{align}

The strain tensor in (3.7) relates to the displacements, $u_i$ , by the geometric equation

(3.9) \begin{align} \epsilon _{\textit{ij}}=\frac {1}{2} \left ( \frac {\partial u_i}{\partial x_j} + \frac {\partial u_j}{\partial x_i} \right )\!. \end{align}

In this way, (3.5), (3.7) and (3.9) solve the displacements of material elements subject to the applied external forces and boundary conditions. The latter will be discussed in § 3.2 specifically for the simulated physical problem.

3.2. Numerical treatments

The numerical simulation of fluid flow is conducted by the in-house software SCP-tran^©, which has been validated in detail and used to solve several problems of compressible two-phase flow (Xiang & Wang Reference Xiang and Wang2017; Wu et al. Reference Wu, Xiang and Wang2018; Xu et al. Reference Xu, Fan, Wu, Wen and Wang2023). The software employs a finite-volume method to discretise governing equations, which involves the fifth-order weighted essentially non-oscillatory scheme for spatial reconstruction (Titarev & Toro Reference Titarev and Toro2004; Coralic & Colonius Reference Coralic and Colonius2014). The Harten-Lax-van Leer contact approximate Riemann solver is employed to calculate the numerical fluxes on the cell edges (Toro Reference Toro2009). As for temporal discretisation, the software utilises the third-order total variation diminishing Runge–Kutta scheme to update the flow field over time (Gottlieb & Shu Reference Gottlieb and Shu1998). Here, uniform grids as well as a constant Courant–Friedrich–Lewis (CFL) number, 0.4, are used in all numerical cases. Considering the results of numerical grid independence tests provided in Appendix A.1, the grid resolution level of 2000 cells per droplet diameter is employed for the fluid simulation.

The calculation of the source terms associated with phase transition is treated as an equivalent relaxation process of chemical potential, which can be decoupled from the hyperbolic operator (Saurel et al. Reference Saurel, Petitpas and Abgrall2008). Here, two processes of phase transition, including the cavitation and condensation, are considered. The condensation process is triggered when the condition $\mu _v\gt \mu _l$ is satisfied, while cavitation is triggered when the condition $p\lt p_{\textit{threshold}}$ is satisfied. Once a phase transition is triggered, the subsequent relaxation calculation will be activated, and the local state variables, including $\alpha _k$ , $\alpha _k\rho _k$ , $p$ and $T$ , will be readjusted under the conservation of mass, momentum and energy, as well as thermal equilibrium in order to achieve $\mu _l=\mu _v$ . For the details of relaxation calculation, please refer to Han et al. (Reference Han, Hantke and Müller2017) and Wu et al. (Reference Wu, Xiang and Wang2018). It should be noted that here we consider two types of cavitation nucleation, each with its own cavitation threshold $p_{\textit{threshold}}$ (Brennen Reference Brennen2013). For the regions away from the solid surface, the bulk homogeneous cavitation is under consideration, where $p_{\textit{threshold}}=p_{\textit{hom}}$ , while for the regions around the solid surface, the near-wall heterogeneous cavitation occurs, corresponding to $p_{\textit{threshold}}=p_{\textit{het}}$ that satisfies (Wu et al. Reference Wu, Liu and Wang2021)

(3.10) \begin{align} (p_s-p_{\textit{het}})^2=(p_s-p_{\textit{hom}})^2\frac {(2-\mathrm{cos}\xi )(1+\mathrm{cos}\xi )^2}{4}. \end{align}

Here $p_s$ denotes the saturated vapour pressure of water and $\xi$ is the static contact angle between liquid water and the solid surface. In practice, the cavitation threshold is related to the impurity and the physical state of water and approximately ranges from $-30$ to $-0.1$ $\mathrm{MPa}$ (Briggs Reference Briggs1950; Caupin & Herbert Reference Caupin and Herbert2006; Brennen Reference Brennen2013; Hong & Son Reference Hong and Son2022). Considering the properties of the rainwater (and seawater for offshore turbines) and the coating surface in the aforementioned physical model, we fix $p_{\textit{hom}}=-10$ $\mathrm{MPa}$ and $\xi =150^\circ$ during the entire simulation. More discussion on the selection of these cavitation-related parameters is provided in Appendix B.

Upon droplet impact, the solid structure is modelled using FE analysis within the commercial ABAQUS FE environment. Both two-dimensional and axisymmetric analyses of the solid structure are performed consistently with the previously described droplet impact simulations. Figure 3 depicts the simulation domain for the FE analyses. In axisymmetric simulations, the blade coating material is modelled by four-node axisymmetric quadrilateral elements with reduced integration (CAX4R; element codes in ABAQUS). Additionally, a layer of one-way infinite elements (CINAX4) is placed at the side and bottom boundaries to prevent wave reflections within the solid structure due to the finite simulation domain (Doagou-Rad & Mishnaevsky Reference Doagou-Rad and Mishnaevsky2020). For the two-dimensional simulations, the aforementioned FE elements are substituted with corresponding plane stress elements (CPS4R and CINPS4, respectively), with the out-of-plane thickness assumed sufficiently small to ensure a two-dimensional representation.

Figure 3.

Schematic illustration of the axisymmetric (or two-dimensional) simulation domain of the solid structure in the FE analysis. The image reflection illustrates the axisymmetric (or two-dimensional symmetric) boundary conditions along the $z$ axis of the simulation domain. The dimensions in the figure are not to scale.

The axisymmetric (or two-dimensional) FE model is used to simulate the behaviour of a W4600 coating with dimensions of length $L_s=10$ mm and height $H_s=1$ mm, over a total simulation time of approximately $5$ $\unicode{x03BC}\mathrm{s}$ . The numerical time step is adaptively determined by the ABAQUS explicit solver, based on local material properties, mesh geometry and ongoing deformations, thereby ensuring CFL stability of the FE analysis throughout the computation (Dassault Systémes Simulia Corporation 2023). The coating in the current FE analysis is uniformly discretised with an element size of $5$ $\unicode{x03BC}\mathrm{ m}$ in all computations, which adheres to the mesh size guidance of less than $2$ % $R_0$ for numerical impact erosion simulations (Zhang et al. Reference Zhang, Zhang, Lv, Li and Guo2019; Doagou-Rad & Mishnaevsky Reference Doagou-Rad and Mishnaevsky2020). Axisymmetric (or two-dimensional symmetric) boundary conditions are applied along the $z$ axis, as shown in figure 3. The material properties of the simulated W4600 coating are detailed in table 1, with the Young’s modulus calculated using (2.1). Given the high material strength and stiffness, we assume a linear elastic model for the coating material. This assumption is based on the fact that a single raindrop, impacting at a velocity of O( $10$ ) $\mathrm{m\,s^{-1}}$ , is unlikely to induce plastic deformation, and is expected to cause only elastic deformation in this high-stiffness material. The numerical grid independence and validation tests for the material part of the current FSI simulation are referred to in the treatments in Appendix A.2 and the authors’ previous work (Hao et al. Reference Hao, Charalambides, Hardalupas, Sergis and Taylor2026a , Reference Hao, Sergis, Taylor, Hardalupas and Charalambidesd ), respectively.

The one-way coupling framework is adaptively implemented with respect to mesh grids on the contact surface between the droplet impact simulation and FE material analysis, applying a time- and space-dependent loading condition to the top surface of the solid structure. This models the impact load of a liquid droplet and subsequent cavitation within the FE analysis. The user-defined loading condition is incorporated using the VDLOAD subroutine in ABAQUS for explicit dynamic analysis (Hao et al. Reference Hao, Charalambides, Hardalupas, Sergis and Taylor2026a ). In the one-way coupling framework, the bottom surface, i.e. solid surface, of the fluid flow simulation remains rigid, whereas the solid in the FE analysis deforms and recovers in response to the impact pressure. The contradiction is addressed by findings in Hao et al. (Reference Hao, Sergis, Taylor, Hardalupas and Charalambides2026d ), which demonstrates that surface motion has a negligible effect on fluid behaviour in one-way coupled FSI systems when the material stiffness exceeds approximately O( $10^{5}$ ) $\mathrm{Pa}$ . Since typical protective coatings of wind turbine blades – including the W4600 coating material considered in this study – exhibit stiffness values in the range of O( $10^9{-}10^{12}$ ) $\mathrm{Pa}$ , the rigid-surface assumption in the fluid flow model is valid. Thus, the use of a one-way coupling approach is justified for the present impact simulations.

4.1. General patterns

For details of the Bernoulli field evolutions during droplet impact, we refer the reader to our earlier studies (Hao et al. Reference Hao, Charalambides, Hardalupas, Sergis and Taylor2026a , Reference Hao, Charalambides, Hardalupas, Sergis and Taylorb ). Here, we move beyond that framework by analysing the complex wave dynamics within the droplet following impact, with a particular focus on the cavitation phenomenon arising from local pressure fluctuations. We aim to elucidate the potential mechanisms by which both bulk homogeneous cavitation and near-wall heterogeneous cavitation are initiated and contribute to surface damage. As a representative case, we consider a two-dimensional droplet with an impact velocity of $U_0=50$ $\mathrm{ms^{-1}}$ . Figure 4 presents the pressure contours at multiple time instants, illustrating the evolution process of internal wave structures and the mechanisms responsible for triggering cavitation.

Figure 4.

Evolution of wave structures and the corresponding cavitation phenomenon inside the droplet ( $U_0=50$ $\mathrm{m\,s^{-1}}$ ), shown by pressure contours. The cavitation zones are presented by the volume fraction isolines of water vapour. Results are shown for the following times: ( $a$ ) 0.450 $\unicode{x03BC}\mathrm{s}$ , ( $b$ ) 0.710 $\unicode{x03BC}\mathrm{s}$ , ( $c$ ) 0.874 $\unicode{x03BC}\mathrm{s}$ , ( $d$ ) 1.038 $\unicode{x03BC}\mathrm{s}$ , ( $e$ ) 1.337 $\unicode{x03BC}\mathrm{s}$ , ( $f$ ) 1.504 $\unicode{x03BC}\mathrm{s}$ , ( $g$ ) 2.283 $\unicode{x03BC}\mathrm{s}$ , ( $h$ ) 2.729 $\unicode{x03BC}\mathrm{s}$ , ( $i$ ) 3.520 $\unicode{x03BC}\mathrm{s}$ .

Figure 4. Long description

Panel A: A heat map showing pressure contours with a water-hammer shock visible. The x-axis ranges from -1.4 to 1.4 and the z-axis ranges from 0 to 2. The color scale ranges from -20.0 to 40.0 MPa, with red indicating higher pressure. Panel B: A heat map showing pressure contours with re-reflected rarefaction waves visible. The x-axis ranges from -1.4 to 1.4 and the z-axis ranges from 0 to 2. The color scale ranges from -20.0 to 40.0 MPa, with red indicating higher pressure. Panel C: A heat map showing pressure contours with a cavitation zone visible. The x-axis ranges from -1.4 to 1.4 and the z-axis ranges from 0 to 2. The color scale ranges from -20.0 to 40.0 MPa, with red indicating higher pressure. Panel D: A heat map showing pressure contours with a collapsing shock and an edge-re-reflected rarefaction wave visible. The x-axis ranges from -1.4 to 1.4 and the z-axis ranges from 0 to 2. The color scale ranges from -6.0 to 12.0 MPa, with red indicating higher pressure. Panel E: A heat map showing pressure contours with near-wall cavitation and a re-edge-re-reflected rarefaction wave visible. The x-axis ranges from -1.4 to 1.4 and the z-axis ranges from 0 to 2. The color scale ranges from -6.0 to 12.0 MPa, with red indicating higher pressure. Panel F: A heat map showing pressure contours with a collapsing shock visible. The x-axis ranges from -1.4 to 1.4 and the z-axis ranges from 0 to 2. The color scale ranges from -6.0 to 12.0 MPa, with red indicating higher pressure. Panel G: A heat map showing pressure contours with near-wall cavitation visible. The x-axis ranges from -1.4 to 1.4 and the z-axis ranges from 0 to 2. The color scale ranges from -3.0 to 3.0 MPa, with red indicating higher pressure. Panel H: A heat map showing pressure contours with near-wall cavitation visible. The x-axis ranges from -1.4 to 1.4 and the z-axis ranges from 0 to 2. The color scale ranges from -3.0 to 3.0 MPa, with red indicating higher pressure. Panel I: A heat map showing pressure contours with a collapsing shock visible. The x-axis ranges from -1.4 to 1.4 and the z-axis ranges from 0 to 2. The color scale ranges from -3.0 to 3.0 MPa, with red indicating higher pressure.

The impingement on the solid wall generates a confined WH shock wave at the bottom of the droplet. As the WH shock propagates upward towards the droplet apex, it is continuously reflected from the internal free surface (figure 4 a). Due to the discrepancy in acoustic impedance between the liquid and air phases, each reflection produces a rarefaction wave. These reflected rarefaction waves primarily consist of two branches, denoted as $\mathrm{Re}\text{-}\mathrm{RW}_{I}$ and $\mathrm{Re}\text{-}\mathrm{RW}_{\textit{II}}$ , as shown in figure 4(b). The two symmetric waves of branch $\mathrm{Re}\text{-}\mathrm{RW}_{\textit{II}}$ converge along the axis of symmetry at approximately $t\sim 0.874$ $\unicode{x03BC}\mathrm{s}$ (figure 4 c), resulting in a substantial local pressure drop and the onset of homogeneous cavitation within the bulk liquid. This phenomenon has been investigated experimentally by Field et al. (Reference Field, Dear and Ogren1989) and analysed numerically by Wu et al. (Reference Wu, Xiang and Wang2018). The subsequent rapid collapse of the cavitation bubbles produces a strong elliptic collapsing shock wave, as noted in figure 4(d). Meanwhile, following their convergence on the central axis, the two symmetric $\mathrm{Re}\text{-}\mathrm{RW}_{\textit{II}}$ waves begin to expand again. As a result, this expansion gives rise to an expanding reflected rarefaction wave (ERRW), which propagates towards the solid wall ahead of the collapsing shock wave – also visible in figure 4(d).

In the fluid simulation, the solid wall is modelled as perfectly rigid, possessing an effectively infinite acoustic impedance ( $Z = \rho c$ ), which greatly exceeds that of liquid water. Therefore, the reflected wave of an ERRW from the solid wall is still a rarefaction wave, denoted as $\mathrm{Re}\text{-}\mathrm{ERRW}$ in figure 4(e). This reflected wave $\mathrm{Re}\text{-}\mathrm{ ERRW}$ creates a low-pressure region near the solid wall, reducing the local pressure below the heterogeneous cavitation threshold $p_{\textit{het}}$ , thereby initiating near-wall heterogeneous cavitation at approximately $t\sim 1.337$ $\unicode{x03BC}\mathrm{s}$ . Following its inception, the near-wall cavitation bubbles are subjected to a pressure recovery primarily driven by the collapsing shock wave emitted from the previously formed bulk cavitation zone. This pressure rise causes the cavitation bubbles to shrink and collapse, producing a series of secondary collapsing shock waves, as demonstrated in figure 4(f).

As previously reported by Wu et al. (Reference Wu, Liu and Wang2021), repeated reflections of collapsing shock waves can lead to recurrent near-wall cavitation – a phenomenon also observed in the present numerical simulation. After approximately $t\sim 1.500$ $\unicode{x03BC}\mathrm{s}$ , the wave structures become increasingly complex, consisting of the $\mathrm{Re}\text{-}\mathrm{ERRW}$ wave, the first and second collapsing waves, and their respective reflections (figure 4 g). Upon interacting with the droplet-surface interface, these shock waves continue to generate low-pressure rarefaction waves that propagate towards the solid wall, ultimately inducing a second near-wall cavitation event at $t\sim 2.729$ $\unicode{x03BC}\mathrm{s}$ (figure 4 h). The accumulated near-wall cavitation bubbles eventually undergo a sustained collapse process beginning around $t\sim 3.200$ $\unicode{x03BC}\mathrm{s}$ , producing multiple collapsing shock waves (figure 4 i). Due to the propagation of these waves, the solid wall continues to be subjected to significant pressure fluctuations after $t\sim 3.200$ $\unicode{x03BC}\mathrm{s}$ , which may contribute to material damage through the mechanism of heterogeneous cavitation erosion.

It is worth noting that, although a perfectly rigid wall is used in the fluid simulation, the relative acoustic impedance between the solid and the fluid is preserved even when the wall is replaced by the compliant coating material considered in the physical model (§ 2). This ensures the relevance and applicability of the simulation results above. Furthermore, as the theoretical analysis of the wave evolution has been thoroughly addressed in previous studies (Wu et al. Reference Wu, Xiang and Wang2018; Xu et al. Reference Xu, Fan, Wu, Wen and Wang2023), the present study will not discuss the details for brevity.

4.2. The effects of cavitation and the impact velocity on wall pressure

In this section we analyse the pressure distribution exerted on the solid wall throughout the droplet impacting and spreading processes. We aim to elucidate the additional erosion potential associated with cavitation phenomena, as well as the influence of impact velocity. To this end, three representative cases are considered: ( $a$ ) $U_0=50$ $\mathrm{m\,s^{-1}}$ with cavitation, ( $b$ ) $U_0=90$ $\mathrm{m\,s^{-1}}$ with cavitation, and ( $c$ ) $U_0=50$ $\mathrm{m\,s^{-1}}$ without cavitation. Figure 5 presents pressure contours and cavitation patterns at three selected time instants for each case. The first instant corresponds to the convergence of the reflected waves $\mathrm{Re}\text{-}\mathrm{RW}_{\textit{II}}$ and the onset of induced homogeneous cavitation. The second and third instants correspond to the first and second collapse events of near-wall cavitation bubbles, respectively. Figure 6 compares the space–time diagrams of the wall pressure distributions for all three cases.

Figure 5.

The effects of cavitation and the impact velocity on the wave evolution inside the droplet in two-dimensional simulations, presented by pressure contours at typical instants. Cavitation zones are denoted in the same way as figure 4. Conditions are ( $a$ ) $U_0=50$ $\mathrm{m\,s^{-1}}$ with cavitation, ( $b$ ) $U_0=90$ $\mathrm{m\,s^{-1}}$ with cavitation, and ( $c$ ) $U_0=50$ $\mathrm{m\,s^{-1}}$ without cavitation. Instants shown are ( $a\,\textrm{i}$ ) 0.874 $\unicode{x03BC}\mathrm{s}$ , ( $b\,\textrm{i}$ ) 0.871 $\unicode{x03BC}\mathrm{s}$ , ( $c\,\textrm{i}$ ) 0.874 $\unicode{x03BC}\mathrm{s}$ , ( $a\,\textrm{ii}$ ) 1.504 $\unicode{x03BC}\mathrm{s}$ , ( $b\,\textrm{ii}$ ) 1.610 $\unicode{x03BC}\mathrm{s}$ , ( $c\,\textrm{ii}$ ) 1.504 $\unicode{x03BC}\mathrm{s}$ , ( $a\,\textrm{iii}$ ) 3.486 $\unicode{x03BC}\mathrm{s}$ , ( $b\,\textrm{iii}$ ) 3.203 $\unicode{x03BC}\mathrm{s}$ and ( $c\,\textrm{iii}$ ) 3.496 $\unicode{x03BC}\mathrm{s}$ .

Figure 6.

Comparison of the space–time diagrams of pressure distributions exerted on the solid wall under three conditions: ( $a$ ) $U_0=50$ $\mathrm{m\,s^{-1}}$ with cavitation, ( $b$ ) $U_0=90$ $\mathrm{m\,s^{-1}}$ with cavitation, and ( $c$ ) $U_0=50$ $\mathrm{m\,s^{-1}}$ without cavitation. The results are from two-dimensional simulations. Regions $(A)$ and $(B)$ are enlarged zones noted by dashed lines in ( $a$ ) and ( $b$ ), respectively.

The comparison between the first two cases highlights the influence of impact velocity on wave evolution and cavitation processes within the droplet. An increase in impact velocity enhances both the strength of the WH shock wave and the velocity of lateral jetting (i.e. wider pressure distributions in $x$ ), as illustrated in figure 6. The transient pressure generated by the WH impact is positively correlated with the impact velocity $U_0$ and can be approximated by the relation $p_{\textit{WH}}=\rho _lU_0(c_l+\chi U_0)$ , following Heymann (Reference Heymann1969), where $\rho _l$ and $c_l$ are the density and sound speed of the undisturbed liquid, respectively, and $\chi =2.0$ is usually taken for liquid water. Consequently, higher impact velocities lead to more pronounced pressure reductions caused by the reflected rarefaction waves (e.g. $\mathrm{Re}\text{-}\mathrm{RW}_{\textit{II}}$ , $\mathrm {Re}\text{-}\mathrm {ERRW}$ and reflected waves of collapsing shocks), thereby generating a greater total volume of both homogeneous and heterogeneous cavitation bubbles. This, in turn, results in more intense bubble collapse events and stronger secondary shock waves, as illustrated in figures 5(a ii,b ii) and 5(a iii,b iii). In the wall pressure distributions shown in figure 6, higher impact velocities not only lead to the more intense and sustained WH pressure before $t\sim 1.250$ $\unicode{x03BC}\mathrm{s}$ , but also more pronounced pressure peaks associated with the first ( $t\sim 1.500$ $\unicode{x03BC}\mathrm{s}$ ) and the second ( $t \gt 3.000$ $\unicode{x03BC}\mathrm{s}$ ) near-wall bubble collapse events. For instance, the maximum wall pressure resulting from the first collapse (near $t\sim 1.500$ $\unicode{x03BC}\mathrm{s}$ ) in case $(b)$ ( $U_0=90$ $\mathrm{m\,s^{-1}}$ ) exceeds $40$ $\mathrm{MPa}$ , whereas the corresponding value in case $(a)$ ( $U_0=50$ $\mathrm{m\,s^{-1}}$ ) is approximately $10$ $\mathrm{MPa}$ . During the second collapse event, although both cases exhibit comparable peak pressures above $50$ $\mathrm{MPa}$ , the spatial pressure distribution differs significantly: in case $(a)$ , the high-pressure regions are more localised and characterised by discrete maxima (figure 6 A), while in case $(b)$ , the pressure is distributed more uniformly across a broader region (figure 6 B). Nonetheless, increasing the impact velocity within the considered range does not alter the fundamental qualitative mechanisms of wave evolution and cavitation formation shown in figure 4.

The comparison between cases with and without cavitation reveals the additional effects of cavitation phenomena on wave propagation and wall pressure. To a certain extent, the collapse of cavitation bubbles modifies the pressure fluctuations that would otherwise result solely from the WH impact. For instance, figures 5(a ii) and 5(c ii) show that the low-pressure region trailing the upward-propagating $\mathrm{Re}\text{-}\mathrm{ERRW}$ wave is disrupted and partially suppressed by disturbances generated from bubble collapse. However, at the relatively low impact velocity of $U_0=50$ $\mathrm{m\,s^{-1}}$ , the first collapse of near-wall cavitation bubbles (around $t\sim 1.500$ $\unicode{x03BC}\mathrm{s}$ ) has a negligible influence on the wall pressure, as seen in figures 6(a) and 6( $c$ ). In contrast, the second collapse event is substantially more intense, producing multiple collapsing shock waves (figure 5 $a\,\textrm{iii}$ ) and resulting in significantly elevated wall pressures (figure 6 A). In this case, local wall pressures can exceed $50$ $\mathrm{MPa}$ at specific locations, whereas in the non-cavitating scenario (case c), the wall pressure remains below $4$ $\mathrm{MPa}$ across the entire wall surface after $t = 3.000$ $\unicode{x03BC}\mathrm{s}$ . These observations clearly demonstrate that cavitation plays a crucial role in amplifying pressure loads on the solid wall. Since a further increase in impact velocity does not fundamentally alter the qualitative cavitation mechanisms, the subsequent analysis of material response and potential erosion is conducted by comparing the cavitating and non-cavitating cases at a fixed impact velocity.

5.1. Rayleigh surface waves

As a liquid droplet impacts the elastic solid, the droplet’s kinetic energy is transferred to the solid, generating waves that propagate through the material from the contact surface. According to stress theory (Lamb Reference Lamb1904; Woods Reference Woods1968), the elastic waves, which consist of the two sources, i.e. longitudinal and transverse components, are initially confined to the loaded area, which is the contact area beneath the droplet, for a transient time $t^*$ , with the radius of disturbance matching the extent of the loaded area. As the growth in the loaded area slows below the sonic speed in the solid, elastic waves detach ahead of the loaded area and propagate at sonic speeds. Among these waves, the longitudinal wave is the first to break away as a head wave (Blowers Reference Blowers1969) ( $X_L$ ), followed by the transverse wave. As the transverse wave detaches from the edge of the loaded area, a Rayleigh wave ( $X_{\kern-0.5pt R}$ ) forms in the transition region between the edge of the loaded area ( $X_{\kern-0.5pt P}$ ) and the transverse wavefront ( $X_T$ ), where surface particles undergo retrograde motion (figure 7 a). This retrograde motion arises from a mismatch in displacements between the inner (loaded) region (I), which remains compressed (Field et al. Reference Field, Dear and Ogren1989), and the transverse wave front. Upon formation, the normal surface displacements at the Rayleigh wave front rises steeply (Miller & Pursey Reference Miller and Pursey1955; Kim et al. Reference Kim, Kim, Kang and Kim2000), causing significant perturbations to the solid surface (Field Reference Field1966). Figure 7(a) illustrates the fully developed profile of the Rayleigh surface wave during liquid droplet impact at $U_0=50$ $\mathrm{m\,s^{-1}}$ , including its formation and propagation stages. While the inner region ( $I$ ) remains compressed (Field et al. Reference Field, Dear and Ogren1989), transverse and longitudinal waves divide the solid surface into the transition region ( $ \textit{II} $ ), where the Rayleigh wave is embedded, and the outer region ( $ \textit{III} $ ) characterised by surface displacements. A detailed (vertical) displacement distribution is shown in figure 7(b), with the maximum surface depression at the centre due to the droplet impact happening at $t_{max}\approx 1.1$ $\unicode{x03BC}\mathrm{s}$ , marked in red.

As shown in figure 7(b), the FE numerical simulation captures the position of the edge of the loaded area ( $X_{\kern-0.5pt P}$ ), corresponding to a change in the surface displacement curvature, which marks the boundary between the inner and transition regions. The position of the Rayleigh wave ( $X_{\kern-0.5pt R}$ ) is also identified by the peak in (positive) surface displacement, labelled in green. Consistent with stress theory (Blowers Reference Blowers1969), the decay between $X_{\kern-0.5pt R}$ and $X_T$ in the outbound of the transition region is steeper than that between $X_{\kern-0.5pt R}$ and $X_{\kern-0.5pt P}$ in the inbound. Moreover, after the formation of the Rayleigh wave, the characteristic slow decay behaviour (O’Connell-Rodwell et al. Reference O’Connell-Rodwell, Arnason and Hart2000) of the wave magnitudes is clearly observed in both time and distance until the end of the simulation at $5.2$ $\unicode{x03BC}\mathrm{s}$ . Thus, the shape and (vertical) displacements of the solid surface (without cavitation) are well characterised both spatially and temporally in the current material analysis. However, we note that the stress theories in literature (Miller & Pursey Reference Miller and Pursey1955; Woods Reference Woods1968; Blowers Reference Blowers1969) typically assume a constant uniform pressure over the solid, whereas in response to a dynamic hydrodynamic loading on the surface, quantitative elastic wave characterisations are more complex.

In contrast, figure 7(c) shows the vertical surface displacements of the solid in response to the droplet impact at the same impact speed but with cavitation. Before cavitation occurs, the solid undergoes the same Rayleigh wave formation process as in the non-cavitation case, with the surface loading primarily characterised by the initiation of WH pressure (Cook Reference Cook1928; Lesser Reference Lesser1995) at the onset of droplet impact. Cavitation, however, occurs later, triggered when the confined WH shock waves propagate to, and are reflected at, the droplet’s free interface (Heymann Reference Heymann1969; Wu et al. Reference Wu, Xiang and Wang2018). During cavitation, we observe a distinct surface depression near the centre of the loaded region ( $I$ ), which leads to a secondary surface depression (denoted by a blue arrow in figure 7 c) shortly after the surface partially recovers from the WH impact. Interestingly, the maximum deflection due to cavitation reaches $9.4$ $\unicode{x03BC}\mathrm{ m}$ , which is comparable to the WH compression of $10.0$ $\unicode{x03BC}\mathrm{ m}$ at the surface centre, but occurs off-centre at $x\approx 0.8R_0$ . Although cavitation has negligible effect extended into the transition ( $ \textit{II} $ ), and hence, the Rayleigh wave, and the outer ( $ \textit{III} $ ) regions, it causes notable additional surface distortions in the loaded (inner) region ( $I$ ). These distortions, particularly the wavy surface curvatures, are comparable to that of the Rayleigh waves on the surfaces and result in peak stresses (Verma et al. Reference Verma, Castro, Jiang and Teuwen2020) (as discussed in subsequent sections) in the solid. According to Bowden & Field (Reference Bowden and Field1964) and Blowers (Reference Blowers1969), such displacement waves could initiate surface pits and, through further interaction, generate ring-shaped surface cracks (Bourne, Obara & Field Reference Bourne, Obara and Field1997), which are a common type of damage in elastic materials.

5.2. Propagation of stress waves within the solid: fatigue insights

Following the propagation of displacement waves both on the surface and through the material, stresses are generated within the solid. Fatigue, driven by stress cycles, is widely recognised as the primary erosion mechanism for impact erosion at subsonic speeds (Ibrahim & Medraj Reference Ibrahim and Medraj2020; Springer Reference Springer1976; Doagou-Rad & Mishnaevsky Reference Doagou-Rad and Mishnaevsky2020), as well as for cavitation erosion (Hobbs Reference Hobbs1966; Preece & Brunton Reference Preece and Brunton1980; Varavka & Kudryakov Reference Varavka and Kudryakov2015), since the water jet induced by bubble collapse generates shock waves similar to those in droplet impacts (Blake & Gibson Reference Blake and Gibson1987; Speirs et al. Reference Speirs, Langley, Pan, Truscott and Thoroddsen2021). In this section we analyse the maximum principal stress in the material, which characterises the fatigue stress state – tensile stresses only – under complex loading conditions. Figure 8(left) illustrates the evolution of maximum principal stress during liquid droplet impact at $U_0=50$ $\mathrm{m\,s^{-1}}$ without cavitation. At the initial impact stage (figure 8 a), the maximum impact stress $S_{\textit{WH}}$ of approximately $25$ $\mathrm{MPa}$ is concentrated at (or near) the surface centre due to the transient WH pressure peak (Hancox & Brunton Reference Hancox and Brunton1966; Field Reference Field1966). Following its generation, stress waves spread through the solid with decreasing magnitude. Similar to the displacement waves, the stress wavefront propagates faster than the growth in the loaded region (denoted by red arrows in figure 8), with its magnitude decaying most slowly near the surface. As a result, we observe a horizontally propagating stress wavefront, labelled $\hat {S}_{\textit{WH}}$ , even after the primary stresses beneath the loaded region have mostly decayed at $t\approx 2.78$ $\unicode{x03BC}\mathrm{s}$ (figure 8 k). In this paper we use the hat symbol to indicate the stress wavefronts. At $t\approx 5.00$ $\unicode{x03BC}\mathrm{s}$ (figure 8 s), the horizontally propagating stress wavefront $\hat {S}_{\textit{WH}}$ still exhibits a notable magnitude of approximately $15$ $\mathrm{MPa}$ . An interesting pattern of the maximum principal stress is that the tensile waves are outboard of the compressive waves. The latter is within the loaded region ( $I$ ), as demonstrated in figure 8(c). This pattern corresponds to opposing motions on either side of the retrograde Rayleigh wave behaviour near the surface (see figure 7 a).

Figure 8.

Comparison of the maximum principal stresses in the two-dimensional solid material between droplet impacts without (a,c,g,i,k,m,o,q,s) and with (b,d,f,h,j,l,n,p,r,t) cavitation at $U_0=50$ $\mathrm{m\,s^{-1}}$ . Each row corresponds to comparable but not identical time instants from the two cases, with approximate simulation times for each row as follows: $t\approx$ $0.59$ $\unicode{x03BC}\mathrm{s}$ (a,b), $0.80$ $\unicode{x03BC}\mathrm{s}$ (c,d), $1.28$ $\unicode{x03BC}\mathrm{s}$ (e,f), $1.56$ $\unicode{x03BC}\mathrm{s}$ (g,h), $1.94$ $\unicode{x03BC}\mathrm{s}$ (i,j), $2.78$ $\unicode{x03BC}\mathrm{s}$ (k,l), $3.78$ $\unicode{x03BC}\mathrm{s}$ (m,n), $4.09$ $\unicode{x03BC}\mathrm{s}$ (o,p), $4.40$ $\unicode{x03BC}\mathrm{s}$ (q,r) and $5.00$ $\unicode{x03BC}\mathrm{s}$ (s,t). The loaded region of each subfigure is denoted by a red line with the edge $X_{\kern-0.5pt P}$ labelled with an arrow.

Figure 8(right) compares the evolution of the maximum principal stress during liquid droplet impact at the same speed but with cavitation. Similar to the non-cavitation case, the solid material experiences the same stress waves before cavitation occurs (figure 8 a–j). However, immediately after cavitation begins (figure 8 l), the material near the surface experiences additional pressure loading from the shock waves induced by the collapse of the cavitation bubbles (see § 4). Shortly thereafter, the cavitation cloud collapses at off-centre positions, creating a local stress concentration in the material as a secondary loading (figure 8 n). Notably, we observe a significant, though localised, stress peak $S_{\textit{ca}v}$ of approximately $25$ $\mathrm{MPa}$ , due to cavitation at $t\approx 4.09$ $\unicode{x03BC}\mathrm{s}$ (figure 8 p), which is nearly equivalent to the stress peak $S_{\textit{WH}}$ caused by the WH pressure. After the generation of the cavitation-induced stresses $S_{\textit{ca}v}$ , they propagate through the solid material and along its surface, labelled as the $\hat {S}_{\textit{ca}v}$ wavefront. Additionally, because cavitation continues to generate bubbles that collapse within the droplet, the resulting stresses $S_{\textit{ca}v}$ persist from the contact surface, creating a large area of high stress in the material that lasts until the end of the simulation at $t\approx 5.00$ $\unicode{x03BC}\mathrm{s}$ (figure 8 t).

Existing computational frameworks for fatigue analysis of rain erosion on wind turbine blades typically adopt one of two approaches. The first approach is based on the material fatigue ( $S{-}N$ ) curve, which characterises the material lifetime as a function of stress level, together with Miner’s rule (Miner Reference Miner1945) for cumulative damage at each stress level. This framework enables prediction of the total fatigue lifetime under realistic meteorological conditions (Springer Reference Springer1976; Slot et al. Reference Slot, Gelinck, Rentrop and van der heide2015). Among these models, the Springer model (Springer Reference Springer1976; DNVGL-RP-0573 2020; Hao et al. Reference Hao, Domenech and Sanchez2026c ) is the most widely adopted, using the peak impact pressure, namely the WH pressure, as the stress metric for fatigue damage accumulation. Within this framework, our results show that cavitation induces a secondary pressure peak comparable in magnitude to the primary WH pressure, highlighting the potentially significant role of cavitation: even a modest overload may accumulate through the log–log dependence of the material $S{-}N$ curve in repeated impacts, leading to a substantially reduced fatigue lifetime. The second approach assumes identical impacts and uses stress or strain metrics within the solid, obtained from numerical simulations, as the basis for fatigue accumulation (Verma et al. Reference Verma, Castro, Jiang and Teuwen2020). This method allows the total fatigue life, expressed in cycles, to be estimated from the material $S{-}N$ curve. From this perspective, the secondary loading induced by cavitation can subject the material to two separate stress events during a single droplet impact, potentially halving the material’s fatigue lifetime. Taken together, these results underscore the significant influence of cavitation on impact-induced fatigue erosion; however, a detailed quantitative assessment of this effect is beyond the scope of the present study.

5.3. Shear erosion

As a liquid droplet impacts the elastic solid, another notable ring-shaped surface erosion pattern forms between the radius of the loaded area and the Rayleigh wave surface cracks due to shear erosion (Lesser & Field Reference Lesser and Field1983; Obreschkow et al. Reference Obreschkow, Dorsaz, Kobel, Bosset, Tinguely, Field and Farhat2011). This impact erosion mechanism is driven by the shearing motion of the solid as the shear wave propagates radially outwards (Brunton Reference Brunton1966), while the central loaded region undergoes elastic compression (Hancox & Brunton Reference Hancox and Brunton1966). As noted by Bourne et al. (Reference Bourne, Obara and Field1997), shear erosion is most intense at the jet root, where lateral sheet ejection occurs (Riboux & Gordillo Reference Riboux and Gordillo2014). It is also more pronounced in ductile materials (for example, PU coatings considered in this study) than in brittle ones (Brunton Reference Brunton1966; Preece & Brunton Reference Preece and Brunton1980). Furthermore, recent studies have highlighted that, for all types of materials, shearing motion can cause significant erosion, especially in multi-layered coating structures, where delamination at the interfaces between layers becomes a severe issue (Mishnaevsky Reference Mishnaevsky2019; Ibrahim & Medraj Reference Ibrahim and Medraj2020). To assess this, shear stress $\tau _{xz}$ is analysed in figure 9, where its evolution characterises the erosion potential of an impacting raindrop on coating materials used in wind turbine blades.

Figure 9.

Comparison of shear stresses in the two-dimensional solid material between droplet impacts without (a,c,g,i,k,m,o,q,s) and with (b,d, f,h, j,l,n,p,r,t) cavitation at $U_0=50$ $\mathrm{m\,s^{-1}}$ . Each row corresponds to the same time instants as in figure 8. The loaded region of each subfigure is denoted by a red line with the edge $X_{\kern-0.5pt P}$ labelled with an arrow.

Figure 9(a,c,e,g,i,k,m,o,q,s) shows the evolution of shear stress $\tau _{xz}$ during liquid droplet impact at $U_0=50$ $\mathrm{m\,s^{-1}}$ without cavitation. Immediately after impact (figure 9 a), positive shear stresses are generated within the loaded region due to the shearing effect from the liquid spreading on the surface (Bourne et al. Reference Bourne, Obara and Field1997). At the same time, negative shear stresses form outside the loaded region due to the mismatch in load distribution on either side of the boundary at $X_{\kern-0.5pt P}$ . Similar to the propagation of principal stresses, the maximum shear stress $\tau _{\textit{WH}}$ is initially concentrated at (or near) the surface centre, driven by the transient WH pressure peak, and spreads horizontally outwards along the surface as the wavefront $\hat {\tau }_{\textit{WH}}$ over time (figure 9 s). Notably, the decay rates of positive and negative shear stresses differ, with negative stresses decaying more slowly.

In comparison, figure 9(b,d,f,h,j,l,n,p,r,t) illustrates the shear stress evolution in response to cavitation-induced loadings. At $t\approx 2.78$ $\unicode{x03BC}\mathrm{s}$ (figure 9 l), shear stresses from cavitation bubble collapse begin to rise and shortly afterwards reach a significant, though localised, stress peak $\tau _{\textit{ca}v}$ (figure 9 p), comparable in magnitudes to $\tau _{\textit{WH}}$ . However, we note a key and unique feature of cavitation loading in comparison to the droplet impact loading, that is the outward-expanding shear stress pattern: negative shear to the left and positive shear to the right, driven by the shock waves emanating from cavitation collapse at off-centre positions. This results in reverse shear patterns in the solid material, as seen in figure 9(n). Moreover, due to the opposite sign of shear stresses $\tau _{\textit{ca}v}$ generated by the continuous collapse of cavitation near the contact surface, and their propagating wavefront $\hat {\tau }_{\textit{ca}v}$ , adjacent shear stresses merge to form a large sheared area with reversed overall shear behaviour. This effect persists until the end of the simulation at $t\approx 5.00$ $\unicode{x03BC}\mathrm{s}$ (figure 9 t). Thus, cavitation introduces a secondary failure zone off-centre, where the material may experience complex shear erosion due to these reversed shearing patterns.

7.1. Cavitation patterns

Figure 13.

Comparison of the cavitation mechanisms in two-dimensional and axisymmetric scenarios, presented by pressure contours at several typical instants. Condition set-up: ( $a$ ) $U_0=50$ $\mathrm{m\,s^{-1}}$ , two-dimensional; ( $b$ ) $U_0=50$ $\mathrm{m\,s^{-1}}$ , axisymmetric; ( $c$ ) $U_0=70$ $\mathrm{m\,s^{-1}}$ , axisymmetric. Instants: ( $a\,\textrm{i}$ ) 0.874 $\unicode{x03BC}\mathrm{s}$ , ( $b\,\textrm{i}$ ) 0.868 $\unicode{x03BC}\mathrm{s}$ , ( $c\,\textrm{i}$ ) 0.880 $\unicode{x03BC}\mathrm{s}$ , ( $a\,\textrm{ii}$ ) 1.038 $\unicode{x03BC}\mathrm{s}$ , ( $b\,\textrm{ii}$ ) 0.979 $\unicode{x03BC}\mathrm{s}$ , ( $c\,\textrm{ii}$ ) 0.984 $\unicode{x03BC}\mathrm{s}$ , ( $a\,\textrm{iii}$ ) 1.337 $\unicode{x03BC}\mathrm{s}$ , ( $b\,\textrm{iii}$ ) 1.861 $\unicode{x03BC}\mathrm{s}$ , ( $c\,\textrm{iii}$ ) 1.874 $\unicode{x03BC}\mathrm{s}$ . Panels ( $a\,\textrm{iv}$ – $c\,\textrm{iv}$ ) are enlarged images of the noted zones by the dashed lines in ( $a\,\textrm{iii}$ – $c\,\textrm{iii}$ ), respectively. The cavitation zones are noted by the volume fraction isolines of water vapour. The $(r,z)$ coordinates are normalised by the droplet’s initial radius $R_0$ .

In this section we analyse the cavitation dynamics in the axisymmetric configuration and compare the results with two-dimensional simulations to elucidate three-dimensional effects. In this way, the potential to obtain physical insight from two-dimensional computations can be evaluated and the possible limitations in predicting the physics of the three-dimensional liquid flow and cavitation in a droplet can be identified. Figure 13 presents pressure contours overlaid with vapour isolines to visualise cavitation zones in three representative scenarios: ( $a$ ) $U_0=50$ $\mathrm{m\,s^{-1}}$ , two-dimensional; ( $b$ ) $U_0=50$ $\mathrm{m\,s^{-1}}$ , axisymmetric; ( $c$ ) $U_0=70$ $\mathrm{m\,s^{-1}}$ , axisymmetric. The results show that the early stage propagation and reflection of the WH shock wave in axisymmetric cases are similar to observations in the two-dimensional configurations. As the WH shock propagates upwards, it is continuously reflected by the droplet surface, generating the reflected rarefaction waves. These rarefaction waves converge, leading to the onset of bulk homogeneous cavitation at approximately $t\sim 0.870$ $\unicode{x03BC}\mathrm{s}$ . During this process, the wave morphology, properties of different wave branches and focal location remain consistent to those in the two-dimensional cases, as clearly demonstrated in figure 13(ai–ci).

However, three-dimensional effects significantly modify the overall pressure distribution and, consequently, the cavitation patterns within the droplet. When a shock wave of given initial intensity propagates and expands in space, the associated high-pressure region behind it decays more rapidly in three-dimensional configurations than in two-dimensional ones. This is primarily due to the additional spatial degree of freedom in three dimensions, which allows the shock wave to spread more freely and reduces the degree of confinement. This phenomenon, commonly referred to as the three-dimensional relieving effect (Anderson Reference Anderson2010), accounts for the attenuated pressure variations observed in axisymmetric simulations. As a result, pressure changes induced by propagating waves – such as the initial WH shock, its reflected rarefaction waves and the subsequent collapsing shock waves – are less pronounced in axisymmetric cases. This mechanism also explains the absence of near-wall cavitation triggered by the reflected expanding rarefaction wave ( $\mathrm{Re}\text{-}\mathrm{ERRW}$ ).

In contrast, an additional three-dimensional mechanism – the intensified convergence effect – can locally amplify pressure variations when shock or rarefaction waves converge along the axis of symmetry. This effect, also reported by Wu et al. (Reference Wu, Liu and Wang2021), arises from the presence of an extra convergence direction in axisymmetric configurations compared with two-dimensional cases. Attributed to this enhanced geometric focusing, the axisymmetric convergence of the rarefaction wave branch $\mathrm{Re}\text{-}\mathrm{RW}_{\textit{II}}$ induces cavitation near the top pole of the droplet (Wu et al. Reference Wu, Liu and Wang2021), as observed in figure 13(bii,c ii). Notably, the resulting top cavitation zone is attached along the axis of symmetry. Reese, Ohl & Rosselló (Reference Reese, Ohl and Rosselló2024) further suggested that such top cavitation may initiate a liquid micro-jet that penetrates or escapes through the droplet’s top surface.

Under the combined influence of the three-dimensional relieving effect and the intensified convergence effect, the patterns of near-wall cavitation differ significantly between two-dimensional and axisymmetric cases. In two-dimensional scenarios, near-wall cavitation is primarily induced by the pronounced pressure drop behind the wave $\mathrm{Re}\text{-}\mathrm{ERRW}$ , as discussed in § 4. The corresponding cavitation zone exhibits a wall-adherent characteristic as shown in figure 13(a iv). However, this mechanism is absent in axisymmetric cases due to the insufficient intensity of the $\mathrm{Re}\text{-}\mathrm{ERRW}$ wave. Instead, in axisymmetric scenarios, near-wall cavitation arises from a different mechanism facilitated by the intensified convergence effect. The expanding shock wave (figures 13 b ii,c ii) undergoes continuous reflection by the curved droplet surface, producing a series of reflected compression and rarefaction waves. Among these, one branch of the reflected rarefaction waves converges near the wall centre, reducing the local pressure below the heterogeneous cavitation threshold and thereby triggering near-wall cavitation, as demonstrated in figure 13(b iii,c iii). As a result, the induced cavitation zone concentrates along the axis of symmetry instead of spreading along the wall. This distinction implies that the intense wall pressure induced by collapsing near-wall bubbles spans a broader area in two-dimensional cases, whereas it is more concentrated at the centre in axisymmetric (three-dimensional) ones. It is also noted that although the axisymmetric cases with $U_0=50$ $\mathrm{m\,s^{-1}}$ and $U_0=70$ $\mathrm{m\,s^{-1}}$ share the same near-wall cavitation mechanism, the wave propagation and reflection processes are still quantitatively affected by the impact velocity. Therefore, the shapes of cavitation zones and the collapse behaviour of the bubbles in the two cases are not identical.

7.2. Material analysis

Before proceeding to the material analyses, we validate the numerical simulation of the axisymmetric solid. Recently, Sun et al. (Reference Sun, Alvarez-Novoa, Andrade, Gutierrez, Gordillo and Cheng2022) measured the temporal evolution of hydrostatic stress distributions within a cross-linked polydimethylsiloxane (PDMS) gel using high-speed stress microscopy with embedded fluorescent polystyrene particle tracers. These tracers enabled visualisation of the deformation of the elastic gel beneath low-speed impacting droplets. The droplets consisted of a 60 wt % aqueous sodium iodide solution, with a fixed diameter of $3.5$ mm and an impact speed of $2.97$ $\mathrm{m\,s^{-1}}$ . Figure 14 compares the hydrostatic stresses in the axisymmetric solid obtained experimentally by Sun et al. (Reference Sun, Alvarez-Novoa, Andrade, Gutierrez, Gordillo and Cheng2022) with those from the present numerical simulations performed at $U_0=50$ $\mathrm{m\,s^{-1}}$ . Cavitation has not yet initiated over the time interval considered and, therefore, does not influence the comparison. Owing to differences in droplet impact conditions and material properties, the hydrostatic stress is normalised by $\rho _lU_0c_s$ to reflect the hydrodynamic loading scale, and time is normalised by $R_0/c_s$ to represent shock-wave propagation within the solid. The numerical results show good qualitative and quantitative agreement with the experiments. In particular, the simulations capture the region of negative hydrostatic stress in front of the spreading contact line, indicating the formation of a Rayleigh wave, as noted by Sun et al. (Reference Sun, Alvarez-Novoa, Andrade, Gutierrez, Gordillo and Cheng2022).

We note that the PDMS gel used by Sun et al. (Reference Sun, Alvarez-Novoa, Andrade, Gutierrez, Gordillo and Cheng2022) is substantially softer (with stiffness of order ${O}(10^5 \mathrm{Pa})$ ) than the 3M PU coating considered in the present study (table 1) (3M 2008). However, experimental visualisation of internal stress fields requires optically transparent materials, which are typically soft. Another example is provided by Yokoyama et al. (Reference Yokoyama, Mitchell, Nassiri, Kinsey, Korkolis and Tagawa2023), who employed integrated photoelasticity to reconstruct three-dimensional stress fields in soft solids with stiffness of order ${O}(10^4 \mathrm{Pa})$ . Furthermore, both techniques (Sun et al. Reference Sun, Alvarez-Novoa, Andrade, Gutierrez, Gordillo and Cheng2022; Yokoyama et al. Reference Yokoyama, Mitchell, Nassiri, Kinsey, Korkolis and Tagawa2023) are restricted to low-speed impacts (order O( $1$ $\mathrm{m\,s^{-1}}$ )), well below the cavitation threshold. To the best of the authors’ knowledge, no experimental study has yet resolved the evolution of internal stress distributions in solids subject to cavitation-induced surface loading. Such measurements are particularly challenging for millimetre-sized droplets, as they require micron spatial resolutions and submillisecond (or even microsecond, when cavitation is involved) temporal resolution (Cheng, Sun & Gordillo Reference Cheng, Sun and Gordillo2022).

Figure 14.

Comparison of hydrostatic stresses in the axisymmetric solid between the experimental measurements of Sun et al. (Reference Sun, Alvarez-Novoa, Andrade, Gutierrez, Gordillo and Cheng2022) (a–e) and the numerical simulations of the present study at $U_0=50$ $\mathrm{m\,s^{-1}}$ (f–j). The impact times shown in each column correspond to the dimensionless times: $tc_s/R_0=0.0$ (a,f), $1.2$ (b,g), $2.4$ (c,h), $3.6$ (d,i) and $4.8$ (e,j). The solid within each panel has dimensionless length $L_0/R_0=1.06$ and height $H_0/R_0=0.29$ ; the hydrostatic stress is normalised by $\rho _lU_0c_s$ .

Following validation, we analyse the coating’s response to impacts with and without cavitation. Figure 15(left) shows the evolution of shear stress $\tau _{rz}$ in the axisymmetric coating material during liquid droplet impact at $U_0=50$ $\mathrm{m\,s^{-1}}$ without cavitation. Similar to the two-dimensional analysis (§ 5.3), positive shear stress $\tau _{\textit{WH}}$ is generated within the loaded region ( $I$ ), with negative shear stress forming outside the loading boundary $X_{\kern-0.5pt P}$ (figure 15 a). The wavefront $\hat {\tau }_{\textit{WH}}$ propagates faster in time than $X_{\kern-0.5pt P}$ (figure 15 v) and evolves into a positive–negative shear pattern (figure 15 s) near the surface, extending beyond the loaded region ( $I$ ). However, the stress magnitudes in the axisymmetric case are generally an order of magnitude smaller than those observed in the two-dimensional simulations (§ 5), primarily due to the diversified shock waves in axisymmetric droplet impact (Wu et al. Reference Wu, Liu and Wang2021). As a result, the impact loadings from reflected shock waves at the droplet’s bottom surface become more pronounced. We observe stresses $\tau _{\textit{ref}}$ induced by these reflection waves, which are depicted in figures 15(g), 15(m) and 15(p) for three distinct reflection waves (both compressive and tensile). These stresses shows notably high magnitudes at the surface centre, with limited radial extension along the surface or in depth. Nevertheless, due to the alternative signs of the reflection waves (Wu et al. Reference Wu, Xiang and Wang2018), they also influence the dominant (for the displayed time duration of figure 15) shear stress around the loading boundary $X_{\kern-0.5pt P}$ , resulting in an intermittent pattern (figures 15 g and 15 j) with an overall reduced magnitude.

Figure 15.

Comparison of shear stresses in the axisymmetric solid material between droplet impacts without cavitation at $U_0=50$ $\mathrm{m\,s^{-1}}$ (a,e,g, f,m,p,s,v), with cavitation at $U_0=50$ $\mathrm{m\,s^{-1}}$ (b,e,h,k,n,q,t,w) and $U_0=70$ $\mathrm{m\,s^{-1}}$ (c, f,i,l,o,r,u,x). Each row corresponds to comparable but not identical time instants from the two cases, with approximate simulation times for each row as follows: $t\approx$ $0.79$ $\unicode{x03BC}\mathrm{s}$ (a,b,c), $1.31$ $\unicode{x03BC}\mathrm{s}$ (d,e,f), $2.02$ $\unicode{x03BC}\mathrm{s}$ (g,h,i), $2.44$ $\unicode{x03BC}\mathrm{s}$ (j,k,l), $2.86$ $\unicode{x03BC}\mathrm{s}$ (m,n,o), $3.09$ $\unicode{x03BC}\mathrm{s}$ (p,q,r), $3.47$ $\unicode{x03BC}\mathrm{s}$ (s,t,u) and $3.93$ $\unicode{x03BC}\mathrm{s}$ (v,w,x). The loaded region of each subfigure is denoted by a red line with the edge $X_{\kern-0.5pt P}$ labelled by an arrow. The contour limits correspond to the same dimensionless stresses $\tau /(\rho _lU_0^2)=-0.6$ (blue) and $0.6$ (red) in each case.

In contrast, at the same impact speed of $U_0=50$ $\mathrm{m\,s^{-1}}$ , with cavitation, the homogeneous cavitation located far from the surface dissipates the shock-wave energy. As a result, the amount and intensity of the shock waves reflected to the droplet’s bottom surface are reduced. The evolution of shear stress $\tau _{rz}$ in the axisymmetric coating material thus forms a complete shear-stressed zone around the loading boundary $X_{\kern-0.5pt P}$ , without the interference from reflected waves (figure 15 w). Instead, upon the arrival of an isolated, strong reflected tensile wave ( $\tau _{\textit{ref}}$ in figure 15 k) at the surface, heterogeneous cavitation occurs immediately afterwards. This leads to an additional localised (positive) shear stress $\tau _{\textit{ca}v}$ induced by the loadings from surface cavity collapses (figures 15 n and 15 q). This positive shear stress follows the negative shear $\tau _{\textit{ref}}$ in propagation, resulting in a negative–positive shear pattern, as shown in figure 15(t). However, we note that, compared with the pattern in figure 9(n), the negative–positive shear pattern in the axisymmetric case is of much smaller magnitudes, with more limited spreading in space, especially in the radial direction. For a higher impact speed ( $U_0=70$ $\mathrm{m\,s^{-1}}$ , figure 15, right-hand side plots), a smoother and more complete shear pattern is observed around the edge $X_{\kern-0.5pt P}$ of the loaded area. This is due to the stronger homogeneous cavitation, which absorbs more shock-wave energy at the focal point of the droplet (Wu et al. Reference Wu, Xiang and Wang2018). Since the heterogeneous cavitation mechanisms remain qualitatively the same at lower and higher impact speeds (see § 7.1), the same negative–positive shear pattern, as seen at $U_0=50$ $\mathrm{m\,s^{-1}}$ , is observed in figure 15(u) and corresponds to surface cavity collapses. In general, the shear motion – and thus the shear erosion potential – in axisymmetric material analysis is still greater in cavitation cases than in non-cavitation cases due to the homogeneous cavities absorbing the energies of reflected shock waves within the droplet; while the effect of heterogeneous cavities is focused at the surface centre due to the geometric converging effect in axisymmetric droplets (Wu et al. Reference Wu, Liu and Wang2021).

Figure 16.

von Mises stress distributions in the axisymmetric material analysis at specific monitoring points: (a) surface centre (point C in figure 3), (b) half-radius off-centre on the surface (point R’ in figure 3), and (c) subsurface location, a quarter-radius beneath the centre (point Z’ in figure 3). The time and stress values are in the dimensionless forms (except time in panel a).

To characterise the influential area of the focused surface cavities at the surface centre, we carry out a detailed analysis of the axisymmetric material response at three monitoring positions, as in § 6. Figure 16(a) shows the temporal distributions of von Mises stress ( $M_s$ ) at the surface centre (point C) of the axisymmetric solid during droplet impact at $U_0=50$ $\mathrm{m\,s^{-1}}$ and $U_0=70$ $\mathrm{m\,s^{-1}}$ , with and without cavitation. Due to the different impact speeds, physical quantities, namely stress magnitudes and wave propagation times (except figure 16(a) at surface centre), are compared in the dimensionless forms in figure 16. Upon impact, von Mises stresses $M_{s,\textit{WH},U_0}$ induced by the confined WH shock wave reach values of $19.5$ $\mathrm{MPa}$ at $U_0=50$ $\mathrm{m\,s^{-1}}$ and $40.6$ $\mathrm{MPa}$ at $U_0=70$ $\mathrm{m\,s^{-1}}$ . After the primary stress peak, wave reflections arriving at the bottom surface of the droplet induces a local stress peak, $M_{s,\textit{ref},50}$ , with a moderate intensity of $14.0$ $\mathrm{MPa}$ . In contrast, when cavitation is considered, the von Mises stresses generated by surface cavity collapses following the reflected waves (see § 7.1) reach secondary peak values of $M_{s,cav,50}=19.7$ $\mathrm{MPa}$ and $M_{s,cav,70}=35.0$ $\mathrm{MPa}$ at the surface centre. These values are comparable to their corresponding primary stress peaks, $M_{s,\textit{WH},U_0}$ , induced by the WH shock wave. We note, in the axisymmetric material analysis, the von Mises stress induced by surface cavity collapse at $U_0=70$ $\mathrm{m\,s^{-1}}$ is near the ultimate tensile strength ( $\textit{UTS}_s$ ) of $37$ $\mathrm{MPa}$ for the W4600 PU coating (table 1). Although impact loadings in practice are likely to be mitigated by several factors – including droplet shape (Zhang et al. Reference Zhang, Zhang, Lv, Li and Guo2019), impact angle (Cimpeanu & Papageorgiou Reference Cimpeanu and Papageorgiou2018), bubble entrapment (Thoroddsen et al. Reference Thoroddsen, Etoh, Takehara, Ootsuka and Hatsuki2005) and the existence of a liquid film on the coating surface (Tropea & Marengo Reference Tropea and Marengo1999) – and these factors influence cavitation dynamics (Xu et al. Reference Xu, Fan, Wu, Wen and Wang2023; Mur et al. Reference Mur, Bußmann, Paula, Adami, Adams, Petkovsek and Ohl2025), we see the potential of cavitation to cause material surface damage upon a single raindrop impact at speeds less than $100$ $\mathrm{m\,s^{-1}}$ .

Figure 16(b) shows the corresponding temporal distributions of von Mises stress ( $M_s$ ) at half-radius off-centre (point R’) on the axisymmetric solid surface for the three cases. At the onset of the impact event, there is a brief delay in stress initiation, with a time duration of $t_{R'}=2t_{Z'}\approx 0.5R_0/c_s$ , calculated based on the wave speed of the propagating wavefront in the solid (6.1). Also, we observe that the stress peak $\hat {M}_{s,\textit{WH},U_0}$ induced by the WH wave occurs at approximately $0.5R_0/U_0$ , which relates to the spreading morphology (Riboux & Gordillo Reference Riboux and Gordillo2016), and hence the speed, of the impacting droplet. In this way, by combining both the impact speed $U_0$ and the sonic speed $c_s$ in the solid, we select the dimensionless time scale $t/(R_0/\sqrt {U_oc_s})$ , which collapses all time series from different impact speeds onto a single master curve, as shown in figure 16(b,c). The important observation, however, is the negligible contribution of cavitation to local stresses at off-centre positions on the surface (figure 16 b). Similarly, at the subsurface (figure 16 c), the cavitation-induced stresses are attenuated to moderate magnitudes, of the same order as $\hat {M}_{s,\textit{ref},U_0}$ induced by the reflected waves in the droplet. This finding is consistent with the observations in figure 15, where surface-cavity-induced stresses have minimal propagation into the solid. In contrast to the results obtained from the two-dimensional analyses in § 6, the diverging effect of impact waves in the axisymmetric droplet is more pronounced (Wu et al. Reference Wu, Liu and Wang2021), leading to cavitation being primarily focused at the surface centre, with limited radial or depth effects along the surface.

The limited influential zone of surface cavities in axisymmetric simulations is further evidenced by the temporal distributions of material displacement at the surface centre (figure 17 a) and the hydrostatic force across the entire surface (figure 17 b). In figure 17(a) we normalise the surface displacements $D$ using the corresponding Bernoulli impact force ratio (i.e. $\rho _lU_0^2/\rho _lU_{\textit{ref}}^2$ ) with a reference speed, $U_{\textit{ref}}$ , of $50$ $\mathrm{m\,s^{-1}}$ , so that the displacement distributions at different impact speeds collapse onto a single master curve, with values comparable to those in figure 12(a). Additionally, we note the identical time scales in figure 17 in contrast to the previous figure 16, as the current macroscopic behaviours are primarily governed by wave propagation in the solid. For figure 17(b), we derive the hydrostatic force by integrating the hydrostatic stresses ( $P_s$ ) over the axisymmetric material surface (see (6.2) for difference), i.e.

(7.1) \begin{align} F_s(t)=\int _0^{2\pi }\int _0^{L_s}P_s(r,t)r\text{d}R\,\text{d}\theta , \end{align}

where $\theta$ is the azimuthal angle in axisymmetric coordinates. Consistent with the trend observed in figures 15 and 16, the magnitudes in the axisymmetric simulations are generally one order of magnitude smaller than those in the two-dimensional simulations. Comparisons of both material displacement and force in figure 17 – between cavitation (under both low- and high-speed cavitation mechanisms) and non-cavitation cases – show indistinguishable contributions from the additional cavitation-induced loadings at the solid surface. Following from the argument in Corollary 1, the negligible effects of cavitation on material displacement and force in figure 17 indicate that the influential zone of surface cavities in the axisymmetric case remains largely confined to the surface centre, with little impact spreading along the surface.

Figure 17.

Distributions of (a) material displacement at the surface centre (point C in figure 3), and (b) hydrostatic force on the surface, in the axisymmetric material analysis. Note the secondary vertical (right) axes, labelled in green, for the case $U_0=70$ $\mathrm{m\,s^{-1}}$ . The displacement (a) and force (b) limits correspond to dimensionless values $D/(U_0^2/U_{\textit{ref}}^2)=-1.2$ to $0.2$ , and $F_s/(\rho _lU_0^2R_0^2)=3.0$ , respectively, in each case.

Overall, compared with the results from the two-dimensional analysis (§§ 5 and 6), the high-fidelity axisymmetric simulations reveal substantially lower magnitudes and more realistic loading conditions during droplet impact (§ 7.1), due to the intensified converging and diverging wave motions in axisymmetric configurations. As a result, the formation of heterogeneous cavitation and subsequent cavitation dynamics, including collapse events, are strongly localised at the surface centre, where the converging waves focus. At this location, we see significant material stresses, as shown in figure 16(a), arising from the localised pressure fields (figure 13 $b\,\textrm{iii}$ , $c\,\textrm{iii}$ ) generated by surface cavity collapses. This localised stress could accelerate fatigue failure in the material subject to duplex stress events in repeated raindrop impacts, or even cause damage on material surfaces from a single raindrop impact. Furthermore, homogeneous cavitation absorbs shock-wave energy in the droplet away from the surface, resulting in a continuous, highly sheared, region in the material (figure 15). This shear may contribute to material failure through shear erosion in ductile materials (Brunton Reference Brunton1966; Preece & Brunton Reference Preece and Brunton1980) or delamination in multi-layered structures (Ibrahim & Medraj Reference Ibrahim and Medraj2020). Fortunately, due to the intensified diverging wave motions in axisymmetric configurations, heterogeneous cavitation remains localised at the surface centre, with rapid decay and minimal influence radially along the surface or at deeper layers (figures 16 and 17). At a quarter-radius beneath the surface centre (point Z’), the intensities of cavitation-induced stresses are of the same order of magnitudes produced by the reflected wavefronts within the droplet (figure 16 c), which have moderate magnitudes in comparison to the primary stress peaks induced by the WH loading.