# Vibration Analysis of Collecting Electrodes by means of the Hybrid Finite Element Method

Authors: I. Adamiec-Wójcik; J. Awrejcewicz; A. Nowak; S. Wojciech
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The paper presents a hybrid finite element method of shell modeling in order to model collecting electrodes of electrostatic precipitators. The method uses the finite element method to reflect elastic features and the rigid finite element method in order to model mass features of the body. A model of dust removal systems of an electrostatic precipitator is presented. The system consists of two beams which are modeled by means of the rigid finite element method and a system of collecting shells modeled by means of the hybrid finite element method. The paper discusses both the procedure of deriving the equations of motion and the results of numerical simulations carried out in order to analyze vibrations of the whole system. Experimental verification of the model is also presented.


---

## Body

## 1. Introduction

Electrostatic precipitators (ESP) are one of the most effective devices used for removing dust from industrial emission gases. The general idea of electrostatic precipitation can be described as charging, collecting, and removing particles from gas flowing through a device. The efficiency of collecting particles can be as high as 99.9% in many cases. The performance and thus the effectiveness of ESP depend on many factors since the whole process is very complex and multidisciplinary, covering fluid dynamics, electrical engineering, mechanical engineering, and chemistry [1]. There are a number of publications dealing with the problems of gas flow and electrical field in ESP [2–5].One of the important factors influencing the effectiveness of ESP is the efficiency of periodic cleaning of the collecting electrodes, which is done by a rapping system (Figure1). Usually there are nine collecting electrodes in one section of the dust-removal system; however, the number of active electrodes can vary between one and nine. The dust is removed from the collecting electrodes by inducing vibrations propagated over the whole section. These vibrations are caused by an axial impact of a beater on a brushing bar (Figure 1). The problem of optimization of the beater has been discussed in [6]. Values of normal and tangent accelerations of vibrations at control points are used to evaluate the process of collecting particles and cleaning the electrodes. These accelerations are an important factor influencing the efficiency of the precipitations since when the vibrations are too small the dust stays on the surface of the electrodes causing corona discharge and a sudden drop in electric field potential. When the vibrations are too large secondary contamination of gases can occur.Figure 1
The dust-removal system.The subject of the paper is a new method called the hybrid finite element method, which has been developed to analyze the vibrations of the collecting electrodes. For this purpose two methods, the finite element method (FEM) and the rigid finite element method (RFEM), are combined. The FEM is used for calculation of the energy of spring-damping deformations, while the RFEM is used to reflect inertial features of the electrodes. The FEM [7] enables us to model complex shapes, and the availability of professional packages such as MSC NASTRAN, Abaqus, and ANSYS is an advantage of the method. In spite of continuous development and modification of the existing commercial software, there are reasons why it is useful to develop specialized models and computer programs. General (universal) software packages have very complex and developed interfaces and they are based on sophisticated mathematical models. The end-user (engineer, designer) often has to spend a long time and much effort in order to learn the package. For small and medium-sized enterprises which design and produce devices of similar structure and activity it is more convenient to use software dedicated to such a type of devices, with a simpler interface adapted to the needs and without many options never used. Moreover, special applications are often more numerically effective and require less advanced hardware. In comparison to commercial packages based on FEM, the software is less universal but more suited to solving specific tasks and fulfilling designers’ requirements [8, 9].The primary division is used in order to define spring deformation energies by means of rectangular shell elements. Since all the elements have the same length along axisx
^, the spring energy of element (
i
,
j
) is independent of i and angle α
j.One of the most evident advantages of the RFEM is the division of the complex structures into rigid elements. This enables us to reflect exactly the mass features of the system. And it makes the discretization intuitive and easy for interpretation. An important feature of this method is a diagonal mass matrix, which considerably simplifies numerical calculations in the case of large systems. Including additional masses (concentrated or distributed) in the system and modeling connections is also straightforward. Until now, the method has been used mainly to model multibody systems with flexible beam-like links. Both large deflections and physical nonlinearities have been taken into account [10–12]. At the same time, the RFEM reflects spring-damping features only approximately and the FEM is much more effective in this respect. This is the reason that in the approach presented the authors intend to combine advantages of both the methods.This paper is the first complete presentation of the hybrid method used for modeling vibrations of shells. The idea of the method lies in dividing not only the upper and bottom beams but also electrodes into finite rigid elements with six degrees of freedom (secondary division of the shells). The rigid generalized coordinates define translational and rotational displacements of the elements with respect to local coordinate systems. The axes of these systems coincide with the main central inertial axes of the elements. This is the reason that mass matrices both of the elements and of the whole system are diagonal. The classical FEM is used in order to reflect elastic features of the electrodes (primary division) and the continuous field of deformation is described by displacements of the nodes (elastic coordinates). Transformation formulae between the elastic and rigid coordinates are defined. Such an approach enables us to describe vibrations by means of displacements of the rigid elements, which are the generalized coordinates of the system. The authors do not know any other research in which this approach is used.In the paper the model of the system is presented. We discuss primary and secondary divisions of the electrodes. Further discretisation of the beams into the rigid finite elements and spring-damping elements and the connection between the beams and the electrodes are described. The theoretical part of the paper concludes with the formulation of the equations of motion of the system. Subsequently, the experimental validation of the method, essential for evaluating its correctness and usefulness, is presented.


## 2. Model of the Dust-Removal System

In the case of the gravity-operated rapping systems, excitation of vibrations is achieved by the axial striking of a beater on the anvil beam, in which an intensive stress-wave is generated.Both the electrodes’ geometrical features and the force impact have an essential influence on tangent and normal accelerations at different points of the electrodes and thus on the effectiveness of the dust-removal process. The dust-removal system, as shown in Figure1, consists of two beams which support a set of collecting electrodes.A collecting electrode is usually connected with the top beam, the suspension bar, by means of one articulated joint, while at the bottom it is connected with the brushing bar using two joints. The model of the system is divided into three subsystems: top and bottom beams and the collecting electrodes.


### 2.1. Model of the SIGMA Electrode

The collecting electrodes are 10- to 16-meter-long shells of complicated shapes (Figure2), which can be considered as a set of long strips with different angles of inclination towards axis y
^. In order to model the electrodes, the hybrid finite element method is proposed. The discretisation procedure is carried out in two steps, which are called primary and secondary division.Figure 2
Primary division of the electrode.


#### 2.1.1. Primary Division

Let us assign a coordinate systemx
^
y
^
z
^ as in Figure 2. Discretization is carried out along axis y
^ into n
y strips and along axis x
^ into n
x segments with the length
(1)
Δ
x
=
L
n
x
,
where L is the length of the electrode. Thus the electrode is divided into
(2)
n
=
n
x
n
y
elements.Let us consider a rectangular element as in Figure3. The following vectors describe displacements of nodes for such an element:
(3)
q
p
′
=
u
p
′
,
v
p
′
,
w
p
′
,
φ
p
′
x
,
φ
p
′
y
,
φ
p
′
z
T
for
p
=
1
,
2
,
3
,
4
,
where u
p
′, v
p
′, w
p
′ are displacements of node p along x
^
′, y
^
′, z
^
′, respectively; φ
p
′
x, φ
p
′
y, φ
p
′
z are respective rotations in node p about axes parallel to x
^
′, y
^
′, z
^
′.Figure 3
A rectangular shell element.Shape functions are defined as(4)
u
′
=
a
1
u
+
a
2
u
x
′
+
a
3
u
y
′
+
a
4
u
x
′
y
′
+
a
5
u
y
′
2
+
a
6
u
x
′
y
′
2
,
v
′
=
a
1
v
+
a
2
v
x
′
+
a
3
v
y
′
+
a
4
v
x
′
y
′
+
a
5
v
x
′
2
+
a
6
v
x
′
2
y
′
,
w
′
=
a
1
w
+
a
2
w
x
′
+
a
3
w
y
′
+
a
4
w
x
′
2
+
a
5
w
x
′
y
′
+
a
6
w
y
′
2
+
a
7
w
x
′
3
+
a
8
w
x
′
2
y
′
+
a
9
w
x
′
y
′
2
+
a
10
w
y
′
3
+
a
11
w
x
′
3
y
′
+
a
12
w
x
′
y
′
3
.It is assumed thatu
p
′, v
p
′, and φ
p
′
z describe shield deformations, while w
p
′, φ
p
′
x, and φ
p
′
ydescribe plate deformations. The angles are defined by the following relations:
(5)
φ
′
x
=
∂
w
′
∂
y
′
,
φ
′
y
=
∂
w
′
∂
x
′
,
φ
′
z
=
1
2
∂
v
′
∂
x
′
-
∂
u
′
∂
y
′
.The elastic strain energy is a sum of energies corresponding to the shield and plate states and can be calculated as [13]
(6)
E
=
1
2
∫
V
σ
(
s
)
T
ε
(
s
)
d
V
+
1
2
∫
V
σ
(
p
)
T
ε
(
p
)
d
V
,
where V
=
Δ
x
Δ
y
h, h is the thickness of the element,
(7)
ε
(
s
)
=
ε
x
s
ε
y
s
γ
x
y
s
=
∂
u
′
∂
x
′
∂
v
′
∂
y
′
∂
u
′
∂
y
′
+
∂
v
′
∂
x
′
,
ε
(
p
)
=
ε
x
p
ε
y
p
γ
x
y
p
=
-
z
∂
2
w
′
∂
x
′
2
∂
2
w
′
∂
y
′
2
2
∂
w
′
∂
x
′
∂
y
′
are vectors of strains,
(8)
σ
(
l
)
=
σ
x
l
σ
y
l
σ
x
y
l
=
D
ε
(
l
)
,
l
∈
s
,
p
are vectors of stresses,
(9)
D
=
E
1
-
υ
2
1
υ
0
υ
1
0
0
0
1
-
υ
2
,
E is Young’s modulus, and υ is the Poisson ratio.After the necessary calculations, the energy of deformation of the element can be expressed in the following form:(10)
E
=
1
2
s
-
T
C
(
s
)
s
-
+
p
-
T
C
(
p
)
p
-
,
where
(11)
s
¯
=
u
1
′
v
1
′
φ
1
′
z
u
2
′
v
2
′
φ
2
′
z
u
3
′
v
3
′
φ
3
′
z
u
4
′
v
4
′
φ
4
′
z
T
,
p
¯
=

w
1
′
φ
1
′
x
φ
1
′
y
w
2
′
φ
2
′
x
φ
2
′
y
w
3
′
φ
3
′
x
φ
2
′
y
w
4
′
φ
4
′
x
φ
2
′
y
T
.The dimension of matricesC
(
s
) and C
(
p
) is 12
×
12. Integrals occurring in elements of matrices can be calculated analytically or by means of Gauss quadratures. Detailed formulae for elements of matrices C
(
s
) and C
(
p
) are presented in [14]. If we assume that the form of these matrices is as follows:
(12)
C
(
l
)
=
C
11
(
l
)
C
12
(
l
)
C
13
(
l
)
C
14
(
l
)
C
21
(
l
)
C
22
(
l
)
C
23
(
l
)
C
24
(
l
)
C
31
(
l
)
C
32
(
l
)
C
33
(
l
)
C
34
(
l
)
C
41
(
l
)
C
42
(
l
)
C
43
(
l
)
C
44
(
l
)
,
l
∈
s
,
p
,
then having taken into account the definition of vectors of generalized coordinates describing plate and shield deformations in (10), the elastic strain energy can be presented in terms of the generalised coordinates as follows:
(13)
E
=
1
2
q
1
′
T
q
2
′
T
q
3
′
T
q
4
′
T
C
11
C
12
C
13
C
14
C
21
C
22
C
23
C
24
C
31
C
32
C
33
C
34
C
41
C
42
C
43
C
44
q
1
′
q
2
′
q
3
′
q
4
′
=
1
2
q
′
T
C
q
′
,
where C is a stiffness matrix with 24
×
24 elements,
(14)
q
′
=
q
1
′
q
2
′
q
3
′
q
4
′
,
and nonzero elements of matrixC from (13) are defined as follows:
(15)
C
i
,
k
q
,
t
=
C
i
,
k
(
s
)
q
,
s
for
q
,
t
=
1,2
C
i
,
k
q
,
6
=
C
i
,
k
(
s
)
q
,
3
for
q
=
1,2
C
i
,
k
q
+
2
,
t
+
2
=
C
i
,
k
(
p
)
q
,
t
for
q
,
t
=
1,2
,
3
C
i
,
k
6
,
q
=
C
i
,
k
(
s
)
3
,
q
for
q
=
1,2
C
i
,
k
6,6
=
C
i
,
k
(
s
)
3,3
.For the assumed values of parametersE
,
υ
,
α
,
h, the stiffness matrices depend only on Δ
x and Δ
y. Having used denotations from Figures 2 and 3 and relation (1), in the case considered, we have C
(
i
,
j
)
=
C
(
j
)
=
C
Δ
x
,
Δ
y
j for i
=
1
,
…
,
n
x,  j
=
1
,
…
,
n
y. Thus, as a result of constant length Δ
x of elements, stiffness matrices for each strip are identical and are calculated only for j
=
1
,
…
,
n
y.The coordinatesx
i
,
j
C
,
y
i
,
j
C
,
z
i
,
j
C of the geometrical center of element (i, j) can be calculated according to the formulas
(16)
x
i
,
j
C
=
i
-
1
Δ
x
+
Δ
x
2
,
y
i
,
j
C
=
y
0
+
∑
s
=
0
j
-
1
b
s
cos
⁡
⁡
α
s
+
b
j
2
cos
⁡
⁡
α
j
,
z
i
,
j
C
=
z
0
+
∑
s
=
0
j
-
1
b
s
sin
⁡
α
s
+
b
j
2
sin
⁡
α
j
,
where y
0, z
0 define coordinates of the left edge of strip number 1 and b
j, α
j are the width and angle of inclination of the jth strip to axis y
^, respectively.Coordinates of nodes of primary elements can be calculated as follows:(17)
r
i
,
j
,
p
=
r
i
,
j
C
+
R
j
r
i
,
j
,
p
′
,
where r
i
,
j
C
=
x
i
,
j
C
y
i
,
j
C
z
i
,
j
C
T is defined in (16),
(18)
R
j
=
1
0
0
0
cos
⁡
α
j
-
sin
α
j
0
sin
α
j
cos
⁡
α
j
,
and r
i
,
j
,
p
′ is the vector of coordinates of the mass centre of the pth segment in the local coordinate system:
(19)
r
i
,
j
,
1
′
=
-
Δ
x
2
-
b
j
2
0
T
,
r
i
,
j
,
2
′
=
Δ
x
2
-
b
j
2
0
T
,
r
i
,
j
,
3
′
=
Δ
x
2
b
j
2
0
T
,
r
i
,
j
,
4
′
=
-
Δ
x
2
b
j
2
0
T
.


#### 2.1.2. Secondary Division of the SIGMA Electrode

