# Dropping Impact Characteristics Analysis of a Cubic Nonlinear Packaging System with a Cantilever Beam Type Elastic Critical Component with Concentrated Tip Mass

Authors: Meng Hao; An-Jun Chen
Journal: Shock and Vibration
 (2015)

---

## Abstract

A mathematical model for a cubic nonlinear packaging system with a cantilever beam type critical component with concentrated tip mass is proposed. The finite element method and the implicit finite difference method together with the Rung-Kutta method are applied to study the dropping impact dynamics of the critical component and the effect of system parameters, such as the value of the concentrated tip mass and the frequency of the main component, is discussed. The results show that the relative displacement and acceleration change remarkably with the length of the cantilever beam, and the maximum internal stress occurs at the joint end of the critical component. With the increase of the value of the concentrated tip mass and/or a higher frequency of the main component, the amplitudes of the responses increase obviously.


---

## Body

## 1. Introduction

Investigation of dynamic models and dropping impact characteristics of the packaging system is essential for the cushioning packaging design. In 1945, the concept of critical component was firstly presented by Mindlin [1], and the dropping impact characteristics of linear and nonlinear packaging systems were studied by established two different dynamic models, single degree of freedom and two degrees of freedom concentrated mass-spring models. Then, the dynamic properties of the two degrees of freedom concentrated mass-spring model with the critical component, such as the damage boundary curve and the three-dimensional shock spectrum, have been discussed by many researchers [2–4]. Wang et al. [5], Bernad et al. [6], and Wu and Chen [7] analyzed the vibration and shock characteristics of the multidegree of freedom packaging system, respectively. Suhir and Burke [8] and Wong et al. [9–11] established dynamic models by using systems with elastic plate part as object and the numerical results showed that the products with elastic parts which are described by concentrated mass-spring models are not accurate for the infinite degree of freedom of elastic parts. Therefore, the nonlinear coupling dynamics models, including the bar, simply supported beam, and cantilever beam type critical components, were suggested by Gao et al. [12–14], the corresponding differential equations of motion were derived and solved by using the explicit finite difference method (EFDM), and the distribution of displacement, acceleration, and internal stress responses of critical components were investigated. However, the cushion packaging of the electromechanical products with elastic cantilever beam with complicated boundary conditions has not attracted enough attention.In this paper, a mathematical model for the packaging system with a cantilever beam type elastic critical component with concentrated tip mass is suggested, and two different numerical methods, the Runge-Kutta method combined with the implicit finite difference method (IFDM) and the finite element method (FEM), are applied to analyze the displacement, acceleration, and internal stress responses. The effects of the value of the concentrated tip mass and the frequency of the main component on the dropping impact responses of the critical component are discussed.


## 2. Modeling, Equations, and Numerical Solution

Assuming the cantilever beam critical component as a uniform elastic body, the main component as a rigid body, the connection between them as a rigid connection and the cushioning material as a cubic nonlinear material, the model of the nonlinear system with critical components under the dropping impact is depicted in Figure1, where m
1 and m
2 denote the mass of main component and critical component, and y
1
(
t
) and y
2
(
x
,
t
) are the displacement response of main component and critical component, respectively. The cantilever beam has radius r and length L. The resilience function of cushioning material is F, F
=
k
0
y
1
+
e
y
1
3, where k
0 and e are the initial elastic constant and the nonlinear constant, respectively. H is the dropping height of system.Figure 1
The dropping model of packaging system with critical component.In order to facilitate the numerical analysis, the coordinate system is established, the static equilibrium position is treated as the origin of coordinate, and the downward direction is regarded as the positive direction. Because the cantilever beam is a slender beam, in the dropping impact process of the system, we will only consider the transverse vibration and bending deformation, neglecting the influence of transverse shearing distortion and the axial deformation. Based on this, equations of motion are expressed with two different methods and solved by FEM and FDM, respectively.


### 2.1. The Dynamic Model of Critical Component

#### 2.1.1. Finite Element Method (FEM)

