Predictive modeling of drop impact force on concave targets

The impact of a raindrop on skin is perceptible, in part not only because of its temperature but also because it imparts a force large enough to be sensed. Impacting drops are ubiquitous and their impact force may be their single most studied dynamic quantity.^1–17 Impact force has been measured with piezoelectrics, accelerometers, and cantilever beams. Despite the number of studies characterizing drop impact force, the surfaces used in studies are flat and horizontal, with few exceptions.^1,4 In a recent review article on liquid drop impact force, Cheng et al.² identify inclined surfaces as an area for future work. In this experimental and computational study, we measure the impact force and film the impacts of falling drops into spherical, concave cups with radii of curvature on the same order as the drop. Drops impact with a Reynolds number $Re=ρDU/μ=2,975−12,800$⁠, where density ρ = 1 g/mL, impact velocity $U≈2gh$⁠, D is the drop diameter, and the assumed drop dynamic viscosity μ = 1 cP. Drops are released from the rest at a height h and fall under gravity g = 9.81 m/s². Our cups are inspired by splash cup plants of Chrysosplenium and Mazus that employ a conical shape to amplify the speed of raindrops to disperse seeds.^18,19 Drops which impact the cups off center attain the greatest dispersal distances at speeds up to five times faster than the oncoming drop. In contrast, our impacting drops remain in their spherical cups. Beyond the impact of drops on concave plants, the force imparted by drops is useful to characterize for sensors, erosion,^20,21 and water capture.^22,23

From a high-level perspective, the force of impact is a direct result of how, and over what timescale, the momentum of a drop is arrested, redirected, and dissipated.^2,8 After first contact on a flat surface, impact force rises quickly and reaches a maximum when the drop equator meets the substrate.² After its peak, impact force decays more gradually than it rose, resulting in an asymmetric impact curve. How temporal impact curves manifest for substrates with curvature on the order of drop diameter is unknown. Impact force arises from momentum redirection and drop deformation, to include breakup, and is governed by drop size and speed, surface tension, viscosity, surface orientation, texture, compliance,^24,25 and advancing contact angle. Not only is this list perhaps incomplete but also stochastic effects further complicate repeatability in instances where exact impact location and drop shape are critical. Many variables governing impact force exist on a continuum and are difficult to precisely control, and, thus, determining a reliable empirical relation through full factorial experimentation is a campaign of intense devotion despite the straightforward nature of trials. Moreover, our experimental system provides the added challenge of limited optical access to the impact zone. Drop dynamics can be measured prior to impact, but the impact itself is obscured from view.

Modern algorithmic tools provide an opportunity to reveal variable relations and dominance in places where tedious experimental campaigns may fail. With a sufficient amount of experimental data, machine learning tools are brilliant in prediction of experimental outcomes, in this case, impact force F, for a vast array of input variables that have not been tested. While machine learning tools have been applied to fluid mechanical systems,^26–28 their application to laboratory experiments in fluid dynamics is relatively limited. A recent example of using machine learning to describe drop experiment data is that of drop impact onto cantilevered fibers.²⁹ Algorithm predictions confirmed scaling relations that suggest the maximum deflection of a fiber impacted on its tip by a falling drop is nearly independent of drop momentum. In another study, machine learning was used to predict the inertial force required to eject drop from flexible, millimetric cantilevers.³⁰ The cantilever drop system exhibits complex dynamics because drops slosh, relocate, and damp vibration.^31,32

We present our experimental and computational approaches in Secs. II and III, respectively. A scaling relation for impact force as a function of system variables is developed in Sec. IV for use in interpreting algorithmic outputs. We discuss outputs and predictions made by our chosen algorithms in Sec. V and conclude our work in Sec. VI.