In the rigid finite element method the flexible body considered is divided into rigid elements reflecting inertial features of the body and spring-damping elements which connect the rigid elements and adopt spring and damping features [9]. For the subsequent analysis damping will be omitted. Since in the hybrid finite element method proposed, spring features of the electrode are included in the elements described above, the method of discretisation into rigid elements is discussed below.Let us consider element (i, j) from the primary division, as in Figure 4.Figure 4
Secondary division of the shell.Coordinate systemx
^, y
^, z
^ is a global coordinate system. Primary element (i, j) is divided into four equal parts called subsegments. Then the rigid finite element (rfe) is formed from one, two, or four subdivisions, which may belong to different primary elements. Thus each segment of the primary element belongs to the following rfes:
(20)
i
,
j
1
=
i
-
1
m
+
j
,
i
,
j
2
=
i
m
+
j
,
i
,
j
3
=
i
m
+
j
+
1
,
i
,
j
4
=
i
-
1
m
+
j
+
1
,
for i
=
1
,
…
,
n
x, j
=
1
,
…
,
n
y, and m
=
n
y
+
1.Figure5 shows an example of primary and secondary division of a rectangular area for n
x
=
3 and n
y
=
4.Figure 5
Division of a rectangular area inton
x
=
3 and n
y
=
4 elements.The number of rigid finite elements is(21)
n
=
n
x
+
1
m
=
n
x
+
1
n
y
+
1
.It can be seen from Figures4 and 5 that rfes are constructed by rigid connections of stiffened subsegments in the following way.(i)
Corner elements:
rfe 1 is subsegment1,1
1 stiffened;
rfem is subsegment 1
,
n
y
4 stiffened;
rfen
x
m is subsegment 1
,
n
x
2 stiffened;
rfen is subsegment n
x
,
n
y
3 stiffened.
(ii)
Other boundary elements:
forj
=
2
,
…
,
n
y,
rfej is formed by connecting and stiffening subsegments 1
,
j
-
1
4 and 1
,
j
1;
rfen
x
m
+
j is formed by connecting and stiffening subsegments n
x
,
j
-
1
3 and n
x
,
j
2;
fori
=
1
,
…
,
n
x
-
1,
rfei
m
+
1 is formed by connecting and stiffening subsegments i
,
1
2 and i
+
1,1
1;
rfe(
i
+
1
)
m is formed by connecting and stiffening subsegments i
,
n
x
3 and i
+
1
,
n
x
4.
(iii)
Internal elements:
fori
=
1
,
…
,
n
x
-
1, j
=
1
,
…
,
n
y
-
1,
rfei
m
+
j
+
1 is formed by connecting and stiffening subsegments i
,
j
3,i
+
1
,
j
4, i
+
1
,
j
+
1
1, and i
,
j
+
1
2.Secondary division of the electrode into rfes is essential for the method discussed since in further considerations the shell is treated as a system ofn rigid bodies with 6 degrees of freedom. A local coordinate system x
^
k
C
y
^
k
C
z
^
k
C, whose axes are central inertial axes of the element, is assigned to each rfe as in Figure 6. Angle β
k defines inclination of the plane of element k with respect to plane x
^
y
^ of coordinate system x
^
y
^
z
^ assigned to the electrode.Figure 6
Rigid finite element.The vector of the generalized coordinates of rfek includes three displacements and three rotations:
(22)
q
k
′
=
Δ
k
′
φ
k
′
=
Δ
x
k
′
,
Δ
y
k
′
,
Δ
z
k
′
,
φ
k
′
x
,
φ
k
′
y
,
φ
k
′
z
T
,
where Δ
k
′
=
Δ
x
k
′
,
Δ
y
k
′
,
Δ
z
′
T and φ
k
′
=
φ
k
′
x
,
φ
k
′
y
,
φ
k
′
z
T.When the vibrations of the collecting electrodes are considered, it can be assumed that anglesφ
k
′
x, φ
k
′
y, φ
k
′
z are small. Thus, when coordinates (
x
′, y
′, z
′
) of a point are known in the local system of the kth rfe, its coordinates in the global coordinate system can be calculated as follows:
(23)
r
=
r
k
rfe
+
R
-
k
r
k
′
C
,
r
k
′
C
=
U
r
′
q
k
′
+
r
′
,
where r
′
=
x
′
,
y
′
,
z
′
T is the coordinate vector in x
^
k
C
,
y
^
k
C
,
z
^
k
C and r
k
rfe is defined in Figure 6, as the coordinate vector of point C
k for a nondeformed rfe (when q
k
′
=
0),
(24)
U
=
1
0
0
0
z
′
-
y
′
0
1
0
-
z
′
0
x
′
0
0
1
y
′
-
x
′
0
,
R
-
k
=
1
0
0
0
cos
⁡
β
k
-
sin
β
k
0
sin
β
k
cos
⁡
β
k
.


#### 2.1.3. Equations of Motion of the SIGMA Electrode

Using the standard RFEM approach, the diagonal mass matrix of thekth rfe can be formulated as follows:
(25)
M
k
=
diag
⁡
m
k
,
m
k
,
m
k
,
I
k
x
,
I
k
y
,
I
k
z
.In the case of SIGMA electrodes it is relatively easy to formulate the algorithm to calculate elements of matrixM
k using the above formula and angles β
k defining the position of central inertial axes of rfes.In view of (20) and (22) the kinetic energy of an rfe can be presented in the form
(26)
T
k
=
1
2
q
˙
k
′
T
M
k
q
˙
k
′
and the kinetic energy of the whole electrode as
(27)
T
=
∑
k
=
1
n
T
k
=
1
2
q
˙
′
T
M
q
˙
′
,
where q
′
=
q
1
′
⋯
q
k
′
⋯
q
n
′
T and M
=
diag
⁡
M
1
,
M
2
,
…
,
M
n.In order to formulate the stiffness matrix for the whole electrode, the displacements of primary elements have to be expressed in terms of generalized coordinates of the rigid finite elements. The vector of generalized coordinates and the stiffness matrix for primary element(
i
,
j
) are defined in (12). Individual segments of primary element (
i
,
j
) belong to various rfes, whose indexes are determined in (20).It results from (20) that node i
,
j
p belongs to rfe:
(28)
k
i
,
j
,
p
=
i
-
1
m
+
j
,
for
p
=
1
,
i
m
+
j
,
for
p
=
2
,
i
m
+
j
+
1
,
for
p
=
3
,
i
-
1
m
+
j
+
1
,
for
p
=
4
.Coordinatesr
i
,
j
,
p in the global coordinate system of the electrode x
^
,
y
^
,
z
^  are defined by (17). Having taken into account (23), node coordinates in local system x
^
k
′
C
,
y
^
k
′
C
,
z
^
k
′
C can be defined as follows:
(29)
r
k
,
p
′
C
=
R
-
k
T
r
i
,
j
,
p
-
r
k
rfe
,
where R
-
k
T
=
R
-
k
-
1, k
=
k
(
i
,
j
,
p
) defined in (28). If rfe k
(
i
,
j
,
p
) has coordinates defined in (22), then coordinates of node i
,
j
p in the global system can be calculated from (23):
(30)
r
k
(
i
,
j
,
p
)
=
r
k
rfe
+
R
-
k
U
r
k
,
p
′
C
q
k
′
+
r
k
,
p
′
C
.In order to define the energy of spring deformation of nodei
,
j
p the node displacements and rotations
(31)
Δ
i
,
j
,
p
′
=
u
i
,
j
,
p
′
,
v
i
,
j
,
p
′
,
w
i
,
j
,
p
′
T
,
φ
i
,
j
,
p
′
=
φ
i
,
j
,
p
′
x
,
φ
i
,
j
,
p
′
y
,
φ
i
,
j
,
p
′
z
T
have to be defined by means of generalized coordinates q
k
′.The displacement vector of nodei
,
j
p belonging to rfe k, as in (28), is defined as follows:
(32)
Δ
i
,
j
,
p
′
=
R
j
T
r
k
i
,
j
,
p
-
r
i
,
j
,
p
.After some transformations, when (30), (29), and (17) are taken into account, the following is obtained:
(33)
Δ
i
,
j
,
p
′
=
R
j
T
R
-
k
U
r
k
,
p
′
C
.Rotation angles can be calculated from the formula(34)
φ
i
,
j
,
p
′
=
R
j
T
R
-
k
S
q
k
′
,
where
(35)
S
=
0
0
0
1
0
0
0
0
0
0
1
0
0
0
0
0
0
1
.
Finally, the following relation is derived:
(36)
q
i
,
j
,
p
′
=
Δ
i
,
j
,
p
′
φ
i
,
j
,
p
′
=
R
j
T
R
-
k
U
r
k
,
p
′
C
R
j
T
R
-
k
S
Δ
k
′
φ
k
′
=
H
k
,
p
q
k
′
.MatrixH
k
,
p enables the coordinates of an rfe to be transformed into the elastic coordinates of a primary element. Thus, the energy of spring deformation of primary element (
i
,
j
) can be written in the following form
(37)
E
i
,
j
=
1
2
∑
p
1
=
1
4
∑
p
2
=
1
4
q
i
,
j
,
p
1
′
C
T
p
1
,
p
2
(
j
)
q
i
,
j
,
p
2
′
,
where C
α
,
β
(
j
) is defined in (12) and (13).Bearing in mind thatk
=
k
(
i
,
j
,
p
), and taking into account (36), the above can be written as follows:
(38)
E
i
,
j
=
1
2
∑
p
1
=
1
4
∑
p
2
=
1
4
q
k
(
i
,
j
,
p
1
)
′
T
H
k
i
,
j
,
p
1
,
p
1
T
C
p
1
,
p
2
(
j
)
H
k
i
,
j
,
p
2
,
p
2
q
k
(
i
,
j
,
p
2
)
′
=
1
2
∑
p
1
=
1
4
∑
p
2
=
1
4
q
k
(
i
,
j
,
p
1
)
′
T
C
p
1
,
p
2
(
i
,
j
)
q
k
(
i
,
j
,
p
2
)
′
,
where C
p
1
,
p
2
(
i
,
j
)
=
H
k
(
i
,
j
,
p
1
)
,
p
1
T
C
p
1
,
p
2
(
j
)
H
k
(
i
,
j
,
p
2
)
,
p
2.This leads to a straightforward method of generating stiffness matrixC using coordinates q
k
′ as generalized coordinates. The appropriate algorithm is as follows:(1)
fork
1
=
1
,
…
,
6
n
;
k
2
=
1
,
…
,
6
n, accept: c
k
1
,
k
2
=
0.
(2)
fori
=
1
,
…
,
n
x; j
=
1
,
…
,
n
y
(i)
forp
=
1
,
…
,
4 compute H
k
(
i
,
j
,
p
)
,
p according to (36)
(ii)
forp
1
=
1
,
…
,
4; p
2
=
1
,
…
,
4
(a)
compute:(39)
k
1
=
k
i
,
j
,
p
1
,
k
2
=
i
,
j
,
p
2
k
1
C
=
6
k
1
-
1
,
k
2
C
=
6
k
2
-
1
(b)
forj
1
=
1
,
…
,
6; j
2
=
1
,
…
,
6 compute:
(40)
c
k
1
C
+
j
1
,
k
2
C
+
j
2
=
c
k
1
C
+
j
1
,
k
2
C
+
j
2
+
c
p
1
,
p
2
(
i
,
j
)
.It is important to note that matrixC is a sparse band matrix and this is used to facilitate computations of dynamics of the collecting electrodes. The potential energy of gravity forces is calculated as follows:
(41)
V
g
=
-
∑
k
=
1
n
x
m
k
g
x
k
C
,
where g is the acceleration of gravity, x
k
C
=
θ
1
r
k
C, and θ
1
=
1
0
0. Bearing in mind (23) we can calculate derivatives of the potential energy
(42)
∂
V
g
∂
q
k
′
=
-
m
k
g
θ
1
R
-
k
U
k
C
.Since the origin of the local coordinate system of rfek is located at the center of mass of the rfe, that is, U
k
C
=
I
0, the derivative of potential energy of gravity forces for the jth strip of the electrode can be presented in the matrix form
(43)
∂
V
g
∂
q
′
=
G
=
G
1
⋯
G
k
⋯
G
n
x
T
,
where G
k
=
-
g
m
k
0
0
0
0
0
T.Equations of motion of a free electrode can be written in the form(44)
M
q
¨
′
+
C
q
′
=
-
∂
V
g
∂
q
′
,
and the eigenvalue problem takes the form
(45)
det
⁡
(
C
-
ω
2
M
)
=
0
.Some important differences between the hybrid finite element method and the finite element method applied to modeling of electrodes are listed in Table1.Table 1
Comparison of features of the finite element method and hybrid finite element method.
Feature
Finite element method
Hybrid element method
Generalized coordinates
Displacements of nodes of elementsq
′
i
,
j
(
p
) (elastic coordinates) (36)
Displacements and rotations of rigid elementsq
′
k (rigid coordinates) (22)
Stiffness matrixC
Interpolation of field displacements and rotations by means of continuous functions
Interpolation of field displacements and rotations by means of continuous functions
Additional transformation of elastic coordinates to rigid coordinates:
q
′
i
,
j
(
p
)
=
H
k
,
p
q
k
′
Mass matrixM
Interpolation of field displacements and rotation velocities by means of continuous functions
M
=
diag
⁡
{
M
1
,
M
2
,
…
,
M
n
}
M
k
=
diag
⁡
{
m
k
,
m
k
,
m
k
,
I
k
x
,
I
k
y
,
I
k
z
}    Inertial moments are taken into account
Additional operations in order to use lumped mass


### 2.2. Models of Beams

