Kurfess T.R. Robotics and Automation Handbook
Подождите немного. Документ загружается.
Flexible Robot Arms 24
-15
The final equations then incorporate first order equations for the joints (one each) and second order
equations for each mass (six each).
˙
θ
J
= D
−1
c
[(−HK
s
H
T
− K
c
)θ
J
+ HK
s
(Z
n
− T
T
Z
0
) + K
c
θ
dJ
]
˙
θ
dJ
(24.31)
d
dt
(
˙
Z
n
) =−J
−1
n
K
s
(Z
n
− T
T
Z
0
− Gθ
J
) (24.32)
d
dt
(
˙
Z
0
) = J
−1
0
TK
s
(Z
n
− T
T
Z
0
− Gθ
J
) + J
−1
0
B
T
F
0
(24.33)
24.2.4 Distributed Models
When a component storeskinetic and elastic potential energy over the same spatially distributed region, the
modeler has twochoices.A larger number of lumped springsand masses can be interspersed to approximate
this distribution, or the distribution can be more fully modeled with partial differential equations or finite
elements. The distribution of lumped elements can be kept to a minimum, but the lumped approximation
becomes increasingly hard to fit to the physical reality. Alternatively, a finite difference approximation
can be made, but then the number of lumped elements increases and the equations become inefficient to
analyze or simulate. In this section we shall examine the alternative of recognizing the distributed mass and
elasticity in the basic modeling of the structural and drive components. As will be seen, approximations
are still necessary to enable analysis and simulation, and the nature and limitations of the approximations
must be understood.
24.2.4.1 Distributed Bar Elements
The simplest element to explain is the uniform bar element, justified when a cross section of material with
a modest aspect ratio is forced along its longitudinal (X) axis, which is the axis along which the cross
section is extruded to form the bar. The assumptions of plane strain and uniform stress over the cross
section permit the simple partial differential equation for the displacement ξ of the mass element of length
dx at location x along the X axis:
ρ
∂
2
ξ
∂t
2
= E
∂
2
ξ
∂x
2
(24.34)
24.2.4.2 Distributed Shaft Elements
Another simple behavior to model is distributed torsional stiffness due to the shear modulus G and
rotational inertia J . Again, the spatial variation of the twist angle φ is along the X axis. The polar moment
of inertia J is typically small for drive shafts relative to driven components such as gears, flywheels, or
beams. For this reason the distributed shaft elements can often be replaced with a lumped rotational spring.
However, the partial differential equation for the distributed model is
µ
∂
2
φ
∂t
2
= GJ
∂
2
φ
∂x
2
(24.35)
24.2.4.3 Distributed Beam Elements in Bending
Beams in bending are a major source of flexible behavior. Also, the kinetic energy and elastic potential
energy is often of equal consequence, especially when higher modes of vibration are considered. The
static model treated previously yielded the stiffness properties of a slender beam. In the distributed case,
a differential mass element is introduced into the partial differential equation. When the slender beam
approximation is augmented in this way with negligible shear deformation and negligible mass moment of
inertia of the differential cross section (not the entire beam), the first and most frequently suitable of model
results, the Bernoulli-Euler beam model. When these two assumptions are not justified, the complexity of
the Timoshenko beam equation (still manageable) must be introduced. The Timoshenko model will not
be examined here as it is typically not necessary for flexible robotic arms.
Copyright © 2005 by CRC Press LLC
24
-16 Robotics and Automation Handbook
For the Bernoulli-Euler beam, deflection w(x) perpendicular to the X axis varies with distance along
X, x and time t. We can consider it possible for the cross section to vary along X also, resulting in the
mass density per unit length m(x) and the area moment of inertia about the neutral axis, I (x), to also
vary. The material properties are assumed constant including the elastic modulus E . The resulting partial
differential equation is
m(x)
∂
2
w(x, t)
∂t
2
+
∂
2
∂x
2
EI(x)
∂
2
w(x, t)
∂x
2
= 0 (24.36)
24.2.4.4 Combined Distributed Effects and Components
By combining distributed effects we mean to combine in a single component the effects of various de-
formation processes: torsion, compression, and bending, for example. When components are attached to
each other, an element of the arm system must interact at its boundary with another element. One beam
may be attached to another or to a payload mass, for example. Both problems will be discussed below.
For small displacements a viable assumption is sometimes to treat each effect on a component as if it
were independent but occurring at the same value of the independent spatial variables. Hence, bending
in two orthogonal planes, axial compression and axial torsion could occur independently in the same
slender arm link. This is clearly an imperfect model, since simultaneous bending in the two directions
can produce torsion which would not be otherwise accounted for. Nevertheless, the combination allows
representative displacements to be considered, and when an element of this type is adjoined to another
element, it transfers energy in all the appropriate modes.
Combining elements or components involves matching boundary conditions on the partial differential
equations with a possible reorientation of the variables due to angle of attachment. The simultaneous
solution of the resulting systems of partial differential equations is extremely difficult in the general
case. Two methods around this obstacle are discussed below: (24.2.5) assume small displacement and
linear behavior and Laplace transform the partial differential equations (Frequency Domain Solutions),
or (24.2.6) convert the partial differential equations into simultaneous ordinary differential equations
(assumed mode solutions).
24.2.5 Frequency Domain Solutions
The treatment here will cover only bending of a planar arm. Bending is the most significant distributed
flexible effect, and the planar arm permits its illustration in a meaningful way. Extensions to more complex
situations of two axis bending, torsion, and compression are also useful. The transform of the differential
equation allows a convenient representation of the dynamics between two points on a structure with a
single independent spatial variable. This representation is known as the transfer matrix representation. An
extensive treatment of this representation is found in Pestel and Leckie [9] and has been applied to arms
as documented in Book et al. [10, 11]. Results that can be obtained from this implementation include
determination of the frequency response, eigenvalues and eigenfunctions (mode shapes), and impulse
response, and through the convolution of the impulse response and response to arbitrary inputs.
The transfer matrix relates variables at two stations along the arm as shown in Figure 24.2. The variables
are referred to as state variables of the partial differential equation but are not to be confused with the
Z
X
Y
w
V
Ψ
M
Station 0
Station 1
FIGURE 24.2 State vector for bending transfer matrix.
Copyright © 2005 by CRC Press LLC
Flexible Robot Arms 24
-17
states of the overall system. Hence for Figure 24.2 we have
z
1
=
−W
M
V
1
=
−displacement
angle
moment
shear force
at station1
z
0
= Bz
1
(24.37)
B is beam transfer matrix.
B =
c
0
−lc
1
ac
2
−alc
3
−β
4
c
3
/lc
0
ac
1
/lac
2
β
4
c
2
/a −β
4
lc
3
/ac
0
−lc
1
−β
4
c
1
/al β
4
c
2
/a −β
4
c
3
/lc
0
(24.38)
where
β
4
=ω
2
l
4
µ/(EI), a = l
2
/(EI)
c
0
=(cosh β + cos β)/2
c
1
=(sinh β + sin β)/(2β)
c
2
=(cosh β − cos β)/(2β
2
)
c
3
=(sinh β − sin β)/(2β
3
)
µ =density/unit length
ω =circular frequency of vibration
E =elastic (Young’s) modulus
I = cross sectional area moment of inertia
The spatial variable has been transformed to the spatial Laplace variable and replaced with its value at
the points of interest. The time Laplace variable remains and is designated as s or as the frequency variable
ω =−js where j =
√
−1. A brief tabulation of important transfer matrices for other planar elements
appears below.
For an angle in the plane of angle ϕ
A =
1/cos ϕ 00 0
0100
0000
m
s
ω
2
sin ϕ tan ϕ 00cosϕ
(24.39)
where m
s
is the sum of all outboard masses from the angle to the end of the arm.
For a rigid mass,
R =
1 −l 00
0100
mω
2
(l − h)(−hl + h
2
)ω
2
m + I ω
2
1 −l
−mω
2
hmω
2
01
(24.40)
where m is the mass of the body, I is the mass moment of inertia about an axis through the center of
mass and perpendicular to the plane of the arm, l is the length of the mass (distance between points of
attachment at stations i and i + 1), and h is the distance to the center of mass from station i.
Copyright © 2005 by CRC Press LLC
24
-18 Robotics and Automation Handbook
station O
station 7
R
1
R
3
R
7
B
2
B
6
A
4
,C
5
FIGURE 24.3 Transfer matrix representation of an arm.
For a controlled joint with the controller transfer function (joint torque/joint angle) = k( j ω)
C =
10 0 0
011/k( j ω)0
00 1 0
00 0 1
(24.41)
When combined together to represent the planar arm with two joints shown in Figure 24.3, the result is a
representation of the entire dynamics wrapped into a matrix product which relates the four state variables
at the two extremes of the arm. The known boundary conditions at these points can be imposed and the
resulting two equations employed to great benefit.
−W
M
V
0
= R
1
B
2
R
3
A
4
C
5
B
6
R
7
−W
M
V
7
(24.42)
In the example of Figure 24.3, the left end is pinned and the right end is free, yielding
z
0
=
−W
M
V
0
=
u
11
u
12
u
13
u
14
u
21
u
22
u
23
u
24
u
31
u
32
u
33
u
34
u
41
u
42
u
43
u
44
−W
M
V
6
(24.43)
u
11
u
12
u
31
u
32
−W
6
=
0
0
24.2.5.1 Eigenvalues and Corresponding Eigenfunctions
After examination of the four equations a set of homogeneous equations can be selected. Such a set
has two solutions, the trivial solution when the system is at rest and the nontrivial solution when the
determinant of the coefficient matrix is zero. Changing the boundary conditions on the arm does not
affect the Equation (24.42) above but does change the selection of the homogeneous equations and hence
Copyright © 2005 by CRC Press LLC
Flexible Robot Arms 24
-19
the frequency determinant.
d(ω) =
u
ik
u
il
u
jk
u
jl
ω=ω
i
= 0 (24.44)
z
i0
= z
j0
= 0 at station 0
z
kn
= 0, z
ln
= 0 at station n
Note that the frequency determinant is a function of the time Laplace variable. It is zero only for ω = ω
i
,
the system eigenvalues. Since the transfer function is transcendental, there will be an infinite number of
solutions which may be real, imaginary, or complex. The solution process is a numerical search over the
complex plane. Considerable care must be exercised to find all the most important eigenvalues.
If the eigenvalues ω
i
are substituted into the transfer matrix for ω, a relationship between the unknown
boundary conditions results. An arbitrary amplification factor remains that can be assigned at will, so, for
example, z
jn
could be given the value 1 and the relative value of the other boundary condition found from
the equation of the form:
u
ik
u
il
u
jk
u
jl
ω=ω
i
z
in
z
jn
=
0
0
(24.45)
Because valid boundary conditions are known, the product transfer matrix equation can be used to find the
states at any station along the arm, and new stations can be inserted along the flexible beams to determine
the mode shapes (eigenfunctions) corresponding to any of the eigenvalues.
24.2.5.2 Frequency Response and Impulse Response
If the boundary conditions are not zero, but are instead sinusoidal forcing functions of a frequency ω, the
analysis proceeds as follows. First, realize that the state variables occur in complementary pairs, and that
if a moment is forced sinusoidally, for example, the angle must also respond sinusoidally and cannot be
independently specified. The shear force for this example or the displacement could be specified to have
nonzero amplitude, although it is seldom helpful to do so. If the equations are rearranged into the forced
boundary variables f
i
and unforced, response boundary variables r
i
, the equations become
z
0
=
f
0
r
0
=
˜
U
11
˜
U
12
˜
U
21
˜
U
22
f
n
r
n
=
˜
Uz
n
(24.46)
and can be solved as
r
0
r
n
=
˜
U
22
˜
U
−1
12
˜
U
21
−
˜
U
22
˜
U
−1
12
˜
U
11
˜
U
−1
12
−
˜
U
−1
12
˜
U
11
f
0
f
n
(24.47)
assuming the inverse exists.
Thefrequencyofexcitationis substituted for ω, and theamplitudeofthe forcingsinusoid(s)is substituted
for the forcing vector. A solution for the response variables can thus be evaluated, and as in the case of the
mode shapes, additional stations can be arbitrarily generated according to the position along the arm.
Since the frequency response can be generated, it is also quite feasible to produce the impulse response
via numerically calculating the inverse transform for any frequency response.
24.2.6 Discretization of the Spatial Domain
It is possible to represent the spatially distributed nature of mass and compliance while retaining the
discrete representation of the system state (as opposed to representing both space and time as continuous
Copyright © 2005 by CRC Press LLC
24
-20 Robotics and Automation Handbook
variables). This is done by describing the continuously variable spatial shapes w(x, t)ofarmdeflection by
an infinite series of fixed basis functions (shapes) φ
i
(x) with time variable amplitudes q
i
(t).
w(x, t) =
∞
i=1
q
i
(t)φ
i
(x) (24.48)
A wise choice of the shapes achieves the necessary accuracy without a large number of basis functions,
which means the infinite series is truncated after a few terms. It is necessary that the basis functions span
the space and that they meet the geometric or essential boundary conditions. Convenient basis functions
include polynomials and unforced solutions to the partial differential equations with assumed boundary
conditions. These later shapes are known as assumed modes and can be particularly efficient choices for
the basis functions. The q
i
become state variables, along with their time derivatives.
The arm joint variables also can be used as state variables. The selection of joint variables and deflection
variables and their corresponding shape functions must be compatible. For example, if the joint angle is
measured from tangent of the preceding beam link to the tangent of the following beam link, called link j,a
clamped boundary condition is appropriate for the first end of link j . In order to permit relative deflection
and rotation of the second end of link j, a free boundary condition should be employed. Solutions of the
partial differential equation for a beam with clamped-free boundary conditions are thus a good choice for
basis shapes. This combination is convenient because the joint variables of the model are the same as a
joint angle sensor would measure.
For other applications of the arm model, it may be preferable to use as joint variables the angle between
lines that connect successive joint axes of the arm. In this case the joint angle accounts for the net deflection
between joints and the deflection variable accounts for the undulations of the shape between joints. In this
case pinned-pinned boundary conditions provide the necessary constraints and are good choices. This
is convenient if one is solving for the inverse dynamics of the flexible arm when moving from point to
point, since the joint angles and deflection variables can be calculated from simple geometry and static
deflections.
Due to space constraints the development of the equations of motion for flexible link manipulators is
abbreviated in this volume. The essence of the method will now be described, and the results will be given in
sufficient detail to allow them to be used. Lagrange’s method will be employed, although other techniques
are preferred by some. More detail on Lagrange’s method applied to structural systems in general may be
found in Meirovitch [12] or as applied specifically to arms with the more complex joint motions included
in Book [13]. An outline of the equation development procedure is now presented.
The kinematics of the flexible arm must first be described. This uses the joint variables and the deflection
variables tolocateeverydifferentialelement of mass, aswellas to expressitsvelocity in terms ofthederivative
of its position and to construct expressions for describing the elastic and gravitational potential energy
of the system as required by Lagrange’s method. Knowing the velocity of each particle of mass in a given
link, integration over the spatial domain of that link allows one to express the kinetic energy of that link.
Similarly, the potential energy of the link can be determined by integrating over the link. We then sum
over all the links to obtain the two scalars, kinetic energy (T) and potential energy (V) of the system as
functions of the joint and deflection variables. The simple conservative form of Lagrange’s equations then
states that
d
dt
∂T
∂
˙
q
i
−
∂T
∂q
i
+
∂V
∂q
i
= F
i
(24.49)
F
i
is the force that does work as q
i
is varied.
These equations hold for every state q
i
, both the joint and the flexible variables. Note that this simplified
form results when T is a function of the variable and its derivative, while V is a function of only the
variable and does not include specifically its derivative. Also note that if clamped boundary conditions are
employed, there is no work performed by the joint actuators when the deflection variables change. The
same is not true when pinned boundary conditions are used and a term for F
i
appears in each equation.
Copyright © 2005 by CRC Press LLC
Flexible Robot Arms 24
-21
Other hints are appropriate for the selection of good basis function shapes. Because we are often
interested in the eigenvalues of the overall system and the subsequent system mode shapes (as opposed to
the assumed mode shapes), one should choose the shapes that most naturally conform to the differential
equation and all interior conditions that form effectively the boundary conditions on the distributed
beams. In mechanics terms, we must choose shapes that meet the geometric boundary conditions and that
are complete. Assumed mode shapes have the advantage of also satisfying the partial differential equation.
If these shapes satisfied the natural boundary conditions (essentially the forces and torques on the beam),
the assumed mode shape would be the true system mode shape. Closely matching the natural boundary
conditions may be worth extra trouble if the most efficient representation is of importance.
To describe the kinematics, we again resort to coordinate transformations as described by 4×4 matrices.
W
j
i
=
x
j
component of O
i
R
j
i
y
j
component of O
i
z
j
component of O
i
000 1
(24.50)
h
i
= W
0
i
h
i
i
= W
i
h
i
i
(24.51)
We will consider transforms to be composed of separate transforms for the joints and undeformed links
and for the link (or joint) deformation
W
j
= W
j−1
E
j−1
A
j
=
W
j−1
A
j
(24.52)
where
A
j
= joint transformation matrix for joint j
E
j−1
=the link transformation matrix for link j − 1 between joints j − 1 and j
W
j−1
=the cumulative transformation from base coordinates to
O
j−1
at the distal end of link j − 1.
Because the axes of our links run in the X direction, we can describe the nominal position to consist of
an X component η and the deflection terms to be described by a series of shapes in X, Y, and Z with
amplitudes δ
ij
h
i
i
(
η
)
=
η
0
0
1
+
m
i
j=i
δ
ij
x
ij
(
η
)
y
ij
(
η
)
z
ij
(
η
)
0
(24.53)
where
x
ij
, y
ij
, z
ij
= the x
i
, y
i
, and z
i
displacement components of mode j of link i’sdeflection, respectively
δ
ij
= the time-varying amplitude of mode j of link i (24.54)
m
i
= the number of modes used to describe the deflection of link i
E
i
=
H
i
+
m
i
j=i
δ
ij
M
ij
(24.55)
where
H
i
=
1001
i
010 0
001 0
000 1
(24.56)
Copyright © 2005 by CRC Press LLC
24
-22 Robotics and Automation Handbook
and
M
ij
=
0 −θ
zi j
θ
yij
x
ij
θ
zi j
0 −θ
xij
y
ij
−θ
yij
θ
xij
0 z
ij
0000
(24.57)
dh
i
dt
=
˙
h
i
=
˙
W
j
i
h
i
+ W
i
i
˙
h
j
(24.58)
˙
W
j
=
˙
W
j
−1
A
j
+
W
j
−1
˙
A
j
(24.59)
and
¨
W
j
=
¨
W
j−1
A
j
+ 2
˙
W
j−1
˙
A
j
+
W
j−1
¨
A
j
(24.60)
where
˙
A
j
=U
j
˙
q
j
(24.61)
¨
A
j
=U
2 j
˙
q
2
j
+U
j
¨
q
j
(24.62)
U
j
= ∂ A
j
/∂q
j
(24.63)
U
2 j
= ∂
2
A
j
/∂q
2
j
(24.64)
q
j
is the joint variable of joint j.
W
j
= W
j
E
j
˙
W
j
=
˙
W
j
E
j
+ W
j
˙
E
j
¨
W
j
=
¨
W
j
E
j
+ 2
˙
W
j
˙
E
j
+ W
j
¨
E
j
˙
E
j
=
m
j
k=1
˙
δ
jk
M
jk
(24.65)
and
¨
E
j
=
m
j
k=1
¨
δ
jk
M
jk
With expressions for the velocity of any point on the slender beam, one can form the kinetic energy of a
differential mass element at that point. Instead of the simpler expression v
T
v the trace operation Tr{vv
T
}
is used to enable later simplifications by exchanging the order of summation and integration.
dk
i
=
1
2
dmTr
˙
W
i
h
i
i
h
i
i
T
˙
W
T
i
+ 2
˙
W
i
h
i
i
˙
h
i
i
T
W
T
i
+ W
i
i
˙
h
i
i
˙
h
i
i
T
W
T
i
(24.66)
where
˙
h
i
i
=
m
j
j=1
˙
δ
ij
[x
ij
y
ij
z
ij
0]
T
(24.67)
The system kinetic energy is found by integrating over each link and summing over all links.
K =
n
i=l
1
i
0
dk
i
(24.68)
Copyright © 2005 by CRC Press LLC
Flexible Robot Arms 24
-23
Careful examination of the expanded form of this equation shows that the following terms appear.
B
li
=
m
1
j=1
m
1
k=1
˙
δ
ij
˙
δ
ik
C
ikj
(24.69)
where
C
ikj
=
1
2
1
1
0
µ[x
ik
y
ik
z
ik
0]
T
[x
ij
y
ij
z
ij
0]dη
(24.70)
isa 4 ×4matrix,it isnonzeroonlyin the 3 ×3 upper-left-handcorner.Itcan also beshown that C
ikj
= C
T
ijk
.
Other terms include
B
2i
=
m
i
j=1
˙
δ
ij
C
ij
+
m
i
k=1
m
i
j=1
δ
ik
˙
δ
ij
C
ikj
(24.71)
where
C
ij
=
1
2
1
i
0
µ[η 001]
T
[x
ij
y
ij
z
ij
0]dη (24.72)
B
3i
= C
i
+
m
i
j=1
δ
ij
C
ik
+C
T
ik
+
m
i
k=1
m
i
j=1
δ
ik
δ
ij
C
ikl
(24.73)
C
i
=
1
2
1
i
0
µ[η 001]
T
[η 001]dη (24.74)
D
ik
= C
ik
+
m
1
1=1
δ
il
C
ilk
(24.75)
G
i
= C
i
+
m
i
k=1
δ
ik
C
ik
+C
T
ik
(24.76)
Lagrange’s equations will require the calculation of derivatives of the kinetic energy
d
dt
(∂ K /∂
˙
q
j
) −∂ K /∂q
j
= 2
n
i=j
Tr
∂W
i
∂q
j
G
i
¨
W
T
i
+
m
i
k=1
¨
δ
ik
D
ik
W
T
i
+ 2
m
i
k=1
˙
δ
ik
D
ik
˙
W
T
i
(24.77)
d
dt
(∂ K /∂
˙
δ
jf
) −∂ K /∂δ
jf
= 2
n
i=j +1
Tr
∂W
i
∂δ
jf
G
i
¨
W
T
i
+
m
i
k=1
¨
δ
ik
D
ik
W
T
i
+ 2
m
i
k=1
˙
δ
ik
D
ik
˙
W
T
i
+Tr
2
˙
W
j
D
jk
+ 2
˙
W
j
m
i
k=1
˙
δ
jk
C
jkf
+ W
j
m
i
k=1
¨
δ
jk
C
jkf
W
T
j
(24.78)
E = Young’s modulus of elasticity of the material.
G = the shear modulus of the material.
I
x
= the polar area moment of inertia of the link’s cross section about the neutral axis.
I
y
, I
z
= the area moment of inertia of the link’s cross section about the y
i
and z
i
axes, respectively.
Here θ
xik
is the angle about the x
i
axis corresponding to the kth mode of link iat the point η.
The potential energy terms in Lagrange’s equations also require integration over each link and summa-
tion over all the links. Integration of the terms results in a “modal stiffness” matrix that is not dependent
Copyright © 2005 by CRC Press LLC
24
-24 Robotics and Automation Handbook
on amplitude.
K
ikl
= K
xikl
+ K
yikl
+ K
zikl
(24.79)
where
K
xikl
=
1
i
0
GI
x
(
η
)
∂θ
xil
∂η
∂θ
xik
∂η
dη
(24.80)
K
yikl
=
1
i
0
EI
y
(
η
)
∂θ
yil
∂η
∂θ
yik
∂η
dη
(24.81)
K
zikl
=
1
i
0
EI
z
(
η
)
∂θ
zil
∂η
∂θ
zik
∂η
dη
(24.82)
Lagrange’s equations then require the calculation of partials of the potential energy
∂V
e
∂q
j
= 0 (24.83)
∂V
e
∂δ
jf
=
m
i
k=1
δ
jk
K
jkf
(24.84)
The values of the coefficients K
jkf
can be determined analytically or numerically, for example, by finite-
element methods.
The potential energy due to gravity requires that a gravity vector be defined
g
T
= [g
x
g
y
g
z
0] (24.85)
This is then used in the partials of gravitational potential energy
∂V
g
∂q
j
=−g
T
n
i=j
∂W
i
∂q
j
r
i
(24.86)
∂V
g
∂δ
jf
=−g
T
n
i=j +1
∂W
i
∂δ
jf
r
i
− g
T
W
j
ε
jf
(24.87)
For j = n,
∂V
g
∂δ
nf
=−g
T
W
n
ε
nf
(24.88)
where
r
i
= M
i
r
ri
+
m
i
k=1
δ
ik
ε
ik
M
i
is the total mass of link i
r
ri
= [r
xi
001]
a vector to the center of gravity from joint i (undeformed)
ε
ik
=
l
i
0
µ[x
ik
y
ik
z
ik
0]
T
dη
Copyright © 2005 by CRC Press LLC