Consider a cantilever beam with concentrated tip mass as illustrated in Figure2, treating the joint between the critical component and the main component as the original point and the length direction of the critical component as the x-axis direction. In the finite element model, the beam is divided into n elements of equal length as shown in Figure 2(a), where element nodes are marked as 1,2
,
…
,
i
,
i
+
1
,
…
,
n
,
n
+
1 in turn, and u
i and u
i
′ are the lateral displacement and the rotation about the cross section of the node i. In arbitrary element i, the length is l
=
L
/
n and each node has two degrees of freedom (Figure 2(b)). Since there are four end displacements (or degree of freedom), a cubic variation in displacement is assumed to approximate the node displacements, in the form [15]:(1)
u
x
=
a
0
+
a
1
x
l
+
a
2
x
l
2
+
a
3
x
l
3
.The four degrees of freedom corresponding to the displacements u
1
(
i
), u
2
(
i
) and the rotations u
1
′
(
i
), u
2
′
(
i
) are given by(2)
u
1
(
i
)
=
u
0
,
u
2
i
=
u
l
,
u
1
′
(
i
)
=
u
′
(
0
)
,
u
2
′
(
i
)
=
u
′
(
l
)
.Write the kinetic and strain energies of the elements as(3a)
T
i
=
1
2
∫
0
l
ρ
A
∂
2
u
∂
t
2
2
d
x
=
1
2
u
˙
i
T
k
u
˙
i
,
(3b)
U
i
=
1
2
∫
0
l
E
I
∂
2
w
∂
x
2
2
d
x
=
1
2
w
i
T
k
w
i
,where E is the elastic modulus, I is the bending moment of inertia, ρ is the mass density, and A is the cross-sectional area of the element. Here u
i is the absolute displacement vector of the i element, w
i
=
u
i
-
y
1
(
t
)
E
∗ (E
∗
=
[
1
0
1
0
]
T) is the relative displacement vector of the i element, and a dot over u
i represents time derivative of u
i. One can obtain, after substituting (1) into (2) and (3a) and (3b), the element stiffness k
i matrix and the element mass m
i [15] matrix as(4)
k
i
=
2
E
I
l
3
6
3
l
-
6
3
l
3
l
2
l
2
-
3
l
2
l
2
-
6
-
3
l
6
-
3
l
3
l
2
l
2
-
3
l
2
l
2
,
m
i
=
ρ
A
l
420
156
22
l
54
-
13
l
22
l
4
l
2
13
l
-
3
l
2
54
13
l
156
-
22
l
-
13
l
-
3
l
2
-
22
l
4
l
2
.Then, the kinetic and strain energies of the critical component are written as(5)
T
=
∑
i
=
1
n
T
i
=
1
2
∑
i
=
1
n
u
˙
i
T
m
i
u
˙
i
,
U
=
∑
i
=
1
n
U
i
=
1
2
∑
i
=
1
n
w
i
T
k
i
w
i
,where u
˙
i
T and w
i
T are transposed matrices of u
˙
i and w
i, respectively.Figure 2
The finite element division of cantilever beam component with tip mass.
(a)
Element division
(b)
Elementiu and w are assumed as the absolute displacement vector and the relative displacement vector of the critical component, respectively, and are given by(6)
u
=
u
1
u
1
′
⋯
u
n
+
1
u
n
+
1
′
T
,
w
=
w
1
w
1
′
⋯
w
n
+
1
w
n
+
1
′
T
.And, the relationship between u
i and u and w
i and w can be expressed as(7)
u
i
=
C
i
u
,
w
i
=
C
i
w
,where C
i
(
4
×
(
2
n
+
2
)
) is the extraction matrix of element i and written as(8)
C
i
=
0
⋯
1
t
h
e
i
c
o
l
u
m
n
0
0
0
⋯
0
0
⋯
0
1
0
0
⋯
0
0
⋯
0
0
1
0
⋯
0
0
⋯
0
0
0
1
⋯
0
.By substituting (7) and C
i into (5), the kinetic and strain energies of the critical component can be rewritten as(9)
T
=
∑
i
=
1
n
T
i
=
1
2
u
˙
T
M
u
˙
,
U
=
∑
i
=
1
n
U
i
=
1
2
w
T
K
w
,where M
(
2
(
n
+
1
)
×
2
(
n
+
1
)
) and K
(
2
(
n
+
1
)
×
2
(
n
+
1
)
) denote the critical component mass matrix and the critical component stiffness matrix and can be given by(10)
M
=
∑
i
=
1
n
C
i
T
m
i
C
i
,
K
=
∑
i
=
1
n
C
i
T
k
i
C
i
.Considering the effect of the concentrated tip mass on the critical component mass matrix M, m
2 is added to the (
2
n
+
1
)
×
(
2
n
+
1
) element of the M.The generalized force vector is the0, since the cantilever beam is not affected by the concentration or distribution force. Neglecting the damping, substituting the mass matrix M, the stiffness matrix K
, and the generalized force vector into Lagrange’s equations, the equation of motion for the critical component can be represented as(11)
M
u
¨
+
K
w
=
0
,where w
=
u
-
y
1
(
t
)
E, E
=
[
1
0
1
0
⋯
1
0
1
0
]
T. Thus, the dynamic response analysis can be transformed into solving the relative displacement (deformation of the critical component) response of the critical component under the inertia force. Substituting u
=
w
+
y
1
(
t
)
E, (11) becomes(12)
M
w
¨
+
K
w
=
-
y
¨
1
(
t
)
M
v
.The initial conditions considered here are given by w
t
=
0
=
[
0
0
⋯
0
]
T and w
¨
t
=
0
=
[
0
0
⋯
0
]
T, and the boundary conditions are specified by w
1
=
0 and w
1
′
=
0. However, (12) cannot be solved since the mass matrix M and the stiffness matrix K are singular matrices with no considering the boundary conditions. Hence, the boundary conditions have to be introduced in the solution. Three methods, the first method that is to eliminate the corresponding rows and columns of matrices of known displacements, the second that is to transform diagonal element into 1, and the third that is to multiply a large number, are usually used to incorporate the boundary conditions, in which the first method has a good applicability to the boundary condition of zero displacement [16]. Based on the first method, the effective mass matrix M
∗ and the effective stiffness matrix K
∗ can be obtained by deleting the first and second rows and columns of M and K. And (12) can be rewritten as(13)
M
∗
w
¨
+
K
∗
w
=
y
¨
1
(
t
)
M
∗
v
.Equation (13), a linear second-order differential equation set with constant coefficient, can be solved theoretically by the Runge-Kutta method. But the efficiency is low, because the dimension of the dynamic equation set is often too high. The effective methods can be divided into two classes: direct integration methods and the mode superposition method, in which direct integration methods mainly include central difference, Houbolt, Wilson, and Newmark methods [15]. Considering the characteristic of equations, the Newmark method is employed in the study to solve (13) and the numerical procedure is implemented in a MATLAB program.In order to investigate natural frequencies and mode shapes of the free vibration, (13) gives the following generalized eigenvalue equation:(14)
K
∗
φ
=
ω
2
M
∗
φ
,where φ represents the amplitudes of the displacements (called the mode shape) and the ω denotes the nature frequency of vibration. Equation (14) can be solved by using the subspace iteration method [16] effectively.According to (13) and equations of motion of main component (mentioned in Section 2.2), the responses of main component and the relative displacement vector w
k
∗ and the relative acceleration vector w
¨
k
∗ of the critical component at any time t
=
k
τ (where τ is the time interval and k is the step number) are obtained. Hence, ignoring the rotations in w
¨
k
∗, the transverse relative displacement matrix w
2
(
n
×
N
) (N is the total step number) is written as(15)
w
2
=
J
w
0
∗
w
1
∗
⋯
w
k
∗
⋯
,where(16)
J
=
1
0
0
0
0
0
⋯
0
0
0
0
1
0
0
0
⋯
0
0
0
0
0
0
1
0
⋯
0
0
⋮
⋮
⋮
⋮
⋮
⋮
⋱
⋮
⋮
0
0
0
0
0
0
⋯
1
0in which J
(
n
×
2
n
) is the extraction matrix of the relative displacement vector.Then, the transverse absolute displacement matrix is given by(17)
u
2
=
w
2
+
y
1
:
,
t
,where u
2
=
[
u
0
u
1
⋯
u
k
⋯
], the relationship between the absolute displacement vector u
k of the critical component, the displacement response of main component y
2, and the screen function u
i
k (mentioned in Section 2.1.2) can be expressed as u
k
(
i
)
=
u
i
k
=
y
2
(
i
l
,
k
τ
), and y
1
(
:
,
t
) is expanding matrix of y
1
(
t
), giving(18)
y
2
(
:
,
t
)
=
y
1
(
0
)
y
1
(
τ
)
⋯
y
1
(
k
τ
)
⋯
y
1
(
0
)
y
1
(
τ
)
⋯
y
1
(
k
τ
)
⋯
⋮
⋮
⋱
⋮
⋮
y
1
(
0
)
y
1
(
τ
)
⋯
y
1
(
k
τ
)
⋯
.At last, the internal stress of the critical component can be obtained by(19)
σ
=
-
E
r
∂
2
y
2
∂
x
2
.


#### 2.1.2. Finite Difference Method (FDM)