The detailed description of the rigid finite element method is presented in [8, 9]. In this paper only the general idea and final equations of motion will be given. More details can be found in [15]. The general idea of discretisation of the top and bottom beams is presented in Figure 7.Figure 7
Connections between the beams and thejth plate and discretisation of beams.The number of elements into which the beams are discretized may be different for both beams. Ifp is the number of electrodes and a is the length of a section of the beam with respect to a single electrode (all electrodes in a rapping system are usually the same) (Figure 2), then in the primary division the length of rigid elements is as follows:
(46)
d
0
s
=
a
L
s
+
a
-
s
2
,
d
i
s
=
a
-
s
,
d
p
s
s
=
a
-
s
2
+
a
R
s
,
mmmmmmlmmmh
for
i
=
1
,
…
,
p
s
-
1
,
where a
-
s
=
a
/
n
-
s, s
∈
t
,
b; t, b denote top or bottom beam, respectively; and p
s
=
n
-
s
p.The secondary division leads to a system of rigid finite elements with the length equal to(47)
l
0
s
=
d
0
s
2
,
l
i
s
=
d
i
-
1
s
+
d
i
s
2
,
l
p
s
+
1
s
=
d
p
s
s
2
,
mmmmmmmmmmml
for
i
=
1
,
…
,
p
s
.We treat both beams as those with a uniform cross-section, whose geometrical parameters exactly reflect mass and stiffness parameters of the real beams.The axes of the local coordinate system assigned to each rfe before deformation are parallel to the axes of the global reference system and are the main central inertial axes of the element. Thus the mass matrix of the element is diagonal and takes the following form:(48)
M
i
s
=
diag
⁡
m
i
s
m
i
s
m
i
s
I
i
,
x
s
I
i
,
y
s
I
i
,
z
s
,
where m
i
s
=
ρ
l
i
s
A
s is the mass of the ith rfe, A
s is the cross-section area, I
x
s
,
I
z
s are inertial moments of the cross-section of the beam, I
i
,
x
s
=
m
i
s
l
i
s
2
/
12
+
I
x
s
/
A
s,  I
i
,
y
s
=
m
i
s
/
A
s
I
x
s
+
I
z
s, and I
i
,
z
s
=
m
i
s
l
i
s
2
/
12
+
I
z
s
/
A
s. The vector of generalized coordinates of rfe i is as follows:
(49)
q
i
s
=
x
~
i
s
y
~
i
s
z
~
i
s
φ
~
i
s
θ
~
i
s
ψ
~
i
s
T
,
where x
~
i
s, y
~
i
s, z
~
i
s are displacements of the center of mass S
i
s of the rfe and φ
~
i
s, θ
~
i
s, ψ
~
i
s are the angles of rotations. The kinetic energy of the ith rfe is calculated according to the formula
(50)
E
i
s
=
1
2
q
˙
i
s
T
M
i
s
q
˙
i
s
.The energy of spring deformation can be presented in the form(51)
V
e
,
i
s
=
1
2
Δ
i
s
T
C
i
s
Δ
i
s
,
where C
i
s
=
diag
⁡
c
i
,
1
s
⋯
c
i
,
6
s, c
i
,
j
s are the stiffness coefficients calculated according to the formulas given in [9], and Δ
i
s is the deformation of spring damping element (sde) i, which depends on q
i
s and q
i
-
1
s.The potential energy of gravity forces of each rfe can be calculated according to the formula(52)
V
s
,
i
g
=
-
m
i
g
g
x
S
i
g
,
where x
S
i
g
=
θ
1
r
S
i
g and θ
1 is defined in (41).Kinetic and potential energies of the beams are the sum of energies of all rfes:(53)
E
s
=
∑
i
=
0
p
s
+
1
E
i
s
,
V
g
s
=
∑
i
=
0
p
s
+
1
V
g
,
i
s
,
and the potential energy of spring deformation is the sum of energies of sdes 0 to p
s:
(54)
V
e
s
=
∑
i
=
0
p
s
V
e
,
i
s
.The vector of generalized coordinates of a beam can be presented as(55)
q
s
=
q
0
s
⋯
q
i
s
⋯
q
p
s
+
1
s
T
.So finally, the following can be written:(56)
d
d
t
∂
E
s
∂
q
˙
s
-
∂
E
s
∂
q
s
+
∂
V
g
s
+
V
e
s
∂
q
s
=
M
s
q
¨
s
+
K
s
q
s
+
G
s
,
where M
s is the diagonal mass matrix, K
s is the stiffness matrix with constant coefficients, and G
s is the vector of gravity forces, s
∈
t
,
b. Thus, the motion of beams is described by the following number of generalized coordinates:
(57)
n
dof
s
=
6
p
s
+
2
.The top beam is connected with the base, and this connection has to be taken into account when deriving the equations of motion.As has been mentioned, the vibrations in the rapping system are excited, when the beater hits the anvil. It is assumed that the impulse of the force is modelled as a central, straight stroke of forceF
=
0
-
S
(
t
)
0
T, and the resulting generalized force follows
(58)
Q
p
b
+
1
(
S
)
=
0
-
S
0
0
0
0
T
.


### 2.3. Modeling Elastic Connections

The collecting electrodes are connected with the suspension and brushing bars and they may be also connected to each other. Such connections can be easily included in the model by means of an elastic connection of two rigid bodies. Let us consider two rigid bodies connected at pointE by means of spring-damping element (sde) as shown in Figure 8.Figure 8
Elastic connectionE between rfes l and r.We assume the following denotations:s
- denotes the system before deformation and s
~ denotes the system rigidly assigned to rfe s, where s
∈
r
,
l.Owing to (22) and (49), the motion of rfes is governed by the vectors
(59)
q
~
s
=
Δ
~
s
φ
~
s
,
s
∈
l
,
r
,
where Δ
~
s
=
x
~
s
y
~
s
z
~
s
T, φ
~
s
=
φ
~
s
θ
~
s
ψ
~
s
T, x
~
s, y
~
s, z
~
s describe the position of the beginning of system {
s
~
} in {
s
-
}, and φ
~
s, θ
~
s, ψ
~
s describe the orientation of system {
s
~
} with respect to {
s
-
}. After deformation, rfes l and r move, and thus the position of point E in the global coordinate system may vary depending on whether point E is considered as a point of rfe l or r:
(60)
r
E
l
=
r
c
l
+
R
l
Δ
~
l
+
R
~
l
r
~
E
l
,
r
E
r
=
r
c
r
+
R
r
Δ
~
r
+
R
~
r
r
~
E
r
,
where r
c
(
s
) is the position vector of the beginning of {
s
-
} in the global system {
}, r
~
E
(
s
) is the vector of coordinates of point E in system {
s
~
}, R
s is the rotation matrix of direction cosines of system {
s
-
} with respect to {
},
(61)
R
~
(
s
)
=
1
-
ψ
~
(
s
)
θ
~
(
s
)
ψ
~
(
s
)
1
-
φ
~
(
s
)
-
θ
~
(
s
)
φ
~
i
(
s
)
1
when the rotation angles are small, and r
~
(
s
) are coordinates of a point in system {
s
~
} rigidly connected with rfe s.Vectors of translational and rotational displacements of sdeE with respect to the global system can be calculated as follows:
(62)
Δ
r
E
=
r
E
(
r
)
-
r
E
(
l
)
=
R
r
x
~
(
r
)
+
S
~
E
(
r
)
φ
~
(
r
)
-
R
l
x
~
(
l
)
+
S
~
E
(
l
)
φ
~
(
l
)
,
Δ
φ
E
=
R
r
φ
~
(
r
)
-
R
l
φ
~
(
l
)
,
where
(63)
S
~
E
(
s
)
=
0
z
~
E
(
s
)
-
y
~
E
(
s
)
-
z
~
E
(
s
)
0
x
~
E
(
s
)
y
~
E
(
s
)
-
x
~
E
(
s
)
0
,
s
∈
r
,
l
.After transformations the following can be obtained:(64)
Δ
~
E
=
Δ
r
~
E
Δ
φ
~
E
=
U
E
(
r
)
q
~
r
-
U
E
l
q
~
l
,
where
(65)
U
E
(
r
)
=
R
E
T
R
r
R
E
T
R
r
S
~
E
(
r
)
0
R
E
T
R
r
,
U
E
(
l
)
=
R
E
T
R
l
R
E
T
R
l
S
~
E
(
l
)
0
R
E
T
R
l
.
It is important to note that the matrices S
~
E
(
l
), S
~
E
(
r
), R
l, and R
r are constant.Finally, the energy of spring deformation of the connection can be presented in the following form:(66)
V
E
=
1
2
Δ
~
E
T
C
E
Δ
~
E
,
where C
E
=
diag
⁡
c
E
x
,
c
E
y
,
c
E
z
,
c
E
φ
,
c
E
θ
,
c
E
ψ. It can be seen from (64) that the energy of spring deformation of the connection depends on generalized coordinates of rfes l and r. Thus the following components will appear in the equations of motion
(67)
∂
V
E
∂
q
~
(
l
)
=
C
l
l
(
E
)
q
~
(
l
)
+
C
l
r
(
E
)
q
~
(
r
)
,
∂
V
E
∂
q
~
(
r
)
=
C
r
l
(
E
)
q
~
(
l
)
+
C
r
r
(
E
)
q
~
(
r
)
,
where
(68)
C
l
l
(
E
)
=
U
E
(
l
)
T
C
E
U
E
(
l
)
,
C
r
r
(
E
)
=
U
E
(
t
)
T
C
E
U
E
(
t
)
,
C
l
r
(
E
)
=
C
r
l
(
E
)
T
=
-
U
E
(
l
)
T
C
E
U
E
(
r
)
.It should be stressed that the element presented is universal. When stiffness coefficients in matrixC
E are assumed to be large enough or equal to zero, any joint (rotary, longitudinal, and spherical) can be modeled.In order to model connections between beams and plates, two joints have been used (as shown in Figure9) and the coefficients of stiffness have been assumed to be c
T
=
1
0
12N/m for translational displacements and as c
R
=
1
0
9 N/rad in the case of rotary displacements.Models of connections between beams and a plate: (a) rotary joint and (b) universal joint.
(a)
(b)


### 2.4. Equations of Motion

The equations of motion for the system of the beams and electrodes can be written in the following form:(69)
M
t
q
¨
t
+
K
t
q
t
=
-
G
t
,
M
j
q
¨
j
+
K
j
q
j
=
-
G
j
,
j
=
1
,
…
,
n
p
,
M
b
q
¨
b
+
K
b
q
b
=
-
G
b
+
f
b
(
S
)
,
where n
p is the number of active electrodes, M
t, K
t, G
t, M
b, K
b, G
b are defined in Section 2.2,   M
j, K
j, G
j are defined in Section 2.1,   j denotes the number of the electrode, f
b
=
0
,
…
,
Q
p
+
1
S
T
T, and Q
p
+
1
(
S
) is defined in (58).It has been assumed that the electrodes are connected, as in Figure10, with the top beam by means of one elastic connection (rotary joint, Figure 9(a)) and with the bottom beam by two elastic connections (all springs are active, Figure 9(b)).Figure 10
Connection of the electrode with the beams.Finally, the equations of motion after necessary calculations take the following form:(70)
M
t
q
¨
t
+
C
11
q
t
+
∑
j
=
1
p
a
C
j
(
1
)
q
j
=
-
G
t
,
M
j
q
¨
j
+
C
j
(
1
)
T
q
t
+
C
(
j
)
q
j
+
C
j
3
T
q
b
=
-
G
j
,
j
=
1
,
…
,
n
p
,
M
b
q
¨
b
+
∑
j
=
1
p
a
C
j
(
3
)
q
j
+
C
33
q
b
=
-
G
b
+
f
b
.The coupling resulting from the connections between the electrodes is reflected by new stiffness matricesC
j
(
1
)and C
j
(
3
), while C
11, C
(
j
), and C
33 are modified matrices K
t, K
j, and K
b supplemented with the connections of the beams with the electrodes, respectively.


## 2.1. Model of the SIGMA Electrode

The collecting electrodes are 10- to 16-meter-long shells of complicated shapes (Figure2), which can be considered as a set of long strips with different angles of inclination towards axis y
^. In order to model the electrodes, the hybrid finite element method is proposed. The discretisation procedure is carried out in two steps, which are called primary and secondary division.Figure 2
Primary division of the electrode.


### 2.1.1. Primary Division

Let us assign a coordinate systemx
^
y
^
z
^ as in Figure 2. Discretization is carried out along axis y
^ into n
y strips and along axis x
^ into n
x segments with the length
(1)
Δ
x
=
L
n
x
,
where L is the length of the electrode. Thus the electrode is divided into
(2)
n
=
n
x
n
y
elements.Let us consider a rectangular element as in Figure3. The following vectors describe displacements of nodes for such an element:
(3)
q
p
′
=
u
p
′
,
v
p
′
,
w
p
′
,
φ
p
′
x
,
φ
p
′
y
,
φ
p
′
z
T
for
p
=
1
,
2
,
3
,
4
,
where u
p
′, v
p
′, w
p
′ are displacements of node p along x
^
′, y
^
′, z
^
′, respectively; φ
p
′
x, φ
p
′
y, φ
p
′
z are respective rotations in node p about axes parallel to x
^
′, y
^
′, z
^
′.Figure 3
A rectangular shell element.Shape functions are defined as(4)
u
′
=
a
1
u
+
a
2
u
x
′
+
a
3
u
y
′
+
a
4
u
x
′
y
′
+
a
5
u
y
′
2
+
a
6
u
x
′
y
′
2
,
v
′
=
a
1
v
+
a
2
v
x
′
+
a
3
v
y
′
+
a
4
v
x
′
y
′
+
a
5
v
x
′
2
+
a
6
v
x
′
2
y
′
,
w
′
=
a
1
w
+
a
2
w
x
′
+
a
3
w
y
′
+
a
4
w
x
′
2
+
a
5
w
x
′
y
′
+
a
6
w
y
′
2
+
a
7
w
x
′
3
+
a
8
w
x
′
2
y
′
+
a
9
w
x
′
y
′
2
+
a
10
w
y
′
3
+
a
11
w
x
′
3
y
′
+
a
12
w
x
′
y
′
3
.It is assumed thatu
p
′, v
p
′, and φ
p
′
z describe shield deformations, while w
p
′, φ
p
′
x, and φ
p
′
ydescribe plate deformations. The angles are defined by the following relations:
(5)
φ
′
x
=
∂
w
′
∂
y
′
,
φ
′
y
=
∂
w
′
∂
x
′
,
φ
′
z
=
1
2
∂
v
′
∂
x
′
-
∂
u
′
∂
y
′
.The elastic strain energy is a sum of energies corresponding to the shield and plate states and can be calculated as [13]
(6)
E
=
1
2
∫
V
σ
(
s
)
T
ε
(
s
)
d
V
+
1
2
∫
V
σ
(
p
)
T
ε
(
p
)
d
V
,
where V
=
Δ
x
Δ
y
h, h is the thickness of the element,
(7)
ε
(
s
)
=
ε
x
s
ε
y
s
γ
x
y
s
=
∂
u
′
∂
x
′
∂
v
′
∂
y
′
∂
u
′
∂
y
′
+
∂
v
′
∂
x
′
,
ε
(
p
)
=
ε
x
p
ε
y
p
γ
x
y
p
=
-
z
∂
2
w
′
∂
x
′
2
∂
2
w
′
∂
y
′
2
2
∂
w
′
∂
x
′
∂
y
′
are vectors of strains,
(8)
σ
(
l
)
=
σ
x
l
σ
y
l
σ
x
y
l
=
D
ε
(
l
)
,
l
∈
s
,
p
are vectors of stresses,
(9)
D
=
E
1
-
υ
2
1
υ
0
υ
1
0
0
0
1
-
υ
2
,
E is Young’s modulus, and υ is the Poisson ratio.After the necessary calculations, the energy of deformation of the element can be expressed in the following form:(10)
E
=
1
2
s
-
T
C
(
s
)
s
-
+
p
-
T
C
(
p
)
p
-
,
where
(11)
s
¯
=
u
1
′
v
1
′
φ
1
′
z
u
2
′
v
2
′
φ
2
′
z
u
3
′
v
3
′
φ
3
′
z
u
4
′
v
4
′
φ
4
′
z
T
,
p
¯
=