Our experimental apparatus used to measure impact force is shown in Fig. 1. Nozzles are positioned at 141 unique heights h ranging from 50 to 494 mm along a 165-cm tall 80/20 1530 extrusion tower affixed to a 3′ × 6′ optical slab (Nexus, Thor Labs, Newton, NJ). Drop heights were selected such that we get a reasonable range in impact velocity but drops reliably strike within cups. Water is pumped to the nozzle by a syringe pump (New Era NE-1010, Farmingdale, NY). Five straight cut needles with outer diameters of 0.5–1.5 mm are used to form drops 2.3–4.4 mm in diameter. Our drop heights do not allow drops to reach their terminal velocity. Drops impact polyacitide (PLA) ‘cups’ manufactured by a three-dimensional printer (Flashforge Creator Pro, City of Industry, CA). A cross-sectional diagram of our cups is shown in Fig. 2(b). The tallest cup post is 20 mm high and has a cup radius of 5 mm. Cups with larger radii r have their bottom positioned 5 mm below the top of a 20 mm high virtual cup. Sketches for each cup showing all pertinent dimensions are given in Fig. S1 of the supplementary material. Cups are adhered to a force gauge (MARK-10 Series 7 M7–012, Copiague, NY) with cyanoacrylate. The force gauge records force readings at 14 000 Hz. A example of the force gauge's digital output can be found in Fig. S2 of the supplementary material. The cup rigidly fixed to the gauge [Fig. 2(a)] behaves as a highly underdamped oscillator and we, thus, take the maximum force value of the gauge output as the measure of impact force F. Impacts are illuminated by LED lights (GS Vitec MultiLED LT, Bad Soden-Salmünster, Germany) and filmed with a high speed camera (Photron AX-100, Tokyo, Japan) at 6400 fps. A mirror set at $45°$ to the lens provides an orthogonal view for quantifying impact eccentricity. Cups are dried with an absorbent cloth between trials.

Videos are analyzed with MATLAB (Mathworks, Natick, MA). From our videos, we measure drop diameter D and the lateral distance from the center of the cup to the center of the drop a. The measurement of a requires both impact views. Offset is limited such that the entirety of the drop remains in the cup, $a<5 mm−D/2$ for r = 5 mm and $a<6 mm−D/2$ for r = 10, 15, and 20 mm. We non-dimensionalize offset as $ε=a/r$ and report the range of values in Table I but do not use this ratio in our algorithm. Two wetting conditions are trialed. The first uses untreated cups, naturally hydrophilic (W = 0), with static and advancing contact angles of $77.5±1.1°$ (N = 3) and $94.7±2.5°$ (N = 3), respectively. For the second, we spray cups with NeverWet to create a hydrophobic condition (W = 1), with static and advancing contact angles of $147.7±5.8°$ (N = 3) and $151.3±3.2°$ (N = 3), respectively. Our experimental parameters are listed in Table I. Our experimental campaign aimed to capture as many trials as time would allow with roughly equal numbers of trials for each cup. Initially, the number of data points for each of four cups ranged from 108 to 132, totaling 476. Data were filtered for eccentricities too large to contain the entire drop, reducing unique data points to 387 with a range of 86–110 per cup. No replicates were performed for two reasons. The first is that true replicates are impossible, owing to deviations in drop size and impact location. The second is that machine learning algorithms are best trained with the largest combinations of unique inputs.³³ 

Experimental parameters listed in Table I are supplied to machine learning algorithms to build a regression model for impact force. Since the output of the prediction model is a continuous variable, F in our case, the prediction problem in this work is a regression problem. A set of base learners is, thus, formed of regressors: Random Forest Regressor (RFR), k-Nearest Neighbor (KNN), the Multi-Layer Perceptron (MLP), and the Gradient Boosting Machine (GBM). These four are well-known learners and frequently used in regression problems. We implement the ensemble using the Scikit-learn library, a free machine learning library for Python programmers.³⁴ 

The performance of any regression problem critically relies upon the selection of the algorithm(s).^29,30,35 In general, the selection of the base learner set³⁵ is guided by data distribution and a bit of trial-and-error. RF and GBM excel with datasets having a high degree of nonlinearity and complex relationships between variables and are truly ensembles within themselves. RFR uses bagging (bootstrap aggregatbootstrap aggregating)³⁶ and GBM using boosting. The key difference between bagging and boosting is that boosting trains its internal base learners to build each decision tree sequentially rather than independently, calculating the average of tree outputs at the end. Each internal base learner within a GBM is dependent on the previously trained base learner. The MLP is a multilayer neural network with full connectivity between the layers, usually trained with backpropagation. Training an MLP requires adjusting model parameters to minimize prediction error. The use of KNN is motivated by small data size and small number of features. For more details on the aforementioned learners, we direct the reader to Orkweha et al.²⁹ 

Ensemble learning combines the outputs of base learners to achieve a more accurate prediction than constituent base learners by lowering error and cross-fitting,³⁷ when individual base learners are selected properly.³⁸ The four aforementioned base learners are combined by a bagging regressor, a combination of bootstrapping and aggregating. The implementation of both bagging and boosting techniques improves the accuracy of our predictive model by providing varying sampling methods. Our ensemble learning implementation method is schematized in Fig. 3. The prediction of each base learner is multiplied by an assigned weight before a summation forms the final prediction.³⁹ 

