Prediction of high-speed debris motion in the framework of time-fractional model: theory and validation

In this work, a newly proposed fractional derivative framework is used for the prediction of high-speed debris motion. The paper focuses on the mathematical formulation of the equation of motion, in which the damping term is generalised using the fractional derivative. The capacity of the proposed approach to predict the motion of debris is justified by the experimental results. Furthermore, the mathematical formulation has been verified by extensive parametric studies on spherical projectiles. The general conclusion is that the elaborated formulation is more reliable compared to the classical approach or, in other words, the fractional viscous damping term (proportional to the fractional velocity of debris) provides a better description of the complexity of the real drag force.


Introduction
Predicting of high-speed debris motion (HSDM) is crucial to determine safety zones in the case of a terrorist attack or the possible action of space debris. Here, we can state that there is a limited amount of research on the prediction of the HSDM. Most works focus on the debris that breaks off the solid obstacle due to the impact of the shock wave -e.g. a work on the statistical description of explosion-produced debris dispersion by van der Voort and Weerheijm [1]; paper by Mébarki et al. [2] on the probability of impact of structural fragments from explosions in industrial facilities; work by Wang et al. [3] in which authors focus on the debris from masonry walls; and the paper that describes the capabilities of the software called KG, which could be used to analyze the risk of explosion in reinforced concrete ammunition magazines [4]. In papers [5] and [6] (both papers by Price et al.), one can find a successful attempt to model the motion of debris using the Euler-Lagrange Point-Particle method. The model given in [5] and consequently developed in [6] is calibrated for spherical particles with the damping coefficient defined as in [7]. The main disadvantage of the mentioned model is that it is limited to spherical particles and requires a sophisticated solver with sufficient computational power. The authors of this paper, based on the experience developed from more than a hundred well-equipped field tests, recognised the need for a fast but flexible phenomenological model of the HSDM, which helps to estimate especially safety zones.
In the previous work of the authors of this article [8], we presented results for the prediction of the debris trajectory using the standard approach. The standard approach to modelling debris trajectories uses the classical form of the equation of motion. However, such a description requires an appropriate model for a damping factor, which may be difficult to obtain because of the irregular geometry of debris, i.e. the damping factor, hence the drag force, is dependent on the fluid flow around an analysed fragment that is sensitive to its shape. Otherwise, the results may be treated only as a rough approximation of a reality. Therefore, in this work, we use a generalised approach using a fractional derivative framework for the prediction of the HSDM. Such a concept allows for calibration of additional (fractional) parameters, which constitutes the new phenomenological model of equations of motion.
The applied fractional modelling (using fractional calculus c.f. seminal books by Miller et al. [9] or Kilbas et al. [10]) follows the trends in modern theoretical mechanics. Different theories which were developed within the fractional framework 1 3 46 Page 2 of 21 include so-called: (i) stress-fractional models (stress-FM) -especially useful for geomechanics -e.g. seminal work by Sumelka where stress-FM was introduced [11]; papers by Sun et al. where stress-FM was applied to model granular soil [12], clay [13] or geomaterials in complex stress state [14]; or work by Lu et al. [15] where stress-FM was devoted to soil modelling; (ii) space-fractional models (s-FCM) -for the description of a strong scale effect for nano-micro bodies -e.g. paper by Sumelka where s-FCM was defined [16]; the works by Patnaik et al. where s-FCM was applied for elastodynamic [17] and variable-order nonlocal elasticity [18]; paper by Pang et al. [19] where dvection-dispersion problem was analysed; or the work by Beda [20] where stability and bifurcation analysis of s-FCM was presented; (iii) time-fractional models (t-FM)-(so-called models with memory) where materials with strong rheological effects are considered -e.g. paper by Failla et al. [21] where t-FM for two-dimensional foundation model was discussed; the work by Di Paola et al. [22] where variableorder fractional viscoelasticity was defined; paper by Sumelka et al. [23] where damage evolution was cover by fractional approach; the work by Kukla et al. [24] where heat conduction in the framework of t-FM was considered; or investigation by Shariyat et al. [25] where viscoelastic annular rotating discs under nonuniform loads were analysed. Herein, the last mentioned modelling technique (time-fractional) will be applied, as the damping term (dependent on the velocity) plays a fundamental role in debris motion analyses, i.e. the description of trajectory, velocity, and acceleration. The benefits of the fractional approach (independently of technical details) lay in the fact that fractional operators are nonlocal, i.e. their definitions are not defined at a point as classical integer order ones but include an interval (surrounding of a point) which defines the range of its action or in other words its nonlocality. Therefore their applicability for modelling is higher compared to integer (local) derivative operators. Furthermore, the fundamental concept of "short memory" proposed in [26] and later developed in [27] (where the basics for fractional damage evolution (FDE) were presented), [28] (where the anisotropic FDE for human brain modelling was shown), [23] (where the implementation of FDE in the finite element method together with 3D AAA model were discussed), which enables to model the memory (nonlocality) as time dependent, will play a crucial role.
The paper is structured as follows. In Sect. 2, we introduce the time-fractional damping equation of motion. Section 3 presents the approach used to solve the equation of motion with time-fractional damping. In Sect. 4, the experimental data acquisition process is described. The effects of an application of fractional derivatives to the classical equation of motion are fully explained in Sect. 5. Moreover, Sect. 5 focusses on the insufficiency of the classical approach for a phenomenological description of debris motion and shows the capability of the fractional derivative approach. Inference is carried out on the basis of experimental results and numerical studies. Finally, Sect. 6 concludes the article.