w
1
′
φ
1
′
x
φ
1
′
y
w
2
′
φ
2
′
x
φ
2
′
y
w
3
′
φ
3
′
x
φ
2
′
y
w
4
′
φ
4
′
x
φ
2
′
y
T
.The dimension of matricesC
(
s
) and C
(
p
) is 12
×
12. Integrals occurring in elements of matrices can be calculated analytically or by means of Gauss quadratures. Detailed formulae for elements of matrices C
(
s
) and C
(
p
) are presented in [14]. If we assume that the form of these matrices is as follows:
(12)
C
(
l
)
=
C
11
(
l
)
C
12
(
l
)
C
13
(
l
)
C
14
(
l
)
C
21
(
l
)
C
22
(
l
)
C
23
(
l
)
C
24
(
l
)
C
31
(
l
)
C
32
(
l
)
C
33
(
l
)
C
34
(
l
)
C
41
(
l
)
C
42
(
l
)
C
43
(
l
)
C
44
(
l
)
,
l
∈
s
,
p
,
then having taken into account the definition of vectors of generalized coordinates describing plate and shield deformations in (10), the elastic strain energy can be presented in terms of the generalised coordinates as follows:
(13)
E
=
1
2
q
1
′
T
q
2
′
T
q
3
′
T
q
4
′
T
C
11
C
12
C
13
C
14
C
21
C
22
C
23
C
24
C
31
C
32
C
33
C
34
C
41
C
42
C
43
C
44
q
1
′
q
2
′
q
3
′
q
4
′
=
1
2
q
′
T
C
q
′
,
where C is a stiffness matrix with 24
×
24 elements,
(14)
q
′
=
q
1
′
q
2
′
q
3
′
q
4
′
,
and nonzero elements of matrixC from (13) are defined as follows:
(15)
C
i
,
k
q
,
t
=
C
i
,
k
(
s
)
q
,
s
for
q
,
t
=
1,2
C
i
,
k
q
,
6
=
C
i
,
k
(
s
)
q
,
3
for
q
=
1,2
C
i
,
k
q
+
2
,
t
+
2
=
C
i
,
k
(
p
)
q
,
t
for
q
,
t
=
1,2
,
3
C
i
,
k
6
,
q
=
C
i
,
k
(
s
)
3
,
q
for
q
=
1,2
C
i
,
k
6,6
=
C
i
,
k
(
s
)
3,3
.For the assumed values of parametersE
,
υ
,
α
,
h, the stiffness matrices depend only on Δ
x and Δ
y. Having used denotations from Figures 2 and 3 and relation (1), in the case considered, we have C
(
i
,
j
)
=
C
(
j
)
=
C
Δ
x
,
Δ
y
j for i
=
1
,
…
,
n
x,  j
=
1
,
…
,
n
y. Thus, as a result of constant length Δ
x of elements, stiffness matrices for each strip are identical and are calculated only for j
=
1
,
…
,
n
y.The coordinatesx
i
,
j
C
,
y
i
,
j
C
,
z
i
,
j
C of the geometrical center of element (i, j) can be calculated according to the formulas
(16)
x
i
,
j
C
=
i
-
1
Δ
x
+
Δ
x
2
,
y
i
,
j
C
=
y
0
+
∑
s
=
0
j
-
1
b
s
cos
⁡
⁡
α
s
+
b
j
2
cos
⁡
⁡
α
j
,
z
i
,
j
C
=
z
0
+
∑
s
=
0
j
-
1
b
s
sin
⁡
α
s
+
b
j
2
sin
⁡
α
j
,
where y
0, z
0 define coordinates of the left edge of strip number 1 and b
j, α
j are the width and angle of inclination of the jth strip to axis y
^, respectively.Coordinates of nodes of primary elements can be calculated as follows:(17)
r
i
,
j
,
p
=
r
i
,
j
C
+
R
j
r
i
,
j
,
p
′
,
where r
i
,
j
C
=
x
i
,
j
C
y
i
,
j
C
z
i
,
j
C
T is defined in (16),
(18)
R
j
=
1
0
0
0
cos
⁡
α
j
-
sin
α
j
0
sin
α
j
cos
⁡
α
j
,
and r
i
,
j
,
p
′ is the vector of coordinates of the mass centre of the pth segment in the local coordinate system:
(19)
r
i
,
j
,
1
′
=
-
Δ
x
2
-
b
j
2
0
T
,
r
i
,
j
,
2
′
=
Δ
x
2
-
b
j
2
0
T
,
r
i
,
j
,
3
′
=
Δ
x
2
b
j
2
0
T
,
r
i
,
j
,
4
′
=
-
Δ
x
2
b
j
2
0
T
.


### 2.1.2. Secondary Division of the SIGMA Electrode

In the rigid finite element method the flexible body considered is divided into rigid elements reflecting inertial features of the body and spring-damping elements which connect the rigid elements and adopt spring and damping features [9]. For the subsequent analysis damping will be omitted. Since in the hybrid finite element method proposed, spring features of the electrode are included in the elements described above, the method of discretisation into rigid elements is discussed below.Let us consider element (i, j) from the primary division, as in Figure 4.Figure 4
Secondary division of the shell.Coordinate systemx
^, y
^, z
^ is a global coordinate system. Primary element (i, j) is divided into four equal parts called subsegments. Then the rigid finite element (rfe) is formed from one, two, or four subdivisions, which may belong to different primary elements. Thus each segment of the primary element belongs to the following rfes:
(20)
i
,
j
1
=
i
-
1
m
+
j
,
i
,
j
2
=
i
m
+
j
,
i
,
j
3
=
i
m
+
j
+
1
,
i
,
j
4
=
i
-
1
m
+
j
+
1
,
for i
=
1
,
…
,
n
x, j
=
1
,
…
,
n
y, and m
=
n
y
+
1.Figure5 shows an example of primary and secondary division of a rectangular area for n
x
=
3 and n
y
=
4.Figure 5
Division of a rectangular area inton
x
=
3 and n
y
=
4 elements.The number of rigid finite elements is(21)
n
=
n
x
+
1
m
=
n
x
+
1
n
y
+
1
.It can be seen from Figures4 and 5 that rfes are constructed by rigid connections of stiffened subsegments in the following way.(i)
Corner elements:
rfe 1 is subsegment1,1
1 stiffened;
rfem is subsegment 1
,
n
y
4 stiffened;
rfen
x
m is subsegment 1
,
n
x
2 stiffened;
rfen is subsegment n
x
,
n
y
3 stiffened.
(ii)
Other boundary elements:
forj
=
2
,
…
,
n
y,
rfej is formed by connecting and stiffening subsegments 1
,
j
-
1
4 and 1
,
j
1;
rfen
x
m
+
j is formed by connecting and stiffening subsegments n
x
,
j
-
1
3 and n
x
,
j
2;
fori
=
1
,
…
,
n
x
-
1,
rfei
m
+
1 is formed by connecting and stiffening subsegments i
,
1
2 and i
+
1,1
1;
rfe(
i
+
1
)
m is formed by connecting and stiffening subsegments i
,
n
x
3 and i
+
1
,
n
x
4.
(iii)
Internal elements:
fori
=
1
,
…
,
n
x
-
1, j
=
1
,
…
,
n
y
-
1,
rfei
m
+
j
+
1 is formed by connecting and stiffening subsegments i
,
j
3,i
+
1
,
j
4, i
+
1
,
j
+
1
1, and i
,
j
+
1
2.Secondary division of the electrode into rfes is essential for the method discussed since in further considerations the shell is treated as a system ofn rigid bodies with 6 degrees of freedom. A local coordinate system x
^
k
C
y
^
k
C
z
^
k
C, whose axes are central inertial axes of the element, is assigned to each rfe as in Figure 6. Angle β
k defines inclination of the plane of element k with respect to plane x
^
y
^ of coordinate system x
^
y
^
z
^ assigned to the electrode.Figure 6
Rigid finite element.The vector of the generalized coordinates of rfek includes three displacements and three rotations:
(22)
q
k
′
=
Δ
k
′
φ
k
′
=
Δ
x
k
′
,
Δ
y
k
′
,
Δ
z
k
′
,
φ
k
′
x
,
φ
k
′
y
,
φ
k
′
z
T
,
where Δ
k
′
=
Δ
x
k
′
,
Δ
y
k
′
,
Δ
z
′
T and φ
k
′
=
φ
k
′
x
,
φ
k
′
y
,
φ
k
′
z
T.When the vibrations of the collecting electrodes are considered, it can be assumed that anglesφ
k
′
x, φ
k
′
y, φ
k
′
z are small. Thus, when coordinates (
x
′, y
′, z
′
) of a point are known in the local system of the kth rfe, its coordinates in the global coordinate system can be calculated as follows:
(23)
r
=
r
k
rfe
+
R
-
k
r
k
′
C
,
r
k
′
C
=
U
r
′
q
k
′
+
r
′
,
where r
′
=
x
′
,
y
′
,
z
′
T is the coordinate vector in x
^
k
C
,
y
^
k
C
,
z
^
k
C and r
k
rfe is defined in Figure 6, as the coordinate vector of point C
k for a nondeformed rfe (when q
k
′
=
0),
(24)
U
=
1
0
0
0
z
′
-
y
′
0
1
0
-
z
′
0
x
′
0
0
1
y
′
-
x
′
0
,
R
-
k
=
1
0
0
0
cos
⁡
β
k
-
sin
β
k
0
sin
β
k
cos
⁡
β
k
.


### 2.1.3. Equations of Motion of the SIGMA Electrode

Using the standard RFEM approach, the diagonal mass matrix of thekth rfe can be formulated as follows:
(25)
M
k
=
diag
⁡
m
k
,
m
k
,
m
k
,
I
k
x
,
I
k
y
,
I
k
z
.In the case of SIGMA electrodes it is relatively easy to formulate the algorithm to calculate elements of matrixM
k using the above formula and angles β
k defining the position of central inertial axes of rfes.In view of (20) and (22) the kinetic energy of an rfe can be presented in the form
(26)
T
k
=
1
2
q
˙
k
′
T
M
k
q
˙
k
′
and the kinetic energy of the whole electrode as
(27)
T
=
∑
k
=
1
n
T
k
=
1
2
q
˙
′
T
M
q
˙
′
,
where q
′
=
q
1
′
⋯
q
k
′
⋯
q
n
′
T and M
=
diag
⁡
M
1
,
M
2
,
…
,
M
n.In order to formulate the stiffness matrix for the whole electrode, the displacements of primary elements have to be expressed in terms of generalized coordinates of the rigid finite elements. The vector of generalized coordinates and the stiffness matrix for primary element(
i
,
j
) are defined in (12). Individual segments of primary element (
i
,
j
) belong to various rfes, whose indexes are determined in (20).It results from (20) that node i
,
j
p belongs to rfe:
(28)
k
i
,
j
,
p
=
i
-
1
m
+
j
,
for
p
=
1
,
i
m
+
j
,
for
p
=
2
,
i
m
+
j
+
1
,
for
p
=
3
,
i
-
1
m
+
j
+
1
,
for
p
=
4
.Coordinatesr
i
,
j
,
p in the global coordinate system of the electrode x
^
,
y
^
,
z
^  are defined by (17). Having taken into account (23), node coordinates in local system x
^
k
′
C
,
y
^
k
′
C
,
z
^
k
′
C can be defined as follows:
(29)
r
k
,
p
′
C
=
R
-
k
T
r
i
,
j
,
p
-
r
k
rfe
,
where R
-
k
T
=
R
-
k
-
1, k
=
k
(
i
,
j
,
p
) defined in (28). If rfe k
(
i
,
j
,
p
) has coordinates defined in (22), then coordinates of node i
,
j
p in the global system can be calculated from (23):
(30)
r
k
(
i
,
j
,
p
)
=
r
k
rfe
+
R
-
k
U
r
k
,
p
′
C
q
k
′
+
r
k
,
p
′
C
.In order to define the energy of spring deformation of nodei
,
j
p the node displacements and rotations
(31)
Δ
i
,
j
,
p
′
=
u
i
,
j
,
p
′
,
v
i
,
j
,
p
′
,
w
i
,
j
,
p
′
T
,
φ
i
,
j
,
p
′
=
φ
i
,
j
,
p
′
x
,
φ
i
,
j
,
p
′
y
,
φ
i
,
j
,
p
′
z
T
have to be defined by means of generalized coordinates q
k
′.The displacement vector of nodei
,
j
p belonging to rfe k, as in (28), is defined as follows:
(32)
Δ
i
,
j
,
p
′
=
R
j
T
r
k
i
,
j
,
p
-
r
i
,
j
,
p
.After some transformations, when (30), (29), and (17) are taken into account, the following is obtained:
(33)
Δ
i
,
j
,
p
′
=
R
j
T
R
-
k
U
r
k
,
p
′
C
.Rotation angles can be calculated from the formula(34)
φ
i
,
j
,
p
′
=
R
j
T
R
-
k
S
q
k
′
,
where
(35)
S
=
0
0
0
1
0
0
0
0
0
0
1
0
0
0
0
0
0
1
.
Finally, the following relation is derived:
(36)
q
i
,
j
,
p
′
=
Δ
i
,
j
,
p
′
φ
i
,
j
,
p
′
=
R
j
T
R
-
k
U
r
k
,
p
′
C
R
j
T
R
-
k
S
Δ
k
′
φ
k
′
=
H
k
,
p
q
k
′
.MatrixH
k
,
p enables the coordinates of an rfe to be transformed into the elastic coordinates of a primary element. Thus, the energy of spring deformation of primary element (
i
,
j
) can be written in the following form
(37)
E
i
,
j
=
1
2
∑
p
1
=
1
4
∑
p
2
=
1
4
q
i
,
j
,
p
1
′
C
T
p
1
,
p
2
(
j
)
q
i
,
j
,
p
2
′
,
where C
α
,
β
(
j
) is defined in (12) and (13).Bearing in mind thatk
=
k
(
i
,
j
,
p
), and taking into account (36), the above can be written as follows:
(38)
E
i
,
j
=
1
2
∑
p
1
=
1
4
∑
p
2
=
1
4
q
k
(
i
,
j
,
p
1
)
′
T
H
k
i
,
j
,
p
1
,
p
1
T
C
p
1
,
p
2
(
j
)
H
k
i
,
j
,
p
2
,
p
2
q
k
(
i
,
j
,
p
2
)
′
=
1
2
∑
p
1
=
1
4
∑
p
2
=
1
4
q
k
(
i
,
j
,
p
1
)
′
T
C
p
1
,
p
2
(
i
,
j
)
q
k
(
i
,
j
,
p
2
)
′
,
where C
p
1
,
p
2
(
i
,
j
)
=
H
k
(
i
,
j
,
p
1
)
,
p
1
T
C
p
1
,
p
2
(
j
)
H
k
(
i
,
j
,
p
2
)
,
p
2.This leads to a straightforward method of generating stiffness matrixC using coordinates q
k
′ as generalized coordinates. The appropriate algorithm is as follows:(1)
fork
1
=
1
,
…
,
6
n
;
k
2
=
1
,
…
,
6
n, accept: c
k
1
,
k
2
=
0.
(2)
fori
=
1
,
…
,
n
x; j
=
1
,
…
,
n
y
(i)
forp
=
1
,
…
,
4 compute H
k
(
i
,
j
,
p
)
,
p according to (36)
(ii)
forp
1
=
1
,
…
,
4; p
2
=
1
,
…
,
4
(a)
compute:(39)
k
1
=
k
i
,
j
,
p
1
,
k
2
=
i
,
j
,
p
2
k
1
C
=
6
k
1
-
1
,
k
2
C
=
6
k
2
-
1
(b)
forj
1
=
1
,
…
,
6; j
2
=
1
,
…
,
6 compute:
(40)
c
k
1
C
+
j
1
,
k
2
C
+
j
2
=
c
k
1
C
+
j
1
,
k
2
C
+
j
2
+
c
p
1
,
p
2
(
i
,
j
)
.It is important to note that matrixC is a sparse band matrix and this is used to facilitate computations of dynamics of the collecting electrodes. The potential energy of gravity forces is calculated as follows:
(41)
V
g
=
-
∑
k
=
1
n
x
m
k
g
x
k
C
,
where g is the acceleration of gravity, x
k
C
=
θ
1
r
k
C, and θ
1
=
1
0
0. Bearing in mind (23) we can calculate derivatives of the potential energy
(42)
∂
V
g
∂
q
k
′
=
-
m
k
g
θ
1
R
-
k
U
k
C
.Since the origin of the local coordinate system of rfek is located at the center of mass of the rfe, that is, U
k
C
=
I
0, the derivative of potential energy of gravity forces for the jth strip of the electrode can be presented in the matrix form
(43)
∂
V
g
∂
q
′
=
G
=
G
1
⋯
G
k
⋯
G
n
x
T
,
where G
k
=
-
g
m
k
0
0
0
0
0
T.Equations of motion of a free electrode can be written in the form(44)
M
q
¨
′
+
C
q
′
=
-
∂
V
g
∂
q
′
,
and the eigenvalue problem takes the form
(45)
det
⁡
(
C
-
ω
2
M
)
=
0
.Some important differences between the hybrid finite element method and the finite element method applied to modeling of electrodes are listed in Table1.Table 1
Comparison of features of the finite element method and hybrid finite element method.
Feature
Finite element method
Hybrid element method
Generalized coordinates
Displacements of nodes of elementsq
′
i
,
j
(
p
) (elastic coordinates) (36)
Displacements and rotations of rigid elementsq
′
k (rigid coordinates) (22)
Stiffness matrixC
Interpolation of field displacements and rotations by means of continuous functions
Interpolation of field displacements and rotations by means of continuous functions
Additional transformation of elastic coordinates to rigid coordinates:
q
′
i
,
j
(
p
)
=
H
k
,
p
q
k
′
Mass matrixM
Interpolation of field displacements and rotation velocities by means of continuous functions
M
=
diag
⁡
{
M
1
,
M
2
,
…
,
M
n
}
M
k
=
diag
⁡
{
m
k
,
m
k
,
m
k
,
I
k
x
,
I
k
y
,
I
k
z
}    Inertial moments are taken into account
Additional operations in order to use lumped mass


