## Abstract

A major advantage of the screw theory is that translations and rotations are treated simultaneously, which can provide greater insight into the vibration phenomena, such as vibration centers and axes. The present study describes how these concepts are extended into beam theory. The stiffness matrix of a beam was derived by incorporating different types of vibrations, such as extension, compression, torsion, and bending, at the same time, which was then used to obtain the equations of motion, including nonstandard forms of boundary conditions. We then presented an analytical method to solve these equations by focusing on two distinct examples, namely the cantilever and robot link. In the first numerical example, the mode shapes of the beam could be regarded as rotations about the vibration centers or axes of the rigid bodies in a discrete system. In the second example, the analytical solutions of mode shapes and natural frequencies of a robot link, for which the revolute joints at both sides are not parallel, were presented to demonstrate the utility of the screw theory. We demonstrated that the screw approach could accurately describe the vibrations of both discrete and continuous systems and that the geometric meaning of the vibration modes of discrete systems can be extended into continuous systems.

## 1 Introduction

The screw theory is often utilized to interpret spatial mechanisms and robotics issues at a high level. The major advantage of this theory is that translations and rotations are considered simultaneously, thus alleviating the complexity of these systems when implementing algebraic and numerical options. From this point of view, the screw theory may be useful for describing vibrations of beams with nonstandard boundary conditions. Previously, Selig and Ding [1,2] applied this theory to derive the equations of static and dynamic beams and identified the natural frequencies and mode shapes for Timoshenko beams with revolute joints at both ends that were not parallel. However, their findings were obtained by using numerical solutions as opposed to implementing an analytical method.

Several research studies have been performed using the screw theory. Ball [3] was the first to describe the oscillation of a rigid body as a repetitive small screw motion. The screw representation has also been used to analyze and design vibration systems. For example, Blanchet [4] introduced the concept of vibration centers and interpreted the mode shape of a planar three-degree-of-freedom (DOF) system as a pure repetitive rotation about the vibration center. That study also showed that the orthocenter of a triangle formed by three vibration centers coincides with its mass center. Some research studies [57] utilized transparent geometric constraints to develop design methods for energy harvesters, which focused on uniform power peaks at equally spaced resonant frequencies to widen the effective bandwidth. However, the nature of spatial vibration systems is so complicated that this approach is often impractical for these systems. One of the ways to alleviate the complexity of spatial vibration systems is to employ the decoupling technique. For instance, decoupling with respect to the plane of symmetry was utilized to design engine mounts, absorbers, and energy harvesters [810]. In addition to decoupling methods, Lee et al. [11] found three geometric constraints that could be applied to practical design methods. One of the constraints was used to design a spatial vibration system with prespecified ratios of peaks and target frequencies [12]. Furthermore, Ryu et al. [13] investigated the geometrical properties of a dual-body system and developed a design method for prespecified ratios among six energy peaks using a planar symmetric dual-body system.

Various analytical methods for nonstandard beams were previously proposed. A few studies have investigated the dynamic behavior of a damaged beam [1418]. For example, Čas et al. [19] considered a three-dimensional two-layer composite beam with interlayer slips. Moreover, Liu and Yang [20] presented an analytical method to describe elastically connected double-beam systems, while Rezaiee-Pajand et al. [21] derived the exact elemental stiffness matrix and calculated an accurate solution for a nonprismatic beam–column.

In the present work, we used the screw approach to describe the vibrations of discrete and continuous systems. Using our approach, analytical solutions for nonstandard beams can be obtained, and our findings reveal that the geometric meaning of the vibration modes of discrete systems can be extended into continuous systems. This paper is organized as follows: Sec. 2 provides a geometrical interpretation of the vibrations of a rigid body. In Sec. 3, the stiffness matrix of the beam expressed in terms of the screw theory is introduced, and the equation of motion (EOM) governing the vibrations of the beams is derived, including the boundary conditions. Analytical solutions of the differential equations for a cantilever are then presented. Section 4 describes the discretization of a continuous system into an equivalent multibody system with nonstandard boundary conditions. In Sec. 5, we demonstrate the geometric meaning of vibration modes by comparing the results of continuous and discrete analysis, and the robot link is selected as the second numerical example to demonstrate the utility of the screw theory. Finally, Sec. 6 presents our conclusions and future research directions.

## 2 Screw Theory of Vibrations

In this section, we describe the screw theory of vibrations. Using the expression, the geometrical properties of single and multibody systems, such as vibration centers and axes, are described via the screw theory.

### 2.1 Single Body.

For a rigid body supported by a number of springs, the EOM for undamped free vibration is given by
$MD¨+KD=0$
(1)
where $M∈R6×6$ and $K∈R6×6$ denote mass and stiffness matrices, respectively, and $D∈R6×1$ is the time-dependent displacement vector. The stiffness matrix can be written as [2224]
$K=[s^1⋯s^n][k10⋱0kn][s^1⋯s^n]T=∑i=1nkis^is^iT$
(2)
where ki is the ith spring constant, and $s^i∈R6×1$ is a line vector indicating the axis of the ith spring. The line vector is represented by $s^i=[siT;riT×siT]T$, where $si$ is a unit direction vector, and $ri$ is a position vector to a point on the axis of the spring.
For a small harmonic displacement, the displacement vector $D$ can be expressed in terms of a screw as
$D=D^ejωt$
(3)
where the time-independent displacement vector $D^(=[δT;ϕT]T)$ is a screw expressed in Plücker's axis coordinates. $δ(∈R3×1)$ is a small translational vector of a point on the rigid body that coincides with the origin of the coordinate frame, and ϕ( ∈ R3×1) is a small angular vector. Substituting Eq. (3) in Eq. (1), we can obtain that
$(K−ω2M)D^=0$
(4)

Equation (4) yields six eigenvectors $D^i(i=1,…,6)$, which correspond to the vibration modes. The vibration modes are generally determined as screws, meaning that the general mode shape is reflected by screw motion in the spatial system. There are three vibration modes in the xy planar system, as shown in Fig. 1: $D^i=[Yi−Xi0001]T$ for i = 1, 2, 3. In this case, the operational mode shape is rotation about the vibration axis $D^i$ (or vibration center (Xi, Yi, 0)).

Fig. 1
Fig. 1
Close modal

### 2.2 Multibody System.

The EOM for undamped free vibration of a dual-body system can be given by
$MD¨+KD=0$
(5)
where $M=[M100M2]∈R12×12$ and $K=[K1−K12−K12K2]∈R12×12$. $Mi$ is the mass matrix, $Ki$ is the stiffness matrix corresponding to all springs attached to the ith body, and $K12$ is one of the springs connected to each other. $D(=[D1TD2T])$ is the 12 × 1 displacement vector, which can be rewritten in terms of a screw chain as [13]
$D=D~ejωt=[D^1D^2]ejωt$
(6)
where $D~$ is the screw chain composed of the two screws $D^1$ and $D^2$. Substituting Eq. (6) into Eq. (5) yields
$(K−ω2M)D~=0$
(7)

From Eq. (7), we can generally obtain 12 eigenvectors (or mode shapes) $D~i(i=1,2,…,12)$. Considering that $D~$ is the screw chain, it can be said that there are 24 screw axes of vibration, and the mode shape can be described by a general vibrational screw motion. In a planar system, there exist 12 vibration centers, and the mode shapes are rotational motions [13].

We consider a n-body vibration system. To obtain the EOM for undamped free vibration, $M$ and $K$ in Eq. (5) become
$M=[M10⋱0Mn]∈R6n×6n$
(8-1)
and
$K=[K1⋯−K1n⋮⋱⋮−K1n⋯Kn]∈R6n×6n$
(8-2)
Similarly, $D$ can be represented by a screw chain composed of n screws as $=D~ejωt$. Using the screw chain, the same form of Eq. (7) is obtained, and thus 6n mode shapes are determined, each consisting of n axes of the screw. In terms of engineering, a system comprising stacked rigid bodies, such as the beam shown in Fig. 2, may be more meaningful than that where all bodies are directly connected to each other. Since only the adjacent bodies are connected by springs, the stiffness matrix can be given by
$K=[K1−K1,20⋯0−K1,2K2−K2,3⋮0−K2,3⋱0⋮Kn−1−Kn−1,n0⋯0−Kn−1,nKn]$
(9)
Fig. 2
Fig. 2
Close modal