Time-fractional damping equation of motion
The translational motion of the debris is considered. In classical analysis, the equation of motion has the form in which the inertia force is related to the product of mass (m) and acceleration ( d 2 ∕dt 2 ), while the damping force is proportional to velocity ( d ∕dt ) with the scalar damping factor c: where the damping factor c is a function of speed, c = c(||d ∕dt||) , and is the vector of the gravitational force, which for gravitational acceleration g is defined in the assumed Cartesian coordinate system as: In the above equations, denotes the vector of coordinates of the position of the debris in the Cartesian coordinate system and t is time.
In this paper, the classical approach that relates the damping force to the velocity is modified. The velocity d ∕dt , which is the first time derivative of the position vector , is replaced by the fractional derivative C t− t D t , where C D is the left-sided Caputo derivative, is an order of derivative and t is a time length scale, which defines in general nonconstant memory, so called "short memory", and therefore, as will be proved in proceeding sections, enables better description of the complexity of the real drag force. The left-sided Caputo derivative of a function f(t) (the function f(t) should have n + 1 continuous bounded derivatives in [t − t , T] for every T > t − t ) is defined as [29]: in which Γ is the Euler gamma function and n is the integer order of a derivative of a scalar function f defined as: where ⌊ ⌋ is an operator that rounds to the nearest integer less than or equal to the number that it acts on, that is, ⌊1.834⌋ = 1 . Such approach, as will be shown in the proceeding sections, considerably improves the description of a drag force, and therefore the precision of modelling debris motion.
In the case of a debris motion analysed, we replace the scalar function f(t) in Eq. (3) with the vector function (t) , which constitutes the solution of our problem. Changing the first derivative component in Eq. (1) by the fractional derivative in Eq. (3), the postulated form of the equation of motion, which takes into account the fractional derivative in the damping term, takes the form: where ĉ is the modified damping factor, which takes into account the inconsistency of units between the classical velocity (the first time derivative of displacement) ( m/s ) and the fractional velocity ( -fractional derivative of displacement) ( m/s ). In this work, the modified damping factor ĉ is defined as a product of the damping factor c and the time-length scale t : The emphasis should be put on the understanding of the parameters related to the fractional derivative in Eq. (5), i.e. the order of the derivative (non-dimensional) and the time length scale t (s in SI units). The parameter can be interpreted as the damping order, whereas t is the motion memory. The parameters and t are postulated to be two additional motion parameters. Both and t need identification. The order of the derivative is assumed to be in the interval: The set of possible values of is presumed based on physical premises. Note that = 0 would represent motion without damping, while for = 2 the damping force would be related to acceleration. Both = 0 and = 2 are excluded. Moreover, for = 1 the classical approach is recovered, that is, Eq. (1).