As for the force analysis, the cantilever beam critical component with concentrated tip mass is depicted as in Figure3, treating the joint between the critical component and the main component as the original point and the length direction of the critical component as the x-axis direction, where y
2
(
x
,
t
) is the displacement distribution function along the x-axial. Consider the representative elemental volume d
x at the x coordinates as illustrated in Figure 3(b) and the force-balance equation is derived by(20)
Q
+
∂
Q
∂
x
d
x
-
Q
=
ρ
A
d
x
∂
2
y
2
∂
x
2
,simplified as(21)
∂
Q
∂
x
=
ρ
A
∂
2
y
2
∂
x
2
,where Q is the cross-sectional shear force of the critical component ρ and A are the mass density and the cross-sectional area, respectively.Figure 3
Force analysis of the critical component.
(a)
(b)Treat the midpoint of the right section of the representative elemental volume as the center of moment. The torque equilibrium equation is established as(22)
(
M
+
∂
M
∂
x
d
x
-
M
)
-
Q
d
x
=
0
,simplified as(23)
Q
=
∂
M
∂
x
,where M
=
-
E
I
∂
2
y
2
/
∂
x
2 is the cross-sectional bending moment of the critical component.Substituting the bending momentM and the shear force Q into (21), the bending vibration equation of the critical component is given by(24)
E
I
∂
4
y
2
∂
x
4
+
ρ
A
∂
2
y
2
∂
t
2
=
0
,where E is the elastic modulus and I is the bending moment of inertia.The initial conditions are considered here as(25)
y
2
(
x
,
0
)
=
0
,
∂
y
2
(
x
,
t
)
∂
t
t
=
0
=
2
g
H
,where g is the acceleration of gravity.The boundary conditions are also specified as(26)
y
2
(
0
,
t
)
=
y
1
,
∂
y
2
(
x
,
t
)
∂
x
x
=
0
=
0
,
∂
2
y
2
(
x
,
t
)
∂
x
2
x
=
L
=
0
,
E
I
∂
3
y
2
(
x
,
t
)
∂
x
3
x
=
L
=
m
2
∂
3
y
2
(
x
,
t
)
∂
t
3
.Solutions of the partial differential equation (24), with nonlinear boundary conditions, are not trivial. Turner et al. [17] analyzed high-frequency response of microscope cantilevers by using the explicit finite difference method (EFDM), which also is presented by Lu et al. [12–14]. A new implicit finite difference scheme is applied here to obtain the numerical approximation, which is useful for boundary conditions (26). Firstly, the solution domain is divided into rectangular net by parallel lines, x
=
x
i
=
i
l (i
=
0,1
,
2
,
…
,
N) and t
=
t
k
=
k
τ (k
=
0,1
,
…), where l
=
L
/
N and τ
>
0. The node coordinate is defined by (
x
i
,
t
k
). Then, assuming u
i
k as the function of the node (
x
i
,
t
k
), (24) is discretized by the implicit difference method as follows:(27)
a
(
u
i
+
2
k
+
1
+
u
i
-
2
k
+
1
)
-
4
(
u
i
+
1
k
+
1
+
u
i
-
1
k
+
1
)
+
6
u
i
k
+
1
l
4
+
u
i
k
+
1
-
2
u
i
k
+
u
i
k
-
1
τ
2
=
0
,simplified as(28)
b
u
i
+
2
k
+
1
-
4
b
u
i
+
1
k
+
1
+
(
6
b
+
0.5
)
u
i
k
+
1
-
4
b
u
i
-
1
k
+
1
+
b
u
i
-
2
k
+
1
=
u
i
k
-
0.5
u
i
k
-
1
,
k
k
k
k
k
k
k
k
k
k
k
k
k
(
i
=
2
,
…
,
N
-
2
)
,where the constants a and b are given by(29)
a
=
ρ
A
E
I
,
b
=
a
τ
2
l
4
.In addition, the truncation error of (28) is O
(
l
4
+
τ
2
). Discretizing versions of the boundary conditions, the “real” nodes, u
1
k
+
1, u
N
-
1
k
+
1, and u
N
k
+
1, are given by(30a)
b
u
3
k
+
1
-
4
b
u
2
k
+
1
+
(
6
b
+
0.5
)
u
1
k
+
1
=
u
1
k
-
0.5
u
1
k
-
1
+
3
b
u
0
k
+
1
,
(30b)
2
b
u
N
k
+
1
+
(
5
b
+
0.5
)
u
N
-
1
k
+
1
-
4
b
u
N
-
2
k
+
1
+
b
u
i
k
+
1
=
u
i
-
1
k
-
0.5
u
N
-
1
k
-
1
,
(30c)
(
b
c
+
2
b
+
0.5
)
u
N
k
+
1
-
4
b
u
N
-
1
k
+
1
+
2
b
u
N
-
2
k
+
1
=
1
+
2
b
c
u
N
k
-
0.5
+
b
c
u
N
k
-
1
,where c
=
2
m
2
l
3
/
E
I
τ
2.An equation of thei
=
0 node is not required since its motion is given directly by the boundary condition in (26). Hence, the approximation of (24) can be written by(31)
A
U
=
B
,where A, U, and B are, respectively, given by(32)
U
=
[
u
1
k
+
1
,
u
2
k
+
1
,
…
,
u
N
-
1
k
+
1
,
u
N
k
+
1
]
T
,
A
=
6
b
+
0.5
-
4
b
b
0
·
·
·
·
-
4
b
6
b
+
0.5
-
4
b
b
0
·
·
·
b
-
4
b
6
b
+
0.5
-
4
b
b
0
·
·
0
·
·
·
·
·
0
·
·
0
b
-
4
b
6
b
+
0.5
-
4
b
b
0
·
·
0
b
-
4
b
6
b
+
0.5
-
4
b
b
·
·
·
0
b
-
4
b
5
b
+
0.5
-
2
b
·
·
·
·
0
2
b
-
4
b
2
b
+
0.5
+
b
c
,
B
=
u
1
k
-
0.5
u
1
k
-
1
+
3
b
u
0
k
+
1
u
2
k
-
0.5
u
2
k
-
1
-
b
u
0
k
+
1
u
3
n
-
0.5
u
3
k
-
1
⋮
u
N
-
3
k
-
0.5
u
N
-
3
k
-
1
u
N
-
2
k
-
0.5
u
N
-
2
k
-
1
u
N
-
1
n
-
0.5
u
N
-
1
k
-
1
1
+
2
b
u
N
k
-
0.5
+
b
c
u
N
k
-
1
.
.Since (31) is a linear system, having N
-
1 equations and N
-
1 unknowns, and the coefficient matrix A is nonsingular matrix, the linear system has a unique solution. Then, the absolute displacements of each node (excepted i
=
0 node) of the critical component can be obtained iteratively.The displacement responses of the critical component relative to the main component are given by(33)
Y
2
(
x
,
t
)
=
y
2
(
x
,
t
)
-
y
1
(
t
)
.And the relative acceleration of the critical component is expressed as(34)
Y
¨
2
(
x
,
t
)
=
y
¨
2
(
x
,
t
)
-
y
¨
1
(
t
)
,where y
¨
2
(
x
,
t
) is the absolute acceleration response of the critical component. At last, the internal stress of the critical component is derived from (19).


### 2.2. The Dynamic Models of Main Component

Considering the dynamic effect of the critical component, the dropping equation of motion of the main component is given by(35)
m
1
y
¨
1
+
F
+
E
I
∂
3
y
2
(
x
,
t
)
∂
x
3
x
=
0
=
0and the initial conditions can be defined as y
1
(
0
)
=
0, y
˙
1
(
0
)
=
2
g
H. However, sometimes the mass of the critical component is much less than the main component, and the effect of the critical component on the main component can be ignored. Then, (35) can be rewritten as(36)
m
1
y
¨
1
+
F
=
0
.We know the numerical solutions of the second-order differential equation by using the Runge-Kutta method that has a good accuracy [18]. So (36) can be solved directly by the Runge-Kutta method. But (35) only can be solved after discretizing the effect formula by the FDM. The effect formula is discretized as(37)
E
I
∂
3
y
2
(
x
,
t
)
∂
x
3
x
=
0
=
E
I
u
3
k
-
2
u
2
k
+
2
u
0
k
-
u
-
1
k
2
l
3
,where u
-
1
k is not actual location on the cantilever beam. This “fictitious” node can be defined by the “real” nodes near the end of the cantilever beam as(38)
u
-
1
k
=
-
u
1
k
+
1
-
2
u
1
k
+
u
1
k
-
1
B
-
u
3
k
-
4
u
2
k
+
6
u
1
k
-
4
u
0
k
,where B
=
E
I
τ
2
/
(
ρ
A
l
4
). Equation (37) is rewritten by coming equation (38), shown as(39)
E
I
∂
3
y
2
(
x
,
t
)
∂
x
3
x
=
0
=
E
I
2
u
3
k
-
6
u
2
k
+
4
u
1
k
+
(
u
1
k
+
1
-
2
u
1
k
+
u
1
k
-
1
)
/
B
2
l
3
.


## 2.1. The Dynamic Model of Critical Component

### 2.1.1. Finite Element Method (FEM)