Assigned weights are determined by the non-negative least squares (NNLS) method such that the final prediction has the lowest possible error given each base learner's prediction accuracy. Regardless of the distribution, the weights always sum to 1. NNLS was first introduced by Lawson and Hanson,⁴⁰ as

where S is the training set such that $S∈Rm×n$ with m features and n number of observations, $α$ is a weight vector, and β is the unknown sample. The NNLS method is implemented after training each base learner with tenfold cross-validation with shuffle which produces ten sets of predictions from each member. k-fold cross-validation is a technique in which the algorithm is trained k times with a fraction $1/k$ of training examples omitted for testing.⁴¹ The weight vector $α$ is distributed among the constituent learning algorithms and predicted values Pr from each base learner are multiplied by a their respective α value and aggregating products to generate a final prediction

Between the moment the drop makes initial contact with the cup and the moment of greatest drop spreading, a maximal impact force F is measured. Due to limitations in viewing angles in a concave, opaque cup, the progression of drop shape and its connection² with temporal force is difficult to obtain optically. The force of an impacting drop arises from the force of impact pressure $Fi$ and the reaction force to deform the drop $Fs$⁠, as shown previously⁴ and we expect both to be of importance to current impacts. The impulse of drop impact, thus, takes the form:

where $U′$ is the vertical velocity of the fluid mass post-collision, to include the velocity of any mass redirected upward by the curved cup, $m=ρπD3/6$ is drop mass, and $Δt∼D/U$ is the impact time. Quantification of $U′$ for every impact scenario is intractable, and so we write $F≈Fi+Fs+mg$⁠.

For horizontal flat surfaces, the incoming drop's velocity contributes to an impact pressure $∼ρU2/2$⁠. Since pressure acts normal to a surface, those on an incline require the normal component of drop velocity $Un$ to scale impact pressure. For the sake of scaling, we assume the reaction of a eccentric impacts in our cups produces forces equivalent to impacting a flat surface with the same inclination as the cup at the impact location. Normal velocity at the point of impact can then be written $Un=U(1−ε2)$⁠. We now define a Weber number that assumes drops release from rest and incorporate eccentricity as

where the surface tension of the liquid is σ. The force arising from momentum is concisely formulated by Soto et al.¹² and which our modification for normal velcoity is written as follows:

Deformation of the drop increases surface energy, and we expect the related force $Fs∼σχ$⁠, where χ is the maximal equivalent diameter of a spreading drop as it impacts. Since we cannot visualize the extent of spreading of our drops, we rely on a relation by Scheller and Bousfield⁴² for impacts on flat, dry surfaces that do not recoil or rebound, $χ/D=0.61(We/Oh)0.166$⁠, where $We=ρU2D/σ$ and the Ohnesorge number $Oh=μ/(ρσD)1/2$⁠. We believe the Scheller and Bousfield relation holds for our dry, curved cups with a modification, but likely with a different scaling constant. We modify their relation by replacing the classical Weber number with $Wen$⁠, such that $χ/D∼(Wen/Oh)0.166$⁠. Now $Fs$ can be written as follows:

Combining Eqs. (3) to (6), we can write an expression for dimensionless force,

which permits $F∼mg$ when h = 0. Unpacking dimensionless groups in Eq. (7) takes the following form:

For the sake of curve-fitting and because Eq. (7) has implicit scaling coefficients in the first two terms on the right-hand side, we may rewrite Eq. (7) in terms of only the parameters which are variable in our experiment

We use Eq. (9) with experiential data and model predictions in Sec. V B.

We impact concave cups with falling water drops to record the maximum impact force recorded by a high-frequency force gauge. Inputs to each trial, h, D, a, r, and W, and output F are passed to the ensemble algorithm described in Sec. III for training and testing. The choice of using h as a parameter over the arguably equivalent U is a matter of measurement error across hundreds of trials. In short, determination of U is prone to more measurement error. Before merging base learners into an ensemble, we resolve the optimum parameters for each base learner by the grid search method. Optimum parameters are shown in Fig. S3 of the supplementary material. The predictive models are trained to predict impact force and employ tenfold cross validation to verify model performance.³⁰ The predicted vs observed impact force is plotted for each base learner in Fig. 4, with root-mean-square-error (RMSE) values printed in Table II. No individual base learner stands out as vastly superior. The ensemble improves RMSE by 15.7% compared to the best base learner, RFs.

