JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

43, 4, pp. 813-824, Warsaw 2005

AN ALTERNATIVE SCHEME FOR DETERMINATION OF

JOINT REACTION FORCES IN HUMAN MULTIBODY

MODELS

Wojciech Blajer

Institute of Applied Mechanics, Technical University of Radom

e-mail: wblajer@poczta.onet.pl

Adam Czaplicki

Department of Biomechanics, External Faculty of Physical Education in Biała Podlaska,

Academy of Physical Education in Warsaw

e-mail: czaplicki@poczta.onet.pl

Multibodymodels are commonly used in the analysis of humanmovements.
The dynamic formulations often useminimal sets of generalized coordinates,
and joint reactions (non-working reactions ofmodel-intrinsic constraints) are
excluded fromevidence.A separatemodeling effort is then required to deter-
mine joint reactions, and the arising numerical procedures are computatio-
nally arduous. In this paper, a novel efficient approach to the determination
of joint reactions is developed, which naturally assists theminimal-form for-
mulations of human body dynamics. The proposed scheme does not involve
matrix inversion, and as such it is well suited for both symbolic manipu-
lations and computer implementations. The method is illustrated with a
seven-segment planarmodel of a human body. Some results from the inverse
dynamics simulation of somersaults on a trampoline are reported.

Key words: human bodymodeling, dynamic analysis, joint reaction forces

1. Introduction

Joint reaction forces play an important role in the dynamics of humanmo-
vements such as walking, jumping and gymnastic exercises; see e.g. Bergman
et al. (2001) andMcNitt-Gray et al. (2001). The question of which loads cross
the joints during routine or sport activities may be of considerable interest
to the clinicians (Wismans et al., 1994). One way to get the information is



814 W.Blajer, A.Czaplicki

to build a mathematical model of the human body dynamics and use it for
the determination of joint reaction forces by means of numerical simulation.
However, dynamical formulations of the humanbody usually useminimal sets
of generalized coordinates, and the joint reactions (non-working reactions of
model-intrinsic constraints) are excluded fromevidence (Eberhard et al., 1999;
Tözeren, 2000; Pandy, 2001; Yamaguchi, 2001). In order to obtain the joint
reactions by using classical multibody codes, described e.g. by Langer et al.
(1987), Nikravesh (1988), Garćıa de Jalón and Bayo (1993), Schiehlen (1997)
and Blajer (2001), a separate modeling effort is required, and the arising nu-
merical procedures may be computationally arduous.

In this paper, a novel approach to the determination of joint reactions, na-
turally assisted with the minimal-form formulation of human body dynamics,
is described. The idea of the scheme is similar to that of Kane and Levinson
(1985), called by them ”bringing noncontributing forces into evidence”. This
extension of Kane’s method, in which the joint reaction forces are originally
eliminated early in the process of deriving dynamic equations, has then been
exploited e.g. by Langer et al. (1987), Lesser (1992), Komistek et al. (1998),
Tisell (2000) andYamaguchi (2001).While inKane’s approach some auxiliary
fictitious generalized speeds are used to identify the noncontributing forces,
here we introduce open-constraint coordinates which express prohibited rela-
tive motions in the joints, in addition to the joint coordinates that describe
relative configurations of adjacent body segments. The followed augmented jo-
int coordinate method yields standard formulae that lead to theminimal-form
dynamic equations, and simultaneously, a pseudo-inverse matrix to the joint
constraint matrix is obtained without much effort. Using the pseudo-inverse,
the formulae for the determination of joint reactions are obtained directly in
a ”resolved” form (nomatrix inversion is required). Thismakes the developed
schemeparticularlywell suited for both symbolicmanipulations and computer
implementations.

2. Augmented joint coordinate method

For simplicity reasons, the method is introduced by using a seven-
segment multibody system shown in Fig.1a, which models the human bo-
dy in planar motion with collateral movements of the lower and upper
extremities (Blajer and Czaplicki, 2001). The position of the nine-degree-of-
freedom system can explicitly be described by k =9 generalized coordinates



An alternative scheme for determination... 815

q = [xH,yH,ϕ1,ϕ2,ϕ3,ϕ4,ϕ5,ϕ6,ϕ7]
>, where xH and yH are the hip co-

