MODELOWANIE INŻYNIERSKIE ISNN 1896-771X 32, s. 45-52, Gliwice 2006
DESIGN OF A BIOMECHANICAL MODEL FOR THE DETERMINATION OF MUSCLE FORCES
IN LOWER EXTREMITIES
WOJCIECH BLAJER
KRZYSZTOF DZIEWIECKI
ZENON MAZUR
Institute of Applied Mechanics, Technical University of Radom
Abstract. Foundations for the design of a human biomechanical model aimed at analyzing planar movements such as vertical or standing long jumps are presented.
The motions in the hip, knee and ankle joints are modeled as enforced by muscle forces applied to the tendon attachment points, and the other joint motions are ac- tuated by torques representing the muscle action. A systematic construction of the related dynamic equations in independent coordinates is developed, enhanced by an effective scheme for determination of reaction forces in the leg joints. The reaction forces from the ground during the support phase can also be obtained.
1. INTRODUCTION AND MOTIVATION
The determination of muscle forces during human movements can play an important role in a deeper understanding of the underlying neural control. It can also be essential for the analysis of internal loads acting on bones and joints. Since the measuring of muscle forces directly within the living beings by means of invasive techniques may be difficult and even injurious, the other possibility is to evaluate the muscle forces following the dynamic analysis based on a biomechanical model and some input data obtained from noninvasive measurements. The reli- ability of results obtained this way is strongly dependent on thoroughness of the mathematical model used.
This paper presents foundations for the design of a biomechanical model aimed at analyzing human movements such as standing long jumps, vertical jumps, jumps down from a height, …, in which both the upper and lower extremities are moving parallel to one another (planar mo- tion). Since the attention is focused on lower limb motion control and loadings, only the mo- tions in the hip, knee and ankle joints are modeled as enforced by muscle forces applied to their tendon attachment points, and the actuation of other (upper body part) joint motions is simpli- fied to torques representing the muscle action. A systematic construction of the related dy- namic equations in independent coordinates is developed, enhanced by an effective scheme for the determination of reaction forces in the leg joints. The reaction forces from the ground dur- ing the support phases of the jumps are also obtained. Using the motion characteristics meas-
ured during the analyzed jumps, the inverse dynamics solution based on the developed mathe- matical model will result in time-variations of the joint torques, which are then shared in the hip, knee and ankle joints on particular muscle forces following an optimization procedure.
Apart from understanding of how the muscle excitations are coordinated during the jumps, the internal loads in the leg joints can also be evaluated.
2. THE BIOMECHANICAL MODEL
x , y
( )
H
y
H
x
H
1
a)
7
5 6
2
4 3 9
8
10
8
21 23 13
18 19 16 15 9 10
b)
τ
τ 4 5
τ9 5
6
τ7
τ8
22 12
17 20 14 11 2
6 7 3 1
Fig. 1. The human biomechanical model
The developed biomechanical model is constructed as a planar kinematic structure seen in Fig. 1a. It consists of N =10 rigid segments branching from the pelvis in the open chain link- ages. The k=12 generalized coordinates that describe the system position with respect to the inertial reference frame are q=[xH yH ϕ K1 ϕ10]T, where x and H y are the hip coordinates H and the angular coordinates ϕ (i i=1 K, ,N ) are measured from the vertical direction.
As said before, the attention of this study is focused on the lower limb motion control and loadings. The actuation in the lower limb joints is then modeled by means of mF =23 muscle forces F K1, ,F23 applied to some tendon attachment points, while the actuation in the other (upper body part) joints is simplified to mτ =5 torques τ K5 , ,τ9 representing the muscle ac- tion in the joints (Fig. 1b). The total vector of m=mF +mτ =28 control inputs is finally
F T
F ]
[ 1 K 23 τ K5 τ9
=
u . Alternatively the control inputs could be represented by m′=9 entries as u′=[τ K1 τ9]T – the resultant torques in all the nine joints (as it was modeled e.g.
in [1,2]). Confronting the nine control inputs with twelve degrees of freedom of the human
model, m′<k, one can easily deduce that in the flying phase we deal with an underactuated system – the mass center of the system moves along a parabola (or vertically) with acceleration g pointed downward and its total angular momentum remains constant, and these motions are not actuated by the internal control inputs. In this contribution we replace the torques τ , 1 τ ,2 τ and 3 τ in the lower limb joints with the muscle forces 4 F K1, ,F23, and a local control redundancy is faced. The distribution of the torques on particular muscle forces is usually achieved by using different optimization procedures [3]. A scheme of this type will be dis- cussed in the sequel.
3. THE INITIAL DYNAMIC EQUATIONS IN ABSOLUTE COORDINATES
In order to derive effectively the dynamic equations of the human model in coordinates q, first the dynamic equations in absolute coordinates p need to be formulated. For the case at hand, the n=3N =30 absolute coordinates are p=[xC1 yC1 θ K1 xC10 yC10 θ10]T, where x , Ci y and Ci θ are the coordinates of the mass center i C and the orientation angle (here we i used θi =ϕi) of the i-th body segment with respect to the inertial frame, i=1 K, ,N. The ini- tial dynamic equations in p are then the constrained Newton-Euler equations in the form
λ C u B f p
M&&= g + − T (1)
where M=diag(m1,m1,JC1,K,m10,m10,JC10) is the generalized mass matrix related to p, m i and J are the mass and mass moment of inertia with respect to Ci C of the i-th segment, i
T g =[0 −m1g 0 L 0 −m10g 0]
f contains the gravitational forces, fu =Bu is the n-vector of generalized control forces, with B being the n×m matrix of distribution of control inputs u in the directions of p, and fc=−CTλ is the n-vector of generalized reaction forces due to the
18 2 ′=
= m
l kinematical joint constraints, with C being the l×n constraint matrix and
T l] [λ L1 λ
=
λ containing the reaction forces in the joints.
l i
C i
a) i
i-th body
i
A
b)
j
B j
F j j
F j
Fig. 2. The muscle forces on the i-th body (a) and angular orientation of the j-th force (b)
The formulation of M and f in Eq. (1) is evident and, as it will be seen in the sequel, C g needs not to be introduced explicitly. The formulation of B is a little more challenging, how- ever. This concerns in particular the first m columns of B related to F F K1, ,F23 (the entries of the last m columns of B are equal to either 0, 1 or –1, appropriately; see [1] for more details). τ
In general many muscle forces can act on the i-th body (Fig. 2a). The influence if the j-th force F (the j-th column of B) can be found in the following four steps. j
1) Detect the two bodies to which the tendons of force F are attached. Only these two j bodies will be affected by F , and only the entries (rows) of the j-th column of B that j correspond to the bodies must be found (the others are in principle equal to zero).
2) Determine the inertial frame coordinates rAj =[xAj yAj]T and rBj =[xBj yBj]T of ten- dons attachment points A and j B (Fig. 2b). Given the segment lengths and coordi-j nates of points A and j B in the local body-fixed reference frames of appropriate j lower limb segments, one can easily determine rAj =rAj(q) and rBj =rBj(q), and then specify the angle α defining the direction of force j F . j
3) Determine the inertial frame coordinates rCi =[xCi yCi]T of mass centers C of the in-i volved lower limb segments. Following the procedure as in Step 2, find rCi =rCi(q). 4) Using basic formulae of mechanics, determine the x and y components of F acting on j
the respective two segments of lower extremities, and the moments of F with respect j to the mass centers of the bodies. Retrieve from this the j-th column of B.
Using the above procedure for j=1 K, ,mF the first m columns of B can be determined. F Then, augmenting them with the last m columns of B defined separately, the whole matrix B τ can be formulated, written symbolically as B(q).
4. DYNAMIC EQUATIONS IN INDEPENDENT COORDINATES
1 4
H
2 2
z
11
1
z
26
z
6 6 2
z
34 3
4
z
4z
6
6 6
Fig. 3. The open-constraint coordinates and the reaction forces in lower extremity joints
The augmented joint coordinate method described in [1,2,4] is applied to both derive the dynamic equations in independent coordinates q and obtain especially useful formulae for the
determination of joint reaction forces. In this method, instead of the constraint equations given explicitly as p=g(q), they are introduced in an augmented form
) , (q z g
p= (2)
where z=[z1 K zl]T are open-constraint coordinates that describe the prohibited relative motions in the joints. Specifically, since z=0, Eq. (2) is virtually equivalent to the explicit constraint equation p=g(q), and the dependence on z is needed only to grasp the prohibited motion directions related to z& – the directions of constraint reactions in the respective joints.
Moreover, the open-constraint coordinates can be introduced only in those joints in which the reaction forces are to be determined. In the case at hand we open thus only the hip, knee and ankle joints, and thus z=[z1 z2 z3 z4 z5 z6]T as seen in Fig. 3. Only l′=6 reaction forces
]T
[λ1 λ2 λ3 λ4 λ5 λ6
′=
λ in the lower extremity joints will thus be determined, where λ′ is a subset of λ introduced in Eq. (1). The explicit forms of relations (2) are not difficult to derive for the present human model, though their number is n=30 and some of them are rather long.
For reasons of compactness they will not be reported here.
By differentiating Eq. (2) with respect to time, and then setting z=0, one obtains z
q E q q D z z
q g q p g
0 0 z
z
&
&
&
&
& = ( ) + ( )
∂ + ∂
∂
= ∂
= =
(3)
while the explicit formulation of constraint equation p=g(q) yields simply p& =D(q)q&. Again, since the maintenance of joint constraints requires z&=0, both the relations are inherently equivalent. The explicit constraint equations at the acceleration level are then
) , ( )
(q q γ q q D
p& && &
& = + (4)
where γ =D& q&. As shown in [1,2,4], the n×k (here 30×12) matrix D is an orthogonal com- plement matrix to the l×n (here 18×30) constraint matrix C introduced in Eq. (1), i.e.
0 C D 0
D
C = ⇔ T T = (5)
The n×l′ (here 30×6) matrix E produced in Eq. (3) has then the features of a pseudo-inverse matrix [5] to the rectangular l′×n matrix C′ being composed of those rows of matrix C which correspond to the open-constraint coordinates z, i.e. C′E=I ⇔ ETC′T =I, where I denotes the l′×l′ identity matrix. On the other hand, ETC′′T =0, where C′′ contains the other rows of C (not contained in C′). Denoted C=[C′T C′′T]T, we have finally
[
I 0]
C 0 E
E I
C ⇔ =
= T T (6)
Starting from the dynamic equations (1) in absolute coordinates p, the dynamic equations in independent coordinates q and the joint reaction forces in the lower extremity joints can be ob- tained effectively using the projection formula in the form
(
M Dq γ f Bu C λ)
0E
D + − − + =
T
T g T
)
( && (7)
The first k =12 components of Eq. (7) lead to the requested dynamic equations in q, i.e.
u B f d q
M&&+ = g + (8)
where M(q)=DTMD is the k×k generalized mass matrix related to q, d(q,q&)=DTMγd and fg =DTfg are the k-vectors of generalized dynamic forces related to q due to the centrifu- gal accelerations and gravitational forces, respectively, and B=DTB is the k×m (here
28
12× ) matrix of distribution of control inputs u in the directions of q. The last l′=6 com- ponents of Eq. (7) lead then to
[
( )]
)
(q,q,q,u E f Bu M Dq γ
λ′ & && = T g + − &&+ (9)
Which offers an effective formula for the determination of joint reaction forces in the lower extremity joints, λ′=[λ1 λ2 λ3 λ4 λ5 λ6]T (Fig. 3). Compared to the more traditional formu- lation λ(q,q,u)=(CM−1C )−1C[M1(fg +Bu)−γ]
& T (refer to [4] for more details), in the pre- sent scheme the joint reactions are obtained in a ‘resolved form (no matrix inversion is in- volved). The scheme (9) does not require the formulation of implicit forms of constraints to produce C, either, which need not to be introduced at all. Finally, in the traditional scheme all
=18
l implicit constraint equations (related to all the system joints) must be introduced to pro- duce the l×n (here 18×30) matrix C(p)→C[g(q)], and all l=18 constraint reaction forces
λ defined in Eq. (1) must be determined simultaneously. In the present scheme we can deter- mine only the chosen joint reactions defined by the introduced open-constraint coordinates z.
The determination of matrix E is not usually concerned with much additional effort.
5. THE INVERSE DYNAMICS SIMULATION
The considered inverse dynamics simulation can be stated as follows: given the motion characteristics qd(t), q&d(t) and q&&d(t) measured during the analyzed motions, the control
)
d(t
u that force the system to complete this specified motion can be determined from the dy- namic equations (8), and then the joint reaction forces λ′d(t) in lower extremities can be found from Eq. (9). While the latter issue is evident, the determination of ud(t) is concerned with two difficulties to overcome.
The first difficulty is concerned with the fact that the ‘flying’ human model is globally un- deractuated – the k-degree-of-freedom system (k=12) is actuated by m′=9 muscle torques
]T
[τ K1 τ9
′=
u applied in the m′ joints, m′<k (see Section 2). Defining the k×m′ matrix B
D
B′= T ′, where B′u′ corresponds to B in Eq. (1) and the u n×m′ matrix B′ is easy to for- mulate, the dynamic equations can be projected into the controlled subspace which leads to
(
1( ))
( ) ( , , )1Bu B q M d f H q u h q q q
M
B′T − = ′T &&+ − − g ⇔ = & && (10)
where H is the m′×m matrix, and h is the m′-vector. On the other hand, the projection into the uncontrolled subspace gives A′(M&q&+d−fg)=0, where A′ is an r×k (r=k−m′=3) orthogonal complement to B′, i.e. A′B′=0, and as such also A′B=0. The projection into the uncontrolled subspace is equivalent to the condition that, during the flying phase, both the
body mass center moves along a parabola (or vertically) and that the total angular momentum of the body remains constant. The r=3 conditions can serve to verify the accuracy of the measured input data and correctness of the mathematical model build.
a) b)
3
1*
*2
* 3
*1
*2
*
Fig. 4. Two modes of foot-ground contact
The situation differs in the support phase. For both the modes of foot-ground contact seen in Fig. 4, r=3 constraints are imposed on the foot. Expressed in q, the constraint equations can be written symbolically as Φ(q)=0.Then, at the velocity and acceleration level they are respectively Cq& =0 and Cq&&−ξ=0, where C(q)=∂Φ ∂q is the r×k constraint matrix and
q C q
ξ(q,&)=−& & is the r-vector of constraint-induced accelerations [4]. The dynamic equations (8) modify then to Mq&&+d=fg +Bu−CTλ*, where λ* =[λ*1 λ*2 λ*3]T are the ground reac- tion forces and moment seen in Fig. 4. Introducing a k×m′ matrix D which specifies the un- constrained subspace – an orthogonal complement to C , i.e. CD=0 ⇔ DTCT =0, the pro- jection into the unconstrained subspace gives
) , , ( ) ( )
(Mq d f H q u h q q q
D u B
DT = T &&+ − g ⇔ = & && (11)
which has the same symbolic and dimensional form as Eq. (10). The projection into the con- strained subspace gives then λ* =(CM−1CT)
(
CM−1(fg −d+Bu)−ξ)
. The computed values of the ground reactions can then be compared with the measured ones.The final question is the distribution of joint torques in the lower extremity joints into the forces F K1, ,F23. Mathematically the problem is associated the solution of m′=9 Eqs. (10) and/or (11) with respect to m=28 unknowns u – the redundant control problem usually solved using different optimization techniques [3]. One possibility of this type is to use a pseudo-inverse to the m′×m matrix H defined in Eqs. (10) and/or (11), i.e. the m×m′ matrix
) 1
( −
=HT HHT
H† such that that HH† =I (identity matrix). Using the pseudo-inverse, ud(t) can be calculated as
(
( )) (
( ), ( ), ( ))
)
(t d t d t d t d t
d H q hq q q
u = † & && (12)
In the case some of the determined muscle forces are negative, they should be set to zero (or small positive values), and the calculations should be repeated until all the muscle forces are either positive or vanish.
6. CONCLUSION
A mathematical model of human body was systematically build, aimed at the inverse dynam- ics simulation of planar movements such as vertical or standing long jumps. A hybrid model of control was used, modeled by means of muscle forces in the lower extremity joints, and simpli- fied to the torques representing the muscle action at the other joints. Two different sets of mo- tion equations were derived, separately for the flying phase and the support phase, that enable for the determination the control (muscle forces/torques) in the joints based on measured mo- tion characteristics. The reaction forces in the lower extremity joints and the ground reactions during the support phase can also be determined. Some guidelines for the solution of the con- trol redundant problem that arises were suggested.
REFERENCES
1. Blajer W., Czaplicki A., „Contact modeling and identification of planar somersaults on the trampoline”, Multibody System Dynamics, 10 (2003), p. 289-312.
2. Blajer W., Czaplicki A., “An alternative scheme for determination of joint reaction forces in human multibody models”, Journal of Theoretical and Applied Mechanics, 43 (2005), p. 813-824
3. Yamaguchi, G.T., Dynamic modeling of musculoskeletal motion: a vectorized approach for biomechanical analysis in three dimensions, Kluwer, Dordrecht, 2001.
4. Blajer W., “On the determination of joint reactions in multibody mechanisms”, Transac- tions of the ASME, Journal of Mechanical Design, 126 (2004), p. 341-350
5. Ben-Israel A., Greville T.N.E., Generalized inverses: theory and applications, Robert E.
Krieger Publishing Company, Huntington, New York, 1980
BUDOWA MODELU BIOMECHANICZNEGO DLA WYZNACZANIA SIŁ MIĘŚNIOWYCH
W KOŃCZYNACH DOLNYCH
Streszczenie. W pracy przedstawiono założenia i sposób budowy modelu matematycznego ciała dla analizy czynności motorycznych typu wyskok do góry czy skok z miejsca do przodu. Ruch w stawach kończyn dolnych wymuszany jest za pomocą sił mięśniowych działających wzdłuż odcinków łączących miejsca ich zaczepów, natomiast ruchy w pozostałych stawach realizowane są za pomocą wypadkowych momentów sił mięśniowych. Zaproponowano systematyczny sposób generowania równań ruchu wymaganych dla rozwiązania zadania symulacji odwrotnej, uzupełnione o efektywny sposób wyznaczania reakcji w stawach kończyn dolnych oraz oddziaływań z podłożem.