## 2.1.1. Primary Division

Let us assign a coordinate systemx
^
y
^
z
^ as in Figure 2. Discretization is carried out along axis y
^ into n
y strips and along axis x
^ into n
x segments with the length
(1)
Δ
x
=
L
n
x
,
where L is the length of the electrode. Thus the electrode is divided into
(2)
n
=
n
x
n
y
elements.Let us consider a rectangular element as in Figure3. The following vectors describe displacements of nodes for such an element:
(3)
q
p
′
=
u
p
′
,
v
p
′
,
w
p
′
,
φ
p
′
x
,
φ
p
′
y
,
φ
p
′
z
T
for
p
=
1
,
2
,
3
,
4
,
where u
p
′, v
p
′, w
p
′ are displacements of node p along x
^
′, y
^
′, z
^
′, respectively; φ
p
′
x, φ
p
′
y, φ
p
′
z are respective rotations in node p about axes parallel to x
^
′, y
^
′, z
^
′.Figure 3
A rectangular shell element.Shape functions are defined as(4)
u
′
=
a
1
u
+
a
2
u
x
′
+
a
3
u
y
′
+
a
4
u
x
′
y
′
+
a
5
u
y
′
2
+
a
6
u
x
′
y
′
2
,
v
′
=
a
1
v
+
a
2
v
x
′
+
a
3
v
y
′
+
a
4
v
x
′
y
′
+
a
5
v
x
′
2
+
a
6
v
x
′
2
y
′
,
w
′
=
a
1
w
+
a
2
w
x
′
+
a
3
w
y
′
+
a
4
w
x
′
2
+
a
5
w
x
′
y
′
+
a
6
w
y
′
2
+
a
7
w
x
′
3
+
a
8
w
x
′
2
y
′
+
a
9
w
x
′
y
′
2
+
a
10
w
y
′
3
+
a
11
w
x
′
3
y
′
+
a
12
w
x
′
y
′
3
.It is assumed thatu
p
′, v
p
′, and φ
p
′
z describe shield deformations, while w
p
′, φ
p
′
x, and φ
p
′
ydescribe plate deformations. The angles are defined by the following relations:
(5)
φ
′
x
=
∂
w
′
∂
y
′
,
φ
′
y
=
∂
w
′
∂
x
′
,
φ
′
z
=
1
2
∂
v
′
∂
x
′
-
∂
u
′
∂
y
′
.The elastic strain energy is a sum of energies corresponding to the shield and plate states and can be calculated as [13]
(6)
E
=
1
2
∫
V
σ
(
s
)
T
ε
(
s
)
d
V
+
1
2
∫
V
σ
(
p
)
T
ε
(
p
)
d
V
,
where V
=
Δ
x
Δ
y
h, h is the thickness of the element,
(7)
ε
(
s
)
=
ε
x
s
ε
y
s
γ
x
y
s
=
∂
u
′
∂
x
′
∂
v
′
∂
y
′
∂
u
′
∂
y
′
+
∂
v
′
∂
x
′
,
ε
(
p
)
=
ε
x
p
ε
y
p
γ
x
y
p
=
-
z
∂
2
w
′
∂
x
′
2
∂
2
w
′
∂
y
′
2
2
∂
w
′
∂
x
′
∂
y
′
are vectors of strains,
(8)
σ
(
l
)
=
σ
x
l
σ
y
l
σ
x
y
l
=
D
ε
(
l
)
,
l
∈
s
,
p
are vectors of stresses,
(9)
D
=
E
1
-
υ
2
1
υ
0
υ
1
0
0
0
1
-
υ
2
,
E is Young’s modulus, and υ is the Poisson ratio.After the necessary calculations, the energy of deformation of the element can be expressed in the following form:(10)
E
=
1
2
s
-
T
C
(
s
)
s
-
+
p
-
T
C
(
p
)
p
-
,
where
(11)
s
¯
=
u
1
′
v
1
′
φ
1
′
z
u
2
′
v
2
′
φ
2
′
z
u
3
′
v
3
′
φ
3
′
z
u
4
′
v
4
′
φ
4
′
z
T
,
p
¯
=

w
1
′
φ
1
′
x
φ
1
′
y
w
2
′
φ
2
′
x
φ
2
′
y
w
3
′
φ
3
′
x
φ
2
′
y
w
4
′
φ
4
′
x
φ
2
′
y
T
.The dimension of matricesC
(
s
) and C
(
p
) is 12
×
12. Integrals occurring in elements of matrices can be calculated analytically or by means of Gauss quadratures. Detailed formulae for elements of matrices C
(
s
) and C
(
p
) are presented in [14]. If we assume that the form of these matrices is as follows:
(12)
C
(
l
)
=
C
11
(
l
)
C
12
(
l
)
C
13
(
l
)
C
14
(
l
)
C
21
(
l
)
C
22
(
l
)
C
23
(
l
)
C
24
(
l
)
C
31
(
l
)
C
32
(
l
)
C
33
(
l
)
C
34
(
l
)
C
41
(
l
)
C
42
(
l
)
C
43
(
l
)
C
44
(
l
)
,
l
∈
s
,
p
,
then having taken into account the definition of vectors of generalized coordinates describing plate and shield deformations in (10), the elastic strain energy can be presented in terms of the generalised coordinates as follows:
(13)
E
=
1
2
q
1
′
T
q
2
′
T
q
3
′
T
q
4
′
T
C
11
C
12
C
13
C
14
C
21
C
22
C
23
C
24
C
31
C
32
C
33
C
34
C
41
C
42
C
43
C
44
q
1
′
q
2
′
q
3
′
q
4
′
=
1
2
q
′
T
C
q
′
,
where C is a stiffness matrix with 24
×
24 elements,
(14)
q
′
=
q
1
′
q
2
′
q
3
′
q
4
′
,
and nonzero elements of matrixC from (13) are defined as follows:
(15)
C
i
,
k
q
,
t
=
C
i
,
k
(
s
)
q
,
s
for
q
,
t
=
1,2
C
i
,
k
q
,
6
=
C
i
,
k
(
s
)
q
,
3
for
q
=
1,2
C
i
,
k
q
+
2
,
t
+
2
=
C
i
,
k
(
p
)
q
,
t
for
q
,
t
=
1,2
,
3
C
i
,
k
6
,
q
=
C
i
,
k
(
s
)
3
,
q
for
q
=
1,2
C
i
,
k
6,6
=
C
i
,
k
(
s
)
3,3
.For the assumed values of parametersE
,
υ
,
α
,
h, the stiffness matrices depend only on Δ
x and Δ
y. Having used denotations from Figures 2 and 3 and relation (1), in the case considered, we have C
(
i
,
j
)
=
C
(
j
)
=
C
Δ
x
,
Δ
y
j for i
=
1
,
…
,
n
x,  j
=
1
,
…
,
n
y. Thus, as a result of constant length Δ
x of elements, stiffness matrices for each strip are identical and are calculated only for j
=
1
,
…
,
n
y.The coordinatesx
i
,
j
C
,
y
i
,
j
C
,
z
i
,
j
C of the geometrical center of element (i, j) can be calculated according to the formulas
(16)
x
i
,
j
C
=
i
-
1
Δ
x
+
Δ
x
2
,
y
i
,
j
C
=
y
0
+
∑
s
=
0
j
-
1
b
s
cos
⁡
⁡
α
s
+
b
j
2
cos
⁡
⁡
α
j
,
z
i
,
j
C
=
z
0
+
∑
s
=
0
j
-
1
b
s
sin
⁡
α
s
+
b
j
2
sin
⁡
α
j
,
where y
0, z
0 define coordinates of the left edge of strip number 1 and b
j, α
j are the width and angle of inclination of the jth strip to axis y
^, respectively.Coordinates of nodes of primary elements can be calculated as follows:(17)
r
i
,
j
,
p
=
r
i
,
j
C
+
R
j
r
i
,
j
,
p
′
,
where r
i
,
j
C
=
x
i
,
j
C
y
i
,
j
C
z
i
,
j
C
T is defined in (16),
(18)
R
j
=
1
0
0
0
cos
⁡
α
j
-
sin
α
j
0
sin
α
j
cos
⁡
α
j
,
and r
i
,
j
,
p
′ is the vector of coordinates of the mass centre of the pth segment in the local coordinate system:
(19)
r
i
,
j
,
1
′
=
-
Δ
x
2
-
b
j
2
0
T
,
r
i
,
j
,
2
′
=
Δ
x
2
-
b
j
2
0
T
,
r
i
,
j
,
3
′
=
Δ
x
2
b
j
2
0
T
,
r
i
,
j
,
4
′
=
-
Δ
x
2
b
j
2
0
T
.


## 2.1.2. Secondary Division of the SIGMA Electrode

In the rigid finite element method the flexible body considered is divided into rigid elements reflecting inertial features of the body and spring-damping elements which connect the rigid elements and adopt spring and damping features [9]. For the subsequent analysis damping will be omitted. Since in the hybrid finite element method proposed, spring features of the electrode are included in the elements described above, the method of discretisation into rigid elements is discussed below.Let us consider element (i, j) from the primary division, as in Figure 4.Figure 4
Secondary division of the shell.Coordinate systemx
^, y
^, z
^ is a global coordinate system. Primary element (i, j) is divided into four equal parts called subsegments. Then the rigid finite element (rfe) is formed from one, two, or four subdivisions, which may belong to different primary elements. Thus each segment of the primary element belongs to the following rfes:
(20)
i
,
j
1
=
i
-
1
m
+
j
,
i
,
j
2
=
i
m
+
j
,
i
,
j
3
=
i
m
+
j
+
1
,
i
,
j
4
=
i
-
1
m
+
j
+
1
,
for i
=
1
,
…
,
n
x, j
=
1
,
…
,
n
y, and m
=
n
y
+
1.Figure5 shows an example of primary and secondary division of a rectangular area for n
x
=
3 and n
y
=
4.Figure 5
Division of a rectangular area inton
x
=
3 and n
y
=
4 elements.The number of rigid finite elements is(21)
n
=
n
x
+
1
m
=
n
x
+
1
n
y
+
1
.It can be seen from Figures4 and 5 that rfes are constructed by rigid connections of stiffened subsegments in the following way.(i)
Corner elements:
rfe 1 is subsegment1,1
1 stiffened;
rfem is subsegment 1
,
n
y
4 stiffened;
rfen
x
m is subsegment 1
,
n
x
2 stiffened;
rfen is subsegment n
x
,
n
y
3 stiffened.
(ii)
Other boundary elements:
forj
=
2
,
…
,
n
y,
rfej is formed by connecting and stiffening subsegments 1
,
j
-
1
4 and 1
,
j
1;
rfen
x
m
+
j is formed by connecting and stiffening subsegments n
x
,
j
-
1
3 and n
x
,
j
2;
fori
=
1
,
…
,
n
x
-
1,
rfei
m
+
1 is formed by connecting and stiffening subsegments i
,
1
2 and i
+
1,1
1;
rfe(
i
+
1
)
m is formed by connecting and stiffening subsegments i
,
n
x
3 and i
+
1
,
n
x
4.
(iii)
Internal elements:
fori
=
1
,
…
,
n
x
-
1, j
=
1
,
…
,
n
y
-
1,
rfei
m
+
j
+
1 is formed by connecting and stiffening subsegments i
,
j
3,i
+
1
,
j
4, i
+
1
,
j
+
1
1, and i
,
j
+
1
2.Secondary division of the electrode into rfes is essential for the method discussed since in further considerations the shell is treated as a system ofn rigid bodies with 6 degrees of freedom. A local coordinate system x
^
k
C
y
^
k
C
z
^
k
C, whose axes are central inertial axes of the element, is assigned to each rfe as in Figure 6. Angle β
k defines inclination of the plane of element k with respect to plane x
^
y
^ of coordinate system x
^
y
^
z
^ assigned to the electrode.Figure 6
Rigid finite element.The vector of the generalized coordinates of rfek includes three displacements and three rotations:
(22)
q
k
′
=
Δ
k
′
φ
k
′
=
Δ
x
k
′
,
Δ
y
k
′
,
Δ
z
k
′
,
φ
k
′
x
,
φ
k
′
y
,
φ
k
′
z
T
,
where Δ
k
′
=
Δ
x
k
′
,
Δ
y
k
′
,
Δ
z
′
T and φ
k
′
=
φ
k
′
x
,
φ
k
′
y
,
φ
k
′
z
T.When the vibrations of the collecting electrodes are considered, it can be assumed that anglesφ
k
′
x, φ
k
′
y, φ
k
′
z are small. Thus, when coordinates (
x
′, y
′, z
′
) of a point are known in the local system of the kth rfe, its coordinates in the global coordinate system can be calculated as follows:
(23)
r
=
r
k
rfe
+
R
-
k
r
k
′
C
,
r
k
′
C
=
U
r
′
q
k
′
+
r
′
,
where r
′
=
x
′
,
y
′
,
z
′
T is the coordinate vector in x
^
k
C
,
y
^
k
C
,
z
^
k
C and r
k
rfe is defined in Figure 6, as the coordinate vector of point C
k for a nondeformed rfe (when q
k
′
=
0),
(24)
U
=
1
0
0
0
z
′
-
y
′
0
1
0
-
z
′
0
x
′
0
0
1
y
′
-
x
′
0
,
R
-
k
=
1
0
0
0
cos
⁡
β
k
-
sin
β
k
0
sin
β
k
cos
⁡
β
k
.