ordinates, and the angular coordinates ϕi (i = 1, . . . ,7) that describe the
angular onfigurations of the seven segments are all measured from the verti-
cal direction. The six control torques due to muscle forces at the joints are
τ = [τ1,τ2,τ3,τ4,τ5,τ6]

>.

Fig. 1. Seven-segmentmultibody model (a), and open-constraint coordinates and
reaction forces in joints of lower extremities (b)



816 W.Blajer, A.Czaplicki

Thederivation of equations ofmotion in q for the above systembelongs to
standardmultibody codes, well described in the literature; see e.g. Kane and
Levinson (1985), Nikravesh (1988), Lesser (1992), Garćıa de Jalón and Bayo
(1993), Schiehlen (1997), andBlajer (2001). The starting point is the dynami-
cal formulation in n=21 absolute coordinates p= [x1,y1,θ1, . . . ,x7,y7,θ7]

>,
where xi, yi and θi are the coordinates of the mass center Ci and the orien-
tation angle (here θi = ϕi) of the ith segment, i = 1, . . . ,7. The absolute
coordinates p are dependent because of the joints, modeled by m=12 kine-
matic constraints on the bodies, k = n−m, called by Langer et al. (1987)
themodel-intrinsic constraints (as distinct from themodel-environment ones).
The constraints express the prohibited relative translations in the joints, and
m local open-constraint coordinates z = [z1, . . . ,zm]

> can be introduced to
describe these prohibited relative motions (illustrated in Fig.1b for the joints
of lower extremities). Since z can be expressed in terms of p, the constraint
equations at the position, velocity and acceleration levels, given in implicit
forms (Schiehlen, 1997; Blajer, 2001), are

z=Φ(p)=0

ż=C(p)ṗ=0 (2.1)

z̈=C(p)p̈−ξ(p, ṗ)=0

where C= ∂Φ/∂p is the m×n constraint matrix of themaximum row-rank,
i.e. constraints (2.1)1 are independent, and the mvector ξ=−Ċṗ involves the
constraint-inducedaccelerations on thebodysegments.Aparticular constraint
equation Φj (j = 1, . . . ,m) depends only on the absolute coordinates of the
adjacent bodies in the joint.
Since the coordinates z describe the prohibited relative translations

of the body segments in the joints, the associated constraint reactions
λ = [λ1, . . . ,λm]

> represent physical joint reaction forces related to z
(Fig.1b). The constrained Newton-Euler equations of motion are then

Mp̈=h−C>λ (2.2)

where M = diag(M(1), . . . ,M(7)) and h = [(h(1))>, . . . ,(h(7))>]> are the
constant generalized mass matrix and the generalized applied forces related
to p, M(i) = diag(mi,mi,JCi), mi and JCi are the mass of the ith segment
and itsmassmoment of inertiawith respect to C1, and h

(i) = [Fxi,Fyi,MCi]
>

are components of the total force applied to the ith segment and the total
torque applied with respect to the mass center Ci.
In the case at hand, h consists of gravitational forces, interactions from

the environment (reactions of themodel-environment constraints) andmuscle
forces (control torques).



An alternative scheme for determination... 817

The standard joint coordinate formulation of multibody system dynamics
(Nikravesh, 1988) is basedon relationsbetween theabsolute coordinates pand
the joint coordinates q, p = g(q), which are the model-intrinsic constraints
given explicitly (Schiehlen, 1997; Blajer, 2001). Implicit constraint equations
(2.1)1 are satisfied identically when expressed in term of q, Φ(g(q))≡0.

In the present augmented joint coordinate method, the explicit constraint
equations are extended by incorporating the open-constraint coordinates z

p= g(q,z) (2.3)

Since z = 0, the above relations are equivalent to the traditional form
p = g(q), and are usually not much more difficult to formulate. In fact, the
dependence on z in (2.3) is needed only to grasp the prohibited motion di-
rections related to ż, called auxiliary fictitious generalized speeds in Kane’s
method (Kane and Levinson, 1985). Namely, by differentiating (2.3) with re-
spect to time and then setting z=0, we obtain

ṗ=
(∂g

∂q

)
∣

∣

∣

z=0
q̇+
(∂g

∂z

)
∣

∣

∣

z=0
ż≡D(q)q̇+E(q)ż (2.4)

