1. Introduction
The impact of droplets on solid surfaces is a common phenomenon in both nature and engineering, seen in processes such as inkjet printing (Derby Reference Derby2010; Lohse Reference Lohse2022), raindrop impact (Abuku et al. Reference Abuku, Janssen, Poesen and Roels2009; Gart et al. Reference Gart, Mates, Megaridis and Jung2015; Gilet & Bourouiba Reference Gilet and Bourouiba2015; Kim et al. Reference Kim, Wu, Esmaili, Dombroskie and Jung2020), droplet freezing (Mishchenko et al. Reference Mishchenko, Hatton, Bahadur, Taylor, Krupenkin and Aizenberg2010; Jung et al. Reference Jung, Tiwari, Doan and Poulikakos2012) and pesticide deposition (Wirth, Storp & Jacobsen Reference Wirth, Storp and Jacobsen1991; Bergeron et al. Reference Bergeron, Bonn, Martin and Vovelle2000). After impact, droplet dynamics can be divided into two stages: spreading and recoiling. In the spreading stage, the kinetic energy is converted into surface energy and dissipated through viscosity, until the droplet reaches its maximum spread. During the recoiling stage, surface tension drives the droplet to retract, converting the surface energy back into kinetic energy, with some energy lost through viscous dissipation. Droplets impacting on a dry rigid substrate can exhibit different modes, including deposition, prompt splashing, corona splashing, receding break-up, partial rebound or complete rebound (Rioboo, Tropea & Marengo Reference Rioboo, Tropea and Marengo2001). The outcome depends on the physical parameters, such as droplet viscosity, impact velocity, surface tension coefficient, solid surface roughness and surface wettability (Josserand & Thoroddsen Reference Josserand and Thoroddsen2016).
Previous research is mainly focused on dynamics of a droplet impacting on rigid surfaces. However, many real-world surfaces, such as leaves, umbrellas, clothing fabrics, are flexible. The effects of substrate flexibility on droplet impact dynamics are still less understood, though the limited studies provide insights. For example, Weisensee et al. (Reference Weisensee, Tian, Miljkovic and King2016) observed that flexible substrates can rebound, returning kinetic energy to the droplet and reducing the contact time by up to 50 %, a phenomenon termed the ‘springboard effect’. Huang et al. (Reference Huang, Song, Wang, Zhao, Liu and Liu2018) found that flexible, superhydrophobic cotton significantly reduces the retraction time and overall contact duration once the impact velocity exceeds a threshold, due to energy stored in the substrate during droplet spreading. Other studies (e.g. Mangili et al. Reference Mangili, Antonini, Marengo and Amirfazli2012; Chen et al. Reference Chen, Bonaccurso, Deng and Zhang2016; Vasileiou et al. Reference Vasileiou, Gerber, Prautzsch, Schutzius and Poulikakos2016; Vasileiou, Schutzius & Poulikakos Reference Vasileiou, Schutzius and Poulikakos2017), mostly using experimental methods, explored similar energy storage and release mechanisms, with flexible substrates enhancing droplet rebound and modifying the conditions required for various impact behaviours.
The role of flexibility in reducing splashing has also been investigated. Pepper, Courbin & Stone (Reference Pepper, Courbin and Stone2008) and Howland et al. (Reference Howland, Antkowiak, Castrejón-Pita, Howison, Oliver, Style and Castrejón-Pita2016) observed that flexible substrates suppress droplet splashing, likely due to the early substrate deformation absorbing energy. Their experiments showed that splashing on the softest substrate requires approximately 70 % more kinetic energy than on a rigid surface. Further, Pegg, Purvis & Korobkin (Reference Pegg, Purvis and Korobkin2018) conducted analytical and numerical studies on the splashing dynamics when a droplet impacts a small flexible plate. They found that the oscillation of the flexible substrate can amplify the splash jet, a phenomenon not observed on rigid surfaces. Similarly, Dressaire et al. (Reference Dressaire, Sauret, Boulogne and Stone2016) investigated droplet impact on a thin flexible fibre and found that the fibre's flexibility affects the threshold capture velocity in a nonlinear way.
Although experimental studies provide valuable insights, numerical simulations offer a more detailed and comprehensive analysis of droplet impact dynamics. However, most existing research on droplet impacts on flexible substrates is experimental, with very few numerical studies available. Xiong, Huang & Lu (Reference Xiong, Huang and Lu2020) conducted two-dimensional (2-D) numerical simulations of droplet impacts on a flexible substrate simply supported at both ends, investigating the effects of Weber numbers within a limited range on droplet impact dynamics. However, these 2-D simulations showed significant discrepancies when compared with actual droplet impact behaviour. In a more recent study, Ma & Huang (Reference Ma and Huang2023) simulated droplet impacts on both 2-D and three-dimensional flexible rectangular substrates, also simply supported at both ends. Their research focused on the influence of Weber number, substrate bending stiffness and solid-to-liquid density ratio on the maximum spreading factor, proposing a scaling law for it. However, they did not analyse other dynamic characteristics influenced by substrate flexibility, such as impact forces. Currently, the dynamic behaviour of droplets impacting a flexible substrate remains unclear, and the underlying mechanisms governing droplet–flexible substrate interactions are still not well understood.
In this study, we numerically simulate a droplet impacting a simply supported flexible circular disk, differing from the rectangular substrates examined by Xiong et al. (Reference Xiong, Huang and Lu2020) and Ma & Huang (Reference Ma and Huang2023). The simulations are conducted within a 2-D axisymmetric framework across a wide range of Weber numbers (We = 1–500). We focus on the effects of the Weber number and substrate bending stiffness on droplet impact dynamics, providing a comprehensive analysis of droplet–flexible substrate interactions. This includes impact modes, maximum spreading factor (βmax), spreading time, scaling laws for βmax, impact forces, energy conversion and more. While the primary analysis focuses on hydrophobic surfaces (contact angle θ = 130°), we also investigate the behaviour of droplets impacting a hydrophilic flexible substrate (θ = 60°) towards the end of this paper.
The remainder of this paper is organised as follows: § 2 describes the computational model, including the physical problem, mathematical formulation, numerical methods and numerical validations. Section 3 presents results and discussions on the dynamics of droplets impacting on a flexible disk. Finally, concluding remarks are provided in § 4.
2. Computational model
2.1. Physical model and mathematical formulation
The schematic of a droplet impacting on a flexible circular plate is shown in figure 1. At the initial moment, a spherical droplet with a diameter D 0, density ρl and viscosity 
$\mu_l$ has a downward impact velocity U 0; directly below the droplet at a distance H (= 0.5D 0), a flat flexible disk with thickness h, diameter L, density ρs and Young's modulus E is simply supported at its edge. The droplet is surrounded by a gas with density ρg and viscosity 
$\mu_g$. In numerical simulations, the gravity is not considered. After impacting, the droplet spreads on the deformed substrate. In this process, the droplet impact gives rise to the deformation and the motion of the substrate, which, in turn, affects the droplet dynamics. As shown in figure 1(b), the spreading factor (β) is defined. At the centre of the substrate, the droplet's thickness and the substrate's deflection are denoted by y and d, respectively. A 2-D axisymmetric model is adopted and the gravity is neglected. The letters r and z represent the radial and axial coordinates, respectively, and the substrate is initially located at z = 0.
Schematic of a droplet impacting on a flexible circular disk: (a) initial configuration; (b) an instant of the droplet spreading on the deformed substrate.
In this system, the fluid flow is governed by the incompressible Navier–Stokes equations,
\begin{gather}{\rho _f}\frac{{\partial {\boldsymbol{v}_f}}}{{\partial t}} + {\rho _f}\textrm{(}{\boldsymbol{v}_f}\boldsymbol{\cdot }\boldsymbol{\nabla }\textrm{)}{\boldsymbol{v}_f} = \boldsymbol{\nabla }\boldsymbol{\cdot }[{-}p\boldsymbol{\mathsf{I}} + \mu \textrm{(}\boldsymbol{\nabla }{\boldsymbol{v}_f} + {\textrm{(}\boldsymbol{\nabla }{\boldsymbol{v}_f}\textrm{)}^\textrm{T}}\textrm{)]} + {\boldsymbol{F}_{st}} + {\rho _f}\boldsymbol{g},\end{gather}
\begin{gather}\boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{v}_f} = \textrm{0,}\end{gather}
where ρf is the fluid density, vf is the fluid velocity vector, 
$\boldsymbol{\mathsf{I}}$ is the unit tensor, 
$\mu_f$ is the fluid dynamic viscosity, p is the pressure, Fst is the surface tension force term of the fluid and t is the time. The flexible substrate is modelled as a shell, the thickness of which is negligible. The deformation and motion of the shell are described by the following structural equations:
\begin{equation}{\rho _s}\frac{{{\partial ^2}{\boldsymbol{u}_s}}}{{\partial {t^2}}} = \boldsymbol{\nabla }\boldsymbol{\cdot }{\textrm{(}\boldsymbol{\mathsf{FS}}\textrm{)}^\textrm{T}} + {\rho _s}\boldsymbol{g},\end{equation}
where ρs is the structure density and us is the displacement vector. In our study, gravity is neglected, i.e. g = 0. The deformation gradient tensor 
$\boldsymbol{\mathsf{F}}$ and the second Piola–Kirchhoff stress tensor 
$\boldsymbol{\mathsf{S}}$ are defined as
\begin{gather}\boldsymbol{\mathsf{F}} = \boldsymbol{\mathsf{I}} + \boldsymbol{\nabla }{\boldsymbol{u}_{s}},\end{gather}
\begin{gather}\boldsymbol{\mathsf{S}} = J{\boldsymbol{\mathsf{F}}^{ - 1}}{\sigma_s}{\boldsymbol{\mathsf{F}}^{ - T}},\end{gather}
respectively, where 
$\sigma_s$ is the Cauchy stress tensor and the volume factor J is given by
\begin{equation}J = \textrm{det}(\boldsymbol{\mathsf{F}}) = \frac{{\textrm{d}V}}{{\textrm{d}{V_0}}},\end{equation}where V 0 is the volume before deformation and V is the volume after deformation. The strain tensor corresponding to the second type of Piola–Kirchhoff stress tensor is the Green–Lagrangian strain tensor:
\begin{equation}\boldsymbol{\mathsf{E}} =
{\textstyle{1 \over 2}}({\boldsymbol{\mathsf{F}}^T}\boldsymbol{\mathsf{F}} - \boldsymbol{\mathsf{I}}).\end{equation}
The flexible disk in this problem is an isotropic linear elastic material, then 
$\boldsymbol{\mathsf{S}}$ and 
$\boldsymbol{\mathsf{E}}$ satisfy Hooke's law:
\begin{equation}\boldsymbol{\mathsf{S}} = \lambda
(\textrm{tr}\,\boldsymbol{\mathsf{E}})\boldsymbol{\mathsf{I}} + 2G\boldsymbol{\mathsf{E}},\end{equation}
where 
$\lambda$ and G are the Lamé coefficient and shear modulus, respectively, which can be calculated from the Young's modulus E and Poisson's ratio v of the material:
\begin{gather}\lambda = \frac{{vE}}{{(1 + v){\kern 1pt} (1 - 2v){\kern 1pt} }},\end{gather}
\begin{gather}G = \frac{E}{{2(1 + v)}}.\end{gather}For the fluid–structure interaction (FSI) boundary, the conditions of velocity compatibility and stress balance must be satisfied:
\begin{gather}\frac{{\partial {\boldsymbol{u}_s}}}{{\partial t}} = {\boldsymbol{v}_f},\end{gather}
\begin{gather}\boldsymbol{n}\boldsymbol{\cdot}{\sigma_s} = \boldsymbol{n}\boldsymbol{\cdot }{\sigma_f},\end{gather}
where the unit vector n is normal to the solid boundary, and fluid stress 
$\sigma_f$ is given by
\begin{equation}{\sigma_f}
={-} p\boldsymbol{\mathsf{I}} + {\mu _f}[\boldsymbol{\nabla
}{\boldsymbol{v}_f} + {(\boldsymbol{\nabla
}{\boldsymbol{v}_f})^\textrm{T}}]{\kern 1pt}
.\end{equation}
In the present study, we choose ρl, D 0 and σ as characteristic quantities. The corresponding reference speed 
${U_{ref}} = \sqrt {\sigma /\textrm{(}{\rho _l}{D_0}\textrm{)}}$. Thus, the following dimensionless governing parameters are given: the Weber number 
$We = {\rho _l}U_0^2{D_0}/\sigma$, the Reynolds number 
$Re = {\rho _l}{U_0}{D_0}/{\mu _l}$, the liquid-to-gas density ratio ρr = ρl/ρg, the liquid-to-gas viscosity ratio 
$\mu_r=\mu_l/\mu_g$, the bending stiffness 
$K = E{h^3}/({\rho _l}U_{ref}^2{L^3})$, the thickness of the flexible disk h* = h/D 0, the disk diameter L* = L/D 0 and the substrate-to-liquid density ratio dr = ρs/ρl. The subscripts s, l and g denote the solid, the liquid and the gas, respectively. The dimensionless time t* = tUref /D 0.
2.2. Numerical method
The numerical simulations are performed using software COMSOL Multiphysics (COMSOL 5.6), which provides a comprehensive platform for simulating FSI phenomena using the finite element method (FEM). In COMSOL Multiphysics, the governing equations for fluid and shell dynamics, i.e. (2.1)–(2.3) are converted into their weak form. The weak form equations are then discretised using FEM, with special attention given to mesh deformation to accurately capture the evolving interface between fluid and structure. The mesh for shell and fluid domains in COMSOL are provided in Appendix A. The fluid region is divided by the triangular elements. In the present model, the thickness of the shell is negligible. The variables of the shell, e.g. displacement vector, are defined on the midplane (z = 0). The shell elements are located on the midplane of the shell. In the simulation, the dynamic mesh deformation is required to account for changes in the geometry of both fluid and shell domains. We can specify mesh refinement criteria and employ adaptive meshing techniques to maintain computational accuracy while capturing complex fluid–structure interactions. Time discretisation is performed using the built-in backward differentiation formula in COMSOL. All variables are fully coupled using the constant Newton method (Deuflhard Reference Deuflhard1974), with a constant damping factor of 0.75. The Yeoh smoothing method (Yeoh Reference Yeoh1993) is employed to constrain the mesh displacement, allowing for the maximum displacement before the mesh is inverted. The MUMPS solver is used for solving the linear systems of the coupled equations during the problem initialisation and transient process. For details regarding finite element discretisation, meshing techniques and other technical aspects, please refer to the COMSOL User's Guide.
In the present simulations, the phase-field method is used to track the gas–liquid interface. In the phase-field method, the surface tension force is calculated based on the following expression:
\begin{gather}{\boldsymbol{F}_{st}} = G\boldsymbol{\nabla }\phi ,\end{gather}
\begin{gather}G = \lambda \left[ { - {\kern 1pt} {\nabla^2}\phi + \frac{{\phi \textrm{(}{\phi^2} - 1\textrm{)}}}{{{\varepsilon^2}}}} \right] = \frac{\lambda }{{{\varepsilon ^2}}}\psi ,\end{gather}where G is the chemical potential, γ is the mobility, λ is the mixing energy density, ε is the interface thickness parameter and ψ is the phase-field auxiliary variable. The diffusion interface is defined by a dimensionless phase-field variable ϕ, ranging from –1 to 1. Here, ϕ = −1 and ϕ = 1 correspond to the gas and liquid, respectively, with densities of ρg and ρl. To track the diffusion interface of two-phase fluids, the following phase-field equation (Cahn–Hilliard equation) is used:
\begin{gather}{\kern 1pt} {\kern 1pt} \frac{{\partial \phi }}{{\partial t}} + {\boldsymbol{v}_f}\boldsymbol{\cdot }\boldsymbol{\nabla }\phi = \boldsymbol{\nabla }\boldsymbol{\cdot }\frac{{\gamma \lambda }}{{{\varepsilon ^2}}}\boldsymbol{\nabla }\psi ,\end{gather}
\begin{gather}\psi ={-} \boldsymbol{\nabla }\boldsymbol{\cdot }{\varepsilon ^2}\boldsymbol{\nabla }\phi + \textrm{(}{\phi ^2} - 1\textrm{)}\phi .\end{gather}Here, σ, λ and ε have the following relationship:
\begin{equation}\sigma = \frac{{2\sqrt 2 \lambda }}{{3\varepsilon }}.\end{equation}In our simulation, we set the interface thickness ε to 0.75 hm, where hm represents the maximum grid size in the region through which the droplet passes. Additionally, γ determines the time scale of diffusion and is closely related to the interface thickness, i.e.
\begin{equation}\gamma = \chi {\varepsilon ^2},\end{equation}whereas χ represents the parameter for adjusting the mobility, which needs to be chosen precisely. If the value is too large, it excessively suppresses the convection term, and if it is too small, it becomes difficult to maintain the correct interface thickness. In our simulation, the value of the mobility adjustment parameter is determined by the following function:
\begin{equation}\chi = \frac{{8{V_{max}}}}{{9\sqrt 2 \sigma }},\end{equation}where Vmax is the maximum velocity in the entire computational domain at the current moment.
The phase-field method treats two-phase flow as a single-phase flow, where the properties in the flow domain such as density ρ and viscosity μ vary based on the volume fractions of different fluid phases,
\begin{gather}\rho = {\rho _g} + \textrm{(}{\rho _l} - {\rho _g}\textrm{)}{V_l},\end{gather}
\begin{gather}\mu = {\mu _g} + \textrm{(}{\mu _l} - {\mu _g}\textrm{)}{V_l}.\end{gather}The volume fractions of the light and heavy fluid phases are related through the phase-field variable ϕ, i.e.
\begin{equation}{V_g} = \frac{{1 - \phi }}{2},\quad {V_\textrm{l}} = \frac{{1 + \phi }}{2}.\end{equation}Open boundary conditions are applied to the upper, lower and lateral sides of the computational domain, while a wetted wall boundary condition is imposed on the substrate, which sets the mass flow rate passing through the wall to zero and specifies the contact angle of the liquid on the wall. This is described by the following equations:
\begin{gather}\boldsymbol{n}\boldsymbol{\cdot }{\varepsilon ^2}\boldsymbol{\nabla }\phi = {\varepsilon ^2}\,\textrm{co}{\textrm{s}^{ - 1}}\,\theta {\bf |}\boldsymbol{\nabla }\phi {\bf |}{\kern 1pt} ,\end{gather}
\begin{gather}\boldsymbol{n}\boldsymbol{\cdot }\frac{{\sigma \lambda }}{{{\varepsilon ^2}}}\boldsymbol{\nabla }\psi = 0,\end{gather}where θ is the static contact angle of the droplet on the disk substrate.
2.3. Numerical validation
Numerical simulations for the typical cases are carried out to quantitatively validate the reliability of the numerical method adopted in the present study. First, we have conducted the simulations for various contact angles within a realistic hysteresis range (i.e. 100° ± 5°) which is consistent with the experiments of Sikalo et al. (Reference Sikalo, Marengo, Tropea and Ganic2002). According to the experiments of Sikalo et al. (Reference Sikalo, Marengo, Tropea and Ganic2002), the key parameters in the simulations are set as follows: initial velocity U 0 = 1.54 m s−1, surface tension coefficient σ = 0.073 N m−1, viscosity 
$\mu_l=1\ \textrm{mPa}\cdot \textrm{s}$, density ρl = 996 kg m−3, diameter D 0 = 2.7 mm. The numerical simulations using the static contact angle model are conducted. According to the experiments of Sikalo et al. (Reference Sikalo, Marengo, Tropea and Ganic2002), three static contact angles (SCA), i.e. 95°, 100° and 105°, are set. In addition, we conducted the numerical simulations using the quasi-dynamic contact angle model, in which the advancing and receding contact angles are 105° and 95°, respectively. Figures 2(a) and 2(b) show the variations of the dimensionless spreading diameter and apex height as a function of dimensionless time, respectively. As shown in figure 2, it can be observed that within the contact angle hysteresis range considered (i.e. 100° ± 5°), the experimental results and the numerical results using the static and quasi-dynamic contact angle models are generally in good agreement. Also, we can see from the numerical results that the discrepancies in the spreading factor and substrate deflection caused by using different contact angle models is not significant. In the present study, we use a static contact angle model for simulations, which is acceptable for some realistic situations where the contact angle hysteresis effect is not significant (Sikalo et al. Reference Sikalo, Marengo, Tropea and Ganic2002).
Time variations of (a) the spreading factor β and (b) the apex height y of the droplet impacting on a rigid wall. The results are obtained from experiments (Sikalo et al. Reference Sikalo, Marengo, Tropea and Ganic2002) and the present simulations using static and quasi-dynamic contact angle models. Three static contact angles (SCA) are set in the static contact angle model, i.e. SCA = 95°, 100° and 105°. In the quasi-dynamic contact angle model, the advancing and receding contact angles are 105° and 95°, respectively, which is denoted as ‘CA= 105°/95°’ in panel (a). The other key parameters in the simulations are set according to the experiments (Sikalo et al. Reference Sikalo, Marengo, Tropea and Ganic2002) as follows: initial velocity U 0 = 1.54 m s−1, surface tension coefficient σ = 0.073 N m−1, viscosity 
$\mu_l=1\ \textrm{mPa}\cdot \textrm{s}$, density ρl = 996 kg m−3 and diameter D 0 = 2.7 mm.
To validate the reliability of FSI calculations involving flexible structures, we conduct numerical simulations of the flow around a rigid cylinder with a flexible strip immersed in a uniform viscous flow. The results are provided in Appendix B. Further, we compared the present results with those of Ma & Huang (Reference Ma and Huang2023). In the simulations of Ma & Huang (Reference Ma and Huang2023), a liquid column impacting on the simply supported flexible plate has been considered; and the key dimensionless parameters are set as follows: the bending stiffness 
${K_B} = EI/({\rho _l}U_{ref}^2{L^3}) = 0.01$, 0.05, the stretching stiffness 
${K_S} = Eh/({\rho _\textrm{l}}U_{ref}^2L) = 100$, the Ohnesorge number 
$Oh = \mu /\sqrt {{\rho _l}{D_\textrm{0}}\sigma } = 0.01$, the liquid-to-gas density ratio ρr = 1000, the liquid-to-gas viscosity ratio 
$\mu_r=50$, the substrate span length L/D 0 = 20, the mass ratio Mr = ρsh/ρlL = 0.01 and the contact angle θ = 170°. Figure 3 shows the maximum spreading factor βmax as a function of We. It is seen that for both rigid and flexible substrate cases, the present results and the results of Ma & Huang (Reference Ma and Huang2023) are generally in good agreement. This illustrates the reliability of our numerical method in solving this multiphase FSI problem. It is noted that the 2-D axisymmetric model with the hydrophobic surfaces (θ = 130°) are considered in the present simulations, which is different from the study of Ma & Huang (Reference Ma and Huang2023) and may cause the data discrepancies.
Maximum spreading factor βmax as a function of We in the present and Ma & Huang's (Reference Ma and Huang2023) simulations.
The grid independence analysis is performed. Figure 4(a,b) shows the time-dependent d and β calculated under three sets of mesh sizes, i.e. Δx = 3D 0/200, 1.5 D 0/200 and D 0/200. Here, the mesh size Δx denotes the maximum mesh size within the refined region at the initial moment. Correspondingly, the numbers of elements in the refined region are 64 560, 224 158 and 480 072, respectively. Since the unrefined domain is composed of air, variations in mesh size in this domain have negligible effects on the droplet–substrate coupling, and thus we limit the maximum size of the rough domain to 0.15D 0 in mesh independence analysis. It is confirmed that Δx/D 0 = 1.5/200 is sufficient to achieve accurate results in the present simulations. Based on our convergence studies with different computational domains, the computational domain R × Z is chosen as 3.5D 0 × 6.5D 0. The domain is large enough so that the blocking effects of the boundaries are not significant.
Grid independence test for We = 50, K = 0.01, dr = 1.25, h* = 0.03, L* = 6, ρr = 1000, 
$\mu_r=100$, v = 0.45 and θ = 130°. (a) Substrate centre deflection and (b) spreading factor.
3. Results and discussion
Here, we present some typical results of a droplet impacting on a flexible substrate. In our study, we focus on the effects of We and K on the droplet impact dynamics, while keeping other parameters constant unless otherwise specified. Specifically, 
$We \in [1,500]$, 
$K \in [0.01,1000]$, Re = 1000, dr = 1.25, h* = 0.03, L* = 6, ρr = 1000, 
$\mu_r=100$, v = 0.45 and θ = 130°. Since the deformation of the substrate can be ignored during impact, the flexible substrate with K = 1000 is regarded as the rigid wall in the following sections.
3.1. Impact mode and regime map
In this section, we present the dynamics of a droplet impacting on a flexible substrate. It is found that compared with the rigid case, the droplet impacting on the flexible substrate has more complex dynamics due to the fluid–structure interaction. Based on a variety of simulations for a wide range of parameters considered here, we have identified five typical modes in terms of the droplet shape evolution during the impact, i.e. the ring-shaped rebound (RSR), the disc-shaped rebound (DSR), the jet-breakup rebound (JBR), the splashing rebound (SR) and edge-breakup rebound (EBR). It is worth noting that since the hydrophobic surface is considered in this section, the droplet bounces completely off the wall in all modes; according to whether the droplets are separated or broken during the rebound process, we can summarise the impact modes into two categories, namely, the intact rebound (IR, e.g. RSR and DSR) and breakup rebound (BR, e.g. JBR, SR or EBR). Note that the splashing phenomenon is considered as a special case of the droplet separation.
Figure 5 shows snapshots of a drop impacting on the flexible substrate for the cases of the typical modes. As shown in figure 5(a), for the case where We = 10 and K = 1000, corresponding to the RSR mode, the droplet spreading factor continues to increase in the spreading stage, while the apex height continues to decrease, the central part of the droplet gradually thins to form a liquid film or even breaks, and a liquid ring is formed when t* = 0.56. Under the action of surface tension, the liquid ring retracts towards the centre. The contraction of the liquid ring leads to the regeneration and the complete rebound of the droplet.
Snapshots of a droplet impacting on a flexible substrate for typical cases: (a) We = 10 and K = 1000; (b) We = 10 and K = 0.01; (c) We = 50 and K = 0.1; (d) We = 500 and K = 1000 and (e) We = 500 and K = 0.1, corresponding to the RSR, DSR, JBR, SR and EBR modes, respectively.
For the case where We = 10 and K = 0.01, corresponding to the DSR mode as shown in figure 5(b), the kinetic energy of the droplet is partially converted into the deformation energy and kinetic energy of the flexible substrate, which leads to a smaller effective We number for droplet impact. Therefore, compared with the case where We = 10 and K = 1000 (figure 5a), the formation of liquid ring is not observed during the droplet spreading stage; instead, the droplet evolves into a disc-like morphology in the late stage of spreading. Compared with the rigid case, when the flexible substrate moves upwards, partial energy of the substrate is converted into the kinetic energy of the droplet, resulting in a liquid film that has an upward velocity (at t* = 0.65), which does not appear in the previous example.
For the case where We = 50 and K = 0.1, corresponding to the JBR mode as shown in figure 5(c), a thin liquid film and thick rim are formed at the centre and the outer edge, respectively. Due to the lager We, spreading will cause the accumulation of the liquid medium at the edges of the droplet. As the substrate rebounds, the edges of the droplet acquire upward velocity, resulting in the conversion of the kinetic energy to the surface energy. This leads to a slight longitudinal stretching and a contraction of the contact line at the edges, partially detaching from the substrate and forming a rising thick rim (at t* = 0.47). After merging, there is a tendency for jetting, resulting in the detachment of a small satellite droplet from the main body of the droplet at the tip (at t* = 0.81). The detached satellite droplet has a greater rising velocity than the main body of the droplet detaching from the substrate.
For the case where We = 500 and K = 1000, corresponding to the SR mode as shown in figure 5(d), after the impact, a non-wetting thin lamella appears at the edge of the droplet (at t* = 0.04). Under the effect of high We, the surface tension is not sufficient to restrain the rapid spreading of the droplet. The capillary instability occurs, making the droplet experience a prompt splashing (at t* = 0.07).
For the case where We = 500 and K = 0.1, corresponding to the EBR mode as shown in figure 5(e), due to the absorption of a portion of the droplet's kinetic energy by the substrate, splashing does not occur although the Weber number is equal to that for the case in figure 5(d) (at t* = 0.08). This phenomenon is similar to the findings of Howland et al. (Reference Howland, Antkowiak, Castrejón-Pita, Howison, Oliver, Style and Castrejón-Pita2016). However, during the subsequent substrate rebound, the thin lamella at the edge acquires an upward velocity, leading to longitudinal stretching, while the contact line retracts. When the substrate reaches its maximum height and starts moving downward, the moving direction of the thin lamella is opposite to that of the central liquid film (at t* = 0.19). This intensifies the longitudinal stretching of the thin lamella, eventually leading to the necking phenomenon and the separation of a liquid ring from the main film. After separation, under the effect of surface tension, the liquid ring begins to contract towards the centre and merges with the liquid film on the rising substrate (at t* = 0.29). The merging process creates a local protrusion on the film during retraction (at t* = 0.40), and the ring-shaped protrusion moves to the centre and finally encloses with an air bubble trapped (at t* = 0.83). The merged droplet longitudinally stretches and undergoes necking due to the velocity differences at different heights. Eventually, the liquid column separates into multiple satellite droplets of varying sizes.
The distribution of the above five modes in a parameter space (We, K) is shown in figure 6. The intact rebound (IR, e.g. RSR or DSR) occurs at 1 ≤ We ≤ 10, while the breakup rebound (BR, e.g. JBR, SR and EBR) occurs at 25 ≤ We ≤ 500.
Five impact modes in a parameter space (K, We) for θ = 130°. Circles and triangles denote the intact rebound (IR) and breakup rebound (BR), respectively. Green circles, blue circles, red triangles, orange triangles and black triangles denote the ring-shaped rebound (RSR), the disc-shaped rebound (DSR), the jet-breakup rebound (JBR), the splashing rebound (SR) and edge-breakup rebound (EBR), respectively.
3.2. Maximum spreading factor and spreading time
Regarding the droplet impact on a solid surface, the maximum spreading factor has always been a research hotspot. In the following, we will focus on the effect of the substrate flexibility on the maximum spreading factor (βmax).
From figure 7(a), it can be observed that as K gradually decreases, βmax exhibits two opposite trends for We < 10 and 10 < We ≤ 250, respectively. For the low-We region (We < 10), as K decreases from 1000 (rigid) to 0.01 (very flexible), βmax shows a weak increasing trend. While for the higher We region (10 < We ≤ 250), βmax first shows a slight increase when K decreases from 1000 to 1, and then decreases significantly when K further decreases from 1 to 0.01. The variation trend of βmax with K for 10 ≤ We ≤ 250, 0.01 ≤ K ≤ 1 is consistent with those reported by Ma & Huang (Reference Ma and Huang2023). The weakened spreading induced by the substrate flexibility is attributed to the fact that the disk absorbs a part of the droplet's kinetic energy for spreading. In contrast, it is observed in figure 7(a) that βmax slightly increases as K decreases for 1 ≤ We ≤ 5 and 0.01 ≤ K ≤ 1. This spreading enhancement due to the substrate flexibility has not been reported by Ma & Huang (Reference Ma and Huang2023) or in other relevant studies. Figure 7(b) shows the βmax as a function of We for the cases of K = 0.01, 3 and 1000, providing a clearer demonstration of the trends in βmax. It is obvious that the spreading enhancement (weakness) occurs for K = 0.01, We < 10 (We > 10). Figure 8(a) demonstrates the time variation curves of the spreading factor β for We = 5 and K = 0.01, 1000. Compared with the rigid case (K = 1000), there is a time lag in the rise of the curve for the case with K = 0.01, due to the wetting hysteresis caused by the downward movement of substrate. As the substrate moves downward, the β is smaller than that for the rigid case in the early stage of droplet spreading. However, as the substrate moves upward, β reaches the maximum value which is larger than that for the rigid case. To understand the changing of the spreading factor, figure 8(d) shows the snapshots from the moment of maximum substrate energy to the moment of maximum spreading factor for the case with We = 5, K = 0.01, contrasted with the snapshot for the rigid case at the certain moment (t* = 0.467) before reaching βmax. At t* = 0.57, due to impacting, the liquid inside the droplet flows from the top to the edge (marked with a blue arrow in figure 8d). As the substrate moves upward, it compresses the bottom of the droplet, and at t* = 0.60 as shown in figure 8(d), both the top and bottom parts of droplet exhibit a tendency to flow to the edge. The flow at the top part becomes weaker, while the flow at the bottom part (marked with a green arrow in figure 8d) becomes stronger as the droplet spreads. After t* = 0.63, the flow to the edge is entirely driven by the upward movement of the substrate. Therefore, the larger βmax at K = 0.01 compared with the rigid case is caused by the upward motion of the substrate in the late stage of droplet spreading. In conclusion, compared with the rigid case, the squeezing effect of substrate on the droplet caused by the substrate's upward movement contributes to the droplet spreading, which is manifested by the extension of spreading time and the increase of the maximum spreading factor.
Variation of the maximum spreading factor βmax (a) with K for different We and (b) with We for different K.
Spreading factor for the cases at (a) We = 5 and K = 0.01 and 1000, (b) We = 100, K = 3 and 1000, and (c) We = 50, K = 0.01 and 1000. The blue dash-dotted line represents the total energy of the substrate. The orange and green arrows indicate the relative increase and reduction of the spreading factor during the substrate's downward and upward motions, respectively. The shaded regions in panels (a), (b) and (c) correspond to the time interval from the moment when the substrate energy reaches its peak to the moment when the droplet achieves maximum spreading. (d) Snapshots of the droplet during the shaded region for K = 0.01 in panel (a). The lower right corner is a snapshot of the droplet at the eve of maximum spread at K = 1000. The blue and green arrows indicate two different flows.
Accordingly, the variations of the spreading factor in figure 8(b,c) can be explained. For figure 8(b), the difference from figure 8(a) is that higher bending stiffness leads to faster oscillation, resulting in the moment of the first downward and upward movement of substrate being significantly earlier than the moment of the maximum spreading. In this case, the squeezing of the substrate on the droplet makes the maximum spreading factor of the droplet increase, but has no obvious effect on the spreading time. For figure 8(c), spreading enhancement caused by the substrate's upward movement is smaller than spreading weakening caused by the substrate's downward movement, Therefore, the net changes in the maximum spreading factor for K = 0.01 are negative, i.e. the maximum spreading factor decreases compared with that for the rigid case.
It is noted that the differences between the present study and the study of Ma & Huang (Reference Ma and Huang2023) are the vibration behaviour of the substrate. In the present study, the substrate vibration is affected by both We and K, which is slightly different from that observed by Ma & Huang (Reference Ma and Huang2023). In the present study, the moment of the first valley of the substrate deflection curve occurs always earlier than the moment of maximum spreading factor. However, the moment of the first valley of the substrate deflection curve can be earlier or later than the moment of maximum spreading factor reported by Ma & Huang (Reference Ma and Huang2023). This difference may be attributed to the different types of the substrates. In our study, a disc-shaped substrate with edge support is used. Additionally, the different parameter scopes considered may also have an influence.
The analysis above clarifies why βmax is greater or smaller than that for the rigid substrate cases. However, it does not fully explain the trends of βmax observed in figure 7(b). To address this point, we quantify the decreasing and increasing effects on the spreading factor during the downward and upward stage of the substrate. The relative spreading weakening caused by the downward movement and the relative spreading enhancement caused by the upward movement, i.e. Δβ 1 (note that Δβ 1 < 0) and Δβ 2, are defined, respectively. These quantities are indicated by the green and orange arrows in figure 8(a–c), respectively. Thus, the net change in spreading factor, Δβ = Δβ 1 + Δβ 2.
Figure 9 shows the changes in the spreading factor as a function of We for K = 0.01 and K = 3. As shown in figure 9(a), as We increases, the spreading weakens, i.e. |Δβ 1|, gradually increases due to the increasing impact energy; however, the spreading enhancement, i.e. Δβ 2, has no significant change. There are possible reasons for these observations. As mentioned above, Δβ 2 is caused by the upward movement of the substrate. Thus, a higher We does not necessarily imply greater Δβ 2 during the upward movement. Whether the droplet continues to spread or retract during the upward stage of the substrate depends on the relative magnitude between droplet kinetic energy and surface energy. However, it is obvious that as We increases, a larger initial impact energy of the droplet will promote the downward movement of the substrate, which results in a larger Δβ 2. The different trends of curves for Δβ 1 and Δβ 2 result in Δβ being positive when We < 10 and negative when We > 10, which corresponds to the spreading enhancement and weakening.
Relative weakening (Δβ 1) and enhancement (Δβ 2) of spreading factor caused by substrate's downward and upward movement, as well as the net effect (Δβ), for two bending stiffness: (a) K = 0.01 and (b) K = 3.
For the cases for K = 3 shown in figure 9(b), both the spreading enhancement and spreading weakening increase as We increases. The spreading enhancement is slightly larger than the spreading weakening, resulting in the positive Δβ which increases with We. This indicates that the spreading enhancement due to the substrate rebound dominates the droplet spreading for the high-We region.
Based on the previous analysis, the upward movement of the substrate during the late spreading stage may prolong the spreading stage. The spreading time, 
${T_{{\beta _{max}}}}$, represents the duration from when the droplet begins to contact the substrate until it achieves maximum spreading. Figure 10 illustrates the variation of 
${T_{{\beta _{max}}}}$ with We and K. Only for We ≤ 25, the spreading time gradually increases with the decreasing K, particularly for K ≤ 10; while for We > 25, 
${T_{{\beta _{max}}}}$ remains almost constant with K.
Effect of K on the spreading time 
${T_{{\beta _{max}}}}$ for different We.
3.3. Scaling law
The maximum spreading factor (βmax) describes the maximum extent of spreading of a liquid droplet upon impact with a solid surface. It is crucial for understanding the interaction between the droplet and the surface, making it an important physical quantity for characterising the dynamics of droplet impact on surfaces. Therefore, how to theoretically predict βmax is a question of widespread interest.
For the case of a liquid droplet impacting a rigid wall, various theoretical prediction models for βmax have been developed based on the conversion relationship among droplet kinetic energy, viscous dissipation and surface energy. For example, in the viscous regime where inertial and viscous forces dominate, research by Chandra & Avedisian (Reference Chandra and Avedisian1991) indicates βmax ∼ Re 1/5; while in the capillary regime where viscosity may be neglected, Collings et al. (Reference Collings, Markworth, McCoy and Saunders1990) obtain another scaling law, βmax ∼ We 1/2, by assuming that all the kinetic energy is transformed into surface energy at the maximal spreading. However, an alternative scaling has been suggested for the same regime. Based on momentum conservation, Clanet et al. (Reference Clanet, Béguin, Richard and Quéré2004) proposed that βmax ∼ We 1/4. In more scenarios, the influences of viscous dissipation or surface tension on spreading cannot be simply excluded, otherwise it will result in the breakdown of the above scaling laws. By interpolating between the We 1/2 and Re 1/5 scaling, Eggers et al. (Reference Eggers, Fontelos, Josserand and Zaleski2010) showed that in the cross-over regime between capillary and viscous regimes, 
${\beta _{max}} \propto R{e^{1/5}}{f_c}(P)$, where fc is a function of the impact parameter P (= We Re −2/5) that varies between zero (capillary regime) and infinity (viscous regime). Further, Laan et al. (Reference Laan, Bruin, Bartolo, Josserand and Bonn2014) constructed the smooth cross-over between these two asymptotes by using a so-called Padé approximant, i.e. 
${\beta _{max}}\,R{e^{ - 1/5}} = {P^{1/2}}/(A + {P^{1/2}})$, where A is a fitting constant. Lee et al. (Reference Lee, Laan, de Bruin, Skantzaris, Skantzaris, Derome, Carmeliet and Bonn2016) extended Laan's approach (Laan et al. Reference Laan, Bruin, Bartolo, Josserand and Bonn2014) by considering wettability. It is found that, after a correction for dynamic wetting using the maximum spreading ratio at zero velocity 
${\beta _{V \to 0}}$, the maximum spreading should still scale with We 1/2, i.e. 
${(\beta _{max}^2 - \beta _{V \to 0}^2)^{1/2}} \sim W{e^{1/2}}$. Building on this conclusion, and considering that Re is constant in this study, a universally applicable prediction model of βmax can be derived by following Laan's approach:
\begin{equation}{(\beta _{max}^2 - \beta _{V \to 0}^2)^{1/2}}\,R{e^{ - 1/5}} = \frac{{W{e^{1/2}}}}{{B + W{e^{1/2}}}},\end{equation}
where B is a fitting coefficient which may vary across different Re cases. However, in the present study, B is considered a constant coefficient due to the fixed value of Re. Also, 
${\beta _{V \to 0}}$ can be calculated from the contact angle (Lee et al. Reference Lee, Laan, de Bruin, Skantzaris, Skantzaris, Derome, Carmeliet and Bonn2016),
\begin{equation}
{\beta _{V \to 0}}=\left\{\begin{array}{@{}ll} {{{\left(
{\dfrac{{4\,s\textrm{i}{\textrm{n}^3}\,{\theta_{V \to
0}}}}{{2 - 3\,\textrm{cos}\,{\theta_{V \to 0}} +
\textrm{co}{\textrm{s}^3}\,{\theta_{V \to 0}}}}}
\right)}^{1/3}}}&{\textrm{if}\;{\theta_{V \to 0}} <
90^\circ \textrm{,}}\\ {{{\left( {\dfrac{1}{{{(2 +
}\textrm{cos}\,{\theta_{V \to
0}})\,\textrm{si}{\textrm{n}^4}({\theta_{V \to
0}}/2)}}{\kern 1pt} }
\right)}^{1/3}}}&{\textrm{if}\;{\theta_{V \to 0}} >
90^\circ .} \end{array}
\right.\end{equation}
For flexible substrate scenario, theoretical predictions of βmax become more challenging because vibration and deformation of substrate will significantly affect the spreading of a droplet. Based on Lee et al.'s (Reference Lee, Laan, de Bruin, Skantzaris, Skantzaris, Derome, Carmeliet and Bonn2016) theoretical model for rigid substrate cases, Ma & Huang (Reference Ma and Huang2023) have proposed a scaling law of βmax for a droplet impacting on a flexible plate simply supported at its both edges, by defining an effective Weber number, i.e.
\begin{equation}W{e_m} = \frac{{We}}{{1 + {\delta _{max}}}}\textrm{,}\end{equation}where δmax is the substrate deflection. Due to momentum conservation, δmax can be predicted as
\begin{equation}{\delta _{max}} = \frac{1}{{2{\rm \pi} }}\frac{m}{{m + {M_s}/3}}\frac{{{U_0}}}{{f{D_0}}},\end{equation}
in which m and M are the mass of the droplet and the flexible substrate, respectively, U 0 is the initial velocity of the droplet, 
$f = {f_0}{[{M_s}/(m + {M_s})]^{{1 / 2}}}$ is the vibration frequency of the system, and 
${f_0} = {\rm \pi}\sqrt {EI/{M_s}{L^3}} /2$ is the natural vibration frequency of a simple-support Euler–Bernoulli beam on both ends. The effective Weber number Wem takes all the factors affecting βmax into account.
Following the prediction model for βmax proposed by Ma & Huang (Reference Ma and Huang2023), the numerical results obtained in the present study have been plotted in the 
$W{e_m} - {\beta _m}$ plane, where 
${\beta _m} = {(\beta _{max}^2 - \beta _0^2)^{1/2}}\,R{e^{ - 1/5}}$ and 
${\beta _0} = {\beta _{V \to 0}}$, as shown in figure 11. In figure 11(a), the effective Weber number Wem is calculated by the theoretical natural frequency of the Euler–Bernoulli beam f 0; the data symbols are spread over the parametric plane for varying K, indicating a poor data collapse. The reason is that the natural frequency of the shell in the present study cannot be accurately predicted with f 0. Since it is difficult to predict the natural frequencies of shells theoretically, we determine the natural frequencies of shells with different stiffnesses K through numerical experiments. As shown in figure 11(b), when the numerically measured natural frequencies of the shells are used to predict Wem by (3.3)–(3.4), the data collapse together. Also, we can see that all data distribute above the solid line plotted by (3.1), indicating that the Wem is slightly underestimated. This is because the maximum deflection δmax is no longer the only characteristic quantity that defines Wem due to the emergence of distinct mechanisms of shell–droplet coupling in the present problem. In the following, according to the revealed physical mechanisms, Wem has been modified to seek a better data collapse.
The 
${\beta _m} = {(\beta _{max}^2 - \beta _0^2)^{1/2}}R{e^{ - 1/5}}$ as a function of Wem for varying stiffness K. (a) Wem is calculated by the theoretical natural frequency of the Euler–Bernoulli beam f 0; (b) Wem is calculated based on the natural frequency of the present shell obtained through numerical experiments. The solid line is plotted by (3.1) and B is 7.6.
From the analysis of figures 8 and 9, the spreading process of a droplet on a flexible substrate can be divided into two stages. The first stage is from when the droplet starts to contact the substrate until the substrate energy reaches the maximum. In this stage, the energy is transferred from the droplet to the substrate which continues to dent downward. In the second stage, from when the substrate energy reaches the peak until the droplet reaches its maximum spreading, the substrate continues to rebound and the substrate energy is transferred to the droplet.
During the first stage, the effective Weber number decreases due to the depression of the substrate and thus the droplet spreading is suppressed; while during the second stage, the effective Weber number increases due to the rebound of the substrate and the droplet spreading is promoted. Therefore, the effect of flexible substrate on the maximum spreading factor depends on the competition of these two mechanisms. Specifically, when |Δβ 1| < Δβ 2 (see figure 8a,b), Δβ = Δβ 1 + Δβ 2 = Δβ 2 − |Δβ 1| > 0, indicating spreading enhancement dominated by the ‘flexible rebound mechanism’; while when |Δβ 1| > Δβ 2 (see figure 8c), Δβ = Δβ 1 + Δβ 2 = Δβ 2 − |Δβ 1| < 0, indicating spreading suppression dominated by ‘flexible compliant mechanism’. It is worth noting that the spreading suppression dominated by the flexible compliant mechanism has been found in previous studies (Vasileiou et al. Reference Vasileiou, Gerber, Prautzsch, Schutzius and Poulikakos2016; Xiong et al. Reference Xiong, Huang and Lu2020; Ma & Huang Reference Ma and Huang2023); however, the spreading enhancement caused by the flexible rebound mechanism has not been reported before.
In the scaling law proposed by Ma & Huang (Reference Ma and Huang2023), the maximum deflection δmax has been used to correct the We, that is, (3.3), which is appropriate for the region where the flexible compliance mechanism dominates. However, in the region where the flexible rebound mechanism dominates, the effective Weber number increases and is related to the rebound displacement of the substrate, δreb, and thus (3.3) is no longer applicable. Instead, δreb has been used to predict the effective Weber number (Wem), i.e. Wem = We/(1 + δreb) for Δβ > 0. Therefore, (3.3) is modified to
\begin{equation}W{e_m} = \frac{{We}}{{1
+ \delta }}\textrm{,}\quad \delta = \left\{
{\begin{array}{@{}ll}
{{\delta_{max}}}&{\textrm{if}\;\Delta \beta > 0\textrm{,}}\\ {{\delta_{reb}}}&{\textrm{if}\;\Delta
\beta < 0\textrm{,}} \end{array}}
\right.\end{equation}
where δreb (note that δreb < 0) and δmax can be obtained from the numerical results. Based on the modified Wem calculated by (3.5), it is shown in figure 12 that 
${\beta _m} = {(\beta _{max}^2 - \beta _0^2)^{1/2}}\,R{e^{ - 1/5}}$ as a function of Wem for varying stiffness K. It is seen that all data are collapsed on the same solid line plotted by (3.1), indicating that the modified Wem calculated by using (3.5) leads to a better data collapse compared with those in figure 11(a,b).
The 
${\beta _m} = {(\beta _{max}^2 - \beta _0^2)^{1/2}}R{e^{ - 1/5}}$ as a function of Wem for varying stiffness K. The modified Wem is calculated by (3.5).
3.4. Impact force
According to the study by Zhang et al. (Reference Zhang, Sanjay, Shi, Zhao, Lv, Feng and Lohse2022) on the impact force for a droplet impacting on a superhydrophobic rigid surface, an inertial shock will lead to the appearance of the first peak of the impact force when the maximum droplet diameter is approximately equal to the initial diameter. The magnitude of this peak is independent of the wettability of the surface. When the droplet retracts, the momentum conservation results in the formation of an upward jet as well as an internal downward jet, which leads to the appearance of the second force peak. However, few studies on the time variation of impact force for flexible substrate situations have been reported. In this section, we focus on the impact force during droplet impact on flexible substrates.
Figure 13(a) makes a comparison between the previous experimental and present numerical results. The parameters in the simulations are identical to those of Zhang et al.'s (Reference Zhang, Sanjay, Shi, Zhao, Lv, Feng and Lohse2022) experiments, except for the contact angle of the substrate. For the same reasons mentioned as in § 2.3, a static contact angle of 130° is used in the simulations. As shown in figure 13(a), the timing and amplitude of the first force peak in the simulations coincide generally with the results from Zhang et al. (Reference Zhang, Sanjay, Shi, Zhao, Lv, Feng and Lohse2022). The impact force reaches a valley near the moment of maximum spreading, before the appearance of the second force peak in the retraction stage of the droplet. The occurrence of the three characteristic moments (corresponding droplet shape snapshots are given in figure 13a) is consistent with the report of Zhang et al. (Reference Zhang, Sanjay, Shi, Zhao, Lv, Feng and Lohse2022). The difference in the contact angle may be responsible for the difference in the valley value of the impact force and the timing of the second peak. The impact force for K = 1000 almost overlaps with that for the rigid case, indicating that the calculation of impact force in our model is reliable.
(a) Validation of impact force on the substrate: the thick solid black line represents the experimental results from Zhang et al. (Reference Zhang, Sanjay, Shi, Zhao, Lv, Feng and Lohse2022), the blue and red solid lines represent the simulation results on the rigid substrate and the substrate with K = 1000, and the red dashed line represents the spreading factor. Evolution of impact force for (b) We = 5 and K = 0.1, (c) We = 50 and K = 0.1 and (d) We = 500 and K = 0.1. In panel (b–d), the blue dashed line represents the deflection curve of the substrate, the thick solid red line represents the impact force on the flexible substrate and the red dashed line represents the impact force on the rigid substrate. Each arrow points to the moment corresponding to the snapshot of the droplet.
Further, we select the cases with K = 0.1, We = 5, 50 and 500 for impact force analysis, which are shown in figures 13(b)–13(d), respectively. For comparison, we calculate the impact forces on a rigid substrate under the same parameters. However, in some cases with K = 1000, the substrate exhibits subtle high-frequency vibrations, which result in severe spikes in the curve of the impact force. Thus, it is not suitable to approximate the substrate with K = 1000 as a rigid substrate for comparison. In this case, we directly compare the impact forces on the flexible substrates with that on the rigid wall.
For the case with We = 5 in figure 13(b), the observed slope changes in the deflection curve in figure 13(b) correspond to the peak and valley of the impact force. Compared with the rigid substrate cases, the impact force curve for the flexible substrate exhibits more frequent fluctuations. During the time period from t* = 0.2 to t* = 0.4, the substrate undergoes the downward movement due to the droplet impact, resulting in a lower impact force on the flexible substrate compared with the rigid wall at the same time. At the moment of two deflection valleys, the impact force reaches two prominent peaks, which are significantly larger than those on the rigid wall.
It can be observed from figure 13(c) that, similar to the case of We = 5, significant peaks of the impact force occur at the two moments corresponding to the valleys of the deflection curve for We = 50. However, in this case, the substrate undergoes upward deflection beyond its initial position during the upward movement stage. As the substrate moves upward, the impact force rapidly decreases, changing from positive to negative, and then reaches a valley near the moment of the deflection peak. This indicates that the effect of the droplet on the substrate gradually shifts from compression to adsorption as the substrate moves upward beyond its initial position, which is a phenomenon not observed in the case of impacting on a rigid wall.
As shown in figure 13(d), when We increases to 500, the substrate exhibits strong oscillations near the equilibrium position at a higher frequency, which leads to alternating peaks and valleys in the load curve, corresponding to the valleys and peaks of the substrate deflection curve, respectively.
3.5. Energy conversion
When a droplet impacts on a hydrophobic rigid wall, the conversion of the kinetic energy of the droplet and its surface energy plays an important role in spreading dynamics of the droplet. However, when the substrate is flexible, it is evident that the energy of the substrate (including the elastic strain energy and the kinetic energy) participates in the energy conversion. In this section, we will discuss the energy conversion in the droplet–substrate coupling. We focus on the interaction between the kinetic energy of the droplet (Ek), the surface energy (Es) and the total energy of the substrate (Et). In the present simulations, Et is calculated by the sum of the substrate kinetic energy and potential energy, 
${E_k} = {\textstyle{1 \over 2}}\int {\rho {v^2}\,\textrm{d}V}$ and 
${E_s} = \sigma \textrm{(}{A_{LV}} - {A_{SL}}\,\textrm{cos}\,\theta \textrm{)}$, where ALV and ASL are the liquid–vapor and liquid–solid contact area, respectively (Du et al. Reference Du, Chamakos, Papathanasiou, Li and Min2021).
The typical processes of the energy interaction are summarised in figure 14. For comparison, figure 14(a) shows the time variations of Ek and Es for K = 1000 (rigid case), where the substrate energy does not participate in the energy interaction. Figures 14(b)–14(d) show the three typical energy interaction processes, respectively, while the corresponding snapshots at key moments are shown in figure 14(e). These energy interaction processes correspond basically to the three types of changes in impact force in figure 13. As shown in figure 14(b), after impact, Ek starts to convert into Es and Et. At t* = b 1, the substrate deflection reaches a peak and, subsequently, the upward motion of the substrate induces an upward velocity of the droplet, resulting in a small increase of droplet kinetic energy (Ek). Near t* = b 3, Es reaches the maximum and the substrate is near its initial position, leading to Et approaching zero. Subsequently, the droplet starts to retract and Es is converted back into Ek. At the moment t* = b 4, Ek reaches another peak. During retraction, the downward momentum within the droplet leads to another impact on the substrate, resulting in the conversion of Ek into Et. At t* = b 5, Et reaches the second peak; meanwhile, the deflection reaches a peak. The occurrences of two peaks of Et (b 2 and b 5) are attributed to the initial impact and the subsequent retracting impact, respectively. In figure 14(c), it is shown that after passing the initial position, the droplet briefly adsorbs the substrate and moves upward for a distance (see figure 13c), which results in the conversion of Ek into Et and Es, and the second peak of Et appears at t* = c 3. The subsequent changes in the energies are similar to those in figure 14(b). In this case, Et exhibits three peaks due to the short absorption of the droplet to the substrate. Figure 14(d) shows the mutual conversion between Ek and Es, which corresponds to the results in figure 13(d).
Time evolution of energy for four typical cases: (a) We = 30 and K = 1000; (b) We = 25 and K = 0.03; (c) We = 50 and K = 0.03 and (d) We = 250 and K = 0.1, and (e) corresponding snapshots at labelled moments in panel (b–d). In panel (a–d), the red solid, blue dash-dotted and black dashed lines represent the droplet's kinetic energy (Ek), the droplet's surface energy (Es) and the total energy of the substrate (Et), respectively.
3.6. Impact dynamics on hydrophilic substrates
To investigate the effect of hydrophilic substrates on the droplet impact dynamics, we set the contact angle to 60°. Figure 15 shows the impact modes in a parameter space (K, We) for θ = 60°. According to whether the droplet completely deposits on the substrate, the impact modes are classified into two main types: deposition (DP) and partial rebound (PR), corresponding to the down-pointing and up-pointing triangles in figure 15, respectively. Based on the droplet shape evolution, we can subdivide the DP mode into complete deposition (CD) and splashing deposition (SD). For the PR mode, we can subdivide it into single-satellite partial rebound (SSPR) and multi-satellites partial rebound (MSPR).
Four impact modes in a parameter space (K, We) for θ = 60°. Black down-pointing triangles, orange down-pointing triangles, red up-pointing triangles and blue up-pointing triangles denote the complete deposition (CD), splashing deposition (SD), single-satellite partial rebound (SSPR) and multi-satellites partial rebound (MSPR) modes, respectively.
Figure 16 presents snapshots for the typical case for each mode. For the case where We = 25 and K = 1000 (see figure 16a), corresponding to the CD mode, the droplet is ultimately deposited on the substrate (at t* = 1.87) after a period of oscillations.
Snapshots of typical cases for four impact modes: (a) We = 25 and K = 1000; (b) We = 25 and K = 0.01; (c) We = 500 and K = 0.05; and (d) We = 500 and K = 1000, corresponding to the CD, SSPR, MSPR and SD modes, respectively.
For the case where We = 25 and K = 0.01 (see figure 16b), corresponding to the SSPR mode, it is shown that the obvious difference from the case in figure 16(a) is that a satellite droplet is separated during the droplet retraction stage, which is attributed to the rebound of the substrate. During the period t* = 1.0–1.38, the droplet is retracting; meanwhile, the droplet apex height is increasing. Also, in the second rebound stage of substrate, the momentum transfers from the substrate to the droplet, leading to a greater velocity in the top part of the droplet. When the substrate moves downward again (t* = 1.38–1.89), the velocities of the top and bottom parts of the droplet are opposite and the liquid column continues to deform vertically. Finally, the liquid column begins to neck and eventually separates a satellite droplet (t* = 1.89–2.04).
For the case where We = 500 and K = 0.05 (see figure 16c), corresponding to the MSPR mode, multiple satellite droplets are separated from the mother droplet. The vibration of the substrate is more intense so that an upright liquid ring is formed directly in the region where the medium is gathered (t* = 0.43). This is followed by the merging of the liquid ring into the central liquid column. The velocity difference between different parts of the elongated liquid column leads to the formation of multiple necks and the subsequent generation of multiple separated satellite droplets (t* = 1.26).
For the case where We = 500 and K = 1000 (see figure 16d), corresponding to the SD mode, the droplet exhibits a prompt splashing in the initial stage of impact (t* = 0.047), then experiences a period of oscillations and ultimately deposits on the substrate (t* = 1.87).
5. Conclusions
We have conducted a series of simulations of droplet impact on a flexible circular disk and have focused on the effects of Weber number (
$We \in [1,500]$), substrate stiffness (
$K \in [0.01,1000]$) and contact angle (θ = 130°, 60°) on the droplet–substrate coupling dynamics. According to whether partial droplet separation occurs during impact process, we have distinguished two categories of the impacting modes for the hydrophobic flexible substrate scenario, namely, the intact rebound (IR, e.g. RSR and DSR) and breakup rebound (BR, e.g. JBR, SR or EBR) modes. The effects of the substrate flexibility on the droplet spreading dynamics have been discussed through the time-variation curves of the maximum spreading factor and spreading time, i.e. βmax and 
${T_{{\beta _{max}}}}$. It has been revealed that the substrate rebound will squeeze the fluid at the bottom of the droplet and induces flows from the bottom to the edges, leading to the enhanced spreading of the droplet. Also, when We ≤ 25 and K ≤ 10, and the bending stiffness is moderate (K ∼ O(1)) or lower, the upward movement of the substrate during the late spreading stage will prolongs the spreading time. Further, the spreading process of a droplet on a flexible substrate has been divided into two stages, during which the droplet spreading is suppressed or promoted due to the flexible compliant and rebound mechanisms, respectively; the effect of a flexible substrate on the maximum spreading factor depends on the competition between the two mechanisms. Based on this, a modified scaling law of βmax has been proposed by introducing the effective Weber number (Wem). The impact force on a flexible substrate is analysed and we summarise three typical impact force variation processes. We have observed that the extreme values in the deflection curves correspond to the peaks and valleys of the impact force curves. The effect of the droplet on the substrate gradually shifts from compression to adsorption as the substrate moves upward beyond its initial position. Furthermore, we have performed the analysis of energy and investigated the interaction among droplet kinetic energy, surface energy and total energy of the substrate, and we summarise three typical energy interaction processes, which correspond basically to the three impact force changes, respectively. Finally, the effect of the hydrophilic substrate on the droplet impact dynamics has been examined. Four typical impact modes, i.e. CD, SD, SSPR and MSPR modes, have been classified in the parameter space (K, We) for θ = 60°.
Funding
This research was supported in part by the National Natural Science Foundation of China (nos 12372263, 12372241) and the Natural Science Foundation of Hubei Province of China (no. 2022CFB131).
Declaration of interests
The authors report no conflict of interest.
Appendix A. The mesh for shell and fluid domains
Figure 17 illustrates the local mesh for both shell and fluid domains. The fluid domain is discretised using triangular meshes. To accurately capture the motions of the droplet and the gas around it, the mesh in the region where the droplet moves is refined. In the refined fluid region, the sizes of most triangular elements are approximately equal, and the number and connectivity of elements remain unchanged during the droplet impact process. Moreover, as shown in figure 17, the red solid circles and the blue solid lines connecting them represent the shell elements, which are located on the midplane of the shell. In the simulation, the dynamic mesh deformation is required to account for changes in the geometry of both fluid and shell domains. We can specify mesh refinement criteria and employ adaptive meshing techniques to maintain computational accuracy while capturing complex fluid–structure interactions.
Mesh for shell and fluid domains. The green solid circles and the green solid lines connecting them represent the shell elements, which is located on the midplane of the shell.
Appendix B. A validation to check the validity of FSI calculations
Figure 18(a) shows the schematic physical model of the fluid–structure interaction benchmark case, which is consistent with Turek & Hron (Reference Turek, Hron, Bungartz and Schäfer2006). The key parameters in the simulation are set as follows: fluid density 
${\rho _l} = 1000\;\textrm{kg}\;{\textrm{m}^{ - 3}}$, dynamic viscosity 
${\mu _l} = 1\;\textrm{Pa} \cdot \textrm{s}$, average inflow velocity 
$U = 2\;\textrm{m}\;{\textrm{s}^{ - 1}}$, cylinder diameter 
$D = 0.1\;\textrm{m}$, density of the flexible strip 
${\rho _s} = 1000\;\textrm{kg}\;{\textrm{m}^{ - 3}}$, Young's modulus 
$E = 5.6\;\textrm{Mpa}$ and Poisson's ratio = 0.4. Under these specified values of parameters, the flexible strip will exhibit a flow-induced vibration. Figures 18(b) and 18(c) shows the time variations of horizontal and vertical positions of strip's trailing edges, respectively. It is evident that the present results are in good agreement with those of Turek & Hron (Reference Turek, Hron, Bungartz and Schäfer2006).
Numerical verification of fluid-structure interaction benchmark case. (a) Schematic diagram of a rigid cylinder with a flexible strip in a flow; (b,c) time variations of (b) horizontal and (c) vertical positions of strip's trailing edges. The solid and dashed lines represent the results from the present simulations and Turk & Hron (2006), respectively. The key parameters in the simulation are set as follows: fluid density 
${\rho _l} = 1000\;\textrm{kg}\;{\textrm{m}^{ - 3}}$, fluid dynamic viscosity 
${\mu _l} = 1\;\textrm{Pa}\mathrm{\ \cdot }\textrm{s}$, average inflow velocity 
$U = 2\;\textrm{m}\;{\textrm{s}^{ - 1}}$, cylinder diameter 
$D = 0.1\;\textrm{m}$, density of the flexible strip 
${\rho _s} = 1000\;\textrm{kg}\;{\textrm{m}^{ - 3}}$, the Young's modulus 
$E = 5.6\;\textrm{Mpa}$ and Poisson's ratio = 0.4.
