...

Kakimpa2012-WindborneDebris.pdf

by user

on
Category: Documents
2

views

Report

Comments

Transcript

Kakimpa2012-WindborneDebris.pdf
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 flight 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 fluid-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 flight modes
associated with the flight 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 fit. Agreement is found to
be good and, based on the findings, an alternative form for the dimensionless flight distance is
presented, which extends the range of the experimental study to longer flight times.
The predictions from the CFD–RBD model are then compared against two quasi-steady analytical
debris flight models. The second model is based on modified force and moment coefficients, which are
informed by the findings from the CFD–RBD model. For plates that have attained a stable, autorotational flight mode, the CFD–RBD and analytical models are in good agreement. Their predictions differ
during the initial stages of flight, 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 defined 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
classified 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) identified 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 definitive
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
flight. 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 coefficients, respectively. This formulation of the debris
flight equations showed that the debris flight 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 ¼ pffiffiffiffiffi ,
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 coefficients, C, of the form
C ¼ C f ðaeff Þ þ C R ðo Þ,
ð4Þ
where C f ðaeff Þ is the fluctuating 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 coefficients
for a static flat 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 flow (Andersen et al., 2005). This
strong plate-wake interaction is especially significant during the
debris launch and initial flight 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 modified quasi-steady torque used in
analytical models differs considerably from the aerodynamic
torque coefficients measured from autorotating flat plates. There
is therefore a need to review the validity of quasi-steady theory in
debris flight 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 flow field 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 findings, therefore, need to be extended to the
flight of low aspect ratio 3D plates in high Reynolds number
flows, 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
flight 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 flow around the shuttle during ascent, and the nonlinear interaction of the foam debris with the flow. The shuttle
foam debris flew for relatively short periods (approximately
0.15 s) compared to windborne debris flight, with limited rotational motion unless perturbed.
This paper presents an improved approach to the numerical
simulation of plate-type windborne debris flight by using a 3D
CFD model sequentially coupled with a six degree of freedom RBD
model to simulate the flight of low aspect ratio plates in wind
flow of Reynolds number 1:3 106 . The paper focuses on platetype debris – such as roofing 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 influence of plate initial orientation, Tachikawa number, mass
moment of inertia, geometric parameters such as thickness-ratio
and aspect-ratio, as well as complex launch flow fields. In the
present paper, Part I, the CFD–RBD predictions are compared
against experimental results where available, with new improved
fit 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 flat, 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-flight 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 fluid relative velocity vector, U is the fluid 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 fluid 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 inflow 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-flight 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 specified 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 findings are presented in
Kakimpa et al. (2010b). The grid independence study was conducted with static and fixed-axis rotating plates, which used the
same spherical inner region as used in the flying 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 difficult
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 fixed-axis
rotation, then this grid was suitable for the flying 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 fine 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 flight. 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 flights 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 fluctuation of body forces on the plate. The complete
flow field around the plate is simulated and coherent flow
structures may be identified (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 define 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-fixed 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 fixed 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-fixed 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 fixed reference frame,
X p Y p Z p , relative to the fixed 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 firstly 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 final 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 fixed 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-fixed 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-fixed reference frame.
A cursory examination of Eq. (14) reveals that if y ¼ 7 p=2,
_ and c
_ can be infinite for finite 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 flight, 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 efficient 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 benefits of
using quaternions are insignificant, 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 final Xp position, and this rotates
the displaced Yt and Zt axes onto their final Yp and Zp positions.
3.2. Review of Euler angle parametrizations
Previous studies of debris flight, 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-fixed
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
fixed angular velocity, xp , is computed from solving the Euler
equation of rotational motion, Eq. (10), in a body-fixed 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
qffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
ð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-fixed 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-fixed
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 efficiency 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 verification 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 fit 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-flight simulations were performed at a Tachikawa
number, K, of 8.3, which is typical for large roofing sheets in a
wind storm. In the first 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 flow 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 figure that there is a large
spread of trajectories, but that these clusters around the experimental fit 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
fit expressions for non-dimensionalised horizontal distance, Kx , and (b) CFD–RBD
predictions for non-dimensionalised horizontal speed, u, together with an experimentally derived fit expression (Lin et al., 2006).
%
Table 1
Polynomial coefficients for the fifth-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 fit 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 flight durations. To
overcome this difficulty, a fifth 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 coefficients 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 fit of Lin et al.
(2006) is shown
pffiffiffiffiffiffiffiffiffiffinffi
u 1e 1:8Kx and su ¼ 0:0814,
ð25Þ
Fig. 4. Nomenclature associated with the flying plates, showing the angle of
attack, az .
where su is the standard deviation of the experimental data. For
long flight 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 flies 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 figure. 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 flight times of t n r 0:8 (approximately 0.6 s),
while the CFD data is derived for plates with longer flight
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 coefficients into a static and an autorotational component according to
Tachikawa (1983). Recent experimental measurements of both
mean and fluctuating forces on a rotating plate (MartinezVazquez et al., 2010) have shown this decomposition to contain
inaccuracies in predicting the fluctuating component of aerodynamic force which Tachikawa (1983) did not originally measure
but rather assumed to be equivalent to the fluctuation of static
force coefficients 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 flat plates that were forced to rotate
at a specified angular speed, using the same spherical inner region
rotating inside a fixed outer domain as was used in the free-flight
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 coefficients
were extracted from the CFD simulations and the mean and rms
of each are plotted in Figs. 6 and 7. Also shown in the figures are
power law and exponential fits 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 identified, 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 fit expressions for the mean and rms force and
moment coefficients, 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 fits for the mean drag, CD, lift, CL, and moment, CM,
coefficients 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
flight 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 findings 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 fits for the rms drag, CD, lift, CL, and moment, CM,
coefficients against non-dimensionalised angular speed, o , for a forced
rotating plate.
%
according to
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
b D Þ2 ðC D Þ2 Þsinð2ap=2Þ,
C D ¼ C D þ 2ððC
C L ¼ k1 C L þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
b L Þ2 ðC L Þ2 Þsinð2aÞ,
2ððC
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
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
coefficients over a complete rotational cycle, calculated using the
fit expressions shown in Figs. 6 and 7.
4.3. Results
In the companion paper, Part II (Kakimpa et al., this issue),
three flight modes for the debris are identified. 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 flutter, transition
or autorotation flight modes. Plates that flutter 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 flutter 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 flutter,
transition and autorotation flight 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
flutter 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 flutter 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 flight,
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 flutter
(left), transitional (centre) and autorotational (right) flight 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 flight. 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 fixed-axis
autorotation and autorotating free flying plates in Section 2.2 of
Part II (Kakimpa et al., this issue).
Although the quasi-steady models based on fixed-axis autorotation theory are able to adequately represent the autorotational mode of flight, they are unable to capture the strongly
unsteady Fluid-Structure Interaction (FSI) involved in the flutter
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 flight (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 coefficients of the
final three cycles of flight 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 coefficients 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 flight. Clearly, this
is to be expected because previous CFD simulations informed the
coefficients used in QS2, but it nonetheless demonstrates that the
QS2 model coefficients are better matched to the CFD. However,
this improved form of the coefficients still does not improve
performance over QS1 in any of the flight modes.
5. Conclusions
The introduction of a robust rigid body model based on
quaternions, coupled with CFD flow simulations, has allowed
the prediction of the free-flight of plate-type debris in a manner
not seen before. CFD–RBD predictions for debris flight have been
found to be in good agreement with available experimental
observations for horizontal displacement and terminal velocity.
Improvements to the lines of best fit used in the prediction of the
mean horizontal displacement have been presented, which do not
rapidly degrade at longer flight 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 flutter mode of flight and hence perform poorly in
these cases. In addition, the Tachikawa assumption about the
fluctuating 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
flight behaviour of the plates under a number of different launch
conditions. This detailed analysis includes a definition of the
various flight 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
financial 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 fluttering
and tumbling plates. Journal of Fluid Mechanics 541, 65–90.
Baker, C.J., 2007. The debris flight 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 flat 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 flat-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 profiles 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 flight. 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 flight using coupled CFD–RBD models. Part II: free and constrained
flight, 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-flight 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 flat plates in uniform flow. Journal of Engineering Mechanics 135-7,
657–668.
Kordi, B., Kopp, G.A., 2009b. ‘‘The debris flight 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 flight
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 fixed 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 flow 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 fluid-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 flows—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 flat plates in uniform flow with application to
wind-generated missiles. Journal of Wind Engineering and Industrial Aerodynamics 14, 443–453.
Wang, Z.J., 2005. Dissecting insect flight. 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.
Fly UP