while the standard formulation p = g(q) yields simply ṗ = D(q)q̇. Again,
since the maintenance of joint constraints assures ż = 0, both relations are
genuinely equivalent.Theexplicit constraint equations at the acceleration level
are then considered in the standard form

p̈=D(q)q̈+γ(q, q̇) (2.5)

where γ = Ḋq̇ is an n-dimensional vector.

As shown in Blajer (2001), the n× k matrix D defined above is an or-
thogonal complement matrix to the constraint matrix C introduced in (2.1)2,
i.e.

CD=0 ⇔ D>C> =0 (2.6)

The said prohibited local motion directions in the joints are then represented
as columns in the n×m matrix E defined in (2.4). An important charac-
teristic of E results from the substitution of (2.4) into (2.1)2, which gives
ż=CDq̇+CEż. Since, according to (2.6), CD=0, it can then be concluded
that

CE= I ⇔ E>C> = I (2.7)

where I denotes the m×m identity matrix. The n×m matrix E produced
in (2.4) has thus features of a pseudo-inverse (Ben-Israel and Greville, 1980)



818 W.Blajer, A.Czaplicki

of the rectangular m×n constraint matrix C. It is worth to note that E is
obtained here symbolically, and the constraintmatrix C is not involved in the
process.

Starting fromdynamic formulation (2.2) in the absolute coordinates p, the
explicit forms of the model-intrinsic constraints, p = g(q) ⇒ ṗ = D(q)q̇ ⇒
p̈=D(q)q̈+γ(q, q̇) are enough for converting the constraint reaction-induced
dynamic equations to a minimal set of constraint reaction-free dynamic equ-
ations in q. Then, the pseudoinversematrix E enables efficient determination
of the ’eliminated’ constraint reactions. The transformation formula is

[

D
>

E
>

]

(

M(Dq̈+γ)−h+C>λ
)

=0 (2.8)

while in the standard formulation of the projection method scheme (Blajer,
2001), E> is replaced by CM−1. The first k equations of (2.8) yield the
requested dynamic equations in q

M(q)q̈+d(q, q̇)=h(q, q̇, t) (2.9)

where M = D>MD is the k × k generalized mass matrix related to q,
d = D>Mγ is the k vector of generalized dynamic forces due to centrifu-
gal and Coriolis accelerations, and h = D>h is the k vector of generalized
applied forces. From the last m equations of (2.8), one obtains

λ(q, q̇, q̈, t)=E>
(

h−M(Dq̈+γ)
)

(2.10)

which offers a novel formula for the determination of joint reaction forces,
distinct from the traditional one obtained after replacing E> with CM−1 in
(2.8), see e.g. Blajer (2001), Garćıa de Jalón and Bayo (1993), and Schiehlen
(1997) for more details, i.e.

λ(q, q̇, t)= (CM−1C>)−1C(M−1h−γ) (2.11)

which is seldom used in biomechanical models.

There are at least three advantages of the new scheme expressed by (2.10)
over the traditional one given by (2.11). Firstly, the joint reactions λ are
directly obtained in a ”resolved” form - no matrix inversion is required as
it is necessary in (2.11). The present scheme is thus particularly well suited
for both symbolic manipulations and computational implementations. Then,
formula (2.10) does not rely on implicit forms (2.1) ofmodel-intrinsic constra-
int equations, which need not to be introduced at all. Both the derivation of



An alternative scheme for determination... 819

minimal-form dynamic equations (2.11) and the determination of joint reac-
tions are now based on the constraint equations given implicitly in augmented
form (2.3). The modification is not usually concerned with much additional
modeling effort. Finally, the present scheme can conveniently be used to deter-
mine only some joint reactions - only respective entries of z canbe introduced,
and a number of columns of E can appropriately be reduced. Such a situation
is illustrated in Fig.1b, where the open-constraint coordinates only in the
joints of lower extremities are involved.

3. Illustration

The proposedmethodwas used to calculate reactions in the joints of lower
extremities of the human planar multibody model shown in Fig.1a, used by
Blajer and Czaplicki (2001) in analysis of somersaults on a trampoline. The
vector h in dynamic equations (2.2) contains gravitational forces,muscle force
torques, and, during the support phase, the interaction from the trampoline
bed.Theway of estimation of this interaction, identification of themodel, and
inverse dynamics procedure applied to solve the inverse dynamics problem is
described in Blajer and Czaplicki (2001).

