Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

56, 1, pp. 239-253, Warsaw 2018
DOI: 10.15632/jtam-pl.56.1.239

DYNAMIC SIMULATION OF A NOVEL “BROOMSTICK” HUMAN

FORWARD FALL MODEL AND FINITE ELEMENT ANALYSIS OF THE

RADIUS UNDER THE IMPACT FORCE DURING FALL

Dariusz Grzelczyk, Paweł Biesiacki, Jerzy Mrozowski, Jan Awrejcewicz

Lodz University of Technology, Department of Automation, Biomechanics and Mechatronics, Łódź, Poland

e-mail: dariusz.grzelczyk@p.lodz.pl; pawel.biesiacki@dokt.p.lodz.pl; jerzy.mrozowski@p.lodz.pl; jan.awrejcewicz@p.lodz.pl

The paper presents dynamic simulation and experimental identification of a human forward
fall model describing the process of “falling like a broomstick” on the outstretched arms.
Themodel implemented inMathematica allows one to estimate time histories of the ground
reaction force in different scenarios of the fall process. These time series are applied as
time-varying load conditions to the numerical analysis of the human radial bone model
created fromthe computed tomographydata.Finally, the obtainednumerical results indicate
that the strain criterion seems to be more useful for estimating the radius fracture site in
comparison to the stress criterion.

Keywords: forward fall, ground reaction force, bone strength, fracture, distal radius

1. Introduction

Falls are common accidents in human daily life. They are severe and inevitable threats during
the walking process of the bipedal movement. Despite all population groups are exposed to
these risk factors, falls are the most serious social health problem among the elderly (Heijnen
and Rietdyk, 2016). As a result, the risk of injuries to hands, torso, head/neck and/or other
parts of the human body is possible. These injuries can eventually lead to death, chronic pain,
disability and/or loss of independence (Robinovitch et al., 2013).

The vast majority of cases of upper extremity injuries occur as a result of a fall with a
direct impact on the extended arms which are exposed to dynamical impact forces (Nevitt and
Cummings, 1993; Palvanen et al., 2000). Distal radius fractures are common in eldery women
with osteoporosis due to their compromised bone density/quality as well as probably due to the
increased risk of falling in this population group.Fractures of the forearmbones represent nearly
20% of all reported fractures worldwide, and themost common type of fracture is the so-called
Colles’ fracture as an injury of distal radius of the forearm (Johnell and Kannis, 2006).

The aforementioned upper extremity injuries may be a result of forward falls, backward
falls or side falls. Colles’ fracture, as an injury of the radius, is a direct result of exceeding the
maximum value of the force allowable for this bone. For instance, distal radius fractures at a
mean force equal to 2260±1010Nwere observedbyFrykmanusing48 cadaveric bones of average
age of 65 years (Frykman, 1967). Spadaro and his co-investigators obtained the mean strength
of a bone fracture at the level of 1640± 980N (Spadaro et al., 1994). Kim and Ashton-Miller
(2009) adopted the value equal to 2400N as a distal radius fracture threshold. Also in one of
the recent papers byBurkhart et al. (2013), the estimated values of force causing fracture of the
studied bones derived from cadavers were approximately equal to 2150N.



240 D. Grzelczyk et al.

2. Literature review of the forward fall models

A forward fall is known as the most frequent type of fall, and more than a half of all falls
among the elderly occur in the forward direction (Nevitt and Cummings, 1993; O’Neill et al.,
1994; Vellas et al., 1998). Recent decades have brought different models related to the impact
of the upper extremities to the ground as a result of a forward or backward fall. One of the
simplest effective mechanical models of the human falling motion is a single-degree-of-freedom
(DoF) linear mass-spring-damper mechanical system subjected to a sudden velocity input or
an impulse force input. Such a model can be readily extended to systems with many DoFs.
Belowwe are focused only on the threemost well-known, widely available and cited forward fall
models.

A fall model proposed by Chiu and Robinovitch (1998) applies to the human forward fall
from a low height on the outstretched and fully extended upper extremities as the worst-case
scenario of such a fall. It is constructed based on a 2-DoFs lumped-parametermechanical system
containing elastic and damping elements responsible for the operation of humanmuscles. David-
son et al. (2006) completed a study with the aim of developing a “risk factor” value in order to
determine whether someone will sustain the radius fracture based on the characteristics of their
fall. Thementioned study used a twomass-spring-dampermodel, referring to the experimental
data obtained from 45 clinical cases of children falling off of playground equipment.