Consider a cantilever beam with concentrated tip mass as illustrated in Figure2, treating the joint between the critical component and the main component as the original point and the length direction of the critical component as the x-axis direction. In the finite element model, the beam is divided into n elements of equal length as shown in Figure 2(a), where element nodes are marked as 1,2
,
…
,
i
,
i
+
1
,
…
,
n
,
n
+
1 in turn, and u
i and u
i
′ are the lateral displacement and the rotation about the cross section of the node i. In arbitrary element i, the length is l
=
L
/
n and each node has two degrees of freedom (Figure 2(b)). Since there are four end displacements (or degree of freedom), a cubic variation in displacement is assumed to approximate the node displacements, in the form [15]:(1)
u
x
=
a
0
+
a
1
x
l
+
a
2
x
l
2
+
a
3
x
l
3
.The four degrees of freedom corresponding to the displacements u
1
(
i
), u
2
(
i
) and the rotations u
1
′
(
i
), u
2
′
(
i
) are given by(2)
u
1
(
i
)
=
u
0
,
u
2
i
=
u
l
,
u
1
′
(
i
)
=
u
′
(
0
)
,
u
2
′
(
i
)
=
u
′
(
l
)
.Write the kinetic and strain energies of the elements as(3a)
T
i
=
1
2
∫
0
l
ρ
A
∂
2
u
∂
t
2
2
d
x
=
1
2
u
˙
i
T
k
u
˙
i
,
(3b)
U
i
=
1
2
∫
0
l
E
I
∂
2
w
∂
x
2
2
d
x
=
1
2
w
i
T
k
w
i
,where E is the elastic modulus, I is the bending moment of inertia, ρ is the mass density, and A is the cross-sectional area of the element. Here u
i is the absolute displacement vector of the i element, w
i
=
u
i
-
y
1
(
t
)
E
∗ (E
∗
=
[
1
0
1
0
]
T) is the relative displacement vector of the i element, and a dot over u
i represents time derivative of u
i. One can obtain, after substituting (1) into (2) and (3a) and (3b), the element stiffness k
i matrix and the element mass m
i [15] matrix as(4)
k
i
=
2
E
I
l
3
6
3
l
-
6
3
l
3
l
2
l
2
-
3
l
2
l
2
-
6
-
3
l
6
-
3
l
3
l
2
l
2
-
3
l
2
l
2
,
m
i
=
ρ
A
l
420
156
22
l
54
-
13
l
22
l
4
l
2
13
l
-
3
l
2
54
13
l
156
-
22
l
-
13
l
-
3
l
2
-
22
l
4
l
2
.Then, the kinetic and strain energies of the critical component are written as(5)
T
=
∑
i
=
1
n
T
i
=
1
2
∑
i
=
1
n
u
˙
i
T
m
i
u
˙
i
,
U
=
∑
i
=
1
n
U
i
=
1
2
∑
i
=
1
n
w
i
T
k
i
w
i
,where u
˙
i
T and w
i
T are transposed matrices of u
˙
i and w
i, respectively.Figure 2
The finite element division of cantilever beam component with tip mass.
(a)
Element division
(b)
Elementiu and w are assumed as the absolute displacement vector and the relative displacement vector of the critical component, respectively, and are given by(6)
u
=
u
1
u
1
′
⋯
u
n
+
1
u
n
+
1
′
T
,
w
=
w
1
w
1
′
⋯
w
n
+
1
w
n
+
1
′
T
.And, the relationship between u
i and u and w
i and w can be expressed as(7)
u
i
=
C
i
u
,
w
i
=
C
i
w
,where C
i
(
4
×
(
2
n
+
2
)
) is the extraction matrix of element i and written as(8)
C
i
=
0
⋯
1
t
h
e
i
c
o
l
u
m
n
0
0
0
⋯
0
0
⋯
0
1
0
0
⋯
0
0
⋯
0
0
1
0
⋯
0
0
⋯
0
0
0
1
⋯
0
.By substituting (7) and C
i into (5), the kinetic and strain energies of the critical component can be rewritten as(9)
T
=
∑
i
=
1
n
T
i
=
1
2
u
˙
T
M
u
˙
,
U
=
∑
i
=
1
n
U
i
=
1
2
w
T
K
w
,where M
(
2
(
n
+
1
)
×
2
(
n
+
1
)
) and K
(
2
(
n
+
1
)
×
2
(
n
+
1
)
) denote the critical component mass matrix and the critical component stiffness matrix and can be given by(10)
M
=
∑
i
=
1
n
C
i
T
m
i
C
i
,
K
=
∑
i
=
1
n
C
i
T
k
i
C
i
.Considering the effect of the concentrated tip mass on the critical component mass matrix M, m
2 is added to the (
2
n
+
1
)
×
(
2
n
+
1
) element of the M.The generalized force vector is the0, since the cantilever beam is not affected by the concentration or distribution force. Neglecting the damping, substituting the mass matrix M, the stiffness matrix K
, and the generalized force vector into Lagrange’s equations, the equation of motion for the critical component can be represented as(11)
M
u
¨
+
K
w
=
0
,where w
=
u
-
y
1
(
t
)
E, E
=
[
1
0
1
0
⋯
1
0
1
0
]
T. Thus, the dynamic response analysis can be transformed into solving the relative displacement (deformation of the critical component) response of the critical component under the inertia force. Substituting u
=
w
+
y
1
(
t
)
E, (11) becomes(12)
M
w
¨
+
K
w
=
-
y
¨
1
(
t
)
M
v
.The initial conditions considered here are given by w
t
=
0
=
[
0
0
⋯
0
]
T and w
¨
t
=
0
=
[
0
0
⋯
0
]
T, and the boundary conditions are specified by w
1
=
0 and w
1
′
=
0. However, (12) cannot be solved since the mass matrix M and the stiffness matrix K are singular matrices with no considering the boundary conditions. Hence, the boundary conditions have to be introduced in the solution. Three methods, the first method that is to eliminate the corresponding rows and columns of matrices of known displacements, the second that is to transform diagonal element into 1, and the third that is to multiply a large number, are usually used to incorporate the boundary conditions, in which the first method has a good applicability to the boundary condition of zero displacement [16]. Based on the first method, the effective mass matrix M
∗ and the effective stiffness matrix K
∗ can be obtained by deleting the first and second rows and columns of M and K. And (12) can be rewritten as(13)
M
∗
w
¨
+
K
∗
w
=
y
¨
1
(
t
)
M
∗
v
.Equation (13), a linear second-order differential equation set with constant coefficient, can be solved theoretically by the Runge-Kutta method. But the efficiency is low, because the dimension of the dynamic equation set is often too high. The effective methods can be divided into two classes: direct integration methods and the mode superposition method, in which direct integration methods mainly include central difference, Houbolt, Wilson, and Newmark methods [15]. Considering the characteristic of equations, the Newmark method is employed in the study to solve (13) and the numerical procedure is implemented in a MATLAB program.In order to investigate natural frequencies and mode shapes of the free vibration, (13) gives the following generalized eigenvalue equation:(14)
K
∗
φ
=
ω
2
M
∗
φ
,where φ represents the amplitudes of the displacements (called the mode shape) and the ω denotes the nature frequency of vibration. Equation (14) can be solved by using the subspace iteration method [16] effectively.According to (13) and equations of motion of main component (mentioned in Section 2.2), the responses of main component and the relative displacement vector w
k
∗ and the relative acceleration vector w
¨
k
∗ of the critical component at any time t
=
k
τ (where τ is the time interval and k is the step number) are obtained. Hence, ignoring the rotations in w
¨
k
∗, the transverse relative displacement matrix w
2
(
n
×
N
) (N is the total step number) is written as(15)
w
2
=
J
w
0
∗
w
1
∗
⋯
w
k
∗
⋯
,where(16)
J
=
1
0
0
0
0
0
⋯
0
0
0
0
1
0
0
0
⋯
0
0
0
0
0
0
1
0
⋯
0
0
⋮
⋮
⋮
⋮
⋮
⋮
⋱
⋮
⋮
0
0
0
0
0
0
⋯
1
0in which J
(
n
×
2
n
) is the extraction matrix of the relative displacement vector.Then, the transverse absolute displacement matrix is given by(17)
u
2
=
w
2
+
y
1
:
,
t
,where u
2
=
[
u
0
u
1
⋯
u
k
⋯
], the relationship between the absolute displacement vector u
k of the critical component, the displacement response of main component y
2, and the screen function u
i
k (mentioned in Section 2.1.2) can be expressed as u
k
(
i
)
=
u
i
k
=
y
2
(
i
l
,
k
τ
), and y
1
(
:
,
t
) is expanding matrix of y
1
(
t
), giving(18)
y
2
(
:
,
t
)
=
y
1
(
0
)
y
1
(
τ
)
⋯
y
1
(
k
τ
)
⋯
y
1
(
0
)
y
1
(
τ
)
⋯
y
1
(
k
τ
)
⋯
⋮
⋮
⋱
⋮
⋮
y
1
(
0
)
y
1
(
τ
)
⋯
y
1
(
k
τ
)
⋯
.At last, the internal stress of the critical component can be obtained by(19)
σ
=
-
E
r
∂
2
y
2
∂
x
2
.