The NNLS problem of Eq. (1) is solved to assign weight to the prediction of each base learner according to Eq. (2). Regression error characteristic (REC) curves are used to visually compare the performance of individual base learners against the ensemble in Fig. 5(a). RECs show the accuracy of a regression model by comparing it against the absolute deviation. Superior predictions come from models which reach 100% accuracy at a smaller value of error tolerance. Simply put, curves which bound more area under the curve (AUC) have less error. Predicted impact force is plotted against measured impact force for the ensemble in Fig. 5(b).

The trained algorithm assigns a parameter “importance score” to each input variable that represents the relative influence, with respect to one another, of each parameter on the recorded impact force. Importance scores are independent of base learners and are found by measuring the decline of model accuracy when shuffling a parameter's value.⁴³ Importance scores are shown graphically in Fig. 6. Our algorithm predicts that the height h of drop release, the driver behind impact velocity, is the single most dominant parameter contributing to impact force. The importance of h is surprising when compared to the importance score of D, which is nearly half. For flat surfaces, $F∼ρgD2h$⁠, and so we might expect the scores for D to be a more important contributor to force in our problem. However, and as we discuss more below, importance scores cannot be read blindly. Importance scores are produced by shuffling each variable and quantifying the error in the prediction vs the measured output, force in this case. In our experiments, we vary D by a factor of 1.9 and h by a factor of 9.9. We, thus, find no surprise that h appears to have such a greater influence on impact force. We would expect these two importance scores to be dramatically different with respect to one another if we modulated D by a factor of 10.

Of lesser but still significant importance is a, a score intimately tied to the magnitude of r. In the limit of very large cup radii, a will have negligible importance. The variables r and W have a small number of discrete values, four and two, respectively. Thus, when a parameter is shuffled, there is a significant chance a parameter value is effectively unchanged. In these cases, the model does not have error contributed by these shuffled parameters and the corresponding importance score is low. Thus, if r and W had a comparable number of discrete values as the others, we would expect their importance scores to be significantly higher. The importance of r and W is assessed by a different means below. Wettability W does not appear in Eq. (9), and the effect of wetting on temporal or maximum impact force has not been characterized in the literature to the authors' knowledge.² 

The nature of our experimental data, purposefully randomized for input into ensemble learning, does not permit the visual validation of Eq. (9). However, using all our experimental data, we run a global optimization scheme to find C₁ and C₂ values which permit Eq. (9) to provide a calculation of impact force with the least amount of error. The lowest RMSE between experimental F/mg and that predicted by Eq. (9) is 50.9 or 38% the mean value of experimental F/mg, with $C1=1.0$ and $C2=224.7$⁠. Here, C₁ is dimensionless and C₂ has units [mm^−0.083]. If we run the optimization scheme on Eq. (9) with $C2=0$⁠, the optimal $C1=1.9$ with $RMSE=58.4$⁠.

We can use Eq. (9) to fit algorithmic predictions in which we vary only one variable at a time. From Eq. (9), we expect F to increase with h, D, and r and decrease with a. Surface tension, which governs the force required to deform a drop, appears in the second term of Eqs. (7) and (8), $∼σ0.92$⁠. While we do not change surface tension in our experiments, we expect impact force to increase with surface tension as drops behave more solid-like.

We employ our trained algorithm to predict impact force in which we can isolate an independent and dependent variable (F). Such an exercise can be performed for any combination of input variables. For simplicity, all invariant inputs are set to their experimental averages, across all trials, reported in Table I. We plot predictions for F/mg across variable D, with h, a, and W = 0 fixed in Fig. 7(a), producing a curve for each experimental cup radius. F/mg generally decreases with increasing D, as expected from Eq. (9). Likewise, F/mg is predicted to increase with cup radius. With offset a fixed at $a¯=2.5$ mm, an increasing cup radius r decreases impact eccentricity ε, thereby increasing $Wen$ and F/mg. A further simplification to Eq. (9) may be done by lumping all constant terms for a particular curve in the predictions of Fig. 7(a),

Equation (10) is fit to prediction points in Fig. 7(a) and shown by smooth curves, with coefficients of determination (⁠$R2$⁠) printed in the legend. We disallow fitting coefficients to be negative and find that for all four curves in Fig. 7(a), C₃ is functionally zero, further indicating that inertial force cannot be considered in isolation, a departure from conventional approaches.^9,11,12,14 With C₄ artificially set to zero, best fit correlation coefficients drop to $0.63−0.69$⁠.

Fixing all inputs according to Table I except for h results in the prediction points shown in Fig. 7(b). Increasing h boosts the kinetic energy of the falling drop and $Wen$⁠, increasing F/mg as expected. The simplification of Eq. (9) for a variable h only takes the following form:

which is fit to data points in Fig. 7(b) and plotted as smooth curves with $R2$ values printed in the legend. As before, C₆ values are dominant, with $O(10−100)$⁠, over C₅ values with O(1). A true test of the dominance of surface tension and viscosity within the deformation term would be to include these as experimental variables and is an area for future work.

Teasing out the influence of a varying impact offset a from Eq. (9) is less tidy than D and h. An expression for F/mg, which is only a function of a, or equivalently r, like that done with Eqs. (10) and (11), creates an equation with six fitting parameters and five terms. We, thus, find fitting such an expression to prediction data of little value. Instead, we apply a second order Savitsky–Golay smoothing filter spanning one-third the prediction domain, as done in previous works.^{29,31,32,44–48} Such a filter reduces prediction noise to provide a better visual representation of the prediction trend. The choice of span is somewhat arbitrary but chosen to capture both local and global curvature. We plot raw and smoothed data for F/mg vs a in Fig. 7(c). For the prediction set associated with each cup radius, a hump in the F/mg curve near a = 2 mm appears. Following the hump, at greater values of a, the force curve flattens. The cause for such behavior is unknown and an area for future work. Perhaps near this offset value, the cup is more effectively capturing both the portion of the normal and tangential (shear) forces imposed by the spreading drop. If so, one would expect the hump for r = 20 mm to be shifted right more than the other cup radii. An alternation in cup shape, from spherical to conical, is expected to dramatically change influence of impact offset. In a conical cup, a drop which does not spread to the cone tip should produce nearly a constant F/mg across a.

The final use of our algorithm is to graphically show the influence of wettability. Contact angle does not appear in Eq. (7) and has been found to have little effect over spreading for impact Reynolds numbers greater than 10 on flat surfaces.^42,49 Maximal spread is also independent of roughness.⁵⁰ It was previously unknown if this effect is minimal for concave targets. Since we use two discrete values of W, we are unable to continuously vary W in predictions to meaningful substance. Thus, we show how F/mg varies with r for W = 0 (hydrophilic) and W = 1 (hydrophobic) in Fig. 7(d), curves which are nearly the same shape as those in Fig. 7(c) rotated $180°$⁠. We predict that the hydrophobic cups experience a greater impact force for all r. Such a prediction is surprising in light of previous works.^42,49,50 We posit that hydrophobicity enables drop mobility and so drops spreading within hydrophobic cups will have a greater bias toward the cup center than a hydrophilic impact, perhaps rolling downhill. Future experiments employing transparent cups where such action can be observed alongside force outputs is a fecund area for future work.

We experimentally impact small millimetric, spherical cups of varying radii with falling drops and record impact force from a high-resolution force gauge. Impacts are filmed with a high-speed camera to record drop diameter and impact distance from the center of the cup. Impact parameters including force, drop height, diameter, impact offset from the cup center, cup radius, and wetting condition are used to train an ensemble machine learning algorithm consisting of four base learners. The ensemble achieves an RMSE of 2.85 mN. We formulate a scaling relation for impact force to compare with algorithmic predictions. Such a comparison reveals that the simple scaling used for impact force on flat surfaces does not well predict impact force in cups. Algorithm predictions show that flatter cups or those with a larger radius of curvature experience a greater force for a fixed experimental offset, and, thus, concave cups experience a lower impact force than flat surfaces. Hydrophobic cups experience a greater impact force than hydrophilic cups, a surprising result considering drop spreading behavior on flat surfaces. Greater experimental data would enhance predictive performance, and future directions of research should include more granularity in the experimental range of surface chemistry, cup radius, and the use of liquids with different density, viscosity, and surface tension. This exploratory study reveals that drop impact into highly curved surfaces is an area fertile for future investigations.

See the supplementary material for the three additional figures.

The authors would like to thank the National Science Foundation, U.S.A. (Nos. CNS-1852002 and CBET-1941341) for support and Gene P. S. Rible for contact angle measurement assistance.

The authors have no conflicts to disclose.

Andrew Keith Dickerson: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Funding acquisition (supporting); Methodology (lead); Project administration (lead); Supervision (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). MD Erfanul Alam: Data curation (equal); Formal analysis (equal); Validation (lead); Visualization (equal); Writing – review & editing (equal). Jacob Buckelew: Investigation (lead); Methodology (supporting); Writing – review & editing (supporting). Nicholas Boyum: Investigation (lead); Methodology (supporting). Damla Turgut: Funding acquisition (lead); Resources (lead); Supervision (equal); Writing – review & editing (equal).

Code, raw videos, and .xls files containing experimental measurements are available in perpetuity via Open Science Framework (https://osf.io/h7s98/).