## 2.1.3. Equations of Motion of the SIGMA Electrode

Using the standard RFEM approach, the diagonal mass matrix of thekth rfe can be formulated as follows:
(25)
M
k
=
diag
⁡
m
k
,
m
k
,
m
k
,
I
k
x
,
I
k
y
,
I
k
z
.In the case of SIGMA electrodes it is relatively easy to formulate the algorithm to calculate elements of matrixM
k using the above formula and angles β
k defining the position of central inertial axes of rfes.In view of (20) and (22) the kinetic energy of an rfe can be presented in the form
(26)
T
k
=
1
2
q
˙
k
′
T
M
k
q
˙
k
′
and the kinetic energy of the whole electrode as
(27)
T
=
∑
k
=
1
n
T
k
=
1
2
q
˙
′
T
M
q
˙
′
,
where q
′
=
q
1
′
⋯
q
k
′
⋯
q
n
′
T and M
=
diag
⁡
M
1
,
M
2
,
…
,
M
n.In order to formulate the stiffness matrix for the whole electrode, the displacements of primary elements have to be expressed in terms of generalized coordinates of the rigid finite elements. The vector of generalized coordinates and the stiffness matrix for primary element(
i
,
j
) are defined in (12). Individual segments of primary element (
i
,
j
) belong to various rfes, whose indexes are determined in (20).It results from (20) that node i
,
j
p belongs to rfe:
(28)
k
i
,
j
,
p
=
i
-
1
m
+
j
,
for
p
=
1
,
i
m
+
j
,
for
p
=
2
,
i
m
+
j
+
1
,
for
p
=
3
,
i
-
1
m
+
j
+
1
,
for
p
=
4
.Coordinatesr
i
,
j
,
p in the global coordinate system of the electrode x
^
,
y
^
,
z
^  are defined by (17). Having taken into account (23), node coordinates in local system x
^
k
′
C
,
y
^
k
′
C
,
z
^
k
′
C can be defined as follows:
(29)
r
k
,
p
′
C
=
R
-
k
T
r
i
,
j
,
p
-
r
k
rfe
,
where R
-
k
T
=
R
-
k
-
1, k
=
k
(
i
,
j
,
p
) defined in (28). If rfe k
(
i
,
j
,
p
) has coordinates defined in (22), then coordinates of node i
,
j
p in the global system can be calculated from (23):
(30)
r
k
(
i
,
j
,
p
)
=
r
k
rfe
+
R
-
k
U
r
k
,
p
′
C
q
k
′
+
r
k
,
p
′
C
.In order to define the energy of spring deformation of nodei
,
j
p the node displacements and rotations
(31)
Δ
i
,
j
,
p
′
=
u
i
,
j
,
p
′
,
v
i
,
j
,
p
′
,
w
i
,
j
,
p
′
T
,
φ
i
,
j
,
p
′
=
φ
i
,
j
,
p
′
x
,
φ
i
,
j
,
p
′
y
,
φ
i
,
j
,
p
′
z
T
have to be defined by means of generalized coordinates q
k
′.The displacement vector of nodei
,
j
p belonging to rfe k, as in (28), is defined as follows:
(32)
Δ
i
,
j
,
p
′
=
R
j
T
r
k
i
,
j
,
p
-
r
i
,
j
,
p
.After some transformations, when (30), (29), and (17) are taken into account, the following is obtained:
(33)
Δ
i
,
j
,
p
′
=
R
j
T
R
-
k
U
r
k
,
p
′
C
.Rotation angles can be calculated from the formula(34)
φ
i
,
j
,
p
′
=
R
j
T
R
-
k
S
q
k
′
,
where
(35)
S
=
0
0
0
1
0
0
0
0
0
0
1
0
0
0
0
0
0
1
.
Finally, the following relation is derived:
(36)
q
i
,
j
,
p
′
=
Δ
i
,
j
,
p
′
φ
i
,
j
,
p
′
=
R
j
T
R
-
k
U
r
k
,
p
′
C
R
j
T
R
-
k
S
Δ
k
′
φ
k
′
=
H
k
,
p
q
k
′
.MatrixH
k
,
p enables the coordinates of an rfe to be transformed into the elastic coordinates of a primary element. Thus, the energy of spring deformation of primary element (
i
,
j
) can be written in the following form
(37)
E
i
,
j
=
1
2
∑
p
1
=
1
4
∑
p
2
=
1
4
q
i
,
j
,
p
1
′
C
T
p
1
,
p
2
(
j
)
q
i
,
j
,
p
2
′
,
where C
α
,
β
(
j
) is defined in (12) and (13).Bearing in mind thatk
=
k
(
i
,
j
,
p
), and taking into account (36), the above can be written as follows:
(38)
E
i
,
j
=
1
2
∑
p
1
=
1
4
∑
p
2
=
1
4
q
k
(
i
,
j
,
p
1
)
′
T
H
k
i
,
j
,
p
1
,
p
1
T
C
p
1
,
p
2
(
j
)
H
k
i
,
j
,
p
2
,
p
2
q
k
(
i
,
j
,
p
2
)
′
=
1
2
∑
p
1
=
1
4
∑
p
2
=
1
4
q
k
(
i
,
j
,
p
1
)
′
T
C
p
1
,
p
2
(
i
,
j
)
q
k
(
i
,
j
,
p
2
)
′
,
where C
p
1
,
p
2
(
i
,
j
)
=
H
k
(
i
,
j
,
p
1
)
,
p
1
T
C
p
1
,
p
2
(
j
)
H
k
(
i
,
j
,
p
2
)
,
p
2.This leads to a straightforward method of generating stiffness matrixC using coordinates q
k
′ as generalized coordinates. The appropriate algorithm is as follows:(1)
fork
1
=
1
,
…
,
6
n
;
k
2
=
1
,
…
,
6
n, accept: c
k
1
,
k
2
=
0.
(2)
fori
=
1
,
…
,
n
x; j
=
1
,
…
,
n
y
(i)
forp
=
1
,
…
,
4 compute H
k
(
i
,
j
,
p
)
,
p according to (36)
(ii)
forp
1
=
1
,
…
,
4; p
2
=
1
,
…
,
4
(a)
compute:(39)
k
1
=
k
i
,
j
,
p
1
,
k
2
=
i
,
j
,
p
2
k
1
C
=
6
k
1
-
1
,
k
2
C
=
6
k
2
-
1
(b)
forj
1
=
1
,
…
,
6; j
2
=
1
,
…
,
6 compute:
(40)
c
k
1
C
+
j
1
,
k
2
C
+
j
2
=
c
k
1
C
+
j
1
,
k
2
C
+
j
2
+
c
p
1
,
p
2
(
i
,
j
)
.It is important to note that matrixC is a sparse band matrix and this is used to facilitate computations of dynamics of the collecting electrodes. The potential energy of gravity forces is calculated as follows:
(41)
V
g
=
-
∑
k
=
1
n
x
m
k
g
x
k
C
,
where g is the acceleration of gravity, x
k
C
=
θ
1
r
k
C, and θ
1
=
1
0
0. Bearing in mind (23) we can calculate derivatives of the potential energy
(42)
∂
V
g
∂
q
k
′
=
-
m
k
g
θ
1
R
-
k
U
k
C
.Since the origin of the local coordinate system of rfek is located at the center of mass of the rfe, that is, U
k
C
=
I
0, the derivative of potential energy of gravity forces for the jth strip of the electrode can be presented in the matrix form
(43)
∂
V
g
∂
q
′
=
G
=
G
1
⋯
G
k
⋯
G
n
x
T
,
where G
k
=
-
g
m
k
0
0
0
0
0
T.Equations of motion of a free electrode can be written in the form(44)
M
q
¨
′
+
C
q
′
=
-
∂
V
g
∂
q
′
,
and the eigenvalue problem takes the form
(45)
det
⁡
(
C
-
ω
2
M
)
=
0
.Some important differences between the hybrid finite element method and the finite element method applied to modeling of electrodes are listed in Table1.Table 1
Comparison of features of the finite element method and hybrid finite element method.
Feature
Finite element method
Hybrid element method
Generalized coordinates
Displacements of nodes of elementsq
′
i
,
j
(
p
) (elastic coordinates) (36)
Displacements and rotations of rigid elementsq
′
k (rigid coordinates) (22)
Stiffness matrixC
Interpolation of field displacements and rotations by means of continuous functions
Interpolation of field displacements and rotations by means of continuous functions
Additional transformation of elastic coordinates to rigid coordinates:
q
′
i
,
j
(
p
)
=
H
k
,
p
q
k
′
Mass matrixM
Interpolation of field displacements and rotation velocities by means of continuous functions
M
=
diag
⁡
{
M
1
,
M
2
,
…
,
M
n
}
M
k
=
diag
⁡
{
m
k
,
m
k
,
m
k
,
I
k
x
,
I
k
y
,
I
k
z
}    Inertial moments are taken into account
Additional operations in order to use lumped mass


## 2.2. Models of Beams

The detailed description of the rigid finite element method is presented in [8, 9]. In this paper only the general idea and final equations of motion will be given. More details can be found in [15]. The general idea of discretisation of the top and bottom beams is presented in Figure 7.Figure 7
Connections between the beams and thejth plate and discretisation of beams.The number of elements into which the beams are discretized may be different for both beams. Ifp is the number of electrodes and a is the length of a section of the beam with respect to a single electrode (all electrodes in a rapping system are usually the same) (Figure 2), then in the primary division the length of rigid elements is as follows:
(46)
d
0
s
=
a
L
s
+
a
-
s
2
,
d
i
s
=
a
-
s
,
d
p
s
s
=
a
-
s
2
+
a
R
s
,
mmmmmmlmmmh
for
i
=
1
,
…
,
p
s
-
1
,
where a
-
s
=
a
/
n
-
s, s
∈
t
,
b; t, b denote top or bottom beam, respectively; and p
s
=
n
-
s
p.The secondary division leads to a system of rigid finite elements with the length equal to(47)
l
0
s
=
d
0
s
2
,
l
i
s
=
d
i
-
1
s
+
d
i
s
2
,
l
p
s
+
1
s
=
d
p
s
s
2
,
mmmmmmmmmmml
for
i
=
1
,
…
,
p
s
.We treat both beams as those with a uniform cross-section, whose geometrical parameters exactly reflect mass and stiffness parameters of the real beams.The axes of the local coordinate system assigned to each rfe before deformation are parallel to the axes of the global reference system and are the main central inertial axes of the element. Thus the mass matrix of the element is diagonal and takes the following form:(48)
M
i
s
=
diag
⁡
m
i
s
m
i
s
m
i
s
I
i
,
x
s
I
i
,
y
s
I
i
,
z
s
,
where m
i
s
=
ρ
l
i
s
A
s is the mass of the ith rfe, A
s is the cross-section area, I
x
s
,
I
z
s are inertial moments of the cross-section of the beam, I
i
,
x
s
=
m
i
s
l
i
s
2
/
12
+
I
x
s
/
A
s,  I
i
,
y
s
=
m
i
s
/
A
s
I
x
s
+
I
z
s, and I
i
,
z
s
=
m
i
s
l
i
s
2
/
12
+
I
z
s
/
A
s. The vector of generalized coordinates of rfe i is as follows:
(49)
q
i
s
=
x
~
i
s
y
~
i
s
z
~
i
s
φ
~
i
s
θ
~
i
s
ψ
~
i
s
T
,
where x
~
i
s, y
~
i
s, z
~
i
s are displacements of the center of mass S
i
s of the rfe and φ
~
i
s, θ
~
i
s, ψ
~
i
s are the angles of rotations. The kinetic energy of the ith rfe is calculated according to the formula
(50)
E
i
s
=
1
2
q
˙
i
s
T
M
i
s
q
˙
i
s
.The energy of spring deformation can be presented in the form(51)
V
e
,
i
s
=
1
2
Δ
i
s
T
C
i
s
Δ
i
s
,
where C
i
s
=
diag
⁡
c
i
,
1
s
⋯
c
i
,
6
s, c
i
,
j
s are the stiffness coefficients calculated according to the formulas given in [9], and Δ
i
s is the deformation of spring damping element (sde) i, which depends on q
i
s and q
i
-
1
s.The potential energy of gravity forces of each rfe can be calculated according to the formula(52)
V
s
,
i
g
=
-
m
i
g
g
x
S
i
g
,
where x
S
i
g
=
θ
1
r
S
i
g and θ
1 is defined in (41).Kinetic and potential energies of the beams are the sum of energies of all rfes:(53)
E
s
=
∑
i
=
0
p
s
+
1
E
i
s
,
V
g
s
=
∑
i
=
0
p
s
+
1
V
g
,
i
s
,
and the potential energy of spring deformation is the sum of energies of sdes 0 to p
s:
(54)
V
e
s
=
∑
i
=
0
p
s
V
e
,
i
s
.The vector of generalized coordinates of a beam can be presented as(55)
q
s
=
q
0
s
⋯
q
i
s
⋯
q
p
s
+
1
s
T
.So finally, the following can be written:(56)
d
d
t
∂
E
s
∂
q
˙
s
-
∂
E
s
∂
q
s
+
∂
V
g
s
+
V
e
s
∂
q
s
=
M
s
q
¨
s
+
K
s
q
s
+
G
s
,
where M
s is the diagonal mass matrix, K
s is the stiffness matrix with constant coefficients, and G
s is the vector of gravity forces, s
∈
t
,
b. Thus, the motion of beams is described by the following number of generalized coordinates:
(57)
n
dof
s
=
6
p
s
+
2
.The top beam is connected with the base, and this connection has to be taken into account when deriving the equations of motion.As has been mentioned, the vibrations in the rapping system are excited, when the beater hits the anvil. It is assumed that the impulse of the force is modelled as a central, straight stroke of forceF
=
0
-
S
(
t
)
0
T, and the resulting generalized force follows
(58)
Q
p
b
+
1
(
S
)
=
0
-
S
0
0
0
0
T
.


## 2.3. Modeling Elastic Connections