### 2.1.2. Finite Difference Method (FDM)

As for the force analysis, the cantilever beam critical component with concentrated tip mass is depicted as in Figure3, treating the joint between the critical component and the main component as the original point and the length direction of the critical component as the x-axis direction, where y
2
(
x
,
t
) is the displacement distribution function along the x-axial. Consider the representative elemental volume d
x at the x coordinates as illustrated in Figure 3(b) and the force-balance equation is derived by(20)
Q
+
∂
Q
∂
x
d
x
-
Q
=
ρ
A
d
x
∂
2
y
2
∂
x
2
,simplified as(21)
∂
Q
∂
x
=
ρ
A
∂
2
y
2
∂
x
2
,where Q is the cross-sectional shear force of the critical component ρ and A are the mass density and the cross-sectional area, respectively.Figure 3
Force analysis of the critical component.
(a)
(b)Treat the midpoint of the right section of the representative elemental volume as the center of moment. The torque equilibrium equation is established as(22)
(
M
+
∂
M
∂
x
d
x
-
M
)
-
Q
d
x
=
0
,simplified as(23)
Q
=
∂
M
∂
x
,where M
=
-
E
I
∂
2
y
2
/
∂
x
2 is the cross-sectional bending moment of the critical component.Substituting the bending momentM and the shear force Q into (21), the bending vibration equation of the critical component is given by(24)
E
I
∂
4
y
2
∂
x
4
+
ρ
A
∂
2
y
2
∂
t
2
=
0
,where E is the elastic modulus and I is the bending moment of inertia.The initial conditions are considered here as(25)
y
2
(
x
,
0
)
=
0
,
∂
y
2
(
x
,
t
)
∂
t
t
=
0
=
2
g
H
,where g is the acceleration of gravity.The boundary conditions are also specified as(26)
y
2
(
0
,
t
)
=
y
1
,
∂
y
2
(
x
,
t
)
∂
x
x
=
0
=
0
,
∂
2
y
2
(
x
,
t
)
∂
x
2
x
=
L
=
0
,
E
I
∂
3
y
2
(
x
,
t
)
∂
x
3
x
=
L
=
m
2
∂
3
y
2
(
x
,
t
)
∂
t
3
.Solutions of the partial differential equation (24), with nonlinear boundary conditions, are not trivial. Turner et al. [17] analyzed high-frequency response of microscope cantilevers by using the explicit finite difference method (EFDM), which also is presented by Lu et al. [12–14]. A new implicit finite difference scheme is applied here to obtain the numerical approximation, which is useful for boundary conditions (26). Firstly, the solution domain is divided into rectangular net by parallel lines, x
=
x
i
=
i
l (i
=
0,1
,
2
,
…
,
N) and t
=
t
k
=
k
τ (k
=
0,1
,
…), where l
=
L
/
N and τ
>
0. The node coordinate is defined by (
x
i
,
t
k
). Then, assuming u
i
k as the function of the node (
x
i
,
t
k
), (24) is discretized by the implicit difference method as follows:(27)
a
(
u
i
+
2
k
+
1
+
u
i
-
2
k
+
1
)
-
4
(
u
i
+
1
k
+
1
+
u
i
-
1
k
+
1
)
+
6
u
i
k
+
1
l
4
+
u
i
k
+
1
-
2
u
i
k
+
u
i
k
-
1
τ
2
=
0
,simplified as(28)
b
u
i
+
2
k
+
1
-
4
b
u
i
+
1
k
+
1
+
(
6
b
+
0.5
)
u
i
k
+
1
-
4
b
u
i
-
1
k
+
1
+
b
u
i
-
2
k
+
1
=
u
i
k
-
0.5
u
i
k
-
1
,
k
k
k
k
k
k
k
k
k
k
k
k
k
(
i
=
2
,
…
,
N
-
2
)
,where the constants a and b are given by(29)
a
=
ρ
A
E
I
,
b
=
a
τ
2
l
4
.In addition, the truncation error of (28) is O
(
l
4
+
τ
2
). Discretizing versions of the boundary conditions, the “real” nodes, u
1
k
+
1, u
N
-
1
k
+
1, and u
N
k
+
1, are given by(30a)
b
u
3
k
+
1
-
4
b
u
2
k
+
1
+
(
6
b
+
0.5
)
u
1
k
+
1
=
u
1
k
-
0.5
u
1
k
-
1
+
3
b
u
0
k
+
1
,
(30b)
2
b
u
N
k
+
1
+
(
5
b
+
0.5
)
u
N
-
1
k
+
1
-
4
b
u
N
-
2
k
+
1
+
b
u
i
k
+
1
=
u
i
-
1
k
-
0.5
u
N
-
1
k
-
1
,
(30c)
(
b
c
+
2
b
+
0.5
)
u
N
k
+
1
-
4
b
u
N
-
1
k
+
1
+
2
b
u
N
-
2
k
+
1
=
1
+
2
b
c
u
N
k
-
0.5
+
b
c
u
N
k
-
1
,where c
=
2
m
2
l
3
/
E
I
τ
2.An equation of thei
=
0 node is not required since its motion is given directly by the boundary condition in (26). Hence, the approximation of (24) can be written by(31)
A
U
=
B
,where A, U, and B are, respectively, given by(32)
U
=
[
u
1
k
+
1
,
u
2
k
+
1
,
…
,
u
N
-
1
k
+
1
,
u
N
k
+
1
]
T
,
A
=
6
b
+
0.5
-
4
b
b
0
·
·
·
·
-
4
b
6
b
+
0.5
-
4
b
b
0
·
·
·
b
-
4
b
6
b
+
0.5
-
4
b
b
0
·
·
0
·
·
·
·
·
0
·
·
0
b
-
4
b
6
b
+
0.5
-
4
b
b
0
·
·
0
b
-
4
b
6
b
+
0.5
-
4
b
b
·
·
·
0
b
-
4
b
5
b
+
0.5
-
2
b
·
·
·
·
0
2
b
-
4
b
2
b
+
0.5
+
b
c
,
B
=
u
1
k
-
0.5
u
1
k
-
1
+
3
b
u
0
k
+
1
u
2
k
-
0.5
u
2
k
-
1
-
b
u
0
k
+
1
u
3
n
-
0.5
u
3
k
-
1
⋮
u
N
-
3
k
-
0.5
u
N
-
3
k
-
1
u
N
-
2
k
-
0.5
u
N
-
2
k
-
1
u
N
-
1
n
-
0.5
u
N
-
1
k
-
1
1
+
2
b
u
N
k
-
0.5
+
b
c
u
N
k
-
1
.
.Since (31) is a linear system, having N
-
1 equations and N
-
1 unknowns, and the coefficient matrix A is nonsingular matrix, the linear system has a unique solution. Then, the absolute displacements of each node (excepted i
=
0 node) of the critical component can be obtained iteratively.The displacement responses of the critical component relative to the main component are given by(33)
Y
2
(
x
,
t
)
=
y
2
(
x
,
t
)
-
y
1
(
t
)
.And the relative acceleration of the critical component is expressed as(34)
Y
¨
2
(
x
,
t
)
=
y
¨
2
(
x
,
t
)
-
y
¨
1
(
t
)
,where y
¨
2
(
x
,
t
) is the absolute acceleration response of the critical component. At last, the internal stress of the critical component is derived from (19).


## 2.1.1. Finite Element Method (FEM)