Since we are interested solely in the joint reactions in lower extremities,
the open-constraint coordinates z= [z1, . . . ,z6]

> are introduced only to joints
H, K and A (Fig.1b). Consequently, denoted g = [(g(1))>, . . . ,(g(7))>]>,
relation (2.3) will take the traditional form p(i) = g(i)(q) for i=1,5,6,7, and
the augmented form p(i) = g(i)(q,z) for i = 2,3,4. Omitting for shortness
the former relations, the latter ones are

x2 =xH +z1+ c2 sinϕ2

y2 = yH +z2− c2 cosϕ2

θ2 =ϕ2

x3 =xH +z1+ l2 sinϕ2+z3+c3 sinϕ3

y3 = yH +z2− l2 cosϕ2+z4− c3 cosϕ3 (3.1)

θ3 =ϕ3

x4 =xH +z1+ l2 sinϕ2+z3+ l3 sinϕ3+z5+ c4 sinϕ4

y4 = yH +z2− l2 cosϕ2+z4− l3 cosϕ3+z6− c4 cosϕ4

θ4 =ϕ4



820 W.Blajer, A.Czaplicki

where l2 and l3 are the total lengths of segments 2 and 3, and c2, c3 and c4
are the distances HC2,KC3 and AC4, respectively, and Ci (i=1, . . . ,7) are
the mass centers of the seven segments as denoted in Fig.1.
Using the relations for i=1, . . . ,7, the 21×9-dimensional matrix D and

the 21-vector γ defined in (2.4) and (2.5) can easily be obtained, and dynamic
equations (2.9) in q canbederived.Theexplicit formof thedynamic equations
is reported in Blajer and Czaplicki (2001). As far as the joint reactions are
concerned, let us denote first E = [(E(1))>, . . . ,(E(7))>]>, where E(i) is the
3×6-dimensionalmatrix related to the ith segment. In the 21×6-dimensional
matrix E, we have then E(1) = E(5) = E(6) = E(7) = 0, and the component
matrices related to the joints of lower extremities (rows 4-12 of E) are

E
(2) =







1 0 0 0 0 0
0 1 0 0 0 0
0 0 0 0 0 0






E
(3) =







1 0 1 0 0 0
0 1 0 1 0 0
0 0 0 0 0 0







(3.2)

E
(4) =







1 0 1 0 1 0
0 1 0 1 0 1
0 0 0 0 0 0







The analytical relations for the joint reactions in H, K and A joints are then

λ1 = Rx−m2(ẍH + c2ϕ̈2 cosϕ2− c2ϕ̇
2
2 sinϕ2)+

− m3(ẍH + l2ϕ̈2 cosϕ2+ c3ϕ̈3 cosϕ3− l2ϕ̇
2
2 sinϕ2− c3ϕ̇

2
3 sinϕ3)+

− m4(ẍH + l2ϕ̈2 cosϕ2+ l3ϕ̈3 cosϕ3+ c4ϕ̈4 cosϕ4− l2ϕ̇
2
2 sinϕ2+

− l3ϕ̇
2
3 sinϕ3− c4ϕ̇

2
4 sinϕ4)

λ2 = Ry − (m2+m3+m4)g−m2(ÿH + c2ϕ̈2 sinϕ2+ c2ϕ̇
2
2 cosϕ2)+

− m3(ÿH + l2ϕ̈2 sinϕ2+c3ϕ̈3 sinϕ3+ l2ϕ̇
2
2 cosϕ2+ c3ϕ̇

2
3 cosϕ3)+

− m4(ÿH + l2ϕ̈2 sinϕ2+ l3ϕ̈3 sinϕ3+ c4ϕ̈4 sinϕ4+ l2ϕ̇
2
2 cosϕ2+

+ l3ϕ̇
2
3 cosϕ3+ c4ϕ̇

2
4 cosϕ4)

λ3 = Rx−m3(ẍH + l2ϕ̈2 cosϕ2+ c3ϕ̈3 cosϕ3− l2ϕ̇
2
2 sinϕ2− c3ϕ̇

2
3 sinϕ3)+