The collecting electrodes are connected with the suspension and brushing bars and they may be also connected to each other. Such connections can be easily included in the model by means of an elastic connection of two rigid bodies. Let us consider two rigid bodies connected at pointE by means of spring-damping element (sde) as shown in Figure 8.Figure 8
Elastic connectionE between rfes l and r.We assume the following denotations:s
- denotes the system before deformation and s
~ denotes the system rigidly assigned to rfe s, where s
∈
r
,
l.Owing to (22) and (49), the motion of rfes is governed by the vectors
(59)
q
~
s
=
Δ
~
s
φ
~
s
,
s
∈
l
,
r
,
where Δ
~
s
=
x
~
s
y
~
s
z
~
s
T, φ
~
s
=
φ
~
s
θ
~
s
ψ
~
s
T, x
~
s, y
~
s, z
~
s describe the position of the beginning of system {
s
~
} in {
s
-
}, and φ
~
s, θ
~
s, ψ
~
s describe the orientation of system {
s
~
} with respect to {
s
-
}. After deformation, rfes l and r move, and thus the position of point E in the global coordinate system may vary depending on whether point E is considered as a point of rfe l or r:
(60)
r
E
l
=
r
c
l
+
R
l
Δ
~
l
+
R
~
l
r
~
E
l
,
r
E
r
=
r
c
r
+
R
r
Δ
~
r
+
R
~
r
r
~
E
r
,
where r
c
(
s
) is the position vector of the beginning of {
s
-
} in the global system {
}, r
~
E
(
s
) is the vector of coordinates of point E in system {
s
~
}, R
s is the rotation matrix of direction cosines of system {
s
-
} with respect to {
},
(61)
R
~
(
s
)
=
1
-
ψ
~
(
s
)
θ
~
(
s
)
ψ
~
(
s
)
1
-
φ
~
(
s
)
-
θ
~
(
s
)
φ
~
i
(
s
)
1
when the rotation angles are small, and r
~
(
s
) are coordinates of a point in system {
s
~
} rigidly connected with rfe s.Vectors of translational and rotational displacements of sdeE with respect to the global system can be calculated as follows:
(62)
Δ
r
E
=
r
E
(
r
)
-
r
E
(
l
)
=
R
r
x
~
(
r
)
+
S
~
E
(
r
)
φ
~
(
r
)
-
R
l
x
~
(
l
)
+
S
~
E
(
l
)
φ
~
(
l
)
,
Δ
φ
E
=
R
r
φ
~
(
r
)
-
R
l
φ
~
(
l
)
,
where
(63)
S
~
E
(
s
)
=
0
z
~
E
(
s
)
-
y
~
E
(
s
)
-
z
~
E
(
s
)
0
x
~
E
(
s
)
y
~
E
(
s
)
-
x
~
E
(
s
)
0
,
s
∈
r
,
l
.After transformations the following can be obtained:(64)
Δ
~
E
=
Δ
r
~
E
Δ
φ
~
E
=
U
E
(
r
)
q
~
r
-
U
E
l
q
~
l
,
where
(65)
U
E
(
r
)
=
R
E
T
R
r
R
E
T
R
r
S
~
E
(
r
)
0
R
E
T
R
r
,
U
E
(
l
)
=
R
E
T
R
l
R
E
T
R
l
S
~
E
(
l
)
0
R
E
T
R
l
.
It is important to note that the matrices S
~
E
(
l
), S
~
E
(
r
), R
l, and R
r are constant.Finally, the energy of spring deformation of the connection can be presented in the following form:(66)
V
E
=
1
2
Δ
~
E
T
C
E
Δ
~
E
,
where C
E
=
diag
⁡
c
E
x
,
c
E
y
,
c
E
z
,
c
E
φ
,
c
E
θ
,
c
E
ψ. It can be seen from (64) that the energy of spring deformation of the connection depends on generalized coordinates of rfes l and r. Thus the following components will appear in the equations of motion
(67)
∂
V
E
∂
q
~
(
l
)
=
C
l
l
(
E
)
q
~
(
l
)
+
C
l
r
(
E
)
q
~
(
r
)
,
∂
V
E
∂
q
~
(
r
)
=
C
r
l
(
E
)
q
~
(
l
)
+
C
r
r
(
E
)
q
~
(
r
)
,
where
(68)
C
l
l
(
E
)
=
U
E
(
l
)
T
C
E
U
E
(
l
)
,
C
r
r
(
E
)
=
U
E
(
t
)
T
C
E
U
E
(
t
)
,
C
l
r
(
E
)
=
C
r
l
(
E
)
T
=
-
U
E
(
l
)
T
C
E
U
E
(
r
)
.It should be stressed that the element presented is universal. When stiffness coefficients in matrixC
E are assumed to be large enough or equal to zero, any joint (rotary, longitudinal, and spherical) can be modeled.In order to model connections between beams and plates, two joints have been used (as shown in Figure9) and the coefficients of stiffness have been assumed to be c
T
=
1
0
12N/m for translational displacements and as c
R
=
1
0
9 N/rad in the case of rotary displacements.Models of connections between beams and a plate: (a) rotary joint and (b) universal joint.
(a)
(b)


## 2.4. Equations of Motion

The equations of motion for the system of the beams and electrodes can be written in the following form:(69)
M
t
q
¨
t
+
K
t
q
t
=
-
G
t
,
M
j
q
¨
j
+
K
j
q
j
=
-
G
j
,
j
=
1
,
…
,
n
p
,
M
b
q
¨
b
+
K
b
q
b
=
-
G
b
+
f
b
(
S
)
,
where n
p is the number of active electrodes, M
t, K
t, G
t, M
b, K
b, G
b are defined in Section 2.2,   M
j, K
j, G
j are defined in Section 2.1,   j denotes the number of the electrode, f
b
=
0
,
…
,
Q
p
+
1
S
T
T, and Q
p
+
1
(
S
) is defined in (58).It has been assumed that the electrodes are connected, as in Figure10, with the top beam by means of one elastic connection (rotary joint, Figure 9(a)) and with the bottom beam by two elastic connections (all springs are active, Figure 9(b)).Figure 10
Connection of the electrode with the beams.Finally, the equations of motion after necessary calculations take the following form:(70)
M
t
q
¨
t
+
C
11
q
t
+
∑
j
=
1
p
a
C
j
(
1
)
q
j
=
-
G
t
,
M
j
q
¨
j
+
C
j
(
1
)
T
q
t
+
C
(
j
)
q
j
+
C
j
3
T
q
b
=
-
G
j
,
j
=
1
,
…
,
n
p
,
M
b
q
¨
b
+
∑
j
=
1
p
a
C
j
(
3
)
q
j
+
C
33
q
b
=
-
G
b
+
f
b
.The coupling resulting from the connections between the electrodes is reflected by new stiffness matricesC
j
(
1
)and C
j
(
3
), while C
11, C
(
j
), and C
33 are modified matrices K
t, K
j, and K
b supplemented with the connections of the beams with the electrodes, respectively.


## 3. Numerical Simulations and Verification of the Method

The model presented has been implemented in Delphi 7.0 on a computer with Microsoft Windows Vista Business 64 Bit. All simulations have been carried out using a personal computer with processor INTEL CORE 2 QUAD Q9550 2.83 GHz LGA775 BOX and 8 MB of RAM.


### 3.1. Free Vibrations

The analysis of free vibrations is carried out for a free rectangular plate with the following dimensions: widtha
=
1 m, length b
=
4 m, and thickness h
=
0.002 m, and the following material parameters: Young modulus E
=
2.06
×
10
11 N/m2, ρ
=
7850 kg/m3, and Poisson number ν
=
0.3. Figure 11 presents the first fifteen frequencies of free vibrations obtained from the hybrid model for a different proportion defined by a coefficient η according to the formula
(71)
η
=
Δ
y
Δ
x
,
where Δ
x
=
a
/
n
x and Δ
y
=
b
/
n
y.Figure 11
Influence of coefficientη on frequencies of free vibrations.It can be seen that each of the analyzed frequencies achieves the maximal value forη
=
1. Thus, for further analysis square elements (Δ
x
=
Δ
y) are used for the discretization of the plate. Results of numerical simulations are compared with the results obtained from Abaqus, where the shell has been discretized using n
x
=
40, n
y
=
10 and n
x
=
60, n
y
=
15 elements.Frequencies of free vibrations obtained from the hybrid model (ω
H
,
i) are compared with results obtained for two models formulated in Abaqus (ω
A
_
L
M
,
i and ω
A
_
C
M
,
i) in Table 2. First six free vibrations are zero as for a rigid body and thus they are not shown in Figure 11.Table 2
Frequencies of free vibrations obtained in the hybrid model and Abaqus.
i
Hybrid model
Abaqus
Lumped mass formulation
Consistent mass formulation
n
x
=
40, n
y
=
10
ω
H
,
i [s−1]
ω
A
_
LM
,
i  ω
i [s−1]
ε
A
_
LM
,
i [%]
ω
A
_
CM
,
i [s−1]
ε
A
_
CM
,
i [%]
ε
A
,
i [%]
1,…,6
0
0
—
0
—
—
7
0.660
0.661
0.13
0.660
0.04
0.15
8
1.608
1.609
0.06
1.608
0.01
0.06
9
1.830
1.838
0.42
1.831
0.06
0.38
10
3.332
3.338
0.18
3.333
0.04
0.15
11
3.605
3.637
0.88
3.609
0.11
0.78
12
5.277
5.299
0.41
5.284
0.13
0.28
13
5.974
6.068
1.55
5.988
0.23
1.34
14
7.535
7.595
0.79
7.556
0.28
0.52
15
8.907
9.124
2.38
8.943
0.40
2.02
16
10.189
10.325
1.32
10.237
0.47
0.86Two types of shell elements, S4R of the first order and S8R of the second order, are used for discretisation in the lumped mass formulation (A
_
L
M) and consistent mass formulation (A
_
C
M), respectively. Quantity ε
A
_
,
i represents percentage error with respect to the Abaqus solutions:
(72)
ε
A
_
,
i
=
ω
A
_
,
i
-
ω
H
,
i
ω
A
_
,
i
·
100
%
,
where index
(73)
A
_
,
i
=
A
_
L
M
,
i
,
for
lumped
mass
formulation
,
A
_
C
M
,
i
,
for
consistent
mass
formulation
,
and quantity ε
A
,
i represents percentage error of the lumped mass formulation with respect to the consistent mass formulation in Abaqus solutions:
(74)
ε
A
,
i
=
ω
A
_
C
M
,
i
-
ω
A
_
L
M
,
i
ω
A
_
C
M
,
i
·
100
%
.Results obtained from the hybrid model differ from those obtained for theA
_
L
M model by no more than 2.4%. Much better compatibility is obtained if comparing the hybrid method with the consistent mass formulation in Abaqus, when the percentage error does not exceed 0.5%. It is important to note that relative error ε
A
,
i between the results obtained for A
_
C
M and A
_
L
M models in Abaqus is larger than 2%. Thus the method proposed enables us to obtain results closer to those from the more exact model A
_
C
M. Values presented in Table 2 prove that the method proposed gives reliable results comparable with those obtained from Abaqus.The influence of numbersn
x and n
y on exactness of the calculations of free vibrations is shown in Figure 12. It can be seen that acceptable exactness is obtained already for n
x
>
20 and n
y
>
6. In order to verify the hybrid method used for complicated shapes, the analysis of free vibrations of the SIGMA electrode was carried out and compared with that obtained using Abaqus and ANSYS packages as well as the strip method [13]. The relative error for the electrode of 10 m has been 7% in comparison with results obtained from ANSYS and larger in comparison with those obtained from Abaqus.Relative error with respect to the number of elements in discretization for (a) tenth (b) thirteenth, and (c) sixteenth frequency, respectively.
(a)
(b)
(c)Some calculations enabling us to analyse the influence of anglesα
j from Figure 2 on frequency of free vibrations of SIGMA plate are also carried out. Figure 13 presents a typical profile of the electrode.Figure 13
Profile of the SIGMA electrode.It is assumed that anglesα
j, which are multiples of 90° (horizontal and perpendicular segments), do not change. Only values of angles α
7
=
α
11
=
4
5
°, α
9
=
α
13
=
-
4
5
°, α
5
=
3
7
°, and α
15
=
-
3
7
° are changed for ±
30
%, while the width of the strips stays constant.Figure14 presents the comparison of frequencies of first 36 free vibrations of a steel plate with the length of 4 m and thickness of 0.0015 m. It can be seen that the change of angles α
j (0.7
α
j and 1.3
α
j) influences the frequencies of free vibrations. The more undulated the profile of the electrode (resulting from the change of angles α
j) the larger the values of frequencies of free vibrations of the SIGMA plate. This conclusion can be an important tip for engineers at the design stage.Figure 14
Free vibrations of the SIGMA plate, comparison for the first 36 frequencies of frr vibrations for different values ofα
j.


### 3.2. Vibration Analysis of the System

