J. Wind Eng. Ind. Aerodyn. 111 (2012) 95–103 Contents lists available at SciVerse ScienceDirect Journal of Wind Engineering and Industrial Aerodynamics journal homepage: www.elsevier.com/locate/jweia An investigation of plate-type windborne debris ﬂight using coupled CFD–RBD models. Part I: Model development and validation B. Kakimpa, D.M. Hargreaves n, J.S. Owen Faculty of Engineering, University of Nottingham, Nottingham NG7 2RD, UK a r t i c l e i n f o abstract Available online 22 August 2012 The development of a coupled computational ﬂuid-dynamics rigid body (CFD–RBD) model is presented. The RBD model deploys rotational quaternions, which are free from the gimbal lock that is associated with Euler rotational matrix. The quaternion model means that the complex 3D spinning ﬂight modes associated with the ﬂight of plate-type windborne debris can be modelled robustly. This paper attempts to determine the accuracy of the CFD–RBD model by comparing the predicted trajectories from a large number of debris simulations with experimentally derived equations of best ﬁt. Agreement is found to be good and, based on the ﬁndings, an alternative form for the dimensionless ﬂight distance is presented, which extends the range of the experimental study to longer ﬂight times. The predictions from the CFD–RBD model are then compared against two quasi-steady analytical debris ﬂight models. The second model is based on modiﬁed force and moment coefﬁcients, which are informed by the ﬁndings from the CFD–RBD model. For plates that have attained a stable, autorotational ﬂight mode, the CFD–RBD and analytical models are in good agreement. Their predictions differ during the initial stages of ﬂight, where the complex non-linear interactions between the plate and its wake are not captured by the analytical models. & 2012 Elsevier Ltd. All rights reserved. Keywords: CFD Rigid body dynamics Debris 1. Introduction Windborne debris can be deﬁned as failed building components, loose items and street furniture that are picked up and carried by the wind during severe storm events. The debris is most commonly classiﬁed into compact, rod and sheet (or plate)-type debris after Wills et al. (2002). Early studies of building performance in severe storms, such as Minor and Beason (1976) identiﬁed windborne debris as a principal cause of building envelope failure and an essential component of the debris damage chain. More recently, damage assessment reports of extreme windstorms in the UK, such as the Birmingham Tornado in 2005 (Marshall and Robinson, 2006) have also revealed a deﬁnitive link between damage to downwind properties and windborne debris generated from upstream structures. Current debris risk and damage models, such as Lin and Vanmarcke (2010), rely on estimates of debris trajectories and impact kinetic energy obtained from analytical models of debris ﬂight. These invariably include two-dimensional (2D) models, as in Tachikawa (1983), Wills et al. (2002), Holmes et al. (2006a), n Corresponding author. Tel.: þ44 115 846 8079; fax: þ 44 115 951 3898. E-mail addresses: [email protected] (B. Kakimpa), [email protected] (D.M. Hargreaves), [email protected] (J.S. Owen). 0167-6105/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jweia.2012.07.008 with no horizontal crosswind motion of the debris, and threedimensional (3D) models, such as Baker (2007), Richards et al. (2008), Kordi and Kopp (2009b), which allow for six degree of freedom (6DOF) debris motion. Tachikawa (1983) proposed a set of non-dimensionalised equations of motion 2 d x ¼ K½ð1uÞ2 þ v 2 ½C D cos bC L sin b, dt ð1Þ 2 d y ¼ 1K½ð1uÞ2 þv 2 ½C D sin b þC L cos b, dt ð2Þ 2 d y ¼ KFr 2L D½ð1uÞ2 þv 2 C M , dt ð3Þ where x ¼ gx=U 2w , y ¼ gy=U 2w , t ¼ gt=U w , u ¼ u=U w , and v ¼ v=U w and where Uw is the mean wind speed, g is the acceleration due to gravity and C D ,C L and CM are the aerodynamic drag, lift and moment coefﬁcients, respectively. This formulation of the debris ﬂight equations showed that the debris ﬂight was controlled by a number of dimensionless parameters. Chief among these was the ratio of aerodynamic force to gravitational force, K¼ rU 2w A 2mg , 96 B. Kakimpa et al. / J. Wind Eng. Ind. Aerodyn. 111 (2012) 95–103 where r is the density of air, A is the area of the plate and m is the mass of the plate. Later, the parameter K was called the Tachikawa number in recognition of this pioneering work (Holmes et al., 2006b). In addition, Tachikawa used a Froude number, Uw FrL ¼ pﬃﬃﬃﬃﬃ , gL and a dimensionless mass moment of inertia parameter, D¼ mL2 , I where L is the characteristic length of the plate and I is the mass moment of inertia. Tachikawa assumed a decomposition of the unsteady aerodynamic force and moment coefﬁcients, C, of the form C ¼ C f ðaeff Þ þ C R ðo Þ, ð4Þ where C f ðaeff Þ is the ﬂuctuating component, which is a function of the plate’s instantaneous angle of attack, aeff , and C R ðo Þ is a mean component, which is a function of the plate’s non-dimensionalised angular velocity, o ¼ oL=U w . Expressions for CR were obtained from experimental measurements of the time-averaged force and torque acting on a rotating plate. Tachikawa assumed C f ðaeff Þ ¼ C s ðaeff ÞC s , ð5Þ where Cs are the instantaneous drag, lift and torque coefﬁcients for a static ﬂat plate and C s is the average of Cs over a complete rotation. The quasi-steady analytical approach however remains limited to simple cases and is incomplete. Studies of free-falling plates have revealed that a quasi-steady empirical treatment of platewake interaction effects is unable to account for cases where the plate interacts strongly with its own wake, or with vortex structures in the surrounding ﬂow (Andersen et al., 2005). This strong plate-wake interaction is especially signiﬁcant during the debris launch and initial ﬂight stages, where very high lift can be generated by these effects. Furthermore, recent experimental measurements on rotating plates (Martinez-Vazquez et al., 2010) have found that the modiﬁed quasi-steady torque used in analytical models differs considerably from the aerodynamic torque coefﬁcients measured from autorotating ﬂat plates. There is therefore a need to review the validity of quasi-steady theory in debris ﬂight modelling. In order to overcome these limitations, it is necessary to develop more complete picture of the aerodynamic forces acting on the plate. The approach adopted in the present work is to employ computational models that allow for the accurate modelling of plate rigid body dynamics (RBD) together with a numerical simulation of the unsteady ﬂow ﬁeld and its non-linear interaction with the plate. The motion of free-falling plates (Andersen et al., 2005; Jin and Xu, 2008) and the dynamics of shuttle ascent foam debris (Murman et al., 2005) have been successfully investigated using such an approach. Jin and Xu (2008) and Andersen et al. (2005), however, focus on low Reynolds number numerical and experimental investigations of the free-fall of 2D high aspect ratio plates. Their ﬁndings, therefore, need to be extended to the ﬂight of low aspect ratio 3D plates in high Reynolds number ﬂows, which have been experimentally observed to predominantly exhibit complex 3D spinning motion (Kordi et al., 2010). Murman et al. (2005) used 3D 6DOF models to investigate the ﬂight of frustum shaped foam debris shed by space shuttles during ascent in a more complete approach. A Computational Fluid Dynamics (CFD) model was coupled with a 3D rigid body dynamics (RBD) model, allowing a direct simulation of the high Mach number ﬂow around the shuttle during ascent, and the nonlinear interaction of the foam debris with the ﬂow. The shuttle foam debris ﬂew for relatively short periods (approximately 0.15 s) compared to windborne debris ﬂight, with limited rotational motion unless perturbed. This paper presents an improved approach to the numerical simulation of plate-type windborne debris ﬂight by using a 3D CFD model sequentially coupled with a six degree of freedom RBD model to simulate the ﬂight of low aspect ratio plates in wind ﬂow of Reynolds number 1:3 106 . The paper focuses on platetype debris – such as rooﬁng sheets, shingles and tiles – which has been found to be the dominant-type in residential settings (NAHB Research Center, 2002). The debris is assumed to be rigid and, for purposes of comparison with existing experimental data and analytical models, simulations have been restricted to simple regular shaped debris although various shapes could be studied using this method. The Euler rotational matrix approach previously used in Kakimpa et al. (2010b) is replaced with a singularity free parametrization of orientation based on rotational quaternions that is valid for all possible plate orientations. In the companion Part II paper (Kakimpa et al., this issue), this coupled CFD–RBD model has been used to investigate the inﬂuence of plate initial orientation, Tachikawa number, mass moment of inertia, geometric parameters such as thickness-ratio and aspect-ratio, as well as complex launch ﬂow ﬁelds. In the present paper, Part I, the CFD–RBD predictions are compared against experimental results where available, with new improved ﬁt expressions for the experimental data proposed. Finally, the results are compared against predictions from analytical models with a discussion of the limitations of quasi-steady theory. 2. CFD modelling A 3.0 kg ﬂat, square plate of side, L, 1.0 m and thickness, t, 0.0254 m is positioned in a domain of size 80L 15L 15L. Fig. 1 shows the computational domain and boundaries. The X-axis is aligned with the alongwind direction, the Y-axis with the vertical and the Z-axis with the crosswind direction. The plate – with its centre of mass initially 4 L from the inlet, 7.55 L from the side walls and 10 L from the bottom wall – is modelled as a rigid wall boundary with a no-slip condition. The plate is held within a spherical inner volume of radius 2 L, which rotates and translates with the plate’s instantaneous rotational and translational velocity as required in the Arbitrary Lagrangian–Eulerian (ALE) method for Fluid–Structure Interaction (Sarrate et al., 2001). The commercial CFD code ANSYS FLUENT (Fluent Inc, 2009) is used to solve the ALE formulation of the unsteady 3D incompressible Navier–Stokes equations using a Finite Volume discretisation r Ur ¼ 0, ð6Þ @U 1 þUr rU þ ðxg UÞ ¼ ðrP þ rsÞ, @t r ð7Þ Ur ¼ Uðxp rÞug , ð8Þ Fig. 1. Computational domain and boundaries of the free-ﬂight simulations. B. Kakimpa et al. / J. Wind Eng. Ind. Aerodyn. 111 (2012) 95–103 Fig. 2. Section through computational grid showing structured hexahedral mesh in the moving inner region surrounding the plate, together with an unstructured tetrahedral outer mesh that is re-meshed at each timestep. where Ur is the ﬂuid relative velocity vector, U is the ﬂuid absolute velocity vector, ug is the translational velocity vector of the moving inner mesh region, xp is the rotational velocity vector of the inner mesh region, r is the position vector of the plate’s centre of mass which corresponds to the origin of the rotating reference frame, P is the ﬂuid pressure, and s is the viscous stress tensor. Values for ug and xp are obtained from a rigid body dynamics (RBD) model described in Section 3. The domain is discretised into a 3D computational grid of approximately 106 cells with a structured mesh in the spherical inner region and an unstructured mesh in the outer region as shown in Fig. 2. FLUENT’s dynamic meshing was used to re-mesh the domain at every timestep in order to accommodate the moving spherical inner region. The inner region was not remeshed in order to preserve the mesh quality in the proximity of the plate. The inlet is modelled as a constant inﬂow velocity boundary, while the outlet is modelled as a constant pressure boundary, and the side boundaries as walls with a free-slip condition. With the free-ﬂight simulations, a mean wind speed, U of 20 m s 1 ðRe ¼ 1:3 106 Þ was imposed at the inlet for most of the simulations, with turbulence intensity and length scale of 1% and 0.02 m respectively, which are typical of low turbulence wind tunnel values (ESDU, 1970). In this study, the variation in wind speed and turbulence intensity in the atmospheric boundary layer has not been considered. A uniform distribution of velocity and turbulence quantities is speciﬁed at the inlet in all cases. Sensitivity studies for turbulence model, discretisation scheme, pressure–velocity coupling scheme, grid size and timestep size were performed and the ﬁndings are presented in Kakimpa et al. (2010b). The grid independence study was conducted with static and ﬁxed-axis rotating plates, which used the same spherical inner region as used in the ﬂying plates simulations discussed here. When using dynamic meshing, the local mesh size around the rotating and translating inner region is constantly changing, making grid independence studies difﬁcult to analyse. The assumption was made that since the inner region, which contains 80% of the cells in the grid, was not remeshed and had passed grid independence studies for static and ﬁxed-axis rotation, then this grid was suitable for the ﬂying plate simulations. As a result of this study, the Realizable k–E model (Shih et al., 1995) has been used. An enhanced wall function that is valid for both ﬁne and coarse near-wall mesh resolutions (Kader, 1981) is required for the near-wall turbulence modelling due to the variations in Reynolds number (if relative wind speed is used as reference velocity) during debris ﬂight. The second order upwind scheme was used for advection terms, with a Rhie–Chow interpolation scheme (Rhie and Chow, 1983) required to interpolate 97 for cell-face pressures from cell-centred pressures in FLUENT’s collocated scheme. The SIMPLE method is used for pressure– velocity coupling and a time-step size of 5 103 s was found to be appropriate for all ﬂights at this Reynolds number. Unlike the quasi-steady models, the CFD simulations do not require any a priori knowledge of the aerodynamic characteristics of the debris and therefore make it possible to apply this model to the study of complex and irregular debris shapes for which no experimental data is readily available to provide a complete aerodynamic characterisation. The pressure and viscous forces acting on the body are easily obtainable at a high spatial resolution during the simulation, making it possible to study the unsteady ﬂuctuation of body forces on the plate. The complete ﬂow ﬁeld around the plate is simulated and coherent ﬂow structures may be identiﬁed (Kakimpa et al., 2010a) allowing a more direct simulation of their effects on plate motion, however this discussion is beyond the scope of the current paper. 3. Rigid body dynamics modelling 3.1. Euler equations Rigid body motion of the plate in real 3D space consists of at most six degrees of freedom—three translational degrees of freedom of a base point (usually the body’s centre of mass) and three independent rotational degrees of freedom about suitably chosen axes. The classical Newton–Euler equations of rigid body dynamics deﬁne a set of six differential equations of motion based on linear and angular momentum conservation principles. These equations provide the complete system of six scalar equations required in order to compute the general six degree of freedom (6DOF) motion of a rigid body m dug ¼ Fg , dt Ip d xp ¼ Mp xp Ip xp , dt ð9Þ ð10Þ where a p subscript indicates that a quantity is expressed in the plate-ﬁxed coordinate system and a g subscript indicates that the quantity is expressed in the global inertial reference frame. Further, m is the mass of the body, I is the mass moment of inertia tensor, u is the translational velocity vector, F is the applied force vector, x is the angular velocity vector and M is the vector of applied torque. Applied forces, F, and torque, M, are a combination of the CFD predicted aerodynamic forces and gravitational forces acting on the plate. The position of the body in three dimensional space is described by specifying its position relative to ﬁxed inertial reference frame represented in Fig. 3 by the Cartesian X g Y g Z g coordinate system. In order to specify the orientation of the body, a set of body-ﬁxed Cartesian axes, X p Y p Z p , with the origin at the body’s centre of mass and corresponding to the body’s principal axes are used. The orientation of the non-inertial body ﬁxed reference frame, X p Y p Z p , relative to the ﬁxed inertial reference frame, X g Y g Z g , is calculated via a translating inertial reference frame, X t Y t Z t , which has its origin at the body’s centre of mass. Fig. 3 illustrates these three reference frames. The Euler angles, f, y and c, result from a series of three rotations on the X t Y t Z t reference frame. The Tait– Bryan ZYX rotation convention is used and this consists of ﬁrstly a rotation of c about the positive Zt-axis, displacing the Xt- and Ytaxes (onto the dotted lines), followed by a rotation of y about the displaced positive Yt-axis, which rotates the displaced Xt-axis onto the ﬁnal Xp-axis. Thirdly, a rotation of f is performed about 98 B. Kakimpa et al. / J. Wind Eng. Ind. Aerodyn. 111 (2012) 95–103 Fig. 3. Illustration of the ﬁxed global inertial reference frame ðX g Y g Z g Þ, the translating inertial reference frame ðX t Y t Z t Þ and the rotating plate-ﬁxed reference frame ðX p Y p Z p Þ with the Euler angles – f, y, c – used to describe the plate’s orientation relative to the translating plate-ﬁxed reference frame. A cursory examination of Eq. (14) reveals that if y ¼ 7 p=2, _ and c _ can be inﬁnite for ﬁnite values of xp . This then f constitutes a mathematical singularity for the Euler angle parametrisation creating what is commonly known as gimbal lock. Stuelpnagel (1964) proves that it is topologically impossible to have a global 3D parametrisation of orientation based on Euler angles that does not have singularities. It is therefore necessary, when using Euler angles to ensure that the motion of the rigid body does not approach the singular orientations. However for chaotic 3D motions such as those in debris ﬂight, this might not be possible to judge a priori. A globally singularity-free parametrisation of orientation is therefore required. In addition, studies by Robinson (1958) proved that Euler angle parametrizations of orientation are less accurate and less computationally efﬁcient than other methods (e.g. Rotational Quaternions), especially when used to integrate incremental changes in orientation over time. This is partly due to the relatively higher computational effort (solving six additional equations) required to enforce the orthogonality constraint on the rotational matrix as well as the accumulating numerical errors in successive calculations, resulting in angular drift. However, when compared with the computational effort required to solve the Navier–Stokes equations, the computational beneﬁts of using quaternions are insigniﬁcant, but they remain a robust and accurate alternative to the Euler angle parametrizations. 3.3. Proposed singularity-free parametrization the Xt-axis, which is now at its ﬁnal Xp position, and this rotates the displaced Yt and Zt axes onto their ﬁnal Yp and Zp positions. 3.2. Review of Euler angle parametrizations Previous studies of debris ﬂight, such as Richards et al. (2008), employ a Euler rotational matrix, R, to transform coordinates from the inertial to the rotated reference frames according to zp ¼ Rðf, y, cÞzg , ð11Þ where zg is a vector in the global inertial reference frame, R is the rotational matrix, and zp is the same vector in the body-ﬁxed reference frame. A rotational matrix must always be orthogonal in order to constitute a pure rotation and its determinant must be equal to þ1 (Robinson, 1958). In addition to the rotational matrix, an Euler angle rates matrix is also required. This matrix, E, relates the rate of change of the Euler angles to the angular velocity of the body _, x ¼ EðWÞW ð12Þ where x is the angular velocity vector, W ¼ ½f, y, c and _ , y_ , c _ . _ ¼ ½f W During rigid body dynamics computations, the angular velocity, x, is not computed in inertial coordinates. Rather, the plate ﬁxed angular velocity, xp , is computed from solving the Euler equation of rotational motion, Eq. (10), in a body-ﬁxed reference frame, X p Y p Z p . As a result of this rotational motion, there will be a change in plate orientation. In order to update the rotational matrix, ðRÞ, it is necessary to compute the updated vector of Euler angles, W, by numerical integration of the Euler angle rates vector _ . This Euler angle rates vector is obtained from the computed W angular velocities via Eq. (13) using the inverse conjugate Euler angle rates matrix, ½E0 1 (Diebel, 2006) _ ¼ ½E0 ðWÞ1 xp , W 2 1 6 ½E ðWÞ1 ¼ 4 0 0 0 ð13Þ sin f tan y cos f tan y 3 cos f sin f sin f sec y cos f sec y 7 5: ð14Þ Due to the limitations in Euler angle parametrization of orientation, alternative parametrizations of orientation have been formulated, among which are the Euler parameters, commonly referred to as rotational or unit quaternions. A quaternion may be represented as a vector in four-dimensional space q ¼ ½q0 ,q1 ,q2 ,q3 T : ð15Þ The quaternion that arises from a rotation, a about an axis n is # " cos 12 a : q¼ ð16Þ n sin 12 a Rotational quaternions are additionally constrained to be unit quaternions, having a unit norm, JqJ, in order for them to represent a pure rotation (Greenwood, 2003). This implies that the algebraic constraint qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ð17Þ JqJ ¼ qT q ¼ q20 þ q21 þ q22 þq23 ¼ 1, must always be enforced in order to prevent any scaling of transformed coordinates, where q ¼ ½q0 ,q1 ,q2 ,q3 T is the adjoint of the quaternion. Given a vector, z in the global inertial reference frame, then z0 , its representation in the body-ﬁxed reference frame, can be obtained using rotational quaternions 0 z0 ¼ q ð18Þ q ¼ Rq ðqÞz, z 2 q20 þ q21 q22 q23 6 Rq ðqÞ ¼ 4 2q1 q2 2q0 q3 2q1 q3 þ 2q0 q2 2q1 q2 þ 2q0 q3 q20 q21 þq22 q23 2q2 q3 2q0 q1 2q1 q3 2q0 q2 3 2q2 q3 þ2q0 q1 7 5, q20 q21 q22 þq23 ð19Þ where Rq ðqÞ is a quaternion based rotational matrix. As with Euler angles, it is necessary to compute the vector of quaternion rates, _ and its mapping to the angular velocity vector x in body-ﬁxed q, coordinates. This is achieved using the inverse conjugate quaternion rates matrix, ½W0 ðqÞT , as _ x ¼ 2WðqÞq, ð20Þ B. Kakimpa et al. / J. Wind Eng. Ind. Aerodyn. 111 (2012) 95–103 q_ ¼ 12½W 0 ðqÞT xp , 2 q1 6 6 q0 ½W0 ðqÞT ¼ 6 6 q3 4 q2 99 ð21Þ q2 q3 q0 q1 q3 15 3 7 q2 7 7: q1 7 5 Expt−Fit CFD−Fit 12 ð22Þ q0 9 Unlike the Euler angle rates matrix, the Quaternion rates matrix is valid for all possible orientations in real 3D space. Because of their accuracy, computational efﬁciency and lack of any singularities, quaternions have become widely applied in rigid body dynamics applications, such as robotics and aerospace (Chou, 1992). The post-correction approach presented in Cline and Pai (2003) is utilized to correct the normality constraint error at each time-step. 6 3 0 0 2 4 6 8 10 Kt∗ 12 14 16 18 20 1.4 4. Model validation 1.2 4.1. Comparison with experiment 1 A preliminary veriﬁcation of the CFD–RBD model against experimental results is described in Kakimpa et al. (2010b). Subsequently, the CFD–RBD model has been rewritten to incorporate the singularity-free approach of Section 3.3. The results from new model have been compared against existing experimental data and ﬁt expressions for the asymptotic velocity and trajectories in Lin et al. (2006). A parametric study was performed to assess the sensitivity of plate-type windborne debris to initial orientation, relative to the wind. Free-ﬂight simulations were performed at a Tachikawa number, K, of 8.3, which is typical for large rooﬁng sheets in a wind storm. In the ﬁrst batch of 36 simulations the plates were held at initial angles of attack, in the vertical XY plane, az , Fig. 4, ranging from 851 to 901 at intervals of 51. The plates were held normal to the ﬂow in the horizontal XZ plane for these simulations. In these simulations there was little crosswind motion even though plates were free to translate and rotate in the crosswind direction. The trajectories, velocity, forces and moments of these simulations are discussed in more detail in Part II. Fig. 5(a) shows the non-dimensionalised horizontal distance, Kxn ¼ Kxg=U 2 , against non-dimensionalised time, Kt n ¼ Ktg=U for the simulations. It is clear from the ﬁgure that there is a large spread of trajectories, but that these clusters around the experimental ﬁt of Lin et al. (2006) Kxn 0:456ðKt n Þ2 0:148ðKt n Þ3 þ 0:024ðKt n Þ5 and sKx ¼ 0:134, 0.8 0.6 0.4 Expt Fit 0.2 0 0 2 4 6 8 10 12 14 Kx∗ Fig. 5. CFD–RBD predicted trajectories for initial angles of attack of 851 r az r 901, showing (a) experimental (Lin et al., 2006) and CFD–RBD based ﬁt expressions for non-dimensionalised horizontal distance, Kx , and (b) CFD–RBD predictions for non-dimensionalised horizontal speed, u, together with an experimentally derived ﬁt expression (Lin et al., 2006). % Table 1 Polynomial coefﬁcients for the ﬁfth-order rational polynomial in Eq. (24). C1 C2 C3 C4 C5 0.0174 0.0666 0.6404 0.0861 0.0030 D0 D1 D2 D3 D4 1.0000 0.9675 0.5874 0.0926 0.0034 n ð23Þ where sKxn is the standard deviation of the trajectory data used by Lin et al. (2006) to produce the ﬁt shown. In their experiments, Lin et al. (2006) varied the mass, dimensions and aspect ratio of the plate as well as the Reynolds number and not just the initial angle of attack, as discussed in the present work. Eq. (23) does not perform well as Ktn becomes large for long ﬂight durations. To overcome this difﬁculty, a ﬁfth order rational polynomial, with better extrapolation properties beyond the range of data used to derive it, is used Kxn ¼ C 1 ðKt n Þ þ C 2 ðKt n Þ2 þ C 3 ðKt n Þ3 þ C 4 ðKt n Þ4 þ C 5 ðKt n Þ5 D0 þ D1 ðKt n Þ þ D2 ðKt n Þ2 þD3 ðKt n Þ3 þ D4 ðKt n Þ4 , ð24Þ where the coefﬁcients Ci and Di are computed as shown in Table 1. Fig. 5(b) shows results for non-dimensionalised horizontal speed, u ¼ u=U, against non-dimensionalised horizontal distance, where U is the mean wind speed. Again the experimental ﬁt of Lin et al. (2006) is shown pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃnﬃ u 1e 1:8Kx and su ¼ 0:0814, ð25Þ Fig. 4. Nomenclature associated with the ﬂying plates, showing the angle of attack, az . where su is the standard deviation of the experimental data. For long ﬂight times, it can be seen that Lin et al. (2006) predict that the debris speed tends towards the mean wind speed. There was, however, a good deal of spread in their data and they were able to observe the over- and under-speeding of the debris shown in B. Kakimpa et al. / J. Wind Eng. Ind. Aerodyn. 111 (2012) 95–103 Fig. 5(b). By over-speeding it is meant that the debris ﬂies faster than the mean wind speed, having converted some of its potential energy to rotational kinetic energy. The resulting autorotation, depending on its sense, will cause over-speeding or underspeeding as seen in the ﬁgure. An explanation of over- and under-speeding is given in Part II (Kakimpa et al., this issue). A larger spread of terminal horizontal wind speeds is observed in the CFD–RBD model results, which give su 0:162, in contrast to experimental observations by Lin et al. (2006), where su 0:0814 was observed. It should also be noted that the expressions of Lin et al. (2006) are derived from a narrow range of data, with typical ﬂight times of t n r 0:8 (approximately 0.6 s), while the CFD data is derived for plates with longer ﬂight durations of 1:6r t n r 2:2. 1.5 Mean force or moment coefficient 100 4.2. Comparison with analytical models 4.2.1. The Quasi-steady 2 (QS2) model One of the fundamental assumptions of existing quasi-steady models, such as QS1, is a decomposition of aerodynamic coefﬁcients into a static and an autorotational component according to Tachikawa (1983). Recent experimental measurements of both mean and ﬂuctuating forces on a rotating plate (MartinezVazquez et al., 2010) have shown this decomposition to contain inaccuracies in predicting the ﬂuctuating component of aerodynamic force which Tachikawa (1983) did not originally measure but rather assumed to be equivalent to the ﬂuctuation of static force coefﬁcients about the mean. In an attempt to provide a more accurate representation of the aerodynamic forces on a rotating plate, an improved quasi-steady force model has been developed and which is called QS2 here. CFD models were created of ﬂat plates that were forced to rotate at a speciﬁed angular speed, using the same spherical inner region rotating inside a ﬁxed outer domain as was used in the free-ﬂight simulations described in Section 2. For a range of angular speeds from 0.0873 to 10.0 rad s 1, the drag, lift and moment coefﬁcients were extracted from the CFD simulations and the mean and rms of each are plotted in Figs. 6 and 7. Also shown in the ﬁgures are power law and exponential ﬁts to the data, the details of which can be found in Kakimpa (2011). The point where the mean nondimensionalised angular speed, o ¼ oL=U w C 0:75, corresponds to the expected autorotational angular speed of the plate, according to the predictions of Iversen (1979). Based on the results of these simulations, three distinct types of motion may be identiﬁed, depending on the angular speed: pre-autorotational with o t 0:75, aurotational with o C0:75 and post-autotational with o \ 0:75. In the pre-autorotational regime, the mean aerodynamic torque (C M in Fig. 6) acting on the plate is positive and serves to accelerate the plate, while in the post-autorotational cases, the aerodynamic torque is negative and opposes the plate rotation, resulting in an aerodynamic damping effect. Using the ﬁt expressions for the mean and rms force and moment coefﬁcients, it is possible for the QS2 model to work at angular speeds above and below from the stable autorotating case. So, the aerodynamic drag, CD, lift, CL, and torque, CM are expressed as functions of both the instantaneous non-dimensionalised rotational speed, o , and the effective angle of attack, a % % % % % 0.5 0 −0.5 −1 0 0.5 1 1.5 2 2.5 Fig. 6. CFD predictions and data ﬁts for the mean drag, CD, lift, CL, and moment, CM, coefﬁcients against non-dimensionalised angular speed, o , for a forced rotating plate. % 2 1.8 rms force or moment coefficient The predictions from CFD–RBD models are compared against those based on quasi-steady numerical solutions to the debris ﬂight equations (Tachikawa, 1983). Two different quasi-steady models have been considered, a recent 2D model proposed by Kordi and Kopp (2009a), hereafter referred to as QS1, and an improved quasi-steady force model based on the ﬁndings on forced rotating plates, hereafter referred to as QS2. 1 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 0 0.5 1 1.5 2 2.5 Fig. 7. CFD predictions and data ﬁts for the rms drag, CD, lift, CL, and moment, CM, coefﬁcients against non-dimensionalised angular speed, o , for a forced rotating plate. % according to qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ b D Þ2 ðC D Þ2 Þsinð2ap=2Þ, C D ¼ C D þ 2ððC C L ¼ k1 C L þ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ b L Þ2 ðC L Þ2 Þsinð2aÞ, 2ððC qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ b M Þ2 ðC M Þ2 Þsinð2aÞ, C M ¼ k1 C M 2ððC where 8 o < k1 ¼ 9o 9 : 0 ð26Þ ð27Þ ð28Þ % % if o a 0 ð29Þ if o ¼ 0 b denote the mean and root mean square of the force The C and C coefﬁcients over a complete rotational cycle, calculated using the ﬁt expressions shown in Figs. 6 and 7. 4.3. Results In the companion paper, Part II (Kakimpa et al., this issue), three ﬂight modes for the debris are identiﬁed. Depending on the B. Kakimpa et al. / J. Wind Eng. Ind. Aerodyn. 111 (2012) 95–103 initial angle of attack, a0 , the plates may enter ﬂutter, transition or autorotation ﬂight modes. Plates that ﬂutter never complete a rotation in one direction before they hit the ground—they simply rotate in one direction and then the other and so on. Autorotating plates immediately enter a stable autorotation state, never reversing their initial sense of rotation. Those plates described as transitional may ﬂutter several times before but eventually entering a state of stable autorotation. The results of the two analytical models, QS1 and QS2, are compared with CFD predictions for a case in each of the ﬂutter, transition and autorotation ﬂight modes. No comparisons are 101 made against the complex 3D spinning mode as the analytical models discussed in this section are limited to 3 DOF motion. The results from four cases are presented in Fig. 8. These include one ﬂutter mode case of ao ¼ 90J , shown in Fig. 8(a, d and g), one transitional mode case of ao ¼ 55J , shown in Fig. 8(b, e and h) and two autorotational mode cases with ao ¼ 730J , shown in Fig. 8(c, f and i). Neither quasi-steady models, QS1 and QS2, agree with the CFD–RBD predictions for the ﬂutter and transition mode plates. In the autorotational cases, shown in Fig. 8(c, f and i), the CFD–RBD and quasi-steady models predict the same mode of ﬂight, 2 2 2 1.5 1.5 1.5 1 1 1 0.5 0.5 0.5 0 0 −0.5 −0.5 −1 −1 −1 −1.5 −1.5 −1.5 −2 −2 −2.5 −2.5 0 −0.5 −2 −2.5 −3 −3.5 −3 −3 −4 −3.5 −3.5 −4.5 0 2 4 6 8 10 12 14 16 0 2 4 6 8 10 12 14 16 1.4 1.4 1.4 1.2 1.2 1.2 1 1 1 0.8 0.8 0.8 0.6 0.6 0.6 0.4 0.4 0.4 0.2 0.2 0.2 0 0 2 4 6 8 10 12 14 16 0 0 2 4 6 8 10 12 14 16 0 0.5 0.5 1 0 0 0.5 −0.5 −0.5 0 −1 −1 −1.5 −1.5 −2 −2 −2.5 −2.5 −3 −3 −3 −3.5 −3.5 −3.5 0 2 4 6 8 10 12 14 16 0 2 4 6 8 10 12 14 16 0 2 4 6 8 10 12 14 16 −0.5 −1 −1.5 −4 0 2 4 6 8 10 12 14 16 −4 −2 −2.5 0 2 4 6 8 10 12 14 16 −4 Fig. 8. CFD–RBD results (CFD) and quasi-steady analytical model (QS) predictions of (a,b,c) rotational speed, (d,e,f) translational speed and (g,h,i) trajectory, for ﬂutter (left), transitional (centre) and autorotational (right) ﬂight modes. 102 B. Kakimpa et al. / J. Wind Eng. Ind. Aerodyn. 111 (2012) 95–103 3.5 2.5 1.5 0.5 0 0 30 60 90 120 150 CFD 180 2.4 0.3 1.6 0.2 0.8 0.1 0 0 −0.8 −0.1 −1.6 −0.2 −2.4 0 CFDavg 30 60 QS1 90 120 150 180 −0.3 QS1avg 0 QS2 30 60 90 120 150 180 QS2avg Fig. 9. A comparison between phase-averaged CFD–RBD aerodynamic forces (CFD) and quasi-steady forces (QS1 and QS2) averaged over consecutive rotational cycles during the stable autorotational stage of ﬂight. Data presented is for a plate with initial az ¼ 301. together with comparable predictions of the debris trajectories. There is a more detailed discussion of the similarities of ﬁxed-axis autorotation and autorotating free ﬂying plates in Section 2.2 of Part II (Kakimpa et al., this issue). Although the quasi-steady models based on ﬁxed-axis autorotation theory are able to adequately represent the autorotational mode of ﬂight, they are unable to capture the strongly unsteady Fluid-Structure Interaction (FSI) involved in the ﬂutter and transitional mode cases. In both modes, a vortex is stably attached to the plate and does not detach until autorotation sets in. Similar limitations of quasi-steady models have been observed in studies involving free-falling plates (Andersen et al., 2005) and hovering insect ﬂight (Wang, 2005) where a similar strong interaction between the plate and its own wake is observed. In these scenarios, CFD–RBD models that directly simulate the complex non-linear interaction involved may offer the best approach for evaluating aerodynamic forces and the dynamic response of the plate. Fig. 9 shows phase averaged aerodynamic coefﬁcients of the ﬁnal three cycles of ﬂight in a plate with an initial angle of attack 301 which would be undergoing stable autorotation. The results are plotted against the effective angle of attack which takes into account the relative horizontal and vertical wind speed. Although comparable values of the average aerodynamic force and torque are obtained for all three models, QS1 predictions for the unsteady force coefﬁcients are different from those obtained from QS2 which gives better agreement with CFD–RBD predictions. This is mainly attributed to the inaccurate force decomposition used in QS1. In QS2, which uses Eqs. (26)–(29), is shown to offer better agreement with CFD for representation of the quasi-steady forces involved in plate autorotation and debris ﬂight. Clearly, this is to be expected because previous CFD simulations informed the coefﬁcients used in QS2, but it nonetheless demonstrates that the QS2 model coefﬁcients are better matched to the CFD. However, this improved form of the coefﬁcients still does not improve performance over QS1 in any of the ﬂight modes. 5. Conclusions The introduction of a robust rigid body model based on quaternions, coupled with CFD ﬂow simulations, has allowed the prediction of the free-ﬂight of plate-type debris in a manner not seen before. CFD–RBD predictions for debris ﬂight have been found to be in good agreement with available experimental observations for horizontal displacement and terminal velocity. Improvements to the lines of best ﬁt used in the prediction of the mean horizontal displacement have been presented, which do not rapidly degrade at longer ﬂight times. CFD–RBD predictions have been contrasted against predictions from a 2D quasi-steady analytical models. Although the 2D quasisteady analytical models are found to agree reasonably well with the CFD results for autorotational mode plates, they do not account for the strongly unsteady FSI involved in the launch stages and ﬂutter mode of ﬂight and hence perform poorly in these cases. In addition, the Tachikawa assumption about the ﬂuctuating component of aerodynamic forces, which is at the core of quasi-steady force models, is found to have some limitations and an alternative quasi-steady force model has been proposed based on the results of CFD forced rotation simulations. Part II (Kakimpa et al., this issue) looks in more detail at the ﬂight behaviour of the plates under a number of different launch conditions. This detailed analysis includes a deﬁnition of the various ﬂight modes discussed herein and the processes that cause over- and under-speeding of the plates. Acknowledgements The authors would like to express their thanks to project partners at the University of Birmingham: Prof. Chris Baker, Drs Mark Sterling, Andrew Quinn and Pedro Martinez–Vasquez. Also, the input of Prof. Peter Richards of the University of Auckland and Dr John Holmes has been invaluable. We also acknowledge the ﬁnancial support of the Engineering and Physical Sciences Research Council in the UK and the Dean of the Faculty of Engineering’s Research Fund at the University of Nottingham. References Andersen, A., Persavento, U., Wang, Z., 2005. Unsteady aerodynamics of ﬂuttering and tumbling plates. Journal of Fluid Mechanics 541, 65–90. Baker, C.J., 2007. The debris ﬂight equations. Journal of Wind Engineering and Industrial Aerodynamics 95, 329–353. Chou, J.C.K., 1992. Quaternion kinematic and dynamic differential equations. IEEE Transactions on Robotics and Automation 8-1, 53–64. Cline, M.B., Pai, D.K., 2003. Post-stabilization for rigid body simulation with contact and constraints. In: Proceedings of the IEEE International Conference on Robotics and Automation. Diebel, J., 2006. Representing Attitude: Euler Angles, Unit Quaternions, and Rotation Vectors. Technical Report, Stanford University, California. ESDU, September 1970. Fluid forces and moments on ﬂat plates. data item 70015. Engineering Science Data Unit, London, UK. Fluent Inc, 2009. FLUENT 12.1 Documentation. Greenwood, D.T., 2003. Advanced Dynamics. Cambridge University Press. Holmes, J., Letchford, C., Lin, N., 2006a. Investigations of plate-type windborne debris. Part II: computed trajectories. Journal of Wind Engineering and Industrial Aerodynamics 94, 21–39. Holmes, J.D., Baker, C.J., Tamura, Y., 2006b. Tachikawa number: a proposal. Journal of Wind Engineering and Industrial Aerodynamics 94, 41–47. Iversen, J.D., 1979. Autorotating ﬂat-plate wings: the effect of the moment of inertia, geometry and Reynolds number. Journal of Fluid Mechanics 92-2, 327–348. B. Kakimpa et al. / J. Wind Eng. Ind. Aerodyn. 111 (2012) 95–103 Jin, C., Xu, K., 2008. Numerical study of the unsteady aerodynamics of freely falling plates. Communications in Computational Physics 3-4, 834–851. Kader, B., 1981. Temperature and concentration proﬁles in fully turbulent boundary layers. International Journal of Heat and Mass Transfer 24-9, 1541–1544. Kakimpa, B., 2011. The Numerical Simulation of Plate-type Windborne Debris Flight. Ph.D. Thesis, Department of Civil Engineering, The University of Nottingham. Kakimpa, B., Hargreaves, D., Owen, J., 2010a. A singularity-free model for 3d windborne debris ﬂight. In: 9th UK Conference on Wind Engineering, University of Bristol, 20–22 September. Kakimpa, B., Hargreaves, D., Owen, J. An investigation of plate-type windborne debris ﬂight using coupled CFD–RBD models. Part II: free and constrained ﬂight, Journal of Wind Engineering Industrial & Aerodynamics, this issue. Kakimpa, B., Hargreaves, D.M., Owen, J.S., Martinez-Vazquez, P., Baker, C.J., Sterling, M., Quinn, A.D., 2010b. CFD modelling of free-ﬂight and auto-rotation of plate type debris. Journal of Wind and Structures 13–2, 169–189. Kordi, B., Kopp, G.A., 2009a. Evaluation of quasi-steady theory applied to windborne ﬂat plates in uniform ﬂow. Journal of Engineering Mechanics 135-7, 657–668. Kordi, B., Kopp, G.A., 2009b. ‘‘The debris ﬂight equations’’ by C.J. Baker. Journal of Wind Engineering and Industrial Aerodynamics 97 (3–4), 151–154. Kordi, B., Traczuk, G., Kopp, G., 2010. Effects of wind direction on the ﬂight trajectories of roof sheathing panels under high winds. Journal of Wind & Structures 13-2, 145–167. Lin, N., Letchford, C.W., Holmes, J.D., 2006. Investigation of plate-type windborne debris: part I. Experiments in wind tunnel and full scale. Journal of Wind Engineering and Industrial Aerodynamics 94, 51–76. Lin, N., Vanmarcke, E., 2010. Windborne debris risk analysis: part I. Introduction and methodology. Journal of Wind and Structures 13-2, 191–206. Marshall, T.P., Robinson, S., 2006. The Birmingham, U.K. Tornado: 28 July 2005. In: 23rd Conference on Severe Local Storms. 103 Martinez-Vazquez, P., Baker, C.J., Sterling, M., Quinn, A., Richards, P., 2010. Aerodynamic forces on ﬁxed and rotating plates. Journal of Wind and Structures 13-2, 127–144. Minor, J.E., Beason, W.L., 1976. Window glass failures in windstorms. ASCE Journal of the Structural Division 102-1, 147–160. Murman, S.M., Aftosmis, M.J., Rogers, S.E., 2005. Characterization of space shuttle ascent debris aerodynamics using CFD methods. In: 43rd AIAA Aerospace Sciences Meeting. NAHB Research Center, 2002. Wind-borne Debris: Impact Resistance of Residential Glazing. Technical Report, U.S. Department of Housing and Urban Development. Rhie, C., Chow, L., 1983. Numerical study of the turbulent ﬂow past an airfoil with trailing edge separation. AIAA Journal 21–11, 1525–1532. Richards, P.J., Williams, N., Laing, B., McCarty, M., Pond, M., 2008. Numerical calculation of the three-dimensional motion of wind-borne debris. Journal of Wind Engineering and Industrial Aerodynamics 96, 2188–2202. Robinson, A.C., December 1958. On the Use of Quaternions in Simulation of Rigidbody Motion. Technical Report, Aeronautical Research Laboratory, Wright Air Development Centre, Ohio, USA. Sarrate, J., Huerta, A., Donea, J., 2001. Arbitrary Lagrangian–Eulerian formulation for ﬂuid-rigid body interaction. Computer Methods in Applied Mechanics and Engineering 190 (24–25), 3171–3188. Shih, T., Liou, W., Shabbir, A., Yang, Z., Zhu, J., 1995. A new k-epsilon eddy-viscosity model for high Reynolds number turbulent ﬂows—model development and validation. Computers & Fluids 24–3, 227–238. Stuelpnagel, J., 1964. On the parametrization of the three-dimensional rotation group. Society for Industrial and Applied Mathematics 6 (4), 422–430. Tachikawa, M., 1983. Trajectories of ﬂat plates in uniform ﬂow with application to wind-generated missiles. Journal of Wind Engineering and Industrial Aerodynamics 14, 443–453. Wang, Z.J., 2005. Dissecting insect ﬂight. Annual Review of Fluid Mechanics 37, 183–210. Wills, J., Lee, B., Wyatt, T., 2002. A model of wind-borne debris damage. Journal of Wind Engineering and Industrial Aerodynamics 90, 555–565.