− m4(ẍH + l2ϕ̈2 cosϕ2+ l3ϕ̈3 cosϕ3+ c4ϕ̈4 cosϕ4− l2ϕ̇
2
2 sinϕ2+

− l3ϕ̇
2
3 sinϕ3− c4ϕ̇

2
4 sinϕ4) (3.3)

λ4 = Ry − (m3+m4)g−m3(ÿH + l2ϕ̈2 sinϕ2+ c3ϕ̈3 sinϕ3+ l2ϕ̇
2
2 cosϕ2+

= c3ϕ̇
2
3 cosϕ3)−m4(ÿH + l2ϕ̈2 sinϕ2+ l3ϕ̈3 sinϕ3+ c4ϕ̈4 sinϕ4+

+ l2ϕ̇
2
2 cosϕ2+ l3ϕ̇

2
3 cosϕ3+c4ϕ̇

2
4 cosϕ4)



An alternative scheme for determination... 821

λ5 = Rx−m4(ẍH + l2ϕ̈2 cosϕ2+ l3ϕ̈3 cosϕ3+ c4ϕ̈4 cosϕ4− l2ϕ̇
2
2 sinϕ2+

− l3ϕ̇
2
3 sinϕ3− c4ϕ̇

2
4 sinϕ4)

λ6 = Ry −m4g−m4(ÿH + l2ϕ̈2 sinϕ2+ l3ϕ̈3 sinϕ3+ c4ϕ̈4 sinϕ4+

+ l2ϕ̇
2
2 cosϕ2+ l3ϕ̇

2
3 cosϕ3+ c4ϕ̇

2
4 cosϕ4)

where m2,m3 and m4 are themasses of segments 2, 3 and 4, Rx and Ry are
the components of the force acting on the feet (segment 4) fromthe trampoline,
and g is the gravity acceleration. It my be worth noting that the above result
is a little different from that possibly obtained by recursive application of
the classical dynamic equilibrium D’Alembert’s principle, starting from the
last (fourth) to the first (second) segment. In the present solution, the joint
reaction forces are obtained in a resolved form,while in the other solution, the
reactions λ3 and λ4 (in K joint) are dependent on λ5 and λ6 (in A joint),
and the reactions λ1 and λ2 (in H joint) are dependent on λ3 and λ4 (in
K joint).

Fig. 2. Reaction forces in joints of lower extremities during back somersault on a
trampoline



822 W.Blajer, A.Czaplicki

A back somersault in the pike position was chosen to present the results
of calculations. By using qd(t), q̇d(t), q̈d(t),Rxd(t) and Ryd(t) frommeasure-
ments previously reported by Blajer and Czaplicki (2001), the joint reactions
were determined from (2.10). Figure 2 shows time variations of the reactions
in lower extremities during a single somersault. Obviously, the joints aremost
loaded at the supporting phase, 0−0.4. One can easily see a significant drop
of the vertical reaction from the ankle to hip level, during this phase.

4. Conclusions

The proposedmethod for the determination of joint reactions is relatively
simple and naturally assists minimal-form formulation (2.9) of human body
dynamics. Computational scheme (2.10) is particularly efficient - no matrix
inversion is required as it is needed in classical code (2.11). As such, scheme
(2.10) is well suited for both symbolicmanipulations and computer implemen-
tations.
The idea of the proposedmethod is not new, it hasmuch in commonwith

Kane’smethodof ”bringingnoncontributing forces into evidence”.Thenovelty
of the present formulation lies in its compactness and clarity. Themethod for
the determination of joint reactions is presented in a systematic and practical
way.
For simplicity reasons, the proposedmethod is introduced and illustrated

in this paper bymeans of a very simple planarmodel of the human body. Ne-
vertheless, the methodology behind is general. The formulation is extendable
to a three dimensional case. A more realistic musculoskeletal model can be
used as well.

References

1. Ben-Israel A., Greville T.N.E., 1980, Generalized Inverses: Theory and
Applications, Robert E. Krieger Publishing Company, Huntington, NewYork

2. Bergman G., Deuretzenbacher G., Heller M., Graichen F., Rohl-
mannA., Strauss J.,DudaG.N., 2001,Hip contact forces and gait patterns
from routine activities, Journal of Biomechanics, 34, 859-871