Odibat's approximation of the Caputo fractional derivative
The calculation of the left-sided Caputo derivative Eq. (3) requires the calculation of the Euler gamma function Γ and the integral Eq. (3). In this work, Odibat's approximation of the Caputo fractional derivative is used [30]. Moreover, the definition of Eq. (3) assumes the known motion history in the time interval from t − t to t, where t is the current instance of time. In the Odibat's procedure, the time interval [t − t , t] is subdivided into s equal intervals as so the time step h is equal to Next, by applying a modification of trapezoidal rule for definite integrals, the approximation of the Caputo fractional derivative of the displacement vector function is then where the parameter n is an integer given by the formula Eq. (4). From Eq. (10) one can state that the Caputo derivative of a certain function is equivalent to a weighted sum of its integer order derivatives values at specified points dependent on the time length scale t and parameter s. Taking into account all possible values of given in Eq. (7), n can be equal to 1 or 2 following the relation:

The equation of motion with fractional derivative damping term for discrete time intervals
The equation of motion Eq. (5) is rewritten at the time instance t = t i+1 : where the Caputo fractional derivative is approximated using Odibat's approximation formula Eq. (10), in which s is an integer related to motion memory, that satisfies the subdivision of the motion history given by Eq. (9). Eq. (10) can be reformulated in the form that separates the quantities related to the motion history ( t ≤ t i ) from the unknown quantities at time t i+1 : where h Γ is defined as: , p 0 is a scalar value calculated from the formula: k is a row vector in which the elements are calculated for k = 1, ..., s − 1: and k is a matrix for which the n-order time derivative is defined as: In this work, n can take only two values: 1 or 2, see Eq. (11), hence, for simplification of notations, starting from Eq. (18) the dot notation is used for integer order derivatives of 1 and 2, i.e. Therefore, for n = 1 , Eq.
Moreover, Odibat's approximation for the fractional derivative with 1 < < 2 ⟹ n = 2 results in the fractional derivative term, defined as in Eq. (19), to be dependent only on the acceleration. This leads to a limited usage of fractional derivatives with 1 < < 2 in this study.

Computation of the first phase of the motion
As long as the left-sided fractional derivative is considered, the calculation of the fractional derivative is based on the motion history related to the time period called motion memory. It is obvious that 'the past' is unknown at the initialisation of the motion. However, the pre-initial state is generated based on the classical initial conditions ( 0 -position, ̇ 0 -velocity, ̈ 0acceleration). The state at time step t −1 = −h is generated for the time interval h calculated based on Eq. (9). The generation of the pre-initial state is done by substitution of known values of initial conditions into the backward finite difference scheme, which gives: At the first two time steps, the memory frame, that is, the time history for which the fractional derivative is calculated, is set to 2h and the number of intervals for the calculation of the fractional derivative m is assumed to be equal to 2. Then, the memory frame increases gradually with the increase of simulation time; see Fig

Time integration
The equation of motion, Eq. (20), is solved numerically using the Newmark method. The reformulated algorithm provided in [31] is used for the calculation. Newmark methods allow for the solution of equations of the following form: where the definitions of matrices ̂ , ̂ , ̂ vary with respect to the value of , see Eq. (20). Below, k is the vector defined in Eq. (16) and k is a matrix, as in Eq. (17). For 0 < ≤ 1: The initial conditions are given as follows: The initial value of ̈ (0) =̈ 0 is unknown. Therefore, an approximation, which does not take into account the fractional derivative, is used: In the Newmark method, the approximation of an unknown function and its time derivatives have the following forms: in which and are called parameters of the Newmark method. For the purpose of this work, the constant average acceleration method is selected (although for more complex motions compared to the one analyzed in this paper, other approaches might be more effective), i.e. = 0.5 , = 0.5 . The time step is chosen on the basis of the convergence analyses.