Consider a cantilever beam with concentrated tip mass as illustrated in Figure2, treating the joint between the critical component and the main component as the original point and the length direction of the critical component as the x-axis direction. In the finite element model, the beam is divided into n elements of equal length as shown in Figure 2(a), where element nodes are marked as 1,2
,
…
,
i
,
i
+
1
,
…
,
n
,
n
+
1 in turn, and u
i and u
i
′ are the lateral displacement and the rotation about the cross section of the node i. In arbitrary element i, the length is l
=
L
/
n and each node has two degrees of freedom (Figure 2(b)). Since there are four end displacements (or degree of freedom), a cubic variation in displacement is assumed to approximate the node displacements, in the form [15]:(1)
u
x
=
a
0
+
a
1
x
l
+
a
2
x
l
2
+
a
3
x
l
3
.The four degrees of freedom corresponding to the displacements u
1
(
i
), u
2
(
i
) and the rotations u
1
′
(
i
), u
2
′
(
i
) are given by(2)
u
1
(
i
)
=
u
0
,
u
2
i
=
u
l
,
u
1
′
(
i
)
=
u
′
(
0
)
,
u
2
′
(
i
)
=
u
′
(
l
)
.Write the kinetic and strain energies of the elements as(3a)
T
i
=
1
2
∫
0
l
ρ
A
∂
2
u
∂
t
2
2
d
x
=
1
2
u
˙
i
T
k
u
˙
i
,
(3b)
U
i
=
1
2
∫
0
l
E
I
∂
2
w
∂
x
2
2
d
x
=
1
2
w
i
T
k
w
i
,where E is the elastic modulus, I is the bending moment of inertia, ρ is the mass density, and A is the cross-sectional area of the element. Here u
i is the absolute displacement vector of the i element, w
i
=
u
i
-
y
1
(
t
)
E
∗ (E
∗
=
[
1
0
1
0
]
T) is the relative displacement vector of the i element, and a dot over u
i represents time derivative of u
i. One can obtain, after substituting (1) into (2) and (3a) and (3b), the element stiffness k
i matrix and the element mass m
i [15] matrix as(4)
k
i
=
2
E
I
l
3
6
3
l
-
6
3
l
3
l
2
l
2
-
3
l
2
l
2
-
6
-
3
l
6
-
3
l
3
l
2
l
2
-
3
l
2
l
2
,
m
i
=
ρ
A
l
420
156
22
l
54
-
13
l
22
l
4
l
2
13
l
-
3
l
2
54
13
l
156
-
22
l
-
13
l
-
3
l
2
-
22
l
4
l
2
.Then, the kinetic and strain energies of the critical component are written as(5)
T
=
∑
i
=
1
n
T
i
=
1
2
∑
i
=
1
n
u
˙
i
T
m
i
u
˙
i
,
U
=
∑
i
=
1
n
U
i
=
1
2
∑
i
=
1
n
w
i
T
k
i
w
i
,where u
˙
i
T and w
i
T are transposed matrices of u
˙
i and w
i, respectively.Figure 2
The finite element division of cantilever beam component with tip mass.
(a)
Element division
(b)
Elementiu and w are assumed as the absolute displacement vector and the relative displacement vector of the critical component, respectively, and are given by(6)
u
=
u
1
u
1
′
⋯
u
n
+
1
u
n
+
1
′
T
,
w
=
w
1
w
1
′
⋯
w
n
+
1
w
n
+
1
′
T
.And, the relationship between u
i and u and w
i and w can be expressed as(7)
u
i
=
C
i
u
,
w
i
=
C
i
w
,where C
i
(
4
×
(
2
n
+
2
)
) is the extraction matrix of element i and written as(8)
C
i
=
0
⋯
1
t
h
e
i
c
o
l
u
m
n
0
0
0
⋯
0
0
⋯
0
1
0
0
⋯
0
0
⋯
0
0
1
0
⋯
0
0
⋯
0
0
0
1
⋯
0
.By substituting (7) and C
i into (5), the kinetic and strain energies of the critical component can be rewritten as(9)
T
=
∑
i
=
1
n
T
i
=
1
2
u
˙
T
M
u
˙
,
U
=
∑
i
=
1
n
U
i
=
1
2
w
T
K
w
,where M
(
2
(
n
+
1
)
×
2
(
n
+
1
)
) and K
(
2
(
n
+
1
)
×
2
(
n
+
1
)
) denote the critical component mass matrix and the critical component stiffness matrix and can be given by(10)
M
=
∑
i
=
1
n
C
i
T
m
i
C
i
,
K
=
∑
i
=
1
n
C
i
T
k
i
C
i
.Considering the effect of the concentrated tip mass on the critical component mass matrix M, m
2 is added to the (
2
n
+
1
)
×
(
2
n
+
1
) element of the M.The generalized force vector is the0, since the cantilever beam is not affected by the concentration or distribution force. Neglecting the damping, substituting the mass matrix M, the stiffness matrix K
, and the generalized force vector into Lagrange’s equations, the equation of motion for the critical component can be represented as(11)
M
u
¨
+
K
w
=
0
,where w
=
u
-
y
1
(
t
)
E, E
=
[
1
0
1
0
⋯
1
0
1
0
]
T. Thus, the dynamic response analysis can be transformed into solving the relative displacement (deformation of the critical component) response of the critical component under the inertia force. Substituting u
=
w
+
y
1
(
t
)
E, (11) becomes(12)
M
w
¨
+
K
w
=
-
y
¨
1
(
t
)
M
v
.The initial conditions considered here are given by w
t
=
0
=
[
0
0
⋯
0
]
T and w
¨
t
=
0
=
[
0
0
⋯
0
]
T, and the boundary conditions are specified by w
1
=
0 and w
1
′
=
0. However, (12) cannot be solved since the mass matrix M and the stiffness matrix K are singular matrices with no considering the boundary conditions. Hence, the boundary conditions have to be introduced in the solution. Three methods, the first method that is to eliminate the corresponding rows and columns of matrices of known displacements, the second that is to transform diagonal element into 1, and the third that is to multiply a large number, are usually used to incorporate the boundary conditions, in which the first method has a good applicability to the boundary condition of zero displacement [16]. Based on the first method, the effective mass matrix M
∗ and the effective stiffness matrix K
∗ can be obtained by deleting the first and second rows and columns of M and K. And (12) can be rewritten as(13)
M
∗
w
¨
+
K
∗
w
=
y
¨
1
(
t
)
M
∗
v
.Equation (13), a linear second-order differential equation set with constant coefficient, can be solved theoretically by the Runge-Kutta method. But the efficiency is low, because the dimension of the dynamic equation set is often too high. The effective methods can be divided into two classes: direct integration methods and the mode superposition method, in which direct integration methods mainly include central difference, Houbolt, Wilson, and Newmark methods [15]. Considering the characteristic of equations, the Newmark method is employed in the study to solve (13) and the numerical procedure is implemented in a MATLAB program.In order to investigate natural frequencies and mode shapes of the free vibration, (13) gives the following generalized eigenvalue equation:(14)
K
∗
φ
=
ω
2
M
∗
φ
,where φ represents the amplitudes of the displacements (called the mode shape) and the ω denotes the nature frequency of vibration. Equation (14) can be solved by using the subspace iteration method [16] effectively.According to (13) and equations of motion of main component (mentioned in Section 2.2), the responses of main component and the relative displacement vector w
k
∗ and the relative acceleration vector w
¨
k
∗ of the critical component at any time t
=
k
τ (where τ is the time interval and k is the step number) are obtained. Hence, ignoring the rotations in w
¨
k
∗, the transverse relative displacement matrix w
2
(
n
×
N
) (N is the total step number) is written as(15)
w
2
=
J
w
0
∗
w
1
∗
⋯
w
k
∗
⋯
,where(16)
J
=
1
0
0
0
0
0
⋯
0
0
0
0
1
0
0
0
⋯
0
0
0
0
0
0
1
0
⋯
0
0
⋮
⋮
⋮
⋮
⋮
⋮
⋱
⋮
⋮
0
0
0
0
0
0
⋯
1
0in which J
(
n
×
2
n
) is the extraction matrix of the relative displacement vector.Then, the transverse absolute displacement matrix is given by(17)
u
2
=
w
2
+
y
1
:
,
t
,where u
2
=
[
u
0
u
1
⋯
u
k
⋯
], the relationship between the absolute displacement vector u
k of the critical component, the displacement response of main component y
2, and the screen function u
i
k (mentioned in Section 2.1.2) can be expressed as u
k
(
i
)
=
u
i
k
=
y
2
(
i
l
,
k
τ
), and y
1
(
:
,
t
) is expanding matrix of y
1
(
t
), giving(18)
y
2
(
:
,
t
)
=
y
1
(
0
)
y
1
(
τ
)
⋯
y
1
(
k
τ
)
⋯
y
1
(
0
)
y
1
(
τ
)
⋯
y
1
(
k
τ
)
⋯
⋮
⋮
⋱
⋮
⋮
y
1
(
0
)
y
1
(
τ
)
⋯
y
1
(
k
τ
)
⋯
.At last, the internal stress of the critical component can be obtained by(19)
σ
=
-
E
r
∂
2
y
2
∂
x
2
.


