行政院國家科學委員會專題研究計畫 成果報告
Timoshenko 平面剛架中含分佈質量的動力分析法 研究成果報告(精簡版)
計 畫 類 別 : 個別型
計 畫 編 號 : NSC 98-2221-E-011-087-
執 行 期 間 : 98 年 08 月 01 日至 99 年 07 月 31 日 執 行 單 位 : 國立臺灣科技大學營建工程系
計 畫 主 持 人 : 蔡相全
計畫參與人員: 碩士班研究生-兼任助理人員:陳俐穎
處 理 方 式 : 本計畫涉及專利或其他智慧財產權,2 年後可公開查詢
中 華 民 國 99 年 06 月 03 日
A Distributed-Mass Approach for Dynamic Analysis of Timoshenko Plane Frames
Hsiang-Chuan Tsai
1Department of Construction Engineering, National Taiwan University of Science and Technology, Taipei, Taiwan.
Abstract
The dynamic stiffness method is the exact method for the dynamic analysis of plane frames using the continuous-coordinate system to consider the effect of mass distribution in beam elements. In this paper, the dynamic stiffness matrices of beam elements are derived from the vibration theory of Timoshenko beams and then assembled to form the dynamic stiffness matrix of the plane frame through the direct stiffness approach. Orthogonal properties of vibration modes in Timoshenko plane frames are theoretically derived, through which the equations of motion in beam elements can be transformed into the decoupled equations of motion in terms of mode amplitudes. The dynamic stiffness matrix may create some null modes where the joints of beam element have null deformation. Unlike the Bernoulli-Euler frames, adding an interior node at the middle of the beam elements cannot normalize all the null modes of flexural poles in the Timoshenko frames. A floating interior node scheme is proposed to eliminate the null modes of flexural poles in the Timoshenko frames.
Keywords: Distributed mass; Dynamic stiffness matrix; Plane frame; Structural dynamics;
Timoshenko beam; Wittrick-Williams algorithm.
1. Introduction
Vibration of a plane frame, which is an assemblage of beam elements, is usually analyzed by a discrete-coordinate system in which the stiffness matrix is derived through the static equilibrium equation of beam elements; the mass matrix may be either formed by lumping the mass at structural joints, called the lumped-mass method, or computed from the static-deformed shape functions of beam elements, called the consistent mass method. These two approximate methods can simulate the dynamic response of lower-frequency modes only, and will create large error on the response of higher-frequency modes [1].
The exact method for the dynamic analysis of plane frames should use the
1
P.O. Box 90-130, Taipei 106, Taiwan. Fax: +8862-27376606. E-mail: [email protected]
continuous-coordinate system to consider the effect of mass distribution in the beam elements, which can be achieved by using the dynamic stiffness method. From the differential equation of beam vibration, the relationships between harmonic forces and displacements at the ends of a beam element can be exactly established. These relations are usually referred to as the dynamic stiffness matrix. After the dynamic stiffness matrices of beam elements are derived, the dynamic stiffness matrix of the whole frame can be assembled by exactly the same procedure used in the static matrix structural analysis. The historical development of using the dynamic stiffness method in structural dynamics has been highlighted by Akesson [2].
There are two types of beam vibration theories. One is Bernoulli-Euler beam theory and the other is Timoshenko beam theory. Timoshenko beam theory considers the effects of shear deformation and rotational inertia, so that it is more accurate in the high-frequency response.
Cheng [3] applied the dynamic stiffness method to analyze Timoshenko plane frames, but did not include the axial deformation effect in the dynamic stiffness matrix, which makes the analysis for the frames having side sway become not robotic. Wang and Lin [4] presented the transfer matrix method to analyze multi-span Timoshenko frames, but required all spans having the same span length and the same column height.
The elements in the dynamic stiffness matrix are nonlinear functions, trigonometric and hyperbolic, of vibration frequency. The dynamic stiffness matrix can be directly employed in solving the steady-state response of the frames subjected to harmonic loadings. To solve the problems of transient response, the natural frequencies and mode shapes of the frames must be first calculated from the dynamic stiffness matrix. The calculation is usually referred to as a nonlinear eigenproblem and is different from the linear eigenproblem encountered in the dynamic analysis of discrete systems. Extending the Sturm sequence property of linear eigenproblems, the Wittrick-Williams algorithm [5,6] can count the number of vibration modes exceeded by a specified frequency from the dynamic stiffness matrix. This ensures that no natural frequencies are missed when using the bisection or other methods [7] to solve the zero determinant equation of the dynamic stiffness matrix. There are several methods to solve the mode shapes from the dynamic stiffness matrix [8,9]. A deflated matrix approach to solve the mode shapes from the dynamic stiffness matrix is proposed by Tsai [10].
The orthogonality of vibration modes is the required condition of using the vibration modes to solve the transient response of the frames. The publications on the orthogonality of vibration modes employed in the dynamic stiffness matrix are few and less rigorous [11].
Tsai [10] derives the mode orthogonality condition of the Bernoulli-Euler plane frames with distributed mass.
In the present paper, the vibration theory of Timoshenko beams is applied to derive the dynamic stiffness matrices of beam elements. Orthogonal properties between vibration modes of different natural frequencies in Timoshenko plane frames considering distributed mass are theoretically derived, through which the equations of motion in beam elements can be
transformed into the decoupled equations of motion in terms of mode amplitudes. Not every vibration mode has distinct natural frequency. Sometimes several modes may have the same natural frequency, which are called the modes of repeated roots. The decoupling process for the modes of repeated roots is derived too.
The natural frequencies of flexural vibration in a single beam of which both ends are restrained from deformation are called the flexural poles of the beam. Similarly, the natural frequencies of axial vibration in a single beam with the deformations at both ends being restrained are called the axial poles of the beam. The dynamic stiffness matrix of beam elements has a characteristic of being infinite at frequencies equal to flexural poles or axial poles of the beam element. In the vibration of frames considering the effect of distributed mass, some vibration modes may have the natural frequency equal to the axial or flexural poles of some beam elements in the frame, called the null modes here. In the null modes, the deformations at all frame joints are null but the deformations in beam elements are not null.
The null mode shapes cannot be directly solved from the dynamic stiffness matrix. Dias and Alves [12] included the degrees of freedom at the joints of beam elements and the constants in the shape functions of beam elements in the eigenproblem, which can eliminate the null modes but has to pay the price of an increased matrix order. For the frame possessing the null modes of axial poles, the mode shape can be calculated by using the approach of nodal force equilibrium [10]. For the Bernoulli-Euler frame possessing the null modes of flexural poles, adding an interior node at the middle of the beam elements owning the flexural poles can normalized the null mode to become a regular mode [9,10]. The approach to normalize the null modes of flexural poles in the Timoshenko frames is studied in the present paper. Unlike the Bernoulli-Euler frames, adding an interior node at the fixed position of the beam elements cannot normalize all the null modes of flexural poles in the Timoshenko frames. A floating interior node scheme is proposed in this paper to eliminate the null modes of flexural poles in the Timoshenko frames.
2. Flexural dynamic stiffness of beam elements
The prismatic beam element shown in Fig. 1 has length L, cross-sectional area A, moment of area I, Young’s modulus E, shear modulus G, shear coefficient k, and density
ρ
. The x axis is the coordinate along the beam axis. Let v denote the transverse displacement of the beam in the y direction, and θ denote the rotational angle of the beam in the counter-clockwise direction. The equations of motion for the transverse free vibration of the Timoshenko beam, which include the effects of rotational inertia and shear deformation, are0 )
( 2
2
2
2 =
∂
− ∂
∂ + ∂
∂
∂
x v kAG x
t
A v θ
ρ
(1)0 )
2 (
2
2
2 =
∂
−∂
∂ +
− ∂
∂
∂
x kAG v
EI x
I t θ θ θ
ρ
(2)with time t. When the beam is vibrating at a specific frequency
ω
, the transverse displacement and rotational angle can be expressed ast
e
ix t
x
v
( , )=ψ
( ) ω (3)t
e
ix t
x φ
ωθ
( , )= ( ) (4)in which
ψ
andφ
are the transverse and rotational shape functions, respectively. Eqs. (1) and (2) becomes0 )
( 2
2
2 + − =
−
d x
d x d kAG d
A ω ψ φ ψ
ρ
(5)0 )
2 (
2
2 − + − =
−
d x
kAG d x
d EI d
I ω φ φ φ ψ
ρ
(6)Differentiate Eq. (6) with respect to x and then bring in
d φ dx
from Eq. (5), which yields 0) 1
( )
( 2 2 2
2 2
4
4
ψ
+ω ρ
+ψ
−ω ρ
−ω ρ ψ
=kAG I EI
A x
d d kAG
EI A I EI
A x
d
d
(7)Differentiate Eq. (5) with respect to x and then bring in
d ψ dx
from Eq. (6), which yields 0) 1
( )
( 2 2 2
2 2
4
4
φ
+ω ρ
+φ
−ω ρ
−ω ρ φ
=kAG I EI
A x
d d kAG
EI A I EI
A x
d
d
(8)Define the following dimensionless parameters
EI
L A
b
=ω
2ρ
;A I r L
1= ;
kAG EI s L
1= (9)
If
brs < 1
, the solutions ofψ
in Eq. (7) andφ
in Eq. (8) have forms as ) sin() cos(
) sinh(
) cosh(
)
( 1 2 3 4
L b x L c
b x L c
b x L c
b x c
x α α β β
ψ
= + + + (10))
~ sin(
)
~ cos(
) sinh(
) ~ cosh(
) ~
( 1 2 3 4
L b x L c
b x L c
b x L c
b x c
x α α β β
φ
= + + + (11)with
2 2 2 2 2
2
4
) (
) 2 (
1
s b r s
r + + − +
−
α
=
(12)2 2 2 2 2
2
4
) (
) 2 (
1
s b r s
r + + − +
β
=
(13)Substitution of Eqs. (10) and (11) into Eq. (7) gives the following relations between the constants
c
i andc~
i2 1
~
c
c L αβ
=
λ
, ~2c
1c L
αβ
=
λ
, ~3c
4c L
αβ
=
μ
, ~4c
3c L
αβ
−
μ
= (14)
with
) ( 2
s
2b
+=
β α
λ
,μ
=α b
(β
2−s
2) (15)The constants can be determined from the transverse displacements and rotational angles at the beam ends,
c
i⎪ ⎪
⎭
⎪ ⎪
⎬
⎫
⎪ ⎪
⎩
⎪ ⎪
⎨
⎧
⎥ ⎥
⎥ ⎥
⎦
⎤
⎢ ⎢
⎢ ⎢
⎣
⎡
−
−
−
−
−
−
=
⎪ ⎪
⎭
⎪ ⎪
⎬
⎫
⎪ ⎪
⎩
⎪ ⎪
⎨
⎧
) (
) (
) 0 (
) 0 ( 1
7 5
1 4
6 7
3 2
7 5
2 4
6 7
3 1
0
4 3 2 1
L L LB
B LB
B
LB B
LB B
LB B
LB B
LB B
LB B
B c c c c
φ ψ φ ψ
λαβ λ
αβ λ
αβ λμ
αβ λ
μαβ μ
αβ μ
αβ λμ
αβ μ
(16)
with
) sin(
) sinh(
) (
)]
cos(
) cosh(
1 [
2
2 20
b b b b
B =
αβ−
α β+
λ−
μ α β) sin(
) sinh(
)]
cos(
) cosh(
1
1 [
b b b b
B
=λ
−α β
−μ α β
) sin(
) sinh(
)]
cos(
) cosh(
1
2 [
b b b b
B
=μ
−α β
+λ α β
) cos(
) sinh(
) sin(
) cosh(
3
b b b b
B
=λ α β
−μ α β
) cos(
) sinh(
) sin(
) cosh(
4
b b b b
B
=μ α β
+λ α β
) sin(
) sinh(
5
b b
B
=λ α
+μ β
) sin(
) sinh(
6
b b
B
=−μ α
+λ β
) cos() cosh(
7
b b
B
=α
−β
The bending moment M and the shear force V in the prismatic Timoshenko beam have the following relations with the transverse displacement and rotational angle
) ( )
(
x EI x
M
=φ
′ ;V
(x
)=kAG
[φ(x
)−ψ′(x
)] (17)Let
M
0,V
0,ψ
0 andφ
0 denote the bending moment, shear force, transverse displacement and rotational angle, respectively, at the beam end ofx = 0
, andM
1,V
1,ψ
1 andφ
1 be the corresponding quantities at the beam end ofx = L
. The positive directions of the end forces and internal forces are shown in Fig. 1, which give the relations) 0
0
M
(M
=− ;V
0 =V
(0) ;ψ
0 =ψ
(0) ;φ
0 =φ
(0) (18) )1
M
(L
M
= ;V
1 =−V
(L
) ;ψ
1 =ψ
(L
) ;φ
1 =φ
(L
) (19) Substituting Eqs. (10) and (11) into Eq. (17) and using the relations in Eqs. (18) and (19), theflexural dynamic stiffness of Timoshenko beams for
brs < 1
is derived as⎪ ⎪
⎭
⎪ ⎪
⎬
⎫
⎪ ⎪
⎩
⎪ ⎪
⎨
⎧
⎥ ⎥
⎥ ⎥
⎦
⎤
⎢ ⎢
⎢ ⎢
⎣
⎡
−
−
−
−
−
−
−
−
=
⎪ ⎪
⎭
⎪ ⎪
⎬
⎫
⎪ ⎪
⎩
⎪ ⎪
⎨
⎧
1 1 0 0
5 2 2 6
2 4
2 1
4 3
6 2 4 5
2 2
4 3
2 1
0 3
1 1
0 0
φ ψ φ ψ
K L LK K
L LK
LK K
LK K
K L LK K
L LK
LK K
LK K
B L
EI
M V M V
(20)
with
)]
cos(
) sinh(
) sin(
) cosh(
)[
( 2 2
2
1
b b b b b
K
=α
+β μ α β
+λ α β
) sin(
) sinh(
)]
( 2
[ )]
cos(
) cosh(
1 )[
( 2 2 2 2 2 2
2
2
b r s b b b s s r b b
K
=αβ
− −α β
+ + −α β
)]
sin(
) sinh(
)[
(
2 22
3
b b b
K =
α+
β λ α+
μ β )]cos(
) )[cosh(
( 2 2
2
4
b b b
K
=αβ α
+β α
−β
)]
cos(
) sinh(
) sin(
) cosh(
)[
(
2 22
5
b b b b b
K =
αβ α+
β λ α β−
μ α β)]
sin(
) sinh(
)[
(
2 22
6
b b b
K =
αβ α+
β−
μ α+
λ βIf
brs > 1
, the solutions ofψ
in Eq. (7) andφ
in Eq. (8) have forms as )sin(
) cos(
) sin(
) cos(
)
( 1 2 3 4
L b x L c
b x L c
b x L c
b x c
x α α β β
ψ
= + + + (21))
~ sin(
)
~ cos(
)
~ sin(
)
~cos(
)
( 1 2 3 4
L b x L c
b x L c
b x L c
b x c
x α α β β
φ
= + + + (22)with
2 2 2 2 2
2
4
) (
) 2 (
1
s b r s
r + − − +
α =
(23)The following relations between the constants
c
i andc~ for
i can be derived from Eq. (7)> 1 brs
2 1
~
c
c L
β α
−
λ
= , ~2
c
1c L
β α
=
λ
, ~3c
4c L
β α
=
μ
, ~4c
3c L
β α
−
μ
= (24)
with
) ( 2
s
2b
− +=
β α
λ
,μ
=α b
(β
2−s
2) (25)The constants
c
i forbrs > 1
become⎪ ⎪
⎭
⎪ ⎪
⎬
⎫
⎪ ⎪
⎩
⎪ ⎪
⎨
⎧
⎥ ⎥
⎥ ⎥
⎦
⎤
⎢ ⎢
⎢ ⎢
⎣
⎡
−
−
−
−
−
−
−
−
=
⎪ ⎪
⎭
⎪ ⎪
⎬
⎫
⎪ ⎪
⎩
⎪ ⎪
⎨
⎧
) (
) (
) 0 (
) 0 ( 1
7 5
1 4
6 7
3 2
7 5
2 4
6 7
3 1
0
4 3 2 1
L L LB
B LB
B
LB B
LB B
LB B
LB B
LB B
LB B
B c c c c
φ ψ φ ψ
β α λ λ
β α λ
β α μ
λ β
α λ
β α μ μ
β α μ
β α μ
λ β
α μ
(26)
with
) sin(
) sin(
) (
)]
cos(
) cos(
1 [
2
2 20
b b b b
B =
αβ−
α β+
λ+
μ α β) sin(
) sin(
)]
cos(
) cos(
1
1 [
b b b b
B
=λ
−α β
+μ α β
) sin(
) sin(
)]
cos(
) cos(
1
2 [
b b b b
B
=μ
−α β
+λ α β
) cos(
) sin(
) sin(
)
3
cos( b b b b
B =
λ α β+
μ α β) cos(
) sin(
) sin(
)
4 cos(
b b b b
B
=μ α β
+λ α β
) sin(
)
5
sin( b b
B =
λ α+
μ β) sin(
)
6
sin( b b
B =
μ α+
λ β ) cos()
7 cos(
b b
B
=α
−β
The flexural dynamic stiffness of Timoshenko beams for has the same form as Eq.
(20) with different definition for
> 1 brs
K
i)]
cos(
) sin(
) sin(
) cos(
)[
( 2 2
2
1
b b b b b
K
=β
−α μ α β
+λ α β
) sin(
) sin(
)]
( 2
[ )]
cos(
) cos(
1 )[
( 2 2 2 2 2 2
2
2
b r s b b b s s r b b
K
=α β
− −α β
+ + −α β
)]
sin(
) sin(
)[
(
2 22
3
b b b
K =
β−
α λ α+
μ β )]cos(
) )[cos(
( 2 2
2
4
b b b
K
=α β β
−α α
−β
)]
cos(
) sin(
) sin(
) cos(
)[
(
2 22
5
b b b b b
K =
αβ β−
α λ α β+
μ α β)]
sin(
) sin(
)[
(
2 22
6
b b b
K =
αβ β−
α μ α+
λ β3. Complete dynamic stiffness of beam elements
In order to be able to robotically form the dynamic stiffness of the frame, the axial stiffness must be included in the dynamic stiffness of beam elements. Let u denote the axial displacement of the beam in the x direction. The equation of motion for the axial free vibration is
2 0
2
2
2 =
∂ ρ ∂
∂ +
− ∂
t A u x
EA u
(27)When the beam is vibrating at a specific frequency
ω
, the axial displacement can be expressed ast
e
ix t
x
u
( , )=χ( ) ω (28)in which χ is the axial shape function. Eq. (27) becomes
2 0
2
2χ +ω ρχ=
E x
d
d
(29)The solution of
χ
in Eq. (29) has a form as )sin(
) cos(
)
( 1 2
L rb x L a
rb x a
x
= +χ (30)
in which the constants
a
i can be determined from the axial displacement at the beam ends,⎭ ⎬
⎫
⎩ ⎨
⎥ ⎧
⎦
⎢ ⎤
⎣
⎡
= −
⎭ ⎬
⎫
⎩ ⎨
⎧
) (
) 0 ( ) sin(
1 ) tan(
1
0 1
2 1
L rb
rb a
a
χ
χ
(31)The axial force N in the prismatic beam has the following relation with the axial displacement,
) ( )
(
x EA x
N
= χ′ (32)Let and denote the axial force and axial displacement, respectively, at the beam end of , and and be the corresponding quantities at the beam end of . The positive directions of the end forces and internal forces in Fig. 1 give the relations
N
0 χ0= 0
x N
1 χ1x = L
) 0
0
N
(N
=− ; χ0 =χ(0) ;N
1=N
(L
) ; χ1 =χ(L
) (33) Substituting Eq. (30) into Eq. (32) and using the relations in Eq. (33), the axial dynamicstiffness is derived as
⎭ ⎬
⎫
⎩ ⎨
⎥ ⎧
⎦
⎢ ⎤
⎣
⎡
−
= −
⎭ ⎬
⎫
⎩ ⎨
⎧
1 0
1 0
) cos(
1
1 )
cos(
)
sin( χ
χ rb rb
rb rb L EA N
N
(34)The complete dynamic stiffness of Timoshenko beams can be obtained by combining Eqs.
(20) and (34),
⎪⎪
⎪⎪
⎭
⎪⎪⎪
⎪
⎬
⎫
⎪⎪
⎪⎪
⎩
⎪⎪⎪
⎪
⎨
⎧
⎥⎥
⎥⎥
⎥⎥
⎥⎥
⎥⎥
⎥⎥
⎥⎥
⎥
⎦
⎤
⎢⎢
⎢⎢
⎢⎢
⎢⎢
⎢⎢
⎢⎢
⎢⎢
⎢
⎣
⎡
−
−
−
−
−
−
−
−
=
⎪⎪
⎪⎪
⎭
⎪⎪⎪
⎪
⎬
⎫
⎪⎪
⎪⎪
⎩
⎪⎪⎪
⎪
⎨
⎧
1 1 1 0 0 0
0 5 2
0 2
0 6 2
0 4
0 2
0 1
0 4
0 3
0 6 2
0 4
0 5 2
0 2
0 4
0 3
0 2
0 1
3
1 1
1 0 0
0
0 0
0 0
0 ) 0
0 tan(
) 0 sin(
0 0
0 0
0 ) 0
0 sin(
) 0 tan(
φ ψ χ φ ψ χ
B K L B LK B
K L B LK
B LK B
K B
LK B
K
rb r
b rb
r b
B K L B LK B
K L B LK
B LK B
K B
LK B
K
rb r
b rb
r b
L EI
M V N M V N
(35)
4. Dynamic stiffness of plane frames
The plane frame shown in Fig. 2 is defined by a global coordinate system (
x
,y
). Letu ,
iv and
iθ
i denote the displacements in the x and y directions and the counter-clockwise rotation at node i, respectively. The total number of the unconstrained degrees of freedom in the joints of the frame is denoted asn
d. The element J has two ends at nodes i and j. LetX ,
iY and
iM denote the forces in the x and y directions and the moment acting at the
i node i of element J, respectively. The transformation between the element local coordinates and the frame global coordinates gives⎪ ⎭
⎪ ⎬
⎫
⎪ ⎩
⎪ ⎨
⎧
θ
⎪ =
⎭
⎪ ⎬
⎫
⎪ ⎩
⎪ ⎨
⎧
φ ψ χ
i i i
v u T
0 0 0
; ⎪
⎭
⎪⎬
⎫
⎪⎩
⎪⎨
⎧ θ
⎪=
⎭
⎪⎬
⎫
⎪⎩
⎪⎨
⎧ φ ψ χ
j j j
v u T
1 1 1
; ⎪
⎭
⎪⎬
⎫
⎪⎩
⎪⎨
⎧
⎪=
⎭
⎪⎬
⎫
⎪⎩
⎪⎨
⎧
i i
i
M Y X
M V N
T
0 0
0
; ⎪
⎭
⎪⎬
⎫
⎪⎩
⎪⎨
⎧
⎪=
⎭
⎪⎬
⎫
⎪⎩
⎪⎨
⎧
j j
j
M Y X
M V N
T
1 1
1
(36)
where
T
is the transformation matrix defined as⎥⎥
⎥
⎦
⎤
⎢⎢
⎢
⎣
⎡
−
−
−
−
−
=
1 0
0
0 ) (
) (
0 ) (
) (
L x x L y y
L y y L x x
i j i
j
i j i
j
T
(37)with (
x
i,y
i) being the coordinates of the joint i, and L being the element length calculated by2
2 ( )
)
(
x
jx
iy
jy
iL
= − + − (38)The dynamic stiffness in Eq. (35) can be transformed to the element stiffness in the global coordinate system by using Eq. (36). The global stiffness matrices of all elements in the frame are then assembled to form the frame stiffness matrix
K
(ω
) by means of the direct stiffness method, which yields0 d
K
(ω
) = (39)in which
K
(ω
) is a matrix and is a function of frequency, and is the displacement vector of the unconstrained degrees of freedom in the frame. Without any external force, the summation of the internal nodal forces leads to the zero vector at the right side. Natural frequencies and shape vectors of vibration modes can be solved from Eq. (39).d
d
n
n
×d
5. Orthogonality of mode shape functions
After the nth mode shape vector of the frame is solved, the axial, transverse and rotational shape functions of element J in the nth mode,
d
nχ
Jn,ψ
Jn andφ
Jn respectively, can be derived from Eqs. (10), (11) and (30), or Eqs. (21) and (22). Similar to Eqs. (1), (2) and (27), the equations of motion for element J vibrating in the nth mode with frequencyω
n are2 0
2
2 − =
−
J Jn J J Jn n J
J
d x
A d E
A ω χ χ
ρ
(40)0 )
( 2
2
2 + − =
−
J J
J J J J J Jn n J
J
d x
d x d G d A k
A ω ψ φ ψ
ρ
(41)0 )
2 (
2
2 − + − =
−
J J J
J J J J Jn J J Jn n J
J
d x
G d A x k
d I d E
I ω φ φ φ ψ
ρ
(42)in which the parameters with the subscript J indicates the corresponding properties of element J. After multiplying Eqs. (40), (41) and (42) by the axial, transverse and rotational shape functions of element J in the mth mode, respectively, and then integrating the summation of the three equations through the element length , the summation of the integrations of all the elements leads to
L
J0 )]
( [
)]
( [
] [
0
2
0
2 0
2
′ =
−
′′ − +
+
− ′′
− ′
′′ + +
∑∫
∑∫
∑∫
J L
J Jn Jn J J J Jn J J Jn n J J Jm
J L
J Jn Jn J J J Jn n J J Jm J
L
J Jn J J Jn n J J Jm
J
J J
dx G
A k I
E I
dx G
A k A
dx A
E A
ψ φ φ
φ ω ρ φ
ψ φ ψ
ω ρ ψ χ
χ ω ρ χ
(43)
which becomes, through integration by part,
{ }
∑
∑∫
∑∫
+
−
−
′ = + ′
− ′
− ′
′ +
− ′
+ +
J
L Jn Jm Jn Jm Jn Jm J
L
J Jn Jm J J Jn Jn Jm Jm J J J Jn Jm J J
J L
J Jn Jm J J Jn Jm Jn Jm J J n
J J
J
M V
N
dx I
E G
A k A
E
dx I
A
0 0
0 2
) (
) )(
(
] )
( [
φ ψ
χ
φ φ ψ
φ ψ φ χ
χ
φ φ ρ ψ ψ χ χ ρ ω
(44)
The right side of Eq. (44) can be expressed in vector form by using the notation similar to Eqs. (18), (19) and (33),
∑
∑ ⎟⎟
⎟ ⎟
⎠
⎞
⎜⎜
⎜ ⎜
⎝
⎛
⎪ ⎭
⎪ ⎬
⎫
⎪ ⎩
⎪ ⎨
⎧
⎪ ⎭
⎪ ⎬
⎫
⎪ ⎩
⎪ ⎨
⎧
⎪ +
⎭
⎪ ⎬
⎫
⎪ ⎩
⎪ ⎨
⎧
⎪ ⎭
⎪ ⎬
⎫
⎪ ⎩
⎪ ⎨
⎧
= +
−
J
Jn Jn
Jn T
Jm Jm Jm
Jn Jn
Jn T
Jm Jm Jm
J
L Jn Jm Jn Jm Jn Jm
M V N
M V N M
V
N
J) (
) (
) (
) (
) (
) (
) (
) (
) (
) (
) (
) ( )
(
0 0
0
0 0 0
1 1
1
1 1 1
0 φ
ψ χ φ
ψ χ φ
ψ
χ (45)
The coordinate transformation in Eq. (36) gives
⎪ ⎭
⎪ ⎬
⎫
⎪ ⎩
⎪ ⎨
⎧
⎪ ⎭
⎪ ⎬
⎫
⎪ ⎩
⎪ ⎨
⎧
⎪ =
⎭
⎪ ⎬
⎫
⎪ ⎩
⎪ ⎨
⎧
⎪ ⎭
⎪ ⎬
⎫
⎪ ⎩
⎪ ⎨
⎧
⎪ =
⎭
⎪ ⎬
⎫
⎪ ⎩
⎪ ⎨
⎧
⎪ ⎭
⎪ ⎬
⎫
⎪ ⎩
⎪ ⎨
⎧
Jn i
Jn i
Jn i T
m i
m i
m i
Jn i
Jn i
Jn i T T
m i
m i
m i
Jn Jn
Jn T
Jm Jm Jm
M Y X v
u
M Y X v
u
M V N
) (
) (
) (
) (
) (
) (
) (
) (
) (
) (
) (
) (
) (
) (
) (
) (
) (
) (
0 0
0
0 0 0
θ θ
φ ψ χ
T
T
(46)which means
0 )
( ) (
) (
) (
) (
) (
) (
) (
) (
) (
) (
) (
) (
) (
) (
) (
) (
) (
0 0
0
0 0 0
1 1
1
1 1 1
=
⎟⎟
⎟ ⎟
⎠
⎞
⎜⎜
⎜ ⎜
⎝
⎛
⎪ ⎭
⎪ ⎬
⎫
⎪ ⎩
⎪ ⎨
⎧
⎪⎭
⎪ ⎬
⎫
⎪⎩
⎪ ⎨
⎧
=
⎟⎟
⎟ ⎟
⎠
⎞
⎜⎜
⎜ ⎜
⎝
⎛
⎪⎭
⎪ ⎬
⎫
⎪⎩
⎪ ⎨
⎧
⎪⎭
⎪ ⎬
⎫
⎪⎩
⎪ ⎨
⎧
⎪⎭ +
⎪ ⎬
⎫
⎪⎩
⎪ ⎨
⎧
⎪⎭
⎪ ⎬
⎫
⎪⎩
⎪ ⎨
⎧ ∑ ∑
∑
i J∈iJn i
Jn i
Jn i T
m i
m i
m i
J
Jn Jn
Jn T
Jm Jm Jm
Jn Jn
Jn T
Jm Jm Jm
M Y X v
u
M V N
M V N
θ φ
ψ χ φ
ψ χ
(47)
where the element
J ∈ i
means the element connecting to node i. Without external force acting at node i, the internal end forces of all the elements connecting to node i are balance, which leads to zero in the above equation. Using Eq. (47), Eq. (44) becomes∑∫
∑∫
′ + ′
− ′
− ′
′ +
′
= +
+
J L
J Jn Jm J J Jn Jn Jm Jm J J J Jn Jm J J
J L
J Jn Jm J J Jn Jm Jn Jm J J n
J
J
dx I
E G
A k A
E
dx I
A
0
0 2
] )
)(
( [
] )
( [
φ φ ψ
φ ψ φ χ
χ
φ φ ρ ψ ψ χ χ ρ ω
(48)
Multiplying the equations of transverse and axial motions of element J vibrating in the
mth mode, respectively, by the transverse and axial shape functions of element J in the nth
mode, the following equation can be established0 )]
( [
)]
( [
] [
0
2
0
2 0
2
′ =
−
′′ − +
+
− ′′
− ′
′′ + +
∑∫
∑∫
∑∫
J L
J Jm Jm J J J Jm J J Jm m J J Jn
J L
J Jm Jm J J J Jm m J J Jn J
L
J Jm J J Jm m J J Jn
J
J J
dx G
A k I
E I
dx G
A k A
dx A
E A
ψ φ φ
φ ω ρ φ
ψ φ ψ
ω ρ ψ χ
χ ω ρ χ
(49)
Using a similar procedure described in the last paragraph, the following equation is derived