Experimental data
In the experiment HSDM is induced by controlled explosion. Measurements were carried out in 2018. In the experiment, the debris (M6 nuts) were placed about 0.9 m above the ground and after an explosion they were recorded by a high-speed camera; see Fig. 2. The spread of debris is limited by horizontal shields, which ensure the movement of debris in one vertical plane. The experimental setup is comprehensively described in a separate work by Sielicki et al. [8]. However, the objective of [8] was different from the objective of this work. Therefore, the processed data from [8] cannot be taken directly as input for current analyses. Hence, the raw data from the experiment were once again processed, as shown in  the following. In this study, none of the components of the experimental setup is important except the motion of the debris themselves.
Trajectories represented by the coordinates x and z are visible in Fig. 3. In Fig. 3, each colour represents a different projectile, and each data point (cross in the graph) represents the position of the projectile at a different time instance. The time plays here a crucial role because the validation requires the comparison of results in four separate domains: x-coordinate vs. time, z-coordinate vs. time, x-velocity vs. time, z-velocity vs. time. The experiment was able to capture only the very beginning of the debris motion. In Fig. 3, at first glance the trajectories seem to be linear. However, the trajectories are curved. Fundamentally, the change of position of a projectile with respect to time is parabolic when damping is neglected. The position-time curve receives even more complicated shapes when considering damping.

Data acquisition and measurement error
The high-speed camera used in the experiment is set up to record a frame each 52.62 s. The frame has a resolution 1280 px x 304 px. In Fig. 2, the line representing 10.8 m has 927 px. The ratio 10.8/927 is used as a scale factor during data acquisition. The measurement error can be calculated directly from a scale factor. Knowing that a pixel is the smallest addressable element in an image, the smallest capturable physical object is about 1.2 cm. The size of the projectiles (M8 nuts) is approximately 1.3 cm. This means that we can estimate the measurement error at about 1.2 cm. However, it should be noted that the distance measurement error of about 1.2 cm does not influence the debris trajectories but produces significant noise in the calculated pseudo-instantaneous velocities. The notion "pseudo-instantaneous" is used to indicate that the velocity is calculated based on measurements, so the time increment is very small, but not infinitesimally small. Data points are collected every 5 frames, which means that the measurement time interval Δt measurement is equal to 263.1 ⋅ 10 −6 s. Marking the distance, which would be measured without any measurement error as L ideal , and the distance measurement error as L , we obtain an expression for a measured velocity: where the velocity measurement error is calculated as L ∕Δt measurement = (1.2 ⋅ 10 −2 )∕(263.1 ⋅ 10 −6 ) = 45.61 m/s . Although the collected data are of the best quality (trajectories can be captured with an accuracy of 1.2 cm), the given velocity measurement error is still significant.
To reduce the noise in an acquired velocity graph, central differences are used for the velocity calculation. The central difference spanning 9 data points is used, which results in 8-th order accuracy of the first derivative approximation. However, this manipulation does not completely solve the problem of noise in measured velocities. Therefore, the exponential regression is used for the velocity approximation. The standard error of the regression is calculated to provide the absolute measure of the distance between the calculated velocities and the regression line. Such acquired experimental data are used hereafter for validation of motion models, see Figs. 4, 5, 6 and 7.

Interpretation of fractional material parameters
Drag force depends on the velocity and spatial orientation of the projectile and is an integral of the stresses induced on the surface of the projectile during its flight. In general, the compressibility of the fluid and the rarefaction effects at the fluidstructure interface are responsible for the distribution of these stresses, hence the drag force. From a physical point of view, compressibility and rarefaction effects depend on Reynolds and Mach numbers [7], which could be practically obtained from the Computational Fluid Dynamics (CFD) analysis. Eventually, detailed calculation of the projectile flight requires coupled motion analysis and CFD modelling to determine a drag force. Furthermore, rotational inertia influences the dissipation of the total energy carried by a moving projectile. As mentioned, the practical application of this research relates to the determination of safety zones in the case of HSDM. Therefore, the usage a fast and robust phenomenological model is crucial in this case. Fundamentally, neither the classical approach nor the incorporation of fractional damping can account for these effects as comprehensively as the coupled motion-CFD model, which, on the other hand, is difficult to be practically used especially due to time consuming calculations. However, the use of fractional damping increases the number of parameters in the equations of motion, improving the ability of the model to predict HSDM, but at the same time, calculations remain fast and robust. According to the classical approach, the equations of motion are of the form given in Eq. (1). Except for the initial conditions, which are evaluated based on the experimental results, the damping factor c is the only unknown parameter of an Eq. (1), in which the drag force is equal to ċ . With the well-known definition of a drag force D : the expression for the damping factor comes naturally: where notion || ⋅ || represents the Euclidean norm of a vector, i.e. ��̇ �� =