3. Blajer W., 2001, A geometrical interpretation and uniform matrix formula-
tion of multibody system dynamics, ZAMM, 81, 247-259



An alternative scheme for determination... 823

4. Blajer W., Czaplicki A., 2001, Modeling and inverse simulation of somer-
saults on the trampoline, Journal of Biomechanics, 34, 1619-1629

5. Eberhard P., Spgele T., Gollhofer A., 1999, Investigations for the dy-
namical analysis of humanmotion,Multibody System Dynamics, 3, 1-20

6. Garćıa de Jalón J., BayoE., 1993,Kinematic and Dynamic Simulation of
Multibody Systems: The Real-time Challenge, Springer, NewYork

7. Kane T.R., Levinson D.A., 1985, Dynamics: Theory and Applications,
McGraw-Hill, New York

8. KomistekR.D., StiehlJ.B.,DennisD.A.,PaxsonR.D., Soutas-Little
R.W., 1998, Mathematical model of the lower extremity joint reaction forces
using Kane’s method of dynamics, Journal of Biomechanics, 31, 185-189

9. Langer F.D., Hemami H., Brown D.B., 1987, Constraint forces in holono-
micmechanical systems,ComputerMethods in AppliedMechanics and Engine-
ering, 62, 255-274

10. LesserM., 1992,Ageometrical interpretationofKane’s equations,Proceedings
of the Royal Society of London,A 436, 69-87

11. McNitt-Gray J.L., Hester D.M.E., Mathiyakom W., Munkasy B.A.,
2001, Mechanical demand and multijoint control during landing depend on
orientation of the body segments relative to the reaction force, Journal of Bio-
mechanics, 34, 1471-1482

12. Nikravesh P.E., 1988, Computer-Aided Analysis of Mechanical Systems,
Price-Hall, Enlewood Cliffs

13. Pandy M.G., 2001, Computermodelling and simulation of humanmovement,
Annual Review of Biomedical Engineering, 3, 245-273

14. Schiehlen W., 1997, Multibody system dynamics: roots and perspectives,
Multibody System Dynamics, 1, 149-188

15. Tisell C., 2000, Integration of Database Technology and Multibody System
Analysis, Doctoral Thesis, Department of Machine Design, Royal Institute o
Technology, Stockholm, Sweden

16. Tözeren A., 2000,Human Body Dynamics: Classical Mechanics and Human
Movement, Springer, NewYork

17. WismansJ.S.H.M.,JanssenE.G.,BeusenbergM.,KoppensW.P.,Lup-
ker H.A., 1994, Injury Biomechanis, Course notes, Faculty ofMechanical En-
gineering, EindhovenUniversity of Technology, Endhoven

18. Yamaguchi G.T., 2001, Dynamic Modeling of Musculoskeletal Motion, Klu-
wer, Norwell, Massachusetts



824 W.Blajer, A.Czaplicki

Alternatywny sposób wyznaczania reakcji w stawach wieloczłonowych

modeli ciała człowieka

Streszczenie

Wieloczłonowe modele ciała człowieka są powszechnie wykorzystywane do ana-
lizy czynności motorycznych. Dla sformułowań dynamiki tych modeli stosowane są
zwykle niezależne współrzędne uogólnione, co powoduje, że reakcje w połączeniach
(idealne reakcje więzów wewnętrznych) są eliminowane na wstępnym etapie mode-
lowania. Dla ich określenia wymagane są dodatkowe procedury modelowania mate-
matycznego, a generowane tą drogą zależności charakteryzują się niską efektywnością
numeryczną. W niniejszej pracy proponowane jest nieco inne podejście do wyzna-
czania reakcji w stawach, w sposób naturalny skojarzone z minimalno-wymiarowym
formułowaniemdynamikiwieloczłonowychmodeli ciała człowieka.Proponowane sfor-
mułowanianiewymagająodwracaniamacierzy, są tym samymefektywne zarównodla
wyprowadzeń symbolicznych, jak i zastosowań numerycznych. Metoda zilustrowana
jest za pomocą siedmioczłonowego płaskiego modelu ciała człowieka. Prezentowane
są wybrane wyniki obliczeń numerycznych odnoszące się do symulacji dynamicznej
odwrotnej sportowca wykonującego salto na trampolinie.

Manuscript received January 05, 2005; accepted for print April 19, 2005