In this paper vibration analysis is carried out for a system of electrodes 16.152 m long, 0.0015 m thick, consisting of 19 strips and the following material parameters: Young modulusE
=
2.06
×
10
11 N/m2, and ρ
=
7850 kg/m3. The course of the force caused by the hammer and inducing the vibrations is presented in Figure 15.Figure 15
Force inducing the vibrations.The scheme of the control system (configuration 1) is presented in Figure16. The system consists of three collecting SIGMA type electrodes (E
1, E
5, and E
9) suspended on a common beam (S
B) and connected at the bottom with a brushing bar (A
B) with the anvil (A) at the end. This system has been used to study the influence of density of discretisation in x-axis direction. The respective computation time has been recorded. Table 3 shows the influence of number n
x, into which the plates are divided, on total calculation time.Table 3
The influence of numbern
x on calculation time.
n
x
200
300
400
500
600
700
800
900
1000
t [s]
275
543
855
919
1240
1453
1599
1996
2246Figure 16
Scheme of the control system (configuration 1).Figures17, 18, and 19 show the influence of number n
x, into which the electrodes are divided, on three components of accelerations at a chosen point of the electrode. Good convergence of the results can be observed.The influence of numbern
x on a
x component of acceleration.
(a)
(b)Figure 18
The influence of numbern
x on a
y component of acceleration.Figure 19
The influence of numbern
x on the a
z component of acceleration.The equations of motion have been integrated using the Newmark method with a constant integration step over time interval0,0.1. Bearing in mind the convergence of the results and a strong dependence of calculation time on number n
x, the results presented below have been obtained for n
x
=
200. Figure 20 presents total accelerations in the middle strip of the electrode at the largest distance from the impact point at 9 different points along the length of the electrode. The origin and the process of wave shifting can be observed.Figure 20
Tangent accelerations at different points.Further, we will compare simulation results with those obtained from experimental measurements. The measurements have been carried out at a special test stand (Figure21) built by a producer of electrostatic precipitators. The system of electrodes was considered during model validation.Figure 21
Test stand.The electrodes were 16 meters long and 1.5 millimeters thick. The configuration of checkpoints corresponds (Figure22) to the arrangement of acceleration sensors on the test stand. The equipment used consists of a 16-channel recorder TEAC lx110, a portable computer with LX Navi and FlexPro software, 5 triaxial vibration ICP sensors, and an ICP force sensor. The signals were recorded with a sampling rate of 24 kHz (per channel) as the response of the system for a single force impulse F
(
t
) applied to the anvil (Figure 15).Figure 22
Arrangements of the acceleration sensors on the electrodes (configuration 2).The authors’ experience shows that one may not directly compare time histories (courses) of accelerations obtained by means of different commercial software packages (ANSYS and Abaqus have been checked) or from numerical simulations and experimental measurements. Therefore, averaged (integral) values of accelerations calculated at points at which the sensors have been placed are used in order to validate the model presented:(75)
R
=
1
T
∫
0
T
a
s
2
d
t
1
/
2
,
where
(76)
a
s
=
a
x
2
+
a
y
2
,
for
s
=
τ
,
a
z
,
for
s
=
υ
.Correspondence of results has been evaluated by means of factors used in computer fluid mechanics [16, 17].The Hit Rate. Consider
(77)
q
=
1
n
p
∑
i
=
1
n
p
N
i
q
,
N
i
q
=
1
,
if
R
c
-
R
e
R
e
<
ε
,
0
,
otherwise
,
where R
c, R
e are calculated according to (75) (superscript c indicates accelerations calculated, while superscript e indicates those obtained from experimental measurements), ε is permissible error, and n
p is the number of points.Factor 2 (FAC2). Consider
(78)
f
=
1
n
p
∑
i
=
1
n
p
N
i
f
,
N
i
f
=
1
,
if
1
2
≤
R
c
R
e
<
2
,
0
,
otherwise
.Figure22 shows arrangements of the measurement points. Due to technical limitations the sensors have been placed below the level of 10 m above the anvil (measurement level 5). The accelerations have not been measured on the two electrodes closest to the beater because their values were larger than 5000 m/s2.Coefficientsf and q, in view of the chosen placement of the sensors, allow us to compare the numerical results and experimental measurements both for the whole system of n electrodes and locally, for example, for a single strip of the electrode. Values of coefficients f (FAC2) and q are presented in Figure 23.FAC2 and hit rateq for normal, tangent, and total accelerations.
(a)
(b)Since maximal values of accelerations are important factors for the producer of electrostatic precipitators, coefficientsq and f are also compared for value R
- defined as follows:
(79)
R
-
=
max
⁡
a
s
c
max
⁡
a
s
e
.The results are shown in Figure24.FAC2 and hit rateq for maximal values of normal, tangent, and total accelerations.
(a)
(b)In both figures hit rateq is presented for different values of ε from (77). In view of the results of validation, the model reflects the main features of the system and it will be of definite benefit for producers of electrostatic precipitators at the design stage in order to predict accelerations induced in the electrodes by the rapping system.


## 3.1. Free Vibrations

The analysis of free vibrations is carried out for a free rectangular plate with the following dimensions: widtha
=
1 m, length b
=
4 m, and thickness h
=
0.002 m, and the following material parameters: Young modulus E
=
2.06
×
10
11 N/m2, ρ
=
7850 kg/m3, and Poisson number ν
=
0.3. Figure 11 presents the first fifteen frequencies of free vibrations obtained from the hybrid model for a different proportion defined by a coefficient η according to the formula
(71)
η
=
Δ
y
Δ
x
,
where Δ
x
=
a
/
n
x and Δ
y
=
b
/
n
y.Figure 11
Influence of coefficientη on frequencies of free vibrations.It can be seen that each of the analyzed frequencies achieves the maximal value forη
=
1. Thus, for further analysis square elements (Δ
x
=
Δ
y) are used for the discretization of the plate. Results of numerical simulations are compared with the results obtained from Abaqus, where the shell has been discretized using n
x
=
40, n
y
=
10 and n
x
=
60, n
y
=
15 elements.Frequencies of free vibrations obtained from the hybrid model (ω
H
,
i) are compared with results obtained for two models formulated in Abaqus (ω
A
_
L
M
,
i and ω
A
_
C
M
,
i) in Table 2. First six free vibrations are zero as for a rigid body and thus they are not shown in Figure 11.Table 2
Frequencies of free vibrations obtained in the hybrid model and Abaqus.
i
Hybrid model
Abaqus
Lumped mass formulation
Consistent mass formulation
n
x
=
40, n
y
=
10
ω
H
,
i [s−1]
ω
A
_
LM
,
i  ω
i [s−1]
ε
A
_
LM
,
i [%]
ω
A
_
CM
,
i [s−1]
ε
A
_
CM
,
i [%]
ε
A
,
i [%]
1,…,6
0
0
—
0
—
—
7
0.660
0.661
0.13
0.660
0.04
0.15
8
1.608
1.609
0.06
1.608
0.01
0.06
9
1.830
1.838
0.42
1.831
0.06
0.38
10
3.332
3.338
0.18
3.333
0.04
0.15
11
3.605
3.637
0.88
3.609
0.11
0.78
12
5.277
5.299
0.41
5.284
0.13
0.28
13
5.974
6.068
1.55
5.988
0.23
1.34
14
7.535
7.595
0.79
7.556
0.28
0.52
15
8.907
9.124
2.38
8.943
0.40
2.02
16
10.189
10.325
1.32
10.237
0.47
0.86Two types of shell elements, S4R of the first order and S8R of the second order, are used for discretisation in the lumped mass formulation (A
_
L
M) and consistent mass formulation (A
_
C
M), respectively. Quantity ε
A
_
,
i represents percentage error with respect to the Abaqus solutions:
(72)
ε
A
_
,
i
=
ω
A
_
,
i
-
ω
H
,
i
ω
A
_
,
i
·
100
%
,
where index
(73)
A
_
,
i
=
A
_
L
M
,
i
,
for
lumped
mass
formulation
,
A
_
C
M
,
i
,
for
consistent
mass
formulation
,
and quantity ε
A
,
i represents percentage error of the lumped mass formulation with respect to the consistent mass formulation in Abaqus solutions:
(74)
ε
A
,
i
=
ω
A
_
C
M
,
i
-
ω
A
_
L
M
,
i
ω
A
_
C
M
,
i
·
100
%
.Results obtained from the hybrid model differ from those obtained for theA
_
L
M model by no more than 2.4%. Much better compatibility is obtained if comparing the hybrid method with the consistent mass formulation in Abaqus, when the percentage error does not exceed 0.5%. It is important to note that relative error ε
A
,
i between the results obtained for A
_
C
M and A
_
L
M models in Abaqus is larger than 2%. Thus the method proposed enables us to obtain results closer to those from the more exact model A
_
C
M. Values presented in Table 2 prove that the method proposed gives reliable results comparable with those obtained from Abaqus.The influence of numbersn
x and n
y on exactness of the calculations of free vibrations is shown in Figure 12. It can be seen that acceptable exactness is obtained already for n
x
>
20 and n
y
>
6. In order to verify the hybrid method used for complicated shapes, the analysis of free vibrations of the SIGMA electrode was carried out and compared with that obtained using Abaqus and ANSYS packages as well as the strip method [13]. The relative error for the electrode of 10 m has been 7% in comparison with results obtained from ANSYS and larger in comparison with those obtained from Abaqus.Relative error with respect to the number of elements in discretization for (a) tenth (b) thirteenth, and (c) sixteenth frequency, respectively.
(a)
(b)
(c)Some calculations enabling us to analyse the influence of anglesα
j from Figure 2 on frequency of free vibrations of SIGMA plate are also carried out. Figure 13 presents a typical profile of the electrode.Figure 13
Profile of the SIGMA electrode.It is assumed that anglesα
j, which are multiples of 90° (horizontal and perpendicular segments), do not change. Only values of angles α
7
=
α
11
=
4
5
°, α
9
=
α
13
=
-
4
5
°, α
5
=
3
7
°, and α
15
=
-
3
7
° are changed for ±
30
%, while the width of the strips stays constant.Figure14 presents the comparison of frequencies of first 36 free vibrations of a steel plate with the length of 4 m and thickness of 0.0015 m. It can be seen that the change of angles α
j (0.7
α
j and 1.3
α
j) influences the frequencies of free vibrations. The more undulated the profile of the electrode (resulting from the change of angles α
j) the larger the values of frequencies of free vibrations of the SIGMA plate. This conclusion can be an important tip for engineers at the design stage.Figure 14
Free vibrations of the SIGMA plate, comparison for the first 36 frequencies of frr vibrations for different values ofα
j.


## 3.2. Vibration Analysis of the System

In this paper vibration analysis is carried out for a system of electrodes 16.152 m long, 0.0015 m thick, consisting of 19 strips and the following material parameters: Young modulusE
=
2.06
×
10
11 N/m2, and ρ
=
7850 kg/m3. The course of the force caused by the hammer and inducing the vibrations is presented in Figure 15.Figure 15
Force inducing the vibrations.The scheme of the control system (configuration 1) is presented in Figure16. The system consists of three collecting SIGMA type electrodes (E
1, E
5, and E
9) suspended on a common beam (S
B) and connected at the bottom with a brushing bar (A
B) with the anvil (A) at the end. This system has been used to study the influence of density of discretisation in x-axis direction. The respective computation time has been recorded. Table 3 shows the influence of number n
x, into which the plates are divided, on total calculation time.Table 3
The influence of numbern
x on calculation time.
n
x
200
300
400
500
600
700
800
900
1000
t [s]
275
543
855
919
1240
1453
1599
1996
2246Figure 16
Scheme of the control system (configuration 1).Figures17, 18, and 19 show the influence of number n
x, into which the electrodes are divided, on three components of accelerations at a chosen point of the electrode. Good convergence of the results can be observed.The influence of numbern
x on a
x component of acceleration.
(a)
(b)Figure 18
The influence of numbern
x on a
y component of acceleration.Figure 19
The influence of numbern
x on the a
z component of acceleration.The equations of motion have been integrated using the Newmark method with a constant integration step over time interval0,0.1. Bearing in mind the convergence of the results and a strong dependence of calculation time on number n
x, the results presented below have been obtained for n
x
=
200. Figure 20 presents total accelerations in the middle strip of the electrode at the largest distance from the impact point at 9 different points along the length of the electrode. The origin and the process of wave shifting can be observed.Figure 20
Tangent accelerations at different points.Further, we will compare simulation results with those obtained from experimental measurements. The measurements have been carried out at a special test stand (Figure21) built by a producer of electrostatic precipitators. The system of electrodes was considered during model validation.Figure 21
Test stand.The electrodes were 16 meters long and 1.5 millimeters thick. The configuration of checkpoints corresponds (Figure22) to the arrangement of acceleration sensors on the test stand. The equipment used consists of a 16-channel recorder TEAC lx110, a portable computer with LX Navi and FlexPro software, 5 triaxial vibration ICP sensors, and an ICP force sensor. The signals were recorded with a sampling rate of 24 kHz (per channel) as the response of the system for a single force impulse F
(
t
) applied to the anvil (Figure 15).Figure 22
Arrangements of the acceleration sensors on the electrodes (configuration 2).The authors’ experience shows that one may not directly compare time histories (courses) of accelerations obtained by means of different commercial software packages (ANSYS and Abaqus have been checked) or from numerical simulations and experimental measurements. Therefore, averaged (integral) values of accelerations calculated at points at which the sensors have been placed are used in order to validate the model presented:(75)
R
=
1
T
∫
0
T
a
s
2
d
t
1
/
2
,
where
(76)
a
s
=
a
x
2
+
a
y
2
,
for
s
=
τ
,
a
z
,
for
s
=
υ
.Correspondence of results has been evaluated by means of factors used in computer fluid mechanics [16, 17].The Hit Rate. Consider
(77)
q
=
1
n
p
∑
i
=
1
n
p
N
i
q
,
N
i
q
=
1
,
if
R
c
-
R
e
R
e
<
ε
,
0
,
otherwise
,
where R
c, R
e are calculated according to (75) (superscript c indicates accelerations calculated, while superscript e indicates those obtained from experimental measurements), ε is permissible error, and n
p is the number of points.Factor 2 (FAC2). Consider
(78)
f
=
1
n
p
∑
i
=
1
n
p
N
i
f
,
N
i
f
=
1
,
if
1
2
≤
R
c
R
e
<
2
,
0
,
otherwise
.Figure22 shows arrangements of the measurement points. Due to technical limitations the sensors have been placed below the level of 10 m above the anvil (measurement level 5). The accelerations have not been measured on the two electrodes closest to the beater because their values were larger than 5000 m/s2.Coefficientsf and q, in view of the chosen placement of the sensors, allow us to compare the numerical results and experimental measurements both for the whole system of n electrodes and locally, for example, for a single strip of the electrode. Values of coefficients f (FAC2) and q are presented in Figure 23.FAC2 and hit rateq for normal, tangent, and total accelerations.
(a)
(b)Since maximal values of accelerations are important factors for the producer of electrostatic precipitators, coefficientsq and f are also compared for value R
- defined as follows:
(79)
R
-
=
max
⁡
a
s
c
max
⁡
a
s
e
.The results are shown in Figure24.FAC2 and hit rateq for maximal values of normal, tangent, and total accelerations.
(a)
(b)In both figures hit rateq is presented for different values of ε from (77). In view of the results of validation, the model reflects the main features of the system and it will be of definite benefit for producers of electrostatic precipitators at the design stage in order to predict accelerations induced in the electrodes by the rapping system.


## 4. Conclusions

Analysis of free vibrations of a single plate shows very good convergence of the method and compatibility of results (the relative error smaller than 0.5%) when compared with those obtained from the commercial software. Thus, the method has been used for modeling and analysis of a system containing thin plates of complicated shapes.In order to validate the model experimental measurements have been carried out on a special test stand consisting of a real system of collecting electrodes used in dry electrostatic precipitators. Consistency factors FAC2 andq have been used for comparative analysis of the results of calculations and measurements. Good correspondence of the results of accelerations of vibrations in three directions has been achieved.The method proposed is intuitive and easy in implementation of complicated systems in which the structure contains shells, beams, and concentrated masses. The authors compared numerical effectiveness of the method presented with some other methods, including commercial software in [18], and it has been shown that due to the shortest calculation time the software developed is an effective tool used by engineers at the design stage.Values of accelerations of vibrations exceed 100 g for the whole surface of the plates. Particularly large vibrations occur at the level of the impact force and accelerations transmitted by the brushing bar exceed 500 g. Since the effectiveness of the dust-removal process depends not only on the values of accelerations but also on their distribution, the analysis of different configurations of the system is very important, especially that the vibrations are quickly suppressed.The software elaborated enables us to analyze the influence of the position of the brushing bar; its present position does not ensure even distribution of vibrations.The overall conclusion is that the process of vibration excitation and wave propagation in the system of electrodes is the result of many factors. This process depends not only on the impact force but also on the physical parameters, geometry, and construction of all the elements that make up this system. In this respect the model presented in this paper is an important novelty since using the testing calculations can predict properties of the future structure as early as at the stage of its design.The possibility of analysis of vibrations for different profiles of the collecting plates with various geometrical parameters makes the software useful for the design process.



---

*Source: 832918-2014-10-30.xml*