If we consider n → ∞, then there are infinitely many natural frequencies and mode shapes in a continuous system such as the cantilever. If the elastic properties of the beam are reflected in the stiffness matrix described by Eq. (9), we may acquire mode shapes similar to those of the beam. Our goal is to show that the dynamic characteristics of this system become equivalent to the ones of the beam as n → ∞ (i.e., each body becomes a plate since the length of the beam remains constant). In the following section, we derive the differential equations of the beam using the screw and solve the equations analytically.

## 3 Vibration of a Continuous System

This section initially derives the stiffness matrix of a beam. Using this stiffness matrix, the EOM for the vibration of a continuous system is obtained and solved analytically.

### 3.1 Stiffness Matrix of a Beam.

A uniform prismatic beam with length L, area of cross section A, Young's modulus E, shear modulus G, and density ρ is modeled by a stack of n identical pieces assumed to be rigid. To track these pieces, we set global and local coordinate frames as follows (Fig. 3):

1. The global frame (OXYZ) is placed at the foot of the beam, the origin O coincides with the centroid of the cross section, and the X- and Y-axes of the frame are aligned with the principal directions of the cross section.

2. The local frame (AiXYZ) is placed in each rigid body, which is aligned with the principal axes of the rigid body.

3. The local frame (BiXYZ) is located in each cut plane, the origin Bi coincides with the centroid of the cut plane, and the X- and Y-axes are aligned with the principal directions.

Fig. 3
Fig. 3
Close modal

Under the assumption that the beam is uniform, it is clear that the Z-axes of all frames are aligned as shown in Fig. 3.

To model the elastic behavior, we considered virtual springs with an initial length of zero on each cut plane, which are connected to the rigid bodies and consist of p sets of three orthogonal springs (Fig. 4). First, in order to determine the spring constants, we suppose that the virtual springs parallel to the Z-axis are stretched by ΔL/n, so the nth rigid body is displaced by ΔL when an axial load FZ is applied, as shown in Fig. 5. In accordance with the Hooke's law, we have
$FZ=pkZΔLn$
(10)
Fig. 4
Fig. 4
Close modal
Fig. 5
Fig. 5
Close modal
From the common normal stress–strain relation of the beam, ΔL can be expressed as
$ΔL=LFZEA$
(11)
From Eqs. (10) and (11), kZ is
$kZ=nEApL=EΔAΔZ$
(12)
where ΔZ( = L/n) is the length of the rigid body and ΔA = A/p.
Referring to Fig. 6, the rigid bodies are dislocated by the application of a shear force FX. Similarly, the shear force can be obtained from the shear stress–strain relationship as follows:
$FX=GAγ≈GA⋅δΔZ$
(13)
where δ is the sliding displacement, and thus, δ can be given by
$δ=FXpkX$
(14)
Fig. 6
Fig. 6
Close modal
Substituting Eq. (14) into Eq. (13) gives
$kX=GΔAΔZ$
(15)
In the same manner, the spring constant kY can be determined as
$kY=GΔAΔZ$
(16)
Consequently, we can derive the stiffness matrix of the beam. Using Eq. (2), the stiffness matrix of the virtual springs on the ith cut plane can be obtained as follows:
$KBi=∑q=1p(kXs^X,qs^X,qT+kYs^Y,qs^Y,qT+kZs^Z,qs^Z,qT)$
(17)
where Bi indicates that the matrix is expressed at BiXYZ. Referring to Fig. 4, $s^X,q$, $s^Y,q$, and $s^Z,q$, which are orthogonal three springs illustrated enlarged, are defined as follows:
$s^X,q=[sX,qT;rqT×sX,qT]T=[100;00−Yq]T$
(18-1)
$s^Y,q=[sY,qT;rqT×sY,qT]T=[010;00Xq]T$
(18-2)
$s^Z,q=[sZ,qT;rqT×sZ,qT]T=[001;Yq−Xq0]T$
(18-3)
where $rq=[XqYq0]T$. Substituting Eqs. (12), (15), (16), and (18) into Eq. (17) gives
$KBi=ΔAΔZ∑q=1p[G000G000E00EYq00−EXq−GYqGXq000−GYq00GXqEYq−EXq0EYq2000EXq2000G(Xq2+Yq2)]$
(19)
As p → ∞ (and thus ΔA → 0), $KBi$ becomes
$KBi=AΔZ×diag(G,G,E,ERX2,ERY2,G(RX2+RY2))$
(20)
where
$RX=∫AY2dAAandRY=∫AX2dAA$

$KBi$ is clearly diagonal since $∫AXdA=∫AYdA=0$.

We may need the stiffness matrix expressed at the global frame. The coordinate transformation from BiXYZ to OXYZ is given by
$KO(Z)=E(Z)TKBiE(Z)$
(21)
where
$E(Z)=[1000Z0010−Z00001000000100000010000001]$
(22)
and Z is the distance from OXYZ to BiXYZ. Substituting Eqs. (20) and (22) into Eq. (21), the stiffness matrix can be obtained as
$KO(Z)=1ΔZK(Z)=AΔZ[G000G000E0−GZ0GZ00000symm.ERX2+GZ2000ERY2+GZ2000G(RX2+RY2)]$
(23)
where $K(Z)$ denotes the stiffness matrix. It must be noted here that $K(Z)$ is the function of the transverse displacement Z. It can be said that $K(Z)$ is the stiffness matrix of the virtual springs on the cross section at a distance Z from the origin O.

### 3.2 EOM for Undamped Free Vibration.

For undamped free vibration, the equilibrium equation can be given by
$MiOD¨i+{Ki+1O(Di−Di+1)+KiO(Di−Di−1)}=0$
(24)
where $Mi∈R6×6$ and $Di∈R6×1$ are the mass matrix and displacement vector of the ith piece, respectively. The mass matrix $MiAi$ can be given by $MAi=ρAΔZ×diag(1,1,1,rgX2,rgY2,rgZ2)$, where ρAΔZ is the mass of the rigid body, and rgX, rgY, and rgZ denote the radius of gyration with respect to the axes of AiXYZ. The radius of gyration can be expressed as follows:
$rgX2=RX2+ΔZ212$
(25-1)
$rgY2=RY2+ΔZ212$
(25-2)
$rgZ2=RX2+RY2$
(25-3)
Then, $MiO$ becomes
$MiO=E(Z)TMAiE(Z)$
(26)
or
$MiO=ΔZρA[1000Z0010−Z000010000−Z0RX2+ΔZ212+Z200Z000RY2+ΔZ212+Z2000000RX2+RY2]$
(27)
Substituting Eqs. (23), (25), and (27) into Eq. (24) and dividing by ΔZ gives
$M(Z)D¨(Z,t)−1ΔZ2{(KD)+−(KD)−}=0$
(28)
where
$M(Z)=MiO/ΔZ$
and
$(KD)±=±K(Z±ΔZ/2)(D(Z±ΔZ,t)−D(Z))$
As ΔZ → 0, we can have the following second-order partial differential equation:
$M(Z)D¨(Z,t)−∂∂Z(K(Z)∂∂ZD(Z,t))=0$
(29)
We again assume that the deflection is small harmonic, i.e., $D(Z,t)=D^(Z)ejωt$. Then,
$ddZ(K(Z)ddZD^(Z))=−ω2M(Z)D^(Z)$
or
$(KD^′)′=−ω2MD^$
(30)
where $D^(Z)(=[δXδYδZ;ϕXϕYϕZ]T)$ is the screw deflection (the prime represents the derivative with respect to Z throughout this paper). It is noted here that the elements of $D^$ are functions of Z. It is also noted that
$M=ρA[100010symm.0010−Z0RX2+Z200Z000RY2+Z2000000RX2+RY2]$
(31)
and
$K=A[G000G000E0−GZ0GZ00000symm.ERX2+GZ2000ERY2+GZ2000G(RX2+RY2)]$
(32)