## 2.1.2. Finite Difference Method (FDM)

As for the force analysis, the cantilever beam critical component with concentrated tip mass is depicted as in Figure3, treating the joint between the critical component and the main component as the original point and the length direction of the critical component as the x-axis direction, where y
2
(
x
,
t
) is the displacement distribution function along the x-axial. Consider the representative elemental volume d
x at the x coordinates as illustrated in Figure 3(b) and the force-balance equation is derived by(20)
Q
+
∂
Q
∂
x
d
x
-
Q
=
ρ
A
d
x
∂
2
y
2
∂
x
2
,simplified as(21)
∂
Q
∂
x
=
ρ
A
∂
2
y
2
∂
x
2
,where Q is the cross-sectional shear force of the critical component ρ and A are the mass density and the cross-sectional area, respectively.Figure 3
Force analysis of the critical component.
(a)
(b)Treat the midpoint of the right section of the representative elemental volume as the center of moment. The torque equilibrium equation is established as(22)
(
M
+
∂
M
∂
x
d
x
-
M
)
-
Q
d
x
=
0
,simplified as(23)
Q
=
∂
M
∂
x
,where M
=
-
E
I
∂
2
y
2
/
∂
x
2 is the cross-sectional bending moment of the critical component.Substituting the bending momentM and the shear force Q into (21), the bending vibration equation of the critical component is given by(24)
E
I
∂
4
y
2
∂
x
4
+
ρ
A
∂
2
y
2
∂
t
2
=
0
,where E is the elastic modulus and I is the bending moment of inertia.The initial conditions are considered here as(25)
y
2
(
x
,
0
)
=
0
,
∂
y
2
(
x
,
t
)
∂
t
t
=
0
=
2
g
H
,where g is the acceleration of gravity.The boundary conditions are also specified as(26)
y
2
(
0
,
t
)
=
y
1
,
∂
y
2
(
x
,
t
)
∂
x
x
=
0
=
0
,
∂
2
y
2
(
x
,
t
)
∂
x
2
x
=
L
=
0
,
E
I
∂
3
y
2
(
x
,
t
)
∂
x
3
x
=
L
=
m
2
∂
3
y
2
(
x
,
t
)
∂
t
3
.Solutions of the partial differential equation (24), with nonlinear boundary conditions, are not trivial. Turner et al. [17] analyzed high-frequency response of microscope cantilevers by using the explicit finite difference method (EFDM), which also is presented by Lu et al. [12–14]. A new implicit finite difference scheme is applied here to obtain the numerical approximation, which is useful for boundary conditions (26). Firstly, the solution domain is divided into rectangular net by parallel lines, x
=
x
i
=
i
l (i
=
0,1
,
2
,
…
,
N) and t
=
t
k
=
k
τ (k
=
0,1
,
…), where l
=
L
/
N and τ
>
0. The node coordinate is defined by (
x
i
,
t
k
). Then, assuming u
i
k as the function of the node (
x
i
,
t
k
), (24) is discretized by the implicit difference method as follows:(27)
a
(
u
i
+
2
k
+
1
+
u
i
-
2
k
+
1
)
-
4
(
u
i
+
1
k
+
1
+
u
i
-
1
k
+
1
)
+
6
u
i
k
+
1
l
4
+
u
i
k
+
1
-
2
u
i
k
+
u
i
k
-
1
τ
2
=
0
,simplified as(28)
b
u
i
+
2
k
+
1
-
4
b
u
i
+
1
k
+
1
+
(
6
b
+
0.5
)
u
i
k
+
1
-
4
b
u
i
-
1
k
+
1
+
b
u
i
-
2
k
+
1
=
u
i
k
-
0.5
u
i
k
-
1
,
k
k
k
k
k
k
k
k
k
k
k
k
k
(
i
=
2
,
…
,
N
-
2
)
,where the constants a and b are given by(29)
a
=
ρ
A
E
I
,
b
=
a
τ
2
l
4
.In addition, the truncation error of (28) is O
(
l
4
+
τ
2
). Discretizing versions of the boundary conditions, the “real” nodes, u
1
k
+
1, u
N
-
1
k
+
1, and u
N
k
+
1, are given by(30a)
b
u
3
k
+
1
-
4
b
u
2
k
+
1
+
(
6
b
+
0.5
)
u
1
k
+
1
=
u
1
k
-
0.5
u
1
k
-
1
+
3
b
u
0
k
+
1
,
(30b)
2
b
u
N
k
+
1
+
(
5
b
+
0.5
)
u
N
-
1
k
+
1
-
4
b
u
N
-
2
k
+
1
+
b
u
i
k
+
1
=
u
i
-
1
k
-
0.5
u
N
-
1
k
-
1
,
(30c)
(
b
c
+
2
b
+
0.5
)
u
N
k
+
1
-
4
b
u
N
-
1
k
+
1
+
2
b
u
N
-
2
k
+
1
=
1
+
2
b
c
u
N
k
-
0.5
+
b
c
u
N
k
-
1
,where c
=
2
m
2
l
3
/
E
I
τ
2.An equation of thei
=
0 node is not required since its motion is given directly by the boundary condition in (26). Hence, the approximation of (24) can be written by(31)
A
U
=
B
,where A, U, and B are, respectively, given by(32)
U
=
[
u
1
k
+
1
,
u
2
k
+
1
,
…
,
u
N
-
1
k
+
1
,
u
N
k
+
1
]
T
,
A
=
6
b
+
0.5
-
4
b
b
0
·
·
·
·
-
4
b
6
b
+
0.5
-
4
b
b
0
·
·
·
b
-
4
b
6
b
+
0.5
-
4
b
b
0
·
·
0
·
·
·
·
·
0
·
·
0
b
-
4
b
6
b
+
0.5
-
4
b
b
0
·
·
0
b
-
4
b
6
b
+
0.5
-
4
b
b
·
·
·
0
b
-
4
b
5
b
+
0.5
-
2
b
·
·
·
·
0
2
b
-
4
b
2
b
+
0.5
+
b
c
,
B
=
u
1
k
-
0.5
u
1
k
-
1
+
3
b
u
0
k
+
1
u
2
k
-
0.5
u
2
k
-
1
-
b
u
0
k
+
1
u
3
n
-
0.5
u
3
k
-
1
⋮
u
N
-
3
k
-
0.5
u
N
-
3
k
-
1
u
N
-
2
k
-
0.5
u
N
-
2
k
-
1
u
N
-
1
n
-
0.5
u
N
-
1
k
-
1
1
+
2
b
u
N
k
-
0.5
+
b
c
u
N
k
-
1
.
.Since (31) is a linear system, having N
-
1 equations and N
-
1 unknowns, and the coefficient matrix A is nonsingular matrix, the linear system has a unique solution. Then, the absolute displacements of each node (excepted i
=
0 node) of the critical component can be obtained iteratively.The displacement responses of the critical component relative to the main component are given by(33)
Y
2
(
x
,
t
)
=
y
2
(
x
,
t
)
-
y
1
(
t
)
.And the relative acceleration of the critical component is expressed as(34)
Y
¨
2
(
x
,
t
)
=
y
¨
2
(
x
,
t
)
-
y
¨
1
(
t
)
,where y
¨
2
(
x
,
t
) is the absolute acceleration response of the critical component. At last, the internal stress of the critical component is derived from (19).


