Abstract
Flying debris is generated in several situations: when a roof is exposed to a storm, when ice accretes on rotating wind turbines, or during inflight aircraft deicing. Four dimensionless parameters play a role in the motion of flying debris. The goal of the present paper is to investigate the relative importance of four dimensionless parameters: the Reynolds number, the Froude number, the Tachikawa number, and the mass moment of inertia parameters. Flying debris trajectories are computed with a fluid-solid interaction model formulated for an incompressible 2D laminar flow. The rigid moving solid effects are modelled in the Navier-Stokes equations using penalization. A VIC scheme is used to solve the flow equations. The aerodynamic forces and moments are used to compute the acceleration and the velocity of the solid. A database of 64 trajectories is built using a two-level full factorial design for the four factors. The dispersion of the plate position at a given horizontal position decreases with the Froude number. Moreover, the Tachikawa number has a significant effect on the median plate position.
1. Introduction
Flying debris, such as ice fragment, can cause serious damage to structures. One way to mitigate the risks associated with their impacts is to use numerical tools to predict their trajectories. However, ice fragment trajectories are stochastic, as many parameters affect the flying path and the flow [1]. In particular, trajectories are sensitive to the size, shape, location, and initial velocity of ice fragments [2]. Consequently, probability distribution maps are built to show the likely path of such fragments, which can lead to the development of risk mitigation solutions. A large number of trajectories are needed to obtain statistically relevant maps in 3D. Mathematical models enable the simulation of multiple trajectories for various initial conditions, which then allows the computation of probability maps around a geometry [3, 4]. The motion of the ice piece is driven mostly by aerodynamic forces. Research works covering wind engineering have established that three dimensionless parameters play a role in the motion of flying debris, apart from the Reynolds number that may play a role in aerodynamic forces [5–7].
For low fidelity models, motion equations are solved with a Runge-Kutta method in a Lagrangian frame of reference attached to the debris. The aerodynamic forces are computed using the static lift, drag, and moment coefficient correlations for the studied debris shape, neglecting the shape effects on the flow field. The two most studied shapes are the sphere [8] and the rectangular flat plate [6]. For the flat plate, the coefficient correlations depend only on the angle of attack, and additional aerodynamic terms must be added to take into account debris rotation [9]. For the sphere, the coefficient correlations depend only on the Reynolds number.
High fidelity models involve a tight coupling of the flow field to the flying debris motion. These models use a computational fluid dynamic (CFD) to solve the Navier-Stokes equations around the moving debris in an Eulerian frame of reference [7, 10], and the forces and moments are computed from the flow solution. The computational costs of such models are higher than those for models uncoupled from the flow field, especially as the flow Reynolds number increases. However, high fidelity models allow the computation of the trajectory of any kind of shape and do not require experimental databases or correlations for aerodynamic forces and moments. Also, high fidelity models take into account Reynolds number effects and lift and drag fluctuations (10% to 20%) caused by vortex shedding [11, 12] in the case of flow separation around debris. Both models are complementary and need research developments to better assess the risk of impacts.
The goal of the present paper is to investigate the relative importance of four dimensionless parameters: the Reynolds number, the Froude number, the Tachikawa number, and the mass moment of inertia parameter. A secondary objective of the paper is to study the effect of launching time on the trajectories. The focus is on the start of the motion, when aerodynamic forces on the ice fragment are high. A trajectory database is built for a 2D rectangular plate at Reynold numbers around 7000, but with other parameter values typical of inflight aircraft deicing. The plate makes an initial angle with the flow field such that flow separation occurs.
The CFD solver uses an immersed boundary method (IBM) to model the flow around the plate geometry [17]. A Cartesian grid is used for the computational domain. The plate boundaries are modelled by adding continuous forcing terms directly to the flow equations using the penalization method [18–20].
The paper is organized as follows. First, the penalized Navier-Stokes equations are presented in Section 2, together with the numerical method based on a Vortex-In-Cell (VIC) scheme. The model used to compute the forces and moments on the plate, along with the appropriate governing equations for fluid-solid interactions, is also described. In Section 3, test cases are presented to validate computations of forces and moments against literature results. In Section 4, a database is created to study the influences of the four parameters and initial angle of attack on the trajectories of a 2D rectangular plate.
2. Mathematical Model and Numerical Method
The fluid-solid interaction model consists of a set of equations for the fluid flow and another set for the solid motion. The effects of the flow on the solid motion are imposed through aerodynamic forces and moments, while those of the solid on the flow are imposed by penalization.
2.1. Fluid Motion
The incompressible laminar Navier-Stokes equations model the flow into the computational domain . A rigid solid , moving at velocity in the domain , is defined with a level set function . In this work, is the signed distance function to , negative inside the solid and positive outside. The boundaries of are located where the sign of the level set function changes from positive to negative. Figure 1 illustrates the signed distance function around the plate. A Cartesian grid is used to discretize the domain .
(a)
(b)
The velocity field is calculated throughout the computational domain, around and inside the solid, as illustrated in Figure 2. The penalization technique is used to impose the solid velocity vector in the area of the computational domain identified by the negative [18].
For a very large penalization parameter, , and the Heavyside function, , that is, if and otherwise, the penalized Navier-Stokes equations in a vorticity formulation are [10]where is the velocity vector, is the kinematic viscosity, is the fluid density, and is the vorticity vector.
The incompressible flow has a divergence-free velocity field. To impose boundary conditions and penalization inside the body, the velocity can be recovered from the vector potential :
The vector potential is a 3D extension of the so-called scalar stream function . This potential vector is imposed to be solenoidal; that is, . By definition, the stream function field is related to by
2.2. Aerodynamic Forces and Moments
The penalization term in (1), , forces the fluid velocity inside the rigid body to be nearly equal to . The local acceleration imposed by the penalization is computed by
The penalization term can never be zero, because an acceleration is always needed to counteract the pressure and viscous forces from the Navier-Stokes equation. Thus,and a large penalization parameter is needed to have . The sensitivity of the solution to the choice of the penalization parameter is discussed in Section 2.4.
The aerodynamic forces and moment are obtained by integration over the computational solid domain :In 2D flows, the moment becomes a scalar.
2.3. Solid Motion
For a fluid-solid interaction the force and moment induce the solid motion:where is the mass of the solid , its moment of inertia, its density, is the characteristic length, the angular velocity, and the gravity vector.
Following [5–7], the dimensionless variables , , and are introduced, being the freestream velocity. In the resulting dimensionless equations of motion, three parameters have significant effects on the solid motion: the Tachikawa number (ratio of the aerodynamic force to the gravity force), the Froude number Fr (or dimensionless plate length ), and the dimensionless mass moment of inertia parameter of the solid:where is the reference area of the solid.
For a flat plate of thickness in 2D motion, , and the mass moment of inertia isThus, and . The dimensionless mass moment of inertia parameter depends solely on the aspect ratio of the plate .
The motion equations in dimensionless form arewith
For the particular case of a solid in motion in gas, most often . Therefore, the solid motion depends on three parameters, separately from the dimensionless forces and moments. Note that, contrary to previous works [5–7], the freestream velocity is used for dimensionless forces and moments instead of the relative velocity.
2.4. Fluid Structure Interaction with Penalized VIC Scheme
The domain is meshed using a uniform fixed Cartesian grid. Equation (1) is solved using a classical VIC scheme [21], detailed in [10]. The grid values of , , and are evaluated at time . Each time step is solved using three substeps as follows:(1)Advection:(2)Stretching and diffusion:(3)Penalization term:
In the substep , the vortex particles located at grid point are displaced with a fourth-order Runge-Kutta time-stepping scheme and remeshed on the grid. In substep , the Laplacian is evaluated on the grid with a second-order accurate standard five-point stencil. Then, the velocity is computed using ( and ). The stream function is computed by solving the linear Poisson equation on the Cartesian grid, using the fast Poisson solver routine of the Intel® Math Kernel Library.
To assess the force acting on the solid, the local acceleration is computed as follows:The aerodynamic forces and moment are then defined by
A first-order explicit Euler scheme is applied for the time integration to update the solid linear and angular velocities. The solid velocities are then used to displace the solid boundary during the time interval . From the displaced solid position, is computed prior to the next time step. The solid boundary is defined by line segments and nodes. For each mesh node, the shortest distance between the node and the displaced line segments becomes the new value. A node is inside the geometry, negative , if a vertical line drawn from the node toward the upper limit of the computational domain crosses an uneven times the solid boundary line segments. The level set function is recalculated at each time step, thus avoiding the possible solid deformations associated with the advection of . In 2D, the computational cost of the recalculation is small compared to the cost of the fluid solution.
The penalization term depends on the values. The case of a 2D cylinder in a square cavity, falling under gravity on a flat plane, is used to study the sensitivity of the solution to the choice of these values. This test case has been presented in detail and used for validation in a previous article [22]. The dimension of the cavity is . The viscosity is 0.01. The density inside and outside the cylinder is, respectively, 1.5 and 1. The cylinder has a radius of 0.125 and is initially located at the point . It accelerates under gravity, set to . Figures 3(a) and 3(b) show the evolution of the cylinder vertical velocity with time. The parameter values range from to with a tenfold increase between values. For values of , a close-up view on velocity values near is provided on Figure 3(b) to see the difference. There is of difference between velocities computed with and . All the calculations in this paper are done with .
(a)
(b)
3. Results and Discussion
3.1. Lescape’s Fluid-Solid Model Validation and Verification
The computations of forces, using the code Lescape (LEvel Set CArtesian PEnalization), on a static body have already been validated in a previous work [23]. In this section, computations of forces and moments for solid motions are validated using an oscillating cylinder and a flapping airfoil.
3.1.1. Inline Oscillating Cylinder in a Fluid at Rest
The interaction of an oscillating circular cylinder with a fluid at rest is a problem well-documented in the literature [13, 14, 24–27]. The two key parameters when dealing with such flows are the Reynolds number and the Keulegan-Carpenter number, , where is the maximum velocity of the cylinder, is the diameter of the cylinder, and is the characteristic frequency of the oscillation. The set-up for this computation corresponds to the LDA experimental and numerical simulations performed by [13].
The motion of the cylinder is described by a simple harmonic oscillation,where is the amplitude of the oscillation. The parameters were set to and . The computational domain is and in the and directions, with the cylinder initially located at the center of the domain. The grid contains cells in the cylinder. Simulations are started in a static fluid field, and the cylinder moves according to (17) until periodic vortex shedding is established (see Figure 4 for an overview of the vortex shedding). Quantitative comparisons with the experimental and numerical data from [13] are considered in Figure 5, which shows the computed velocity profiles at four locations and three different phase angles. The numerical results from [13] were obtained by solving the Navier-Stokes equations for the velocity and pressure with a 2D finite volume method imposing an oscillatory flow around a cylinder at rest. The maximum difference between the numerical results is below for both and velocity, which indicates a good agreement, taking into account the differences in the mathematical models and the numerical methods. In the paper by [13], the comparison between the numerical and the experimental data is considered to be very satisfactory, even if the experimental error is not given. It is expected that the velocity distributions should be symmetric for and antisymmetric for . The dimensional drag force divided by the cylinder diameter over one period is also in good agreement with literature results (see Figure 6).
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(e)
(f)
3.1.2. Validation of Force Calculations through an Oscillating Airfoil
An oscillating NACA airfoil simultaneously experiencing pitching and heaving motions is modelled. The pitching axis is located along the airfoil chord at position . The airfoil motion, described by [15], is defined as follows: is the pitching amplitude, and the phase difference is set to 90°. The heaving velocity is then given by
Based on the imposed motion and the upstream flow conditions, the airfoil experiences an effective angle of attack and an effective upstream velocity defined bywhere the freestream velocity far upstream of the oscillating airfoil is .
To validate our simulations, a regime corresponding to the parameters , , , , and was computed. A view of the motion is sketched in Figure 7.
The computational domain is ; the boundary conditions are set to Dirichlet on the left, top, and bottom, and a Neumann boundary condition is set on the right. The numerical results, obtained using and , are then compared to the force predictions presented by [15, 16].
The force coefficients , , and the moment coefficients , are presented in Figure 8. The drag, lift, and moment coefficients are defined, respectively, aswhere the -axis is aligned with and the -axis is perpendicular to the far field velocity vector. Considering the Lescape’s drag coefficient (Figure 8(a)), the time evolution of the drag coefficient coincides with the results of Kinsey et al. and Campobasso et al. The numerical error for is below 3% for the maximum and minimum amplitudes. Lescape slightly overpredicts (by around 10%) during the decreasing phase of the drag. The evolution of the lift (Figure 8(b)) coincides with literature results. Lescape is in phase advance compared to Kinsey et al. and Campobasso et al., and the maximum time lag is below 5% just before the mid cycle . The error in amplitude is globally below 5% during the cycle, with the maximum error occurring at the end of the cycle, and being about 10%. The pitching moment (Figure 8(c)) is the one with the least agreement with the literature results, with a maximal error of 20% at and , the remaining error being below 10%. It must be noticed that, for pitching moments, Lescape shows a tendency to underestimate maximum (positive values) and overestimate minimum (negative values). The evolution in time of follows the literature results; it increases when literature results increase and decreases in accordance with published results (meaning that, in the situation of fluid-solid interaction, the Lescape code will predict the rotation angle and its sign to an acceptable extent). For , a small time lag at a maximum of 3% is observed. , , and numerical results are satisfactory in the perspective of trajectory computations.
(a)
(b)
(c)
3.2. Plate Trajectories Database
For the parametric study, a sharp corner plate is used, as in Figure 2. During these simulations, the rotation axis of the rectangular plate coincides with its barycenter. The effects of , , , and on flat plate trajectories are studied in order to identify the important parameters at the start of the trajectory computation. A full factorial design with four factors and two levels [28] is used to build the test matrix.
The simulations are performed using and a domain of , with the plate’s center initially located at , and with starting angles of . For each factor, minimal and maximal values are selected around a base value: , , , and . Except for the Reynolds number, these are the typical values for shed ice pieces around an aircraft in a deicing situation [29]. For , two values of the aspect ratio are used, namely, and . This is only a variation, but this is a significant thickness change. Table 1 gives the 16 run factors. The 16 runs are executed for two plate starting angles and two release times. Thus, the database consists of trajectories.
The starting angle creates a separated flow around the plate. Vortices are shed periodically at low frequencies. The flow is first solved around a static plate, which goes on until the periodic flow is well established. The plate is released into the flow either at time s or at time s. The dimensionless position of the plate when it crosses the line will be studied.
The computed trajectories are illustrated in Figures 9 and 10. To help visualize the individual trajectories, the results are separated into four. For a positive plate starting angle, these results are presented in Figure 9, with a further separation between the two shedding times. Similarly, the results for a negative plate starting angle are shown in Figure 10. The positions and are plotted against each other. Since the initial angle of the plate is not null, vortex shedding is created behind the plate. The influence of vortex shedding is not negligible for a particular trajectory. For different starting times, the plate positions at the end of the trajectory calculations are different, as comparison between (a) and (b) shows. However, the results at different launching time will be grouped together to study the effects of other parameters. It can be seen that, in general, the plate accelerates faster in the direction than in the direction. The gravity effects curve the trajectory. Here, the plate motion is mostly driven by aerodynamic forces and the starting angle has a dominant effect.
(a)
(b)
(a)
(b)
Figure 11 presents an example of the trajectories obtained for a plate with an initial positive angle. The dimensionless positions of the center of gravity for the run parameters with an initial positive angle of attacks and a release at s is shown. The computed average trajectory of the runs is also plotted, together with the upper and lower limits for a confidence interval. There is a probability that the mean position of the plate gravity center at a given position lies between the upper and lower limits. The confidence interval is computed using the sample variance and a Student probability distribution [30]. As can be seen, many trajectories are clearly outside the confidence interval. The mean trajectory gives a representative idea of the plate motion, but all the trajectory must be plotted to take into account the stochastic nature of the motion.
The plate geometry, corresponding to run , is drawn at a constant time interval in Figure 12. For this run, the plate is rotating only slightly. The horizontal acceleration of the plate is visible as the spacing between its locations increases with the position. The mean trajectory of the center of gravity is also plotted with the corresponding confidence interval. For run , the average trajectory is representative of the plate motion.
The trajectory footprint probability maps used for the 3D trajectories analysis by [3] become histograms for 2D trajectories. At , a line is drawn in the direction and split into intervals. The heights of the histogram bar are computed by dividing the count of plate passage into an interval by , which represents the total number of runs [28]. The histograms associated with the two initial plate angular orientations are shown in Figure 13. In this figure, the position is on the horizontal axis, and two histograms are plotted, one for each plate orientation. The length of each bar gives the probability of passage of the plate gravity center in an interval of . The initial plate orientation has a clear effect on the histograms, with the positive angle giving positive values, and the negative angle giving negative values. This is expected since an initial positive angle creates a positive lift coefficient, and a negative angle creates a negative lift coefficient. The two histograms are placed almost symmetrically around the line, indicating that gravity is not the dominant force of the problem at the beginning of the motion.
The initial plate orientation is the parameter that has the most effect on the trajectories. To study the effects of the other parameters, the database is split into two sets of runs: set one for the starting positive angle, and set two for the starting negative angle. Box plots are used to compare parameter effects [28]. Figure 14 shows the relation between the box plot and the probability histogram for the results of the flat plate with a starting negative angle. The central mark is the median of the runs. The edges of the box are , the 25th percentile, and , the 75th percentile. The whiskers extend to and , or the most extreme data points not considered outliers. The points for which are smaller than or greater than are considered outliers. One outlier is visible in Figure 14.
For the trajectories computed with a starting angle of (positive orientation) (Figure 15), the box plots show that the parameter has a significant effect on the median value. Also, the median value for is lower than the 25th percentile for the distribution, which confirms that gravity forces are more important for lower values. For the two other dimensionless parameters, the median does not change significantly between minimum and maximum values. For example, for , the median is between the edges of the box for distribution, while that for is between the edges of the box for distribution. The result dispersion is also changed by values and increases with increasing values. The whiskers extend from to for and from to for . According to the box plot, lower values will make the plate fly higher and will reduce the dispersion. As is only present in (12), the result dispersion is caused by the plate rotation. We should recall that is the dimensionless length of the plate; . It means that longer dimensionless length plates fly more easily and in a more predictable fashion.
The box plots for a plate at a starting angle of (negative orientation) are shown in Figure 16. The number once again has the most significant effects on the plate location at . The standard deviation increases noticeably with an increase of the number. Decreasing the decreases the median position of the plate, as was the case for the plate with a starting positive angle. Decreasing and has almost no effects on the final center of gravity position.
This database parametric study focuses on the start of the trajectory. To give an order of magnitude, in typical ice trajectory experiments [29], the ice piece travels a dimensionless horizontal distance of , three times the horizontal distance in this database. The major limitation of this study is the 2D computational domain that precludes any comparisons against 3D experimental results. The computed 2D results should be equivalent to 3D results with a plate of infinite span and different from the results obtained with a plate of short span extent. Short span plates are typical of most engineering applications and experiments. Another limitation is the small number of runs used for the statistical analysis. Although runs are used to build the database, the effect of the plate starting angle on the final position is strong, and runs must be separated into two databases. Trajectory depends on launch time and multiple launch times should be used to populate a larger database. Nevertheless, the database size is comparable to the database size used by [7] for parametric study.
The numerical errors on aerodynamic coefficients observed during the oscillating airfoil validation could also have an impact on the database results. The aerodynamic moment coefficient oscillates typically between −0.6 and 0.6 during one trajectory calculation. The number of oscillation cycles is usually above three. Assuming the airfoil validation errors apply for the moving plate, the moment coefficient maximum values are underestimated by and the minimum values overestimated by . For one cycle, the numerical errors have probably neutral effects. Also, it is worth mentioning that, for Lagrangian methods, often used for trajectory calculations, of errors on static aerodynamic coefficients are usually assumed, together with an unknown uncertainty on the damping function used to model the dynamic caused by rotational velocity [5]. An additional limitation specific to aeronautic is due to the values of the Re, which are one level of magnitude lower than those for aerodynamic applications. This may not be a critical limitation as the static experimental lift, drag, and momentum correlation for a rectangular plate do not show any dependence on Re [31]. However, Lagrangian trajectory calculations use corrections to the lift, drag, and momentum that depend on Reynolds. Actual CFD results show some number effects.
4. Conclusions
A database of trajectories for a 2D plate released in a moderate Reynolds number flow has been built using CFD. Apart from the plate starting angles and launch times, four parameters influence the trajectories: the Reynolds number, the Froude number, the Tachikawa number, and the dimensionless mass moment of inertia parameter. Both the Froude and the Tachikawa numbers have significant effects on the trajectories. It appears that increasing the Froude number increases the dispersion of the trajectories in a vertical plan. Lowering the Tachikawa numbers decreases the median value of the dimensionless vertical position. Launch times change significantly the trajectory, as indicated by scatter plots of the dimensionless vertical position as a function of horizontal position, due to the unsteady separated flow around the static plate.
The CFD tool used to build the database solves the penalized Navier-Stokes equations. The numerical method for the flow is based on a VIC scheme. The aerodynamic forces and moments are computed by integrating the penalty term and are used to compute the velocity from the solid motion equation. The velocity is imposed on the solid area of the flow, identified by a level set function. The forces and moments are validated against literature results for an inline oscillating cylinder and a flapping airfoil.
The mean trajectory with a confidence interval gives a representative idea of the plate motion as long as the two starting angles are processed separately. The histograms of passage probabilities at a vertical position show the strong effect of the starting plate angle. The box plot representation allows a comparison of the effects of the four parameters on a single graph.
In the future, it could be interesting to carry out the same study with a higher Reynolds number more representative of the aerodynamic flow in deicing situations. Also, the same kind of study could be conducted in 3D to see if the same conclusions stand. Progress in high performance computing is however sorely needed if large enough databases are to be built.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
(i) Numerical experiments presented in this paper were carried out using MCIA (Mésocentre de Calcul Intensif Aquitain) facilities, the PLaFRIM experimental testbed, and the supercomputer Mammouth Parallèle 2 from Sherbrooke University, managed by Calcul Québec and Compute Canada. The operation of this supercomputer is funded by the Canada Foundation for Innovation (CFI), the ministère de l’economie, de la science et de l’innovation du Québec (MESI), and the Fonds de recherche du Québec-nature et technologies (FRQ-NT). PLaFRIM cluster is being developed under the Inria PlaFRIM development action with support from LABRI and IMB and other entities: Conseil Régional d Aquitaine, FeDER, Université de Bordeaux, and CNRS; see https://www.plafrim.fr/fr/accueil/. (ii) This study has been carried out with financial support from the French State, managed by the French National Research Agency (ANR) in the frame of the “Investments for the future” Programme IdEx Bordeaux-CPU (ANR-10-IDEX-03-02).