Since we assumed that the beam is prismatic, i.e., dA/dZ = 0, the term A of Eqs. (31) and (32) can be canceled in Eq. (30). However, if dA/dZ ≠ 0, e.g., a tapered beam, then the area A must be considered.

As is well known, Eq. (30) is similar to the Sturm–Liouville problem (SLP). Considering that $M$ and $K$ are positive definite, this may be assumed to be a regular problem (although this is not proved in the present paper). One of the important properties of the regular SLP is that there is an infinite but countable number of ω2, suggesting the presence of infinitely many natural frequencies and corresponding mode shapes.

### 3.3 Boundary Conditions.

Before solving Eq. (30), it is essential to specify the boundary conditions. We specify the three types of joints, namely, fixed, free, and revolute, which are listed in Table 1. First, the fixed type supports prevent all movements in any direction, and the screw deflection is obviously $D^(0)=0$ (fixed foot) or $D^(L)=0$ (fixed end). Second, when the end of the beam is free, there are no reaction forces and moments. Then, we can get $D^′(L)=0$. Although we do not provide an analytical description of the static beams in this paper, the boundary conditions can be derived from the statics of the beam. Similar to Eq. (30), static deflection can be given by
$(KD^′)′=−w^$
(33)
where $w^$ is the applied wrench per length. For the free end of the beam, Eq. (33) becomes $KD^′=0$ at Z = L because the beam does not experience any force and moment at its end. Then, it is clear that $D^′(L)=0$ since $K$ is positive definite. Lastly, we consider a beam supported by revolute joints, as shown in Fig. 7. The rigid element attached to the joint can be freely rotated about the joint axis, but it is clamped in all other directions. Also, there will be no torque reaction about the axis. For the revolute joints that are parallel or aligned with the X-axis, the boundary conditions can be written as
$D^(Z=0)={000;100}T$
(34-1)
$ϕX′(Z=0)=0$
(34-2)
$D^(Z=L)={0L0;100}T$
(34-3)
$ϕX′(Z=L)=0$
(34-4)
Fig. 7
Fig. 7
Close modal
Table 1

Boundary conditions of fixed, free, and revolute joints