## 2.2. The Dynamic Models of Main Component

Considering the dynamic effect of the critical component, the dropping equation of motion of the main component is given by(35)
m
1
y
¨
1
+
F
+
E
I
∂
3
y
2
(
x
,
t
)
∂
x
3
x
=
0
=
0and the initial conditions can be defined as y
1
(
0
)
=
0, y
˙
1
(
0
)
=
2
g
H. However, sometimes the mass of the critical component is much less than the main component, and the effect of the critical component on the main component can be ignored. Then, (35) can be rewritten as(36)
m
1
y
¨
1
+
F
=
0
.We know the numerical solutions of the second-order differential equation by using the Runge-Kutta method that has a good accuracy [18]. So (36) can be solved directly by the Runge-Kutta method. But (35) only can be solved after discretizing the effect formula by the FDM. The effect formula is discretized as(37)
E
I
∂
3
y
2
(
x
,
t
)
∂
x
3
x
=
0
=
E
I
u
3
k
-
2
u
2
k
+
2
u
0
k
-
u
-
1
k
2
l
3
,where u
-
1
k is not actual location on the cantilever beam. This “fictitious” node can be defined by the “real” nodes near the end of the cantilever beam as(38)
u
-
1
k
=
-
u
1
k
+
1
-
2
u
1
k
+
u
1
k
-
1
B
-
u
3
k
-
4
u
2
k
+
6
u
1
k
-
4
u
0
k
,where B
=
E
I
τ
2
/
(
ρ
A
l
4
). Equation (37) is rewritten by coming equation (38), shown as(39)
E
I
∂
3
y
2
(
x
,
t
)
∂
x
3
x
=
0
=
E
I
2
u
3
k
-
6
u
2
k
+
4
u
1
k
+
(
u
1
k
+
1
-
2
u
1
k
+
u
1
k
-
1
)
/
B
2
l
3
.


## 3. Numerical Example

In the present study, a numerical example is cited to illustrate the procedure, and the effects of some parameters are investigated. To validate the method, a small electromechanical product with the mass of the main componentm
1
=
10 kg, the mass of the concentrated tip mass m
2
=
0.02 kg, and the characteristics of the critical component as follows, cross-sectional radius r
=
0.003 m, length L
=
0.1 m, density ρ
=
7850 kg/m3, elastic modulus E
=
200 GPa, and elastic limit σ
e
=
180 MPa, is analyzed. The allowable stress of the critical component is 150 Mpa, since the material safety coefficient is defined as 1.2. And the resilience function of cushioning material is F, F
=
k
0
y
1
+
e
y
1
3, where k
0
=
100 N/cm, e
=
72 N/cm3 [19]. The whole packaging system drops from the height H
=
0.6 m.The dropping impact responses, such as the absolute and relative displacements, the absolute and relative accelerations, and the internal stress and their extremums, are obtained from the coupled equations of (13), (17), (19), and (35), as shown in Figure 4 and Table 1, respectively. From Figures 4(a), 4(b), and 4(c), we can see that the absolute displacement and acceleration do not obviously change along the cantilever beam. Figures 4(d) and 4(e) show that the relative displacement and acceleration vary significantly. Maximums of the displacement and acceleration appear at the free end of the cantilever beam, while the maximum internal stress occurs at the joint end.Table 1
The extremum values of dropping impact responses.
Absolute displacement(mm)
Relative acceleration(mm)
Absolute displacement(m/s2)
Relative acceleration (m/s2)
Internal stress(Mpa)
Extremum
41.03
0.41
−554.05
−17.06
80.76
Time(s)
0.0158
0.0157
0.0157
0.0150
0.0157
Location
Free end
Free end
Free end
Free end
Joint endFigure 4
The dropping impact responses of the critical component.
(a)
The absolute displacement
(b)
The absolute acceleration
(c)
The internal stress
(d)
The relative displacement
(e)
The relative accelerationDiscussing the difference between responses obtained by the FEM and the IFDM, the relative displacement and acceleration of the free end and the internal stress of the joint end of the cantilever beam are depicted in Figure5. Compared with responses obtained by the FEM, the relative errors of the relative displacement and acceleration and the internal stress are 7.32%, 5.22%, and 3.42%, respectively, and maximums occur at the similar time. Therefore, the FEM and the IFDM have good adaptability for solving these problems.Figure 5
The dropping impact responses of the critical component by using two different methods.
(a)
The relative displacement
(b)
The relative acceleration
(c)
The internal stressConsidering the effect of the critical component on the main component, Figure6 shows the difference between the relative displacement and acceleration of the free end and the internal stress of the joint end with various values of concentrated tip mass (0.01, 0.02, 0.03, and 0.04 kg). Clearly, with the value of the concentrated tip mass increasing, maximums of the relative displacement and acceleration and the internal stress move up. In addition, maximums occur at the similar time.Figure 6
The dropping impact responses of the critical component with different values of concentrated tip mass.
(a)
The relative displacement
(b)
The relative acceleration
(c)
The internal stressIn order to investigate the effect of the frequency of main component (defined by (36)), the packaging system with the same initial conditions and characteristics, mentioned in this part firstly, is considered. The nonlinear constant e of the cushioning material changes in a large range. Then, the corresponding frequency can be solved by the variational iteration method [20], as shown in Table 2. And the five frequencies of the critical component can be obtained by (14), as shown in Table 3. The variation of the relative displacement and acceleration of the free end and the internal stress of the joint end with different frequencies is presented in Figure 7. It can be seen that amplitudes of the relative displacement and acceleration increase significantly, when the frequency of the main component is closer to the first frequency. And when the frequency is 26.91 s−1, the internal stress has exceeded the allowable stress. Therefore, the cushioning material needs to be redesigned.Table 2
The vibration frequencies of the main component.
Nonlinear constante (N/cm2)
72
216
500
1000
Frequency (s−1)
14.32
18.55
22.72
26.91Table 3
The five frequencies of the critical component.
Order
1
2
3
4
5
Frequency (s−1)
195.63
1967.7
6144.4
12689
21613Figure 7
The dropping impact responses of the critical component with different frequencies.
(a)
The relative displacement
(b)
The relative acceleration
(c)
The internal stress


## 4. Conclusion

A finite element method and a finite difference method for the dropping impact analysis of the nonlinear packaging system with a cantilever beam type critical component with concentrated tip mass have been studied in this paper, having good accuracy needed for comparisons with each other. Considering the dropping impact responses of the critical component, the absolute displacement and acceleration do not obviously change along the cantilever beam. However, maximums of the relative displacement and acceleration appear at the free end of the cantilever beam, while the maximum internal stress occurs at the joint end. Therefore the structure of elastic critical component cannot be ignored.The effects of the concentrated tip mass and the frequency of main component on the response also have been discussed. Increasing the value of the concentrated tip mass and/or the frequency of the main component can lead to the remarkable rise of amplitudes of the displacement, acceleration, and the internal stress, especially when the frequency of the main component is closer to the first frequency of the critical component. Therefore, the selected cushioning material should make the frequency of the main component far away from the nature frequency of the critical component. The results may lead to a thorough understanding of the damage mechanism of packaged product and design of cushioning packaging.



---

*Source: 602984-2015-07-27.xml*