√̇
⋅̇ , C D is the drag coefficient, is the density of the fluid (air) and A p is the projectile area projected onto the plane perpendicular to the direction of movement. Here, the drag coefficient C D , as well as fluid density , are assumed to be constant. For the given projectiles (M8 nuts used in the experiment), the value of C D , which corresponds to the experimental data, has been found to be equal to 2.0, for the air density taken as 1.2041 kg ⋅ m −3 .
data points that represent the experimental results reveal a significant scatter for the horizontal speed, and the vertical speed. This is due to the limited accuracy of the measurement, which was explained in Sect. 4.2. However, the regression curves could be plotted, which allows one to recognise the trend interpreted as the physical behaviour. Both the classical approach and the fractional derivative approach perform similarly, indicating that they could both represent the phenomena. Furthermore, in this section, the benefits of using the fractional derivative approach are shown, for example, in the explanation of Figs. 8 and 9 (the validation flowchart is shown in Appendix A and the raw data from the field experiments is given in Appendix B).
The calculated values from the classical formulation, presented in Figs. 4, 5, 6 and 7, are given for two extreme values of the projected frontal areas of the projectile: A p−min = 105.5 mm 2 and A p−max = 129.3 mm 2 . The distance and elevation curves are not significantly affected by the value of A p in the range that the experiment covers; however, the differences in the speed-time relationship, especially the one that represents the horizontal speed, are noticeably influenced by the assumed frontal area of the projectile. Note that regardless of the frontal area value, the speed vs. time curves calculated using the classical approach are in between the dashed lines, which represent the standard error of the regression. The area marked by the yellow colour represents the effect of the projectile obliquity angle, which can be recovered by applying the time-fractional derivative damping term. Technically, it is obtained by ranging the parameter from 0 to 1.
The fractional derivative solutions presented, yellow areas in Figs. 4, 5, 6 and 7, are developed based on the minimum frontal area of the projectile, and hence for the case of = 1 , the dotted black line is recovered. Again, the curves that represent the relationship between speed and time are particularly affected by changes in the parameters. The experimental results that the validation study refers to are given for time periods around 0.01 -0.02 s. Standard convergence studies as for finite difference calculations have shown that for these It is shown (Figs. 4, 5, 6 and 7) that manipulation of fractional derivative parameters gives promising potential to adjust the calculated motion to experimental results. Due to the fact that the time-fractional approach is nonlocal, the impact of the fractional derivative parameters on the projectile motion accumulates during the process of motion, which can represent the influence of compressibility, rarefaction effects as well as the lift. The classical approach does not have sufficient number of parameters to do the same.
Hence, the non-local time-fractional derivative approach introduces in a homogenised way complex phenomena connected with local effects between air and moving projectile. As will be presented in Sect. 5.3 the effect of fractional modelling is increasing with respect to the advancement of projectile motion.