Z$D^$$D^′$
Fixedfoot (Z = 0)$D^=0$
end (Z = L)$D^=0$
Freefoot (Z = 0)$D^′=0$
end (Z = L)$D^′=0$
Revolute joint (X-axis)foot (Z = 0)$D^={000;100}T$$ϕX′=0$
end (Z = L)$D^={0L0;100}T$$ϕX′=0$
Revolute joint (Y-axis)foot (Z = 0)$D^(Z=0)={000;010}T$$ϕY′=0$
end (Z = L)$D^(Z=L)={−L00;010}T$$ϕY′=0$
Z$D^$$D^′$
Fixedfoot (Z = 0)$D^=0$
end (Z = L)$D^=0$
Freefoot (Z = 0)$D^′=0$
end (Z = L)$D^′=0$
Revolute joint (X-axis)foot (Z = 0)$D^={000;100}T$$ϕX′=0$
end (Z = L)$D^={0L0;100}T$$ϕX′=0$
Revolute joint (Y-axis)foot (Z = 0)$D^(Z=0)={000;010}T$$ϕY′=0$
end (Z = L)$D^(Z=L)={−L00;010}T$$ϕY′=0$
Equations (34-1) and (34-3) express the line vectors (represented in Plücker's axis coordinates) aligned with the axes of the joints. Equations (34-2) and (34-4) indicate that there is no torque about the joint axes. In the same manner, the boundary conditions for the revolute joints with their axes parallel or aligned with the Y-axis can be given by
$D^(Z=0)={000;010}T$
(35-1)
$ϕY′(Z=0)=0$
(35-2)
$D^(Z=L)={−L00;010}T(35−3)$
(35 -3)
$ϕY′(Z=L)=0(35−4)$
(35 - 4)

Besides the aforementioned supports, we can also obtain the boundary conditions for any support using the reciprocal relation between wrenches and instantaneous velocity [2].

### 3.4 Analytical Solutions.

In this subsection, Eq. (30) is solved for the cantilever. We again assume a uniform and prismatic beam with length L, area of cross section A, Young's modulus E, shear modulus G, and density ρ. The global coordinate frame OXYZ is set as illustrated in Fig. 3.

After substituting Eqs. (31) and (32) into Eq. (30), Eq. (30) can be decoupled into four sets of second-order differential equations as follows:
$EδZ″=−ω2ρδZ$
(36)
$GϕZ″=−ω2ρϕZ$
(37)
$G(δX+ZϕY)″−GϕY′=−ω2ρ(δX+ZϕY)$
(38-1)
$EϕY″+GRY2(δX+ZϕY)′−GRY2ϕY=−ω2ρϕY$
(38-2)
$G(δY−ZϕX)″+GϕX′=−ω2ρ(δY−ZϕX)$
(39-1)
$EϕX″−GRX2(δY−ZϕX)′−GRX2ϕX=−ω2ρϕX$
(39-2)

From these equations, we can see that they are decoupled into axial vibrations, rotational vibrations with respect to the Z-axis, and bending vibration about the principal directions of the cross section. Consequently, this indicates that the natural frequencies (ω′s) are determined from each decoupled equation, and the corresponding mode shapes are axial and rotational vibrations about the Z-axis and bending in the principal directions of the cross section.

Homogeneous solutions of Eqs. (36) and (37) can be easily obtained as follows:
$δZ=sin((2p−1)π2⋅ZL)andωp=1LEρ⋅2p−12π$
(40)
$ϕZ=sin((2p−1)π2⋅ZL)andωp=1LGρ⋅2p−12π$
(41)

It should be noted here that ωp can be determined by the boundary conditions: δX(0) = 0, $δX′(L)=0$, ϕZ(0) = 0, and $ϕZ′(L)=0$.

To solve Eq. (38), we can rewrite this equation as
$(δX+ZϕY)″=ϕY′−ω2ρG(δX+ZϕY)$
(42-1)
and
$ϕY″=−GERY2(δX+ZϕY)′+1E(GRY2−ω2ρ)ϕY$
(42-2)
Then, the second-order differential equation can be reduced to a first-order one by
$u′=Pu$
(43)
where $u={δX+ZϕYϕY(δX+ZϕY)′ϕ′y}T$ and
$P=[00100001−ρω2G00101E(GRY2−ρω2)−GERY20]$
(44)
Suppose that λi and $ui$ for i = 1, …, 4 are the eigenvalues and eigenvectors of $P$, respectively. Then, the solution of $u$ can be expressed as
$u=∑q=14aqeλqZuq$
(45)
where ai are the coefficients, and uji are the jth elements of $ui$. Application of the boundary conditions (δX(0) = 0, $δX′(L)=0$, ϕY(0) = 0, $δY′(L)=0$) to Eq. (45) yields
$∑q=14aqu1q=0$
(46-1)
$∑q=14aqu2q=0$
(46-2)
$∑q=14aqu3qeλqL=(δX+ZϕY)′(L)=ϕY(L)$
(46-3)
$∑q=14aqu4qeλqL=0$
(46-4)
It should be noted here that the derivative of (δX + Y) at Z = L is
$(δX+ZϕY)′(L)=δX′(L)+LϕY′(L)+ϕY(L)=ϕY(L)$
Since $ϕY(L)=∑q=14aqu2qeλqL$, we can have
$U{a1a2a3a4}={0000}$
(47)
where
$U={u11u12u13u14u12u22u23u24(u31−u21)eλ1L(u32−u22)eλ2L(u33−u23)eλ3L(u34−u24)eλ4Lu41eλ1Lu42eλ2Lu43eλ3Lu44eλ4L}$
(48)

For a nontrivial solution, the following characteristic equation must be zero: $detU=0$. Recalling that λi and $ui$ are the eigenvalues and the corresponding eigenvectors of $P$, respectively, it can be said that $detU$ depends on ω. The determinant of $U$ is not always zero, but there may exist ω so that $detU=0$, suggesting that ω is the natural frequency for $detU=0$. As previously mentioned, when considering the regular SLP, there are infinitely many values of ω that satisfy $detU=0$. If the above is iterated, then we can find the natural frequencies ωp for p = 1, 2, …, n, …, ∞. Once ωp is found, the mode shapes depending on ai, λi, and $vi$ can be determined.

Next, Eq. (39) can be solved in a similar manner to Eq. (38). Equation (39) can be reduced to
$v′=Qv$
(49)
where $v={δY−ZϕXϕX(δY−ZϕX)′ϕ′X}T$ and
$Q=[00100001−ρω2G00−101E(GRX2−ρω2)GERX20]$
(50)
Supposing that μi and vi for i = 1, …, 4 are the eigenvalues and eigenvectors of $Q$, respectively, then $v$ can be written as
$v=∑q=14bqeμqZvq$
(51)
where bi are coefficients. The boundary conditions (δY(0) = 0, $δY′(L)=0$, ϕX(0) = 0, $δX′(L)=0$) provide the characteristic equation: $detV=0$, where
$V={v11v12v13v14v21v22v23v24(v31+v21)eμ1L(v32+v22)eμ2L(v33+v23)eμ3L(v34+v24)eμ4Lv41eμ1Lv42eμ2Lv43eμ3Lv44eμ4L}$
(52)
and vji is the jth element of $vi$. Again, we need to specify ω and compute $detV$. For $detV=0$, the natural frequencies and the corresponding mode shapes can be determined.

There may be a degenerate case that both $detU$ and $detV$ become zero for the same ω. Since this is a linear vibration system, the solutions can be expressed by linear combinations of vibrational modes with the same natural frequency. In other words, the solutions can be a mixture of some modes with the same ω. The degeneracy sometimes occurs in symmetrical beams such as a cylindrical beam.

Until now, we have described how to solve the differential equations. As mentioned above, $U$ and $V$ are dependent on the boundary conditions, as opposed to $P$ and $Q$ which are not. This means that the natural frequencies and the corresponding mode shapes are affected by the boundary conditions. Therefore, it is important to construct $U$ and $V$ using the boundary conditions and find ω through the iterative method presented. However, it should be mentioned that the natural frequencies for bending vibrations cannot be explicitly expressed.

## 4 Discrete Analysis of Beam

This section introduces a method to discretize a continuous system with the aforementioned boundary conditions into an equivalent multibody system. We again model a beam as a stack of n rigid bodies, and then the EOM is given by Eq. (5). As already mentioned, in order for a multibody system to have the same dynamic characters of the beam, the mass properties and elasticity are reflected into the mass and stiffness matrices $M$ and $K$. The mass matrix of the ith rigid body $Mi$ can be easily obtained by the coordinate transformation in Eq. (26):
$Mi=E((i−12)ΔZ)TMAiE((i−12)ΔZ)$
(53)
where ΔZ = L/n. However, we should consider not only elasticity but also boundary conditions when determining the stiffness matrix. We assume that the boundary conditions are given by constraints on both ends of the beam. Then, they affect obviously only $K1$ and $Kn$. So, we can use Eqs. (20) and (21) to determine the stiffness matrices except for $K1$ and $Kn$ as follows:
$Ki=E((i−1)ΔZ)TKBi−1E((i−1)ΔZ)+E(iΔZ)TKBiE(iΔZ)fori=2,…,n−1$
(54)
and
$Ki−1,i=E((i−1)ΔZ)TKBi−1E((i−1)ΔZ)$
(55)

For $K1$ and $Kn$, it is necessary to modify the EOM according to the type of supports. The following subsections explain how to apply boundary conditions for each case.

### 4.1 Fixed and Free Ends.

For a beam fixed at the foot, the first rigid body does not move, i.e., the displacement vector is $D1=0$. As is well known, in such cases, the rows and columns including both $M1$ and $K1$ are canceled as follows:
$M=[M20⋱0Mn]$
and
$K=[K2−K2,30⋯0−K2,3K3−K3,4⋮0−K3,4⋱0⋮Kn−1−Kn−1,n0⋯0−Kn−1,nKn]$

For a beam fixed at end, the rows and columns including $Mn$ and $Kn$ are vanished. When a beam is fixed at both ends, we should erase the rows and columns containing all of $M1$, $K1$, $Mn$ and $Kn$.

When the foot of a beam is free, it is clear that $K1=K1,2$ since the first rigid body is connected with the second one. For free end, $Kn=Kn−1,n$.

Now that we have determined both $M$ and $K$ for fixed and free ends, eigenanalysis can be performed to obtain natural frequencies and vibration modes after assuming small harmonic displacements.

### 4.2 General Boundary Conditions.

In this subsection, we aim to discuss supports with p-DOF. While the method in the previous subsection can cover 0- and 6-DOF supports, they are insufficient for describing supports with n-DOF such as revolute, cylindrical, and spherical joints. Therefore, we introduce a more general approach that can cover arbitrary supports in this subsection.

Suppose that the supports have p- and q-DOF, respectively. Then, the displacement vectors of the first and last rigid bodies can be expressed as
$D1=∑i=1paiA^i=[A^1⋯A^p]a=Aa$
(56-1)
and
$Dn=∑i=1qbiB^i=[B^1⋯B^q]b=Bb$
(56-2)
where $a(=[a1⋯ap]T∈Rp×1)$ and $b(=[b1⋯bq]T∈Rq×1)$ are time-dependent vectors. The columns of $A(∈R6×p)$ and $B(∈R6×q)$ are linearly independent displacement vectors of the first and last rigid bodies, respectively. Substituting Eq. (56) into Eq. (5) yields the following equations:
$M1Aa¨+K1Aa−K1,2D2=0$
(57-1)
$M2(n−1)D¨2(n−1)+K2(n−1)D2(n−1)−{K1,2Aa0(6n−24)×1Kn−1,nBb}=0$
(57-2)
$MnBb¨+KnBb−Kn−1,nDn−1=0$
(57-3)
where
$D2(n−1)={D2T⋯Dn−1T}T∈R(6n−12)×1$
$M2(n−1)=[M20⋱0Mn−1]∈R(6n−12)×(6n−12)$
and
$K2(n−1)=[K2−K2,30⋯0−K2,3K3−K3,4⋮0−K3,4⋱0⋮Kn−2−Kn−2,n−10⋯0−Kn−2,n−1Kn−1]∈R(6n−12)×(6n−12)$
Meanwhile, $K1$ represents the stiffness matrix of all the springs attached to the first rigid body. So, it can be written as $K1=K0,n+K1,2$, where $K0,n$ is the stiffness concerning the joint at the foot. The joint cannot resist to displacement belonging to the column space of $A$, i.e., $K0,nA=0$, and therefore
$K1A=K1,2A$
(58)
Similarly, we can have
$KnB=Kn−1,nB$
(59)
Substituting Eqs. (58) and (59) into Eq. (57) and premultiplying the left side of Eqs. (57-1) and (57-3) by $AT$ and $BT$, respectively, yield
$[m~1]ϕ¨1+[k~1]ϕ1−ATK1,2D^2=0$
(60-1)
$[m~n]ϕ¨1+[k~n]ϕn−BTKn−1,nD^n−1=0$
(60-2)
where $[m~1]=ATM1A$, $[m~n]=BTMnB$, $[k~1]=ATK1,2A$, and $[k~n]=BT−n−1,nB$. Using Eqs. (57-2) and (60), we can reconstruct the EOM including the displacement vector, mass, and stiffness matrices as follows:
$MD¨+KD=0$
where
$D={aTD2T⋯Dn−1TbT}T$
(61-1)
$M=[[m~1]0M2(n−1)0[m~n]]$
(61-2)
and
$K=[[k~1]−ATK1,20⋯0−K1,2AK2−K2,3⋮0−K2,3⋱0⋮Kn−1−Kn−1,nB0⋯0−BTKn−1,n[k~n]]$
(61-3)

We have derived the EOM governing vibrations of a multibody system equivalent to a beam for each type of support. The EOM is determined by the DOF of the rigid bodies at both ends, and the natural frequencies and vibration modes can be obtained through eigenanalysis. Furthermore, this method can be applied to a wider range of cases, considering its ability to handle boundary conditions of other bodies. This suggests that it can be also effectively utilized for the discrete analysis of plates.

## 5 Numerical Examples

This section introduces two numerical examples, namely the cantilever and the robot link. First, for the cantilever, we explain the approach employed in this paper and compare the accuracy of this method with numerical solutions (or discrete system analysis). Second, we provide an analytic solution for the robot link (similar to the example described in Ref. [2]) where the revolute joint axes are not parallel and show the discretization into an equivalent multibody system.

### 5.1 Cantilever Beam.

In this example, we consider a uniform and prismatic cantilever with a rectangular cross section (Fig. 8). The parameters of the beam are given in Table 2. Using these parameters, Eqs. (40) and (41) can provide the natural frequencies for axial and rotational vibrations about the Z-axis, as listed in Table 3.

Fig. 8
Fig. 8
Close modal
Table 2

System parameters of the cantilever

 Material $ρ=8050(kg/m3)$⁠, $E=200×109(N/m2)$⁠, $G=7.6923×1010(N/m2)$ Geometry h = 0.2 (m) (height), w = 0.1 (m) (width), L = 1.5 (m), RX = 0.0289 (m2), RY = 0.0577 (m2)
 Material $ρ=8050(kg/m3)$⁠, $E=200×109(N/m2)$⁠, $G=7.6923×1010(N/m2)$ Geometry h = 0.2 (m) (height), w = 0.1 (m) (width), L = 1.5 (m), RX = 0.0289 (m2), RY = 0.0577 (m2)
Table 3

Natural frequencies ωp (rad/s) for axial (δZ) and rotational (ϕZ) vibrations

Typesω1ω2ω3ω4ω5
δZ5219.7015,659.1026,098.5136,537.9146,977.31
ϕZ3237.129711.3616,185.6122,659.8529,134.09
Typesω1ω2ω3ω4ω5
δZ5219.7015,659.1026,098.5136,537.9146,977.31
ϕZ3237.129711.3616,185.6122,659.8529,134.09
In order to solve for bending vibrations, we substitute the values of these parameters into Eqs. (44) and (50)
$P=[00100001−1.0465×10−7ω20010115.38−4.025×10−8ω2−115.380]$
(62)
and
$Q=[00100001−1.0465×10−7ω200−10461.54−4.025×10−8ω2461.540]$
(63)

We then compute the eigenvalues and the corresponding eigenvectors of $P$ and $Q$, respectively, iteratively for 0 < ω < 10,000[rad/s] with a step size of 0.01. For $detU=0$ or $detV=0$, the solutions are listed in Table 4. It is noted that $U$ and $V$ are described by Eqs. (48) and (52), respectively. Now, we can line up the natural frequencies and the corresponding mode shapes for all vibrations, as illustrated in Fig. 9. As shown, our results are similar to those of the classical beam theory.

Fig. 9
Fig. 9
Close modal
Table 4

Natural frequencies ωp (rad/s) for bending vibrations

ω1ω2ω3
$δX$, $ϕY(P)$444.252603.726676.63
$δY$, $ϕX(Q)$224.161379.683758.78
ω1ω2ω3
$δX$, $ϕY(P)$444.252603.726676.63
$δY$, $ϕX(Q)$224.161379.683758.78
Consequently, we performed a discrete analysis to evaluate the application of the screw approach to a continuous vibration system. From the boundary conditions of the cantilever, the mass and stiffness matrices can be given by
$M=[M20⋱0Mn]$
and
$K=[K2−K2,30⋯0−K2,3K3−K3,4⋮0−K3,4⋱0⋮Kn−1−Kn−1,n0⋯0−Kn−1,nKn−1,n]$

$Mi$, $Ki$, and $Ki−1,i$ were calculated by using Eqs. (53)(55). To perform eigenanalysis, we again assume small vibrations, i.e., $D=D~ejωt$, where $D~={D^2TD^3T⋯D^nT}T$. $D^i$ is the vibrational axis of the ith rigid body, and the eigenequation is $−ω2MD~+KD~=0$, which was solved in matlab for n = 5, 10, 20, and 30.

Our results for the natural frequencies are listed in Table 5. It can be seen that the numerical solutions approach exact solutions as n increases. For bending vibration (except ω5), the error rates between the exact and numerical results decrease from approximately 27% (n = 5) to 3% (n = 30). In addition, the error rate for the axial vibration (ω5) decreases from 11% to 2%. The error rates are calculated by $(E−N)/E×100(%)$, where E and N denote the exact and numerical solutions, respectively. These results suggest that the error rates depend on the type of vibration and the number of rigid bodies n. In contrast, the error rates are independent of high and low frequencies. It is generally known that a larger error occurs at a higher frequency. In this sense, a beam modeled as a stack of rigid bodies and the concept of virtual springs may be utilized for the analysis of beams with more complex geometry.

Table 5

Exact and numerical solutions of the natural frequencies

nω1ω2ω3ω4ω5ω6
Exact solutions224.16444.251379.682603.723237.123758.78
Numerical solutions (error rate)5279.68 (24.8%)552.98 (24.4%)1759.19 (27.5%)3262.97 (25.3%)3578.57 (10.5%)4802.90 (27.8%)
10248.92 (11.0%)492.83 (10.9%)1538.20 (11.5%)2883.30 (10.7%)3403.62 (5.1%)4195.63 (11.6%)
20235.90 (5.2%)467.31 (5.2%)1452.68 (5.3%)2733.36 (5.0%)3319.22 (2.5%)3956.33 (5.3%)
30231.86 (3.4%)459.37 (3.4%)1427.09 (3.4%)2688.08 (3.2%)3291.60 (1.7%)3886.24 (3.4%)
nω1ω2ω3ω4ω5ω6
Exact solutions224.16444.251379.682603.723237.123758.78
Numerical solutions (error rate)5279.68 (24.8%)552.98 (24.4%)1759.19 (27.5%)3262.97 (25.3%)3578.57 (10.5%)4802.90 (27.8%)
10248.92 (11.0%)492.83 (10.9%)1538.20 (11.5%)2883.30 (10.7%)3403.62 (5.1%)4195.63 (11.6%)
20235.90 (5.2%)467.31 (5.2%)1452.68 (5.3%)2733.36 (5.0%)3319.22 (2.5%)3956.33 (5.3%)
30231.86 (3.4%)459.37 (3.4%)1427.09 (3.4%)2688.08 (3.2%)3291.60 (1.7%)3886.24 (3.4%)

As listed in Table 6, the vibration axes of the first, third, and sixth modes lie on the ZX plane, whereas those of the second and fourth modes lie on the YZ plane, and the vibration axes of the fifth modes are aligned with the Z-axis. All mode shapes are shown in Fig. 10, and our results are almost the same as the ones obtained by the screw approach. Referring to Fig. 11, the mode shapes can be thought of as rotations of rigid bodies about the respective vibration axes, and this thought is the most important concept in the screw theory of a linear vibration, such as the vibration center and axis [4]. Therefore, our results suggest that this approach can be extended to describe the vibration of two connected rigid bodies [13] but also to a n-body vibration system.

Fig. 10
Fig. 10
Close modal
Fig. 11
Fig. 11
Close modal
Table 6

Numerical solutions obtained by the discrete analysis for n = 30

Form of vectorResults of vibration axes
1st mode$D^i=ϕX{0Z0;100}T$(parallel to the X-axis)$D^5=0.07{00.120;100}T$
$D^10=0.14{00.230;100}T$
$D^20=0.21{00.380;100}T$
$D^30=0.21{00.430;100}T$
2nd mode$D^i=ϕY{−Z00;010}T$(parallel to the Y-axis)$D^5=0.07{−0.1100;010}T$
$D^10=0.14{−0.2200;010}T$
$D^20=0.21{−0.3700;010}T$
$D^30=0.22{−0.4100;010}T$
3rd mode$D^i=ϕX{0Z0;100}T$(parallel to the X-axis)$D^5=−0.09{00.100;100}T$
$D^10=−0.08{00.060;100}T$
$D^20=0.12{01.250;100}T$
$D^30=0.21{01.190;100}T$
4th mode$D^i=ϕY{−Z00;010}T$(parallel to the Y-axis)$D^5=−0.08{−0.0700;010}T$
$D^10=−0.08{0.0200;010}T$
$D^20=0.12{−1.2500;010}T$
$D^30=0.21{−1.1800;010}T$
5th mode$D^i=ϕZ{000;001}T$(aligned with the Z-axis)$D^5=0.06{000;001}T$
$D^10=0.12{000;001}T$
$D^20=0.22{000;001}T$
$D^30=0.26{000;001}T$
6th mode$D^i=ϕX{0Z0;100}T$(parallel to the X-axis)$D^5=0.11{00.070;100}T$
$D^10=−0.02{02.090;100}T$
$D^20=−0.06{00.520;100}T$
$D^30=0.22{01.310;100}T$
Form of vectorResults of vibration axes
1st mode$D^i=ϕX{0Z0;100}T$(parallel to the X-axis)$D^5=0.07{00.120;100}T$
$D^10=0.14{00.230;100}T$
$D^20=0.21{00.380;100}T$
$D^30=0.21{00.430;100}T$
2nd mode$D^i=ϕY{−Z00;010}T$(parallel to the Y-axis)$D^5=0.07{−0.1100;010}T$
$D^10=0.14{−0.2200;010}T$
$D^20=0.21{−0.3700;010}T$
$D^30=0.22{−0.4100;010}T$
3rd mode$D^i=ϕX{0Z0;100}T$(parallel to the X-axis)$D^5=−0.09{00.100;100}T$
$D^10=−0.08{00.060;100}T$
$D^20=0.12{01.250;100}T$
$D^30=0.21{01.190;100}T$
4th mode$D^i=ϕY{−Z00;010}T$(parallel to the Y-axis)$D^5=−0.08{−0.0700;010}T$
$D^10=−0.08{0.0200;010}T$
$D^20=0.12{−1.2500;010}T$
$D^30=0.21{−1.1800;010}T$
5th mode$D^i=ϕZ{000;001}T$(aligned with the Z-axis)$D^5=0.06{000;001}T$
$D^10=0.12{000;001}T$
$D^20=0.22{000;001}T$
$D^30=0.26{000;001}T$
6th mode$D^i=ϕX{0Z0;100}T$(parallel to the X-axis)$D^5=0.11{00.070;100}T$
$D^10=−0.02{02.090;100}T$
$D^20=−0.06{00.520;100}T$
$D^30=0.22{01.310;100}T$

In this subsection, a robot link is used as an example similar to the one proposed in Ref. [2] to demonstrate the utility of the screw theory. In general, the classical beam theory does not deal with hinges that are not parallel at both ends, such as the robot link. However, the boundary conditions only affect $U$ and $V$, as shown in the methodology presented in this paper, and we can therefore easily obtain analytical solutions for the robot link.

The robot link consists of a beam with a square cross section. The hinge at the foot is aligned with the X-axis, and the other hinge is parallel to the Y-axis. The system parameters are given in Table 7. Referring to Table 1, the boundary conditions are given by
$D^={000;100}TandϕX′=0atZ=0$
(64)
$D^={−L00;010}TandϕY′=0atZ=L$
(65)
Table 7

System parameters of the robot link

 Material $ρ=8050(kg/m3)$⁠, $E=200×109(N/m2)$⁠, $G=7.6923×1010(N/m2)$ Geometry $h=w=0.1(m)$⁠, $L=1.5(m)$⁠, $RX=0.0289(m2)$⁠, $RY=0.0289(m2)$
 Material $ρ=8050(kg/m3)$⁠, $E=200×109(N/m2)$⁠, $G=7.6923×1010(N/m2)$ Geometry $h=w=0.1(m)$⁠, $L=1.5(m)$⁠, $RX=0.0289(m2)$⁠, $RY=0.0289(m2)$
The joints at both ends do not allow the translation along Z-axis and the rotation about Z-axis. It implies that the solutions for axial and rotational vibrations are equivalent to the case fixed at both ends. In the fixed case, Eqs. (36) and (37) can be used to get the natural frequencies of axial and rotational vibrations ωa and ωr as follows:
$ωa=10439.40,20878.81,31318.21[rad/s]$
$ωr=6474.24,12948.49,19422.73[rad/s]$
To solve for bending vibration, we compute $P$ and $Q$
$P=[00100001−1.0465×10−7ω20010461.54−4.025×10−8ω2−461.540]$
(66)
and
$Q=[00100001−1.0465×10−7ω200−10461.54−4.025×10−8ω2461.540]$
(67)

To construct $U$ and $V$, the boundary conditions are arranged as

$δX=ϕY=0(atZ=0)$
(68-1)
$δX+LϕY=ϕY′=0(atZ=L)$
(68-2)
$δY=ϕX′=0(atZ=0)$
(69-1)
$δY=ϕX=0(atZ=L)$
(69-2)
By introducing the eigenvalues λi and eigenvectors $ui$ of $P$, the boundary conditions of Eq. (68) can be rewritten as
$[u11u12u13u14u21u22u23u24u11eλ1Lu12eλ2Lu13eλ3Lu14eλ4Lu41eλ1Lu42eλ2Lu43eλ3Lu44eλ4L]{a1a2a3a4}={0000}$
(70)
and thus, $U$ is
$U=[u11u12u13u14u21u22u23u24u11eλ1Lu12eλ2Lu13eλ3Lu14eλ4Lu41eλ1Lu42eλ2Lu43eλ3Lu44eλ4L]$
(71)
Similarly, Eq. (69) gives $V$
$V=[v11v12v13v14v41v42v43v44v11eμ1Lv12eμ2Lv13eμ3Lv14eμ4Lv21eμ1Lv22eμ2Lv23eμ3Lv24eμ4L]$
(72)

By iteration, we can find the six degenerate solutions for 0 < ω < 20 000[rad/s] as listed in Table 8. That is, $detU=detV=0$ for the solutions. Considering the symmetry of the robot link, we can say that the degeneracy is natural. Therefore, two vibration modes occur with the same natural frequency. It means that mode shape can be expressed by linear combinations of the two vibration modes, and therefore, there are infinitely many mode shapes for a natural frequency. Figure 12 illustrates mode shapes where |ϕX(Z = 0)| = |ϕY(Z = L)| among numerous vibration modes for the first and second natural frequencies, respectively. The thick lines depict the axes of revolute joints. This figure indicates that the robot link does not have pure bending modes. In terms of finding the analytic solutions for the robot link, the proposed screw approach can be useful for the vibration analysis of a beam that experiences an arbitrary bending.

Fig. 12
Fig. 12
Close modal
Table 8

Analytical solutions of natural frequency for the robot link

ω1ω2ω3ω4ω5ω6
971.343069.506186.0810,146.8614,781.8419,943.06
ω1ω2ω3ω4ω5ω6
971.343069.506186.0810,146.8614,781.8419,943.06
In what follows, we performed a multibody analysis for the robot link. The boundary conditions of Eqs. (64) and (65) yield
$D1=a{000;100}T$
(73-1)
and
$Dn=b{−L00;010}T$
(73-2)
Obviously, the revolute joints do not produce forces and moments with respect to $D^1$ and $D^n$, respectively, i.e., $K0,1D^1=0$ and $K0,nD^n=0$. It can be used to follow the reconstruction of the mass and stiffness matrices in Eqs. (57)(61) as follows:
$M=[m~10M2(n−1)0m~n]$
and
$K=[k~1−D^1TK1,20⋯0−K1,2D^1K2−K2,3⋮0−K2,3⋱0⋮Kn−1−Kn−1,nD^n0⋯0−D^nTKn−1,nk~n]$
where $m~i=D^iTMiD^i$, $k~1=D^1TK1,2D^1$, and $k~n=D^nTKn−1,nD^n$. The mass and stiffness matrices can be calculated by using the system parameters (Table 7) and the number of bodies n. The displacement vector is defined as $D={aD2T⋯Dn−1Tb}T$. For small harmonic displacements, the vector can be written as $D=D~ejωt$ where $D~={a^D^2T⋯D^n−1Tb^}T$. Once again, the eigenanalysis was conducted for n = 5, 10, 20, and 30.

Table 9 shows the results of the natural frequencies ωi with degeneracy as well as the nondegenerate ones. The vibration modes for n = 30 are listed in Table 10. Observation of Tables 9 and 10 reveals that the degenerate vibration modes are two linearly independent vectors and that the nondegenerate ones are axial or rotational vibrations when the robot link is fixed at both ends, i.e., a = b = 0. In addition, the error rates decrease as n increases, as shown in Table 10. Consequently, we can say that the numerical solutions converge to exact solutions.

Table 9

Numerical solutions of the natural frequencies

nω1,2ω3,4ω5,6ω7ω8ω9
51235.91 (27.1%)3964.45 (29.1%)7828.93 (26.6%)7886.4012,716.4414,572.17
n$ω1,2$$ω3,4$$ω5,6$$ω7$$ω8,9$$ω10$
101082.99 (11.5%)3435.11 (11.9%)6932.27 (12.1%)7157.14 (10.5%)11,347.25 (11.8%)11,540.54 (10.5%)
201023.08 (5.3%)3235.22 (5.4%)6520.89 (5.4%)6807.23 (5.14%)10,689.36 (5.3%)10,976.33 (5.14%)
301005.01 (3.5%)3176.54 (3.5%)6401.46 (3.5%)6694.21 (3.4%)10,496.16 (3.4%)10,794.10 (3.4%)
nω1,2ω3,4ω5,6ω7ω8ω9
51235.91 (27.1%)3964.45 (29.1%)7828.93 (26.6%)7886.4012,716.4414,572.17
n$ω1,2$$ω3,4$$ω5,6$$ω7$$ω8,9$$ω10$
101082.99 (11.5%)3435.11 (11.9%)6932.27 (12.1%)7157.14 (10.5%)11,347.25 (11.8%)11,540.54 (10.5%)
201023.08 (5.3%)3235.22 (5.4%)6520.89 (5.4%)6807.23 (5.14%)10,689.36 (5.3%)10,976.33 (5.14%)
301005.01 (3.5%)3176.54 (3.5%)6401.46 (3.5%)6694.21 (3.4%)10,496.16 (3.4%)10,794.10 (3.4%)
Table 10

Numerical solutions obtained by the discrete analysis for n = 30

$a^$$D^16$$b^$
ω1,2$D~1$$0.23$${0−0.130−0.05400}T$$0$
$D~2$$0$${0.0240000.070}T$$0.20$
ω3,4$D~3$$0.22$${0−0.140−0.2100}T$$0$
$D~4$$0$${0.16000−0.180}T$$0.21$
ω5,6$D~5$$0.23$${0−0.130−0.05400}T$$0$
$D~6$$0$${−0.060000.120}T$$0.21$
ω7$D~7$$0$${00000−0.26}T$$0$
ω8,9$D~8$$0.23$${0−0.130−0.05400}T$$0$
$D~9$$0$${−0.130000.160}T$$0.21$
ω10$D~10$$0$${00−0.26000}T$$0$
$a^$$D^16$$b^$
ω1,2$D~1$$0.23$${0−0.130−0.05400}T$$0$
$D~2$$0$${0.0240000.070}T$$0.20$
ω3,4$D~3$$0.22$${0−0.140−0.2100}T$$0$
$D~4$$0$${0.16000−0.180}T$$0.21$
ω5,6$D~5$$0.23$${0−0.130−0.05400}T$$0$
$D~6$$0$${−0.060000.120}T$$0.21$
ω7$D~7$$0$${00000−0.26}T$$0$
ω8,9$D~8$$0.23$${0−0.130−0.05400}T$$0$
$D~9$$0$${−0.130000.160}T$$0.21$
ω10$D~10$$0$${00−0.26000}T$$0$

Figure 13 illustrates the mode shapes of robot links. The vibration modes with degeneracy can be expressed by linear combinations of two vectors: $D~=αD~i+βD~i+1$, and thus, there are infinitely many mode shapes with a natural frequency. The mode shapes illustrated in Fig. 13 were selected among the mode shapes so that |ϕX(Z = 0)| = |ϕY(Z = L)|. In this case, the vibration axis of the ith piece $D^i$ is not a line but a screw. That is, the mode shape is a repetitive screw motion. It suggests that the mode shape can be considered as screw motions as well as rotations. This is a key insight provided by the screw theory.

Fig. 13
Fig. 13
Close modal

So far, we have interpreted the dynamic characteristics of the robot link both analytically and numerically. Through both approaches, we were able to obtain nondegenerate and degenerate solutions and also confirmed that the numerical solutions converge to the analytical solutions. In the future work, we will apply the approaches to problems with comply geometry and plate theory.

## 6 Conclusions

The present study demonstrates that the screw theory could accurately describe the vibrations of discrete and continuous systems. We initially modeled a beam as a stack of rigid bodies connected by virtual springs to each other, and then derived the EOM governing different types of vibrations, such as extension, compression, torsion, and bending, and introduced a method to solve these equations analytically for each boundary condition.

Our study makes three major contributions. First, the stiffness matrix of the beam was derived by using virtual springs acting as normal and shear stresses. This representation could be very useful for deformation problems of continuous systems, such as in ones evaluating the dynamics and statics of the beam. In addition, this approach also enabled us to easily obtain the stiffness matrix of a multibody system by coordinate transformation. Second, it was shown that a continuous system could be approximated as a stack of rigid bodies connected by virtual springs to each other, and the effectiveness of this approach was confirmed by a comparison between analytical and numerical results. Our findings also revealed that the mode shapes could be thought of as rotations (or generally screw motions) of rigid bodies about the vibration axes. Lastly, we analytically solved the differential equations for nonstandard forms of the classical beam theory. We showed that analytical solutions could be easily obtained by the proposed iterative method, provided that the boundary conditions were clearly specified.

Throughout this paper, we only considered deformations in one direction (Z-axis). As a result, the stiffness matrix of the beam was derived, and the mode shapes could be understood as rotations of rigid bodies. However, these concepts could be extended into deformations in two directions, such as the plate theory. It is expected that the stiffness matrix can still be used, and the mode shapes of the plate could be considered as screw motions (or rotations) of small rigid plates. In the future, we intend to investigate the applicability of this approach by considering the screw theory.

## Acknowledgment

This research was supported by a grant from R&D Program (Development of core technology for evacuation control in the accident case and passenger safety, PK2302A2) of the Korea Railroad Research Institute and the Korea Agency for Infrastructure Technology Advancement (KAIA) grant funded by the Ministry of Land, Infrastructure and Transport (Grant RS-2023-00238018, Development of technology for cognizing, predicting and responding to high-risk disasters for deep railway).

## Funding Data

• Korea Railroad Research Institute (KRRI) and Korea Agency for Infrastructure Technology Advancement (KAIA).

## Conflict of Interest

There are no conflicts of interest. This article does not include research in which human participants were involved. Informed consent was obtained for all individuals. Documentation provided upon request. This article does not include any research in which animal participants were involved.

## Data Availability Statement

The authors attest that all data for this study are included in the paper.

## Nomenclature

$s$ =

unit direction vector of $s^i$

A =

area of cross section

E =

Young's modulus

G =

shear modulus

L =

length of a beam

$D$ =

time-dependent displacement vector

$K$ =

stiffness matrix

$M$ =

mass matrix

$s^i$ =

line vector expressed in Plücker's axis coordinates

$D~$ =

screw chain

$D^$ =

time-independent displacement vector (or vibration axes)

$ri$ =

position vector of $s^i$

$Ki$ =

stiffness matrix of all springs attached to ith rigid body

$Ki,j$ =

stiffness matrix of springs connecting ith and jth bodies

$Mi$ =

mass matrix of ith rigid body

OXYZ =

global frame

$D^(Z)$ =

screw deflection

AiXYZ =

local frame placed in ith rigid body

BiXYZ =

local frame located in ith cut plane

$δ$ =

small translational vector

ρ =

density

ϕ =

small angular vector

ωp =

pth natural frequency

## References

1.
Selig
,
J.
, and
Ding
,
X.
,
2001
, “
A Screw Theory of Static Beams
,”
Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems
,
Maui, HI
,
Oct. 29–Nov. 3
, pp.
312
317
.
2.
Selig
,
J.
, and
Ding
,
X.
,
2009
, “
A Screw Theory of Timoshenko Beams
,”
ASME J. Appl. Mech
,
76
(
3
), p.
031003
.
3.
Ball
,
R. S.
,
1900
,
A Treatise on the Theory of Screws
,
Cambridge University Press
,
Cambridge, MA
.
4.
Blanchet
,
P.
,
1998
, “
Linear Vibration Analysis Using Screw Theory
,”
Ph.D. thesis
,
Georgia Institute Technology
,
Atlanta, GA
.
5.
Kim
,
J. W.
,
Lee
,
S.
, and
Choi
,
Y. J.
,
2015
, “
Design Method of Planar Vibration System for Specified Ratio of Energy Peaks
,”
J. Sound Vib.
,
344
, pp.
363
376
.
6.
Kim
,
H. S.
,
Kim
,
J. W.
,
Park
,
S.-B.
, and
Choi
,
Y. J.
,
2017
, “
Design of Serial Linkage-Type Vibration Energy Harvester With Three Resonant Frequencies
,”
Smart Mater. Struct.
,
26
(
11
), p.
115030
.
7.
Kim
,
H. S.
,
Ryu
,
W.
,
Park
,
S.-B.
, and
Choi
,
Y. J.
,
2019
, “
3-Degree-of-Freedom Electromagnetic Vibration Energy Harvester With Serially Connected Leaf Hinge Joints
,”
J. Intell. Mater. Syst. Struct.
,
30
(
2
), pp.
308
322
.
8.
Hong
,
M. B.
, and
Choi
,
Y. J.
,
2011
, “
The Conditions for Mode Decoupling of a Linear Vibration System With Diagonalizable Stiffness Matrices via Screw Theory
,”
Proc. Inst. Mech. Eng. C, J. Mech. Eng. Sci.
,
225
(
1
), pp.
101
111
.
9.
Jang
,
S. J.
, and
Choi
,
Y. J.
,
2007
, “
Geometrical Design Method of Multi-Degree-of-Freedom Dynamic Vibration Absorbers
,”
J. Sound Vib.
,
303
(
1–2
), pp.
343
356
.
10.
Park
,
S.-B.
,
Jang
,
S.-J.
,
Kim
,
I.-H.
, and
Choi
,
Y. J.
,
2017
, “
Broadband Vibration Energy Harvester Utilizing Three Out-of-Plane Modes of One Vibrating Body
,”
Smart Mater. Struct.
,
26
(
10
), p.
105049
.
11.
Lee
,
Y. G.
,
Ryu
,
W.
,
Choi
,
H. S.
, and
Choi
,
Y. J.
,
2020
, “
The Conditions for a Linear Vibration System to Have Only Pure Rotation Modes
,”
IEEE Access
,
8
, pp.
75860
75873
.
12.
Lee
,
Y. G.
, and
Choi
,
Y. J.
,
2020
, “
Design of a Spatial Vibration System for Prescribed Ratio of Energy Peaks
,”
IEEE Access
,
8
, pp.
104371
104385
.
13.
Ryu
,
W.
,
Lee
,
Y. G.
,
Jeon
,
C. G.
, and
Choi
,
Y. J.
,
2020
, “
Design for Prespecified Ratios Between Six Energy Peaks Using a Planar Symmetric Dual-Body Vibration System
,”
IEEE Access
,
8
, pp.
117817
117832
.
14.
Lele
,
S. P.
, and
Maiti
,
S. K.
,
2002
, “
Modeling of Transverse Vibration of Short Beams for Crack Detection and Measurement of Crack Extension
,”
J. Sound Vib.
,
257
(
3
), pp.
559
583
.
15.
Lin
,
H.-P.
,
2004
, “
Direct and Inverse Methods on Free Vibration Analysis of Simply Supported Beams With a Crack
,”
Eng. Struct.
,
26
(
4
), pp.
427
436
.
16.
Whalen
,
T. M.
,
2008
, “
The Behavior of Higher Order Mode Shape Derivatives in Damaged, Beam-Like Structures
,”
J. Sound Vib.
,
309
(
3–5
), pp.
426
464
.
17.
Shafiei
,
M.
, and
Khaji
,
N.
,
2011
, “
Analytical Solutions for Free and Forced Vibrations of a Multiple Cracked Timoshenko Beam Subject to a Concentrated Moving Load
,”
Acta Mech.
,
221
(
1–2
), pp.
79
97
.
18.
Liu
,
S.
,
Zhang
,
L.
,
Chen
,
Z.
,
Zhou
,
J.
, and
Zhu
,
C.
,
2016
, “
Mode-Specific Damage Identification Method for Reinforced Concrete Beams: Concept, Theory and Experiments
,”
Constr. Build. Mater.
,
124
, pp.
1090
1099
.
19.
Čas
,
B.
,
Planinc
,
I.
, and
Schnabl
,
S.
,
2018
, “
Analytical Solution of Three-Dimensional Two-Layer Composite Beam With Interlayer Slips
,”
Engi. Struct.
,
173
, pp.
269
282
.
20.
Liu
,
S.
, and
Yang
,
B.
,
2019
, “
A Closed-Form Analytical Solution Method for Vibration Analysis of Elastically Connected Double-Beam Systems
,”
Compos. Struct.
,
212
, pp.
598
608
.
21.
Rezaiee-Pajand
,
M.
,
Masoodi
,
A. R.
, and
Bambaeechee
,
M.
,
2019
, “
Tapered Beam-Column Analysis by Analytical Solution
,”
Proc. Inst. Civ. Eng. Struct. Build.
,
172
(
11
), pp.
789
804
.
22.
Griffis
,
M.
, and
Duffy
,
J.
,
1991
, “
Kinestatic Control: A Novel Theory for Simultaneously Regulating Force and Displacement
,”
J. Mech. Design
,
113
(
4
), pp.
508
515
.
23.
Ciblak
,
N.
, and
Lipkin
,
H.
,
1998
, “
Synthesis of Stiffnesses by Springs
,”
Proceedings of the ASME Design Engineering Technical Conferences
,
Atlanta, GA
,
Sept. 13–16
, pp.
1
10
.
24.
Huang
,
S.
, and
Schimmels
,
J. M.
,
2000
, “
The Bounds and Realization of Spatial Compliances Achieved With Simple Serial Elastic Mechanisms
,”
IEEE Trans. Robot. Autom.
,
16
(
1
), pp.
99
103
.