DeGoede and Ashton-Miller (2003) employed Adams software to devolop a half-body sym-
metric model of the human forward fall consisting of five segments (legs, torso with head and
neck, upper arm, forearm and hand). In that model, the movement of the hand and the pivot
point in the ankle toe are fixed, and friction is not taken into account.Moreover, thewrist, elbow
and hip joints are adopted as flat joints without friction, whereas the shoulder joint connecting
the arm to the body viamassless blades is modelled as a linear frictionless sliding joint with the
movement axis perpendicular to the longitudinal body axis. The main goal of the authors was
to study the possibility of injury in olderwomen.Therefore, to determine themodel parameters,
baseline height and weight equal to 1.63m and 62kg, respectively, were used (Kroemer et al.,
1997).

KimandAshton-Miller (2009) proposed another flatmodel of the forward fall as a twoDoFs
system constructed based on a mechanical model of a double pendulum rotating freely around
a pivot corresponding to the ankles of the lower human extremities. To provide amathematical
model, the mechanical system was reduced to a system of linear translational movements with
2-DoFs with spring-damper elements responsible for attenuation action of the human muscles.
Numerical simulation for the fall height of 1.5m showed that the maximum impact force was
doubled from1250Nto a value of 2610N, dependingon the fall scenario. In thisway, the authors
showed that for the same fall conditions (the same faller falling from the same height), too rapid
movement of the arms in comparison to the rest of the body could cause the occurence of an
excessive force acting on the upper extremity more than 2350N, which caused a fracture of the
distal radius.

To conclude, it can be stated that, in general, the considered falling process was usually
mathematically modelled based on flat and linear mechanical systems consisting of two rigid
bodies with masses moved by transverse motion connected by linear spring-damper elements.
Thesemodels were usually presented as a set of the second-order ordinary differential equations
of motion (Chiu andRobinovitch, 1998), but also the appropriate equations were written in the
state space (Kim and Ashton-Miller, 2009). On the other hand, in the case of a more complex
mechanicalmodel, numerical simulationswereperformedonlyusing commercialAdams software
(DeGoede and Ashton-Miller, 2003). In the present paper, we propose our own mathematical
model of the human forward fall on the outstretched arms. In addition, unlike previous models
met in the literature, our model takes into account a direct influence of different human speed



Dynamic simulation of a novel “broomstick” human forward fall model... 241

just before the trip over an obstacle and starting the falling process, which affects the final force
of impact of the upper extremity on the ground.

3. The proposed “broomstick” forward fall model

3.1. Mathematical modelling

Thehuman “like a broomstick” forward fall on the outstretched arms is schematically shown
in Fig. 1a. Equations ofmotion of the analysed systemhave been obtained by theNewton-Euler
method, and Free Body Diagrams of the system are shown in Fig. 1b.

Fig. 1. Mechanical modelling of the forward fall process: (a) the proposed forward fall biomechanical
model as a planarmechanical systemwith 2-DoFs embedded in the Cartesian coordinate system;

(b) Free Body Diagrams of the system

Themodel consists of two rigid bodies, 1 and 2, withmassesm1,m2 andmoments of inertia
I1, I2 around centres of gravity of the bodies, respectively. Body 1 corresponds to the human
torso with legs, neck and head, whereas body 2 corresponds to the human upper extremities
(including arms, forearms andhands).Theparametersa1 anda2 represent thedistances between
the gravity centres of individual bodies and their rotation axes, whereas l1 and l2 denote the
total lengths of thementionedbodies, respectively.Theangle θ1(t) denotes theangle between the
horizontal x axis and the longitudinal axis of body 1.The angle θ2(t) is the anglemeasured from
the axis of body 1 to the axis of body 2. The groundwith non-linear contact law is characterised
by stiffness and damping coefficients ky and by, respectively.

In this model, we take the vectors θ1(t)= [0,0,θ1(t)]
T, θ2(t)= [0,0,θ2(t)]

T of the angles in
the joints j1 and j2, vectors rC1(t), rC2(t) of displacements of gravity centres of the first and
the second body, as well as the vectors l1(t), l2(t) of displacements of the joint j2 and the end
of body 2 which can impact to the ground, respectively, in the following form

rC1(t)= [x1(t),y1(t),0]
T = [a1cosθ1(t),a1 sinθ1(t),0]

T

rC2(t)= [x2(t),y2(t),0]
T = [l1cosθ1(t)+a2cosα(t), l1 sinθ1(t)−a2 sinα(t),0]

T

l1(t)= [l1cosθ1(t), l1 sinθ1(t),0]
T

l2(t)= [l1cosθ1(t)+ l2cosα(t), l1 sinθ1(t)− l2 sinα(t),0]
T

(3.1)

whereα(t) = 180◦−θ1(t)−θ2(t). The forcesQ1 = [0,−m1g,0]
T andQ2 = [0,−m2g,0]

T are the
gravity forces acting on the gravity centres of bodies 1 and 2, respectively, with the acceleration



242 D. Grzelczyk et al.

of gravity g (g = 9.81m/s2). The force R(t) = [Rx(t),Ry(t),0]
T is the reaction force in the

joint j1. The unknown force in the joint j2 and presented in the Free Body Diagrams (see
Fig. 1b) is denoted as P(t) = [Px(t),Py(t),0]

T. Finally, the force F(t)= [Fx(t),Fy(t),0]
T is the

ground reaction force acting on body 2 at the moment of its impact to the ground.
Let M1(t) = [0,0,0]

T and M2(t) = [0,0,M2(t)]
T denote the torques generated in the joints

j1 and j2, respectively. As a result, for the two considered free bodies presented in Fig. 1b, we
can write the following equations of motion in the vector form

m1r̈C1(t)=R(t)+Q1+P(t)

I1θ̈1(t)=M1(t)−M2(t)−rC1(t)×R(t)+ [l1(t)−rC1(t)]×P(t)

m2r̈C2(t)=−P(t)+Q2+F(t)

I2θ̈2(t)=M2(t)− [l1(t)−rC2(t)]×P(t)+ [l2(t)−rC2(t)]×F(t)

(3.2)

Equations (3.2) can be reduced to the scalar form

I1θ̈1(t)=−M2(t)+a1Rx(t)sinθ1(t)−a1Ry(t)cosθ1(t)

+(l1−a1)Py(t)cosθ1(t)− (l1−a1)Px(t)sinθ1(t)

I2θ̈2(t)=M2(t)+a2Py(t)cosα(t)+a2Px(t)sinα(t)

+(l2−a2)Fy(t)cosα(t)+(l2−a2)Fx(t)sinα(t)

(3.3)

where

Px(t)=Fx(t)−m2ẍ2(t)

Py(t)=Fy(t)−m2ÿ2(t)−m2g

Rx(t)=m1ẍ1(t)−Px(t)=m1ẍ1(t)+m2ẍ2(t)−Fx(t)

Ry(t)=m1ÿ1(t)+m1g−Py(t)=m1ÿ1(t)+m2ÿ2(t)+m1g+m2g−Fy(t)

(3.4)

At the time of tripping over an obstacle, a human instinctively and quickly pulls his arms
to the front in order to arrest and/or to absorb the fall. This process is described by the time
histories of the angle θ2(t), which can be estimated based on the kinematics of the falling
process registered with the use of a digital camera. Based on the second equation of (3.3) and
the function θ2(t) obtained in this way, we can calculate the torqueM2(t) generated by arms in
the joint shoulder during the falling process. Taking the torque M2(t) in the first Eq. of (3.3),
one obtains eventually the following equation of motion around the joint j1

I1θ̈1(t)+ I2θ̈2(t)= a1Rx(t)sinθ1(t)−a1Ry(t)cosθ1(t)+(l1−a1)Py(t)cosθ1(t)

− (l1−a1)Px(t)sinθ1(t)+a2Py(t)cosα(t)+a2Px(t)sinα(t)

+(l2−a2)Fy(t)cosα(t)+(l2−a2)Fx(t)sinα(t)

(3.5)

with the function θ1(t) as a solution to this equation of motion.

3.2. Ground reaction force

To predict the vertical ground reaction force (GRF) Fy(t), a non-linear model of impact at
the wrist-ground interface has been employed in the form

Fy(t)= ky|y(t)|
3(1− byẏ(t))J(−y(t)) (3.6)

where parameters ky, by denote the ground stiffness and damping coefficients in the vertical
direction, respectively, y(t) = l1 sinθ1(t)− l2 sinα(t), and the function J(−y(t)) is the step
function defined as

J(−y(t)) =

{

1 if y(t)< 0

0 if y(t)­ 0
(3.7)



Dynamic simulation of a novel “broomstick” human forward fall model... 243

For instance, this formulation of GRF had been used previously to model the heel strike in
running (Gerritsen et al., 1995) and tomodel GRF at the hand-ground interface (DeGoede and
Ashton-Miller, 2003). In our model, we have not included the horizontal force component Fx(t)
which is required to prevent the hand from sliding. However, the total GRF is less than 5%
greater than the vertical GRF (DeGoede and Ashton-Miller, 2003).

3.3. Initial conditions

The considered model has been reduced to a single DoF model with the function θ1(t) as
the solution and the following initial conditions: θ1(0)-initial angular position, and θ̇1(0)-initial
angular velocity. At thebeginningof the trip (t=0), ahuman is usually in the standingposition,
and thus we take θ1(0)= 90

◦ (see Fig. 1a). The initial angular velocity θ̇1(0) is estimated based
on the transverse human speed during walking just before the moment of the trip. For this
reason, we assume that the linear speed of the human gait decreases from v0 to 0, and the
rotational velocity θ̇1(t) increases from 0 to θ̇1(0) at the same time. Taking into account the
principle of conservation ofmomentum, rotarymotion of the humanbody around the pivot axis
is governed as

I
∆θ̇1(t)

∆t
= rF ⇒ ∆θ̇1(t)=−

(m1+m2)v0r

I
= θ̇1(0) (3.8)

where F∆t is a force impulse causing rotation of the human body around the rotation axis
placed in feet (ankles), I is themoment of inertia of the human body around the pivot axis and
r is the distance between the humangravity centre and the pivot axiswith the upper extremities
adjusted along the body (typical position of the human body during walking). On the contrary
to the previous investigations met in the literature, the assumption adopted in this paper allows
one also to study kinematic and dynamic parameters during the falling process for different
walking speeds of the faller just before the moment of the trip over an obstacle.

4. Identification of the fall model parameters

4.1. GrabCAD model of the human body

To determine the human body parameters required in the proposed fall model, the full 3D
scanned human body model based on GrabCAD (2016) (see Fig. 2) as well as the appropria-
te AutoDesk Inventor commands have been used. The mass of the whole body is calculated
assuming that the average density equals to ρ = 1050kg/m3, and the model parameters are:
a1 = 1.01m, a2 = 0.22m, l1 = 1.4m, l2 = 0.53m, r = 1.03m, m1 = 61.0kg, m2 = 7.4kg,
I1 =9.9kg·m

2, I2 =0.232kg·m
2, I =82.86kg·m2.

4.2. Kinematics of the fall process

The kinematic analysis of the faller between the moment of tripping over an obstacle and
hitting his hands to the groundwere observed using anOptitrack optoelectronicmotion analysis
system installed at theDepartment of Automation, Biomechanics andMechatronics, LodzUni-
versity of Technology, Lodz, Poland. The front panel of the computer program used to operate
this systemwith the location of the individualmarkers on the faller’s body is presented inFig. 3,
where the process of forward fall from the standing position to the ground is shown. Before the
experimenal test, the volunteer (one of the authors of this paper) has been instructed to safely
fall “like a broomstick”. The fall of the human body was performed using a certificated soft
insurance mattress with the full size 2m×2m×0.3m, and thus there was no a direct threat of
injury to the upper extremities and other body parts.



244 D. Grzelczyk et al.

Fig. 2. Three-dimensional scanned human body model used for the estimation of its geometrical and
mechanical properties (GrabCAD, 2016)

Fig. 3. Forward fall from the standing position to the soft mattress observed by Optitrack system using
37 passive reflective markers distributed on the faller’s body

Using Optitrack system we observed that the faller extends his arms from initial position
θ20 ≈ 5

◦ to the final position θ2max ≈ 80
◦ at the moment of the impact to the ground. That is

why, in all further presented investigations we use the time history of angle θ2(t) described by
the formula

θ2(t)=

{

θ20+(θ2max−θ20)sin
2(λt) if t¬T

θ2max if t>T
(4.1)

whereT denotes duration of the fall (the timebetween tripping andhitting the ground),whereas
λ corresponds to the speed of movements of the faller’s arms. The parameters λ and T strongly
depend on the walking speed of the faller before the trip. The averaged time history of the
angle θ2(t) obtained from the experiment for v0 =1.5m/s and its approximation are presented
in Fig. 4.

4.3. Hand-ground contact parameters

Figure 5a shows the schematics (the first frame of our animation created inMathematica) of
the initial configuration of the subject falling in the forward direction from the initial shoulder
height of 1m (DeGoede and Ashton-Miller, 2003). The red line represents the torso with legs,
the blue line represents upper extremities, whereas the green circle is the head of the faller.



Dynamic simulation of a novel “broomstick” human forward fall model... 245

Fig. 4. Time history of the angle θ2(t) (red points and red curve) obtained from the experiment and its
approximation by an analytical smooth function (green curve). In this experiment, the faller started to

fall from not a fully standing position (θ1(0)< 90
◦)

The presented initial position was obtained for the following initial conditions: θ1(0) = 45
◦,

θ̇1(0) = 0 and θ2(t) = 83.5
◦ = const. During numerical experiments we tested different values

of the parameters ky and by to obtain GRF which corresponds (both from a qualitative and
quantitative point of view) to the GRFs presented by DeGoede and Ashton-Miller (2003). In
that paper, the authors tested five healthy youngmale volunteers aged between 22 and 28 years
with the average body mass of 72± 7kg and the overall height of 173± 3cm (DeGoede and
Ashton-Miller, 2002). Finally, the best degree of fit was obtained for ky = 50000N/m

3 and
by =0.6s/m (the obtained results are shown in Fig. 5b).

Fig. 5. Initial configuration of the faller’s body, i.e. the fall from the shoulder height of 1m (a), and the
time history of GRF obtained numerically for the proposed forward fall model for ky =50000N/m

3

and by =0.6s/m, which corresponds to the GRF obtained experimentally for a representative fall from
the shoulder height of 1m (b) (DeGoede and Ashton-Miller, 2003)

5. Numerical simulation of the forward fall model

5.1. Parameters and relationships used in simulation

Experimental observations confirmed that with an increase in the walking speed v0, the
human instinctively faster pulls his arms in the forward direction and the value of thementioned
parameterλ also increases. This is why in all our numerical simulationswe use the time histories



246 D. Grzelczyk et al.

of the angle θ2(t) governed by Eq. (4.1) with different parameters θ2max and λ, depending on
the angle φArm (the angle between the arm and the vertical axis of the Cartesian coordinate
system at themoment of the impact to the ground) and the speed v0 of the humanwalking just
before the trip. Table 1 presents the dependence of the angle θ2max on the angle φArm, while
Table 2 presents the dependence of parameters λ and T on the speed v0 of the faller before the
trip. The parameters of the proposed fall model (identified in Section 4) allow one to conduct
broader analysis for different parameters and different fall scenarios.

Table 1.Values of the angle θ2max corresponding to different values of the angle φArm

φArm [
◦] 0 5 10 15 20 25 30

θ2max [
◦] 67.75 72.85 78.11 83.55 89.16 94.93 100.86

Table 2.Values of parameters T and λ corresponding to different values of walking speed v0

v0 [m/s] 0.5 1.0 1.5 2.0 2.5 3.0

T [s] 1.0 0.7 0.58 0.5 0.4 0.35

λ [1/s] 1.58 2.23 2.70 3.15 3.92 4.50

5.2. Time histories of ground reaction forces

Figure 6 shows the time history of GRF obtained for a fall from the full standing position
(θ1(0) = 90

◦), v0 =1.5m/s and φArm =15
◦, which corresponds to the most common falls met

in real situations. In this case, the maximum value of GRF reaches 2400N, which corresponds
to the distal radius fracture threshold (Kim and Ashton-Miller, 2009). Thus, we demonstrated
that the taken typical and themost common falling conditions are enough to cause compressive
fracture of the distal radius. The investigated falling process has been also observed by using
the animation created inMathematica (see Fig. 7).

Fig. 6. Time histories of GRF from the standing position and v0 =1.5m/s

Fig. 7. Animation snapshots of the faller’s body plotted in regular time intervals in different phases of
the fall from the full standing position and walking speed v0 =1.5m/s (obtained in

Mathematica software)



Dynamic simulation of a novel “broomstick” human forward fall model... 247

Figure 8 illustrates times histories of the forceFy(t) acting on a single hand for v0 =1.5m/s
and for different values of the angleφArm. For all threepresented cases, the plotted timehistories
of GRF are similar. However, it can be noticed that themaximumvalue of GRF increased from
2374N for φArm =0 to 2486N for φArm =30

◦.

Fig. 8. GRFs for v0 =1.5m/s and different values of angle φArm: (a) time histories Fy(t),
(b) maximum values Fymax

The results in Fig. 9 show the influence of different values of the velocity v0 on time histories
of the forceFy(t). For larger values of the parameter v0, the duration of the fall is smaller while
themaximumvalue ofGRF is greater. The value ofGRF increases from2246N for v0 =0.5m/s
to 2534N for v0 =2.0m/s. It means that for a smaller speed v0 (i.e. v0 < 1.5m/s), the GRF is
less than the distal radius fracture threshold, whereas for larger values of v0 (i.e. v0 > 1.5m/s),
this threshold is exceeded.

Fig. 9. GRFs for φArm =15
◦ and different values of the velocity v0: (a) time histories Fy(t),
(b) maximum values Fymax

Figure 10 presents maximum values of GRF as a function of the velocity v0 for different
values of the angle φArm. The maximum value of GRF increases with the increasing speed v0
as well as the angle φArm. For the lowest presented value of v0 equal to 0.5m/s, the maximum
value of GRF increases from 2217N for φArm = 0 to 2335N for φArm = 30

◦. For the largest
value of v0 equal to 3.0m/s, the maximum value of GRF increases from 2898N for φArm = 0
to 3000N for φArm = 30

◦. If one considers the case φArm = 0, the maximum value of GRF
increases from 2217N for v0 = 0.5m/s to 2898N for v0 = 3.0m/s. In the case of φArm = 30

◦,
the maximum value of GRF increases from 2335N for v0 =0.5m/s to 3000N for v0 =3.0m/s.



248 D. Grzelczyk et al.

For small values of v0 (v0 ¬ 1m/s), the maximum value of GRF is less than the distal radius
fracture threshold, regardless of the angleφArm. In the case of v0 =1.5m/s, themaximumvalue
of GRF is less than the mentioned fracture threshold when 0 ¬ φArm ¬ 10

◦, whereas GRF
exceeds this threshold for φArm > 15

◦. For v0 ­ 2.0m/s, the maximum of GRF exceeds the
distal radius fracture threshold regardless of the angle φArm. Concluding, for a small walking
speed before the falling process, the value of the angle φArm does not have a great influence.
At a higher speed of human gait, the distal radius fracture threshold is usually exceeded, finally
leading to injuries and/or fractures of the upper extremities. In the case of v0 = 1.5m/s, the
correct configuration of the human body (the correct value of the angle φArm) at themoment of
the impact to the ground can reduce the maximum value of GRF and, as a result, can protect
the faller from potential injuries or fractures of the upper extremity bones.

Fig. 10. Maximum values Fymax of GRF as a function of the velocity v0 for different values of
the angle φArm

6. Finite element analysis of the radial bone

6.1. Geometry and material properties of the radial bone model

TheDICOMdata used in this paper come fromcadavers of 35-years-oldmanwith aheight of
1.73m and weight of 75kg. These data have been obtained using a Siemens 64 Slice computed
tomography (CT) Scanner in the Department of Forensic Medicine, Jagiellonian University
Medical College, Cracow, Poland. The DICOM file of all upper extremity bones composed of
the total number of 274 slices with the slice thickness equal to 1.5mm, pixel size equal to
0.977mm and resolution 512×512 was imported to Mimics. During the next steps, the radius
was separated, the computer model of this bone was obtained, and the computational mesh of
the radius was corrected to avoid further numerical errors. As a result, the mesh of the radius
was reduced to 3444 surface elements, and the correctness of the constructed computational
mesh was verified using a Fix Wizard function. Finally, a realistic 3D FE model of the radius
consisting of 15751 FEs was obtained and used for the strength analysis in Ansys.We used the
SOLID185FE-shaped tetrahedron element.Theassumed in this paper isotropy of thematerial is
not ideal, but theCTdata provided only scalar information, and the determination of principal
material directions had to be inferred (Neuert et al., 2013).

Material inhomogeneity of the radiuswasmodelled inMimics based on theCT images. First,
we tested three different density-elasticity relationships describing human long bones to obtain
the correct range of Young’s modulus of the radius based on the CT images and the Hounsfield



Dynamic simulation of a novel “broomstick” human forward fall model... 249

Unit scale. Finally, mechanical properties of the considered radial bone were calculated based
on the density-elasticity relationships proposed by Rho et al. (1995)

ρ=1.067HU+131 E=0.004ρ2.01 (6.1)

whereρ expressed in [kg/m3] representsdensityof thebone,HUis thenondimensionalHounsfield
Unit,E expressed in [MPa] denotesYoung’smodulus,whereas ν is Poisson’s ratio. The obtained
rangeofYoung’smodulusof the radial bone is close to the rangeof the appropriatevalues for this
bone presented in the literature. Young’s moduli of the radius vary in the range 600-12000MPa
(see Fig. 11), whereas Poisson’s ratio equals 0.3 for all FEs.

Fig. 11. Histogramand cross-sections of the radius showing quantitive distributions of FEs characterised
by different Hounsfield Units and spatial distribution of inhomogeneousmaterial properties

6.2. Load and boundary conditions

Numerical results presented inFig. 6were applied as the load conditions of the developedFE
upper extremitymodel analysed inAnsys. These time histories of GRFs correspond to theGRF
for the standing position and thewalking speed v0 =1.5m/s. As far as the boundary conditions
applied in the analysed numerical FE model of the radius were concerned, all six spatial DoFs
in the proximal radius (region of the elbow joint)were fixed.The angle between the longitudinal
axis of the radius and axis of the gravity field was equal to φArm, while the GRF was applied
in the region of the radial neck in the vertical direction (direction of the gravity field). In all
the investigated cases, we considered the same time-varying load condition (see Fig. 6), but the
analysis was conducted for different directions of theGRF acting on the radius (i.e. for different
angles φArm).

6.3. Stress and strain analyses

To predict bone fracture sites, we decided to measure the maximum von Mises stresses as
the main criterion first. Moreover, we also measured maximal strains of the bones as another
criterion leading to the additional information useful for a better understanding of the upper
extremity behaviour under the applied load conditions, i.e. the assessment of the radius fracture
site and failure load. The determination of themoment of timewhen themaximum stresses and
strains were the highest enabled finding the fracture site and to investigate the state of the rest
of the bones. The mentioned results are presented in Figs. 12 and 13 in the form of von Mises
stress and strain distributions.

Values of both the maximal von Mises stress and bone strains significantly depend on the
value of the angle φArm. Namely, with an- increase in the angle φArm, von Mises stresses and
strains also increase. It should be noted that in the case of axial compression of the bone
(φArm =0), there are no bending torques and the bone itself is more resistant in such a configu-
ration in comparison to the nonaxial compression. In the case of nonaxialGRF, bending torques
are significant and, therefore, there are considerable values of both stresses and strains. For all
the considered cases, the maximum values of stress and strain occur for φArm =30

◦.



250 D. Grzelczyk et al.

Fig. 12. Stress distributions in the radius for different angles φArm

Fig. 13. Strain distributions in the radius for different angles φArm

The results presented in Figs. 12 and 13 show also other regularity. Namely, the maximum
stresses occur on themedial side (diaphysis) of the radius bonewhile themaximumstrains occur
in the distal region of this bone. As it is known, Colles’ fracture is the most common type of
injury related to the forward fall on the outstretched upper extremities. Therefore, the presented
results indicate that the strain criterion can be more useful for estimating the radius fracture
site (the maximum strains are concentrated in the distal radius, see Fig. 13).

Table 3.Maximum von Mises stress andmaximum strains for three different angles φArm

Radius Time [s]

φArm =0
von Mises stress [MPa] 65.54 0.122

Strain [–] 0.0154 0.126

φArm =15
◦
von Mises stress [MPa] 78.88 0.132

Strain [–] 0.0186 0.132

φArm =30
◦
von Mises stress [MPa] 148.58 0.124

Strain [–] 0.0363 0.124

In contrast topreviousnumerouspapersdealingusuallywith resting conditions inour studies
we carried out also a transient state analysis in Ansys. Values of maximum stresses and strains
of the radius bone for different angles φArm are presented in Table 3. The maximum values of



Dynamic simulation of a novel “broomstick” human forward fall model... 251

the vonMises stress and strains do not occur at themoment of themaximumpeak of GRF but
later, i.e. at the time about 0.12-0.13s.

7. Conclusions

The forward fall model proposed in this paper enables one to estimate the vertical ground reac-
tion forces acting on the hands in various scenarios of the human falling process. The obtained
numerical simulations fit with other results presented in the literature (Kim andAshton-Miller,
2009), both from a qualitative and quantitative point of view. Moreover, the simulations show
that the parameters describing the human body and parameters modelling biomechanical pro-
perties between the palmar cartilages and the groundhave an essential influence on the obtained
results.

It should be emphasized that the developed model has some limitations. First of all, the
movement of the shoulder grid with respect of the torso and stiffness/damping properties of the
shoulder joint have not been implemented in our model. Moreover, the horizontal component
of the ground reaction force has not been considered. Nevertheless, the mentioned limitations
may be of interest for our future study. On the other hand, our further modifications and
improvements of the proposed model can be oriented on the increase of the number of DoFs of
themechanical model describing the human body aswell as taking into consideration the rotary
stiffness and damping properties in each of the human joints.

The choice of a linear material law is justified by small displacements of the radius bone
observed during compressive tests presented in the literature (Bosisio et al., 2007). Although
all bones have been modelled as a linear elastic isotropic material, in our opinion, the applied
simplification should not significantly affect the ability of the proposedmodel to predict fracture
sites and failure load of the radius boneunder loads resulting froma forward fall. The variability
of Poisson’s ratio has not been evaluated in our investigations.

The performed numerical investigations with the time history of theGRF that occur during
the falling process in the forward direction on the outstretched arms are sufficient to determi-
ne potential fracture sites and the obtained results agree with numerical/experimental results
presented in the literature (Edwards and Troy, 2012). In addition, the obtained results indicate
that the angle φArm and, consequently, the direction of load applied to the radius have a strong
impact on the fracture strength of this bone. It means that falls from a standing position on the
outstretched arms generate the value of GRF which can exceed the mean human distal radius
fracture threshold.Moreover, we have also shown that themaximal strain criterion seems to be
more useful for the estimation of the fracture site than the appropriate von Mises stress crite-
rion. Althoughwe obtained various numerical results, unfortunately wewere unable to compare
these results with own experimental studies from a quantitative point of view. However, the
obtained numerical results show that our model provides a realistic estimation of radius bone
strain, fracture strength and fracture side estimation under various loading scenarios simulating
a forward fall.

Ethical Approval

This article does not contain any studies performed on animals. The presented experimental studies
have been performed using one of the author of this paper (Paweł Biesiacki) without any other human
participants.

Acknowledgments

The work has been partially supported by the National Science Centre of Poland under the grant

OPUS 9 No. 2015/17/B/ST8/01700 for years 2016-2018.



252 D. Grzelczyk et al.

References

1. Bosisio M.R, Talmant M., Skalli W., Laugier P., Mitton D., 2007, Apparent Young’s
modulus of human radius using inverse finite-element method, Journal of Biomechanics, 40,
2022-2028

2. Burkhart T.A., Andrews D.M., Dunning C.E., 2013, Multivariate injury risk criteria and
injury probability scores for fractures to the distal radius, Journal of Biomechanics, 46, 973-978

3. Chiu J., Robinovitch S.N., 1998, Prediction of upper extremity impact forces during falls on
the outstretched hand, Journal of Biomechanics, 31, 1169-1176

4. DavidsonP.L., ChalmersD.J., Stephenson S.C., 2006,Prediction of distal radius fracture in
children, using biomechanical impactmodel and case-control data on playground free falls, Journal
of Biomechanics, 39, 503-509

5. DeGoedeK.M.,Ashton-Miller J.A., 2002,Fall arrest strategy affects peak hand impact force
in a forward fall, Journal of Biomechanics, 35, 843-848

6. DeGoede K.M., Ashton-Miller J.A., 2003, Biomechanical simulations of forward fall arrests:
effects of upper extremity arrest strategy, gender and aging-related declines in muscle strength,
Journal of Biomechanics, 36, 413-420

7. Edwards W.B., Troy K.L., 2012, Finite element prediction of surface strain and fracture
strength at the distal radius,Medical Engineering and Physics, 34, 290-298

8. FrykmanG., 1967,Fractureof thedistal radius including sequelae shoulder-hand-finger syndrome,
disturbance in the distal radio-ulnar joint and impairment of nerve function, Acta Orthopaedica
Scandinavica, 108, 1-153

9. Gerritsen K.G.M., van den Bogert A.J., Nigg B.M., 1995, Direct dynamics simulation of
the impact phase in heel-toe running, Journal of Biomechanics, 28, 661-668

10. GrabCAD, 2016, Open CAD library (https://grabcad.com/library)

11. Heijnen M.J.H., Rietdyk S., 2016, Falls in young adults: Perceived causes and environmental
factors assessed with a daily online survey,Human Movement Science, 46, 86-95

12. JohnellO.,Kannis J.A., 2006,Anestimateof theworldwideprevalenceanddisabilityassociated
with osteoporotic fractures,Osteoporosis International, 17, 1726-1733

13. Kim K.-J., Ashton-Miller J.A., 2009, Segmental dynamics of forward fall arrests: A system
identification approach,Clinical Biomechanics, 24, 348-354

14. Kroemer K.H.E.,KroemerH.J., Kroemer-Elbert K.E., 1997,Engineering Physiology: Ba-
ses of Human Factors/Ergonomics, Van Nostrand Reinhold, NewYork

15. Neuert M.A.C., Austman R.L., Dunning C.E., 2013, The comparison of density-elastic mo-
dulus equations for the distal ulna at multiple forearm positions: a finite element study, Acta of
Bioengineering and Biomechanics, 15, 37-43

16. Nevitt M.C., Cummings S.R., 1993, Type of fall and risk of hip and wrist fractures: the study
of osteoporotic fractures, Journal of the American Geriatrics Society, 41, 11, 1226-1234

17. O’Neill T.W., Varlow J., Silman A.J., Reeve J., Reid D.M., Todd C., Woolf A.D.,
1994,Age and sex influences on fall characteristics,Annals of the Rheumatic Diseases, 53, 773-775

18. Palvanen M., Kannus P., Parkkari J., Pitkajarvi T., Pasanen M., Vuori I., Jarvinen
M., 2000, The injurymechanisms of osteoporotic upper extremity fractures among older adults: a
controlled study of 287 consecutive patients and their 108 controls,Osteoporosis International, 11,
822-831

19. Rho J.Y., Hobatho M.C., Ashman R.B., 1995, Relations of mechanical properties to density
and CT number in human bone,Medical Engineering and Physics, 17, 347-355



Dynamic simulation of a novel “broomstick” human forward fall model... 253

20. Robinovitch S.N., Feldman F., Yang Y., Schonnop R., Leung P.M., Sarraf T., Sims-
GouldJ., LoughinM., 2013,Video capture of the circumstancesof falls in elderlypeople residing
in long-term care: an observational study,The Lancet, 381, 47-54

21. Spadaro J.A., Werner F.W., Brenner R.A., Fortino M.D., Fay L.A., Edwards W.T.,
1994, Cortical and trabecular bone contribute strength to the osteopenic distal radius, Journal of
Orthopaedic Research, 12, 211-218

22. Vellas B.J., Wayne S.J., Garry P.J., Baumgartner R.N., 1998, A two-year longitudinal
study of falls in 482 community-dwelling elderly adults, The Journals of Gerontology, Series A,
Biological Sciences and Medical Sciences, 53, 264-274

Manuscript received June 7, 2017; accepted for print September 21, 2017