Analysis of projectile trajectory and velocity
The numerical examples are constructed based on the information arising from the experimental data and calculations, which have been graphically summarised in Figs. 4, 5, 6 and 7. The initial speed of a projectile is taken equal to 500 m/s, which is a working estimate of the speed of M8 nut 4 m from the blast source. The drag coefficient C D is assumed to be constant and equal to 2.0. The area of the projectile projected on the plane perpendicular to the moving direction A p is taken as the minimum possible value, which, for M8 nut is equal to 105.5 mm 2 .
Taking the minimum frontal area as the reference is justified by two reasons: (A) the minimum frontal area corresponds to the minimum drag, which gives conservative results; (B) flying projectiles tend to stabilise their motion to be parallel to the inertia axis corresponding to the lowest moment of inertia, which corresponds to a small or the smallest possible frontal area. With these conditions, the following example shows how the order of the fractional derivative and a time length scale t influence the trajectory and speed of the projectile in motion calculated based on Eq. (20) 1 .
In Fig. 8, a change in trajectory is given due to alteration of the order of the fractional derivative for constant values of the time length scale. The calculation is made for varied from 0.0001 to 1.0 with the step equal to 0.01. Resulting trajectories are given in grey colour. The red and upper black curves represent the cases in which trajectories have been obtained using the classical equation of motion, i.e., = 1.0 , respectively, for the case of the maximum and minimum frontal area of the M8 nut. The intermediate black curves are the trajectories for from 0.0001 to 1.0 with step 0.1. Each black line is labelled with the value of to which it corresponds. It is shown that for a given time length scale, an appropriate choice of the derivative order modifies the trajectory in a way that the resulting trajectories range between the curves representing maximum and minimum drag, which is along with physical interpretation discussed in Sect. 5.1.
Furthermore, the effect of the change in order on the fractional derivative of the trajectory depends on the value of the time-length scale t . In Fig. 9, the ranges of trajectories obtained for different values of t are given.The higher the value of the time length scale, the wider the range of trajectories that could be calibrated using the order of the derivative . Interpreting the parameter t as motion memory, it can be said that, taking into account information about the motion for a longer time, it is possible to cover a wider variety of possible trajectory curves (cf. Sect. 5.1).
The value of the time-length scale has its limit, which is equal to the total time of motion at the current time step. It is visible in Fig. 10, where an analysis is carried out on the order of fractional derivative = 0.85 and the influence of the time length scale on the trajectory is shown. Furthermore, considering the curves representing t = 1 − 10 ms in Fig. 10, it could be observed that decreasing the value of the time length scale, at some point, results in a convergent solution, i.e. further decrease in t does not influence the solution. Therefore, a solution for a given ∶ 0 < < 1 may be calibrated using the time length scale t , however, such manipulation has its lower and upper bounds.
As mentioned in the Introduction, the velocity of debris is of special interest in the analyses of casualty and injury risk from being hit by debris. The correlation between velocity and death probability can be explicitly given as in [32] or by reference to the kinetic energy [33,34]. Regardless of a correlation, safety zones are always defined based on the relation between the velocity and distance of the flying projectiles. With this respect, continuing a consideration about the given numerical example, the changes in the projectile's velocity with respect to the processed distance are shown in Figs. 11 and 12: first, due to change of the order of the fractional derivative ; second, due to change of the time length scale t . The notions given in the subsequent figures correspond to the notions used before; therefore, no redundant explanations are given.
It is shown that it is possible to find a particular time-length scale for which altering the order of the fractional derivative from 0 to 1 covers the range of velocities between the cases calculated using the classical approach for minimum and maximum drag; see Fig. 11. On the other hand, increasing the time length scale (motion memory), Fig. 12, reduces the velocitydistance curves, which means that the projectiles tend to experience higher deceleration. This corresponds to the findings of the shape of the trajectory curves shown in Figs. 9 and 10.
The value of the time length scale influences the process of a deceleration of projectiles as shown in Fig. 13. In Fig. 13, the velocity histories resulting from the fractional derivative approach are related to the minimum drag case calculated using the classical approach. For comparison, the curves related to maximum drag and intermediate drag, i.e., A p = 0.5(A p−min + A p−max ) , are given. Therefore, one more time, the physical interpretation shown in Sect. 5.1, is captured.
Furthermore, in the classical approach, the relation between the velocity history for minimum drag and the velocity history for cases with higher drag is proportional (black, red, and  Fig. 13), while modifying the order of fractional derivative and the time length scale makes it possible to deflect this relation in a way that increases the drag at lower velocities. This capability of the fractional approach is even more promising when the results are set together with experimental measurements of drag coefficients, which is explained in the following section, where the consideration are carried out for equivalent spherical projectiles due to three reasons: (1) there exist experimental data on the relation between the Reynolds number and the drag coefficient for such geometry, (2) rotation does not affect the drag, hence, (3) the Reynolds number can be precisely calculated at each instance of time.

The equivalent drag coefficient analysis
The term equivalent drag coefficient is developed for the purpose of this article. It aims to make it possible to compare the drag resulting from the fractional derivative approach with the drag coefficient, as it is understood in the classical approach. Recall that in the classical approach, Eq. (1), the drag force drag is defined as a product of the damping factor c and the velocity vector ̇ : On the other hand, in the fractional equation of motion, Eq. (5), the role of the damping force is played by a component related to the fractional derivative, so: However, from the equilibrium condition, both in the classical approach and in the fractional derivative approach, the damping force must be equal to: where is defined in Eq. (2). Using Eq. (45) for the evaluation of the drag force is more convenient in the post-processing of data, especially in the case of fractional derivative equations of motion. Knowing the value of the drag force, it is possible to calculate a relation for the (variable) drag coefficient, which gives an equivalence between the fractional solution and the classical approach, namely: where is the density of the fluid (1.2041 kg/m 3 for air) and A ps is the frontal area of a spherical projectile. A vector resultant from Eq. (46) has its components in the direction of the Cartesian coordinate system. To obtain a scalar value of the drag coefficient in the direction parallel to the projectile movement C eq D , the vector eq D must be projected on the direction given by the velocity vector ̇ , hence: in which the nominator represents the dot product of the vectors additional model for the evolution of the drag coefficient. It is postulated that the calibration of the equations of motion by the two scalar parameters and t is more convenient than the manipulation of an additional function of the drag coefficient. During the motion of each sphere, an equivalent drag coefficient, defined in Eq. (47), and the Reynolds number in flow past a sphere are calculated, giving the possibility of validating the results with the experimental data given in [7] and [35]. C eq D vs. Re relationships are given in Figs. 14, 15, 16 and 17. Note that at the beginning the moving sphere has the maximum speed, so the Reynolds number is the highest. Then, the projectile slows down to the point where the gravitational force starts to dominate and the sphere slightly speeds up. Therefore, the graphs in Figs. 14, 15, 16 and 17 could be analysed from right to left, analogously to the way the projectile fly.
The shape of C eq D vs. Re graphs is correlated with the known experimental values of C eq D for the sphere, which means that due to the change in the range of Re for different diameters of spheres, the resulted shape of C eq D vs. Re relations have to be adjusted. It is done by modifying the time length scale in a way to fit the position of the peak value of C eq D and then magnifying the effect by selecting appropriate orders of the fractional derivative .
It has been found that, as a rule of thumb, the value of the time length scale t is responsible for the position and shape of the peak at C eq D vs. Re graph, whereas the order of the fractional derivative magnifies the value of C eq D . Furthermore, as the time length scale increases, the peek visible in the C eq D vs. Re graphs becomes more smooth.
In Figs. 14, 15, 16 and 17 it is shown that the range of Reynolds numbers of moving spheres decreases with increasing sphere diameter. However, for the expected velocities of spherical debris, the maximum experienced Reynolds number does not significantly exceed the value of 300,000 at which a sudden drop of C eq D is encountered during experiments. In such cases, adjusting the parameters of the fractional derivatives provides a better fit for an equivalent drag coefficient. In conclusion, by comparison to the results for spherical particles, it is shown that the fractional derivative approach, with an appropriate calibration, can provide sufficiently precise results. However,  [7,35] the parameters given in the fractional derivative approach do not necessarily limit the model to the cases of spherical particles. The fractional derivative approach can be successfully used in estimating HSDM and its modelling capacity is greater compared to the classical case.

Conclusions
The work can be concluded by 5 points: • Using fractional derivatives for the damping term of the equations of motion of debris allows us to make the motion model more flexible; i.e., it gives two more parameters that may be modified to capture the real motion of debris. • Adjusting the two parameters of the fractional derivative: the order of the derivative and the time length scale t allows for a better fit of the debris motion curves, i.e. position vs. time and velocity vs. time, without incorporation of additional, sometimes complicated, functions of the drag coefficient, as in the classical approach. • The fractional derivatives approach has no limitations with respect to the classical approach, since the classical approach is a special case of the fractional derivatives approach for = 1.0.  [7,35] influences the range (distance) debris flight. However, in detail, the influence of these parameters on the motion of debris is more complex. • The way in which the combination of and t influences the motion of debris may be explained using the (variable) equivalent drag coefficient, i.e., modifying and t results in a motion, which in the classical approach would be modelled with a non-linear drag coefficient.
In summary, the fractional derivative approach (non-local in time) is more reliable compared to the classical approach (local in time). The effectiveness of the fractional derivative framework has been verified on the basis of selected experimental results.

Appendix B. Position -time data from experiments on PBIED
See Tables 1, 2, 3 and 4.