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.

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 [25].

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 (Figure 1). 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.

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 axis , the spring energy of element is independent of and angle .

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 [1012]. 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 Figure 1, 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 (Figure 2), which can be considered as a set of long strips with different angles of inclination towards axis . 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.

2.1.1. Primary Division

Let us assign a coordinate system as in Figure 2. Discretization is carried out along axis into strips and along axis into segments with the length where is the length of the electrode. Thus the electrode is divided into elements.

Let us consider a rectangular element as in Figure 3. The following vectors describe displacements of nodes for such an element: where , , are displacements of node along , , , respectively; , , are respective rotations in node about axes parallel to , , .

Shape functions are defined as

It is assumed that , , and describe shield deformations, while , , and describe plate deformations. The angles are defined by the following relations:

The elastic strain energy is a sum of energies corresponding to the shield and plate states and can be calculated as [13] where , is the thickness of the element, are vectors of strains, are vectors of stresses, 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: where

The dimension of matrices and is . Integrals occurring in elements of matrices can be calculated analytically or by means of Gauss quadratures. Detailed formulae for elements of matrices and are presented in [14]. If we assume that the form of these matrices is as follows: 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: where is a stiffness matrix with elements, and nonzero elements of matrix C from (13) are defined as follows:

For the assumed values of parameters , the stiffness matrices depend only on and . Having used denotations from Figures 2 and 3 and relation (1), in the case considered, we have for ,  . Thus, as a result of constant length of elements, stiffness matrices for each strip are identical and are calculated only for .

The coordinates of the geometrical center of element (, ) can be calculated according to the formulas where , define coordinates of the left edge of strip number 1 and , are the width and angle of inclination of the th strip to axis , respectively.

Coordinates of nodes of primary elements can be calculated as follows: where is defined in (16), and is the vector of coordinates of the mass centre of the th segment in the local coordinate system:

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 (, ) from the primary division, as in Figure 4.

Coordinate system , , is a global coordinate system. Primary element (, ) 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: for , , and .

Figure 5 shows an example of primary and secondary division of a rectangular area for and .

The number of rigid finite elements is

It can be seen from Figures 4 and 5 that rfes are constructed by rigid connections of stiffened subsegments in the following way.(i)Corner elements:rfe 1 is subsegment stiffened;rfe is subsegment stiffened;rfe is subsegment stiffened;rfe is subsegment stiffened.(ii) Other boundary elements:for ,rfe is formed by connecting and stiffening subsegments and ;rfe is formed by connecting and stiffening subsegments and ;for ,rfe is formed by connecting and stiffening subsegments and ;rfe is formed by connecting and stiffening subsegments and .(iii)Internal elements:for , ,rfe is formed by connecting and stiffening subsegments ,, , and .

Secondary division of the electrode into rfes is essential for the method discussed since in further considerations the shell is treated as a system of rigid bodies with 6 degrees of freedom. A local coordinate system , whose axes are central inertial axes of the element, is assigned to each rfe as in Figure 6. Angle defines inclination of the plane of element with respect to plane of coordinate system assigned to the electrode.

The vector of the generalized coordinates of rfe includes three displacements and three rotations: where and .

When the vibrations of the collecting electrodes are considered, it can be assumed that angles , , are small. Thus, when coordinates , , of a point are known in the local system of the th rfe, its coordinates in the global coordinate system can be calculated as follows: where is the coordinate vector in and is defined in Figure 6, as the coordinate vector of point for a nondeformed rfe (when ),

2.1.3. Equations of Motion of the SIGMA Electrode

Using the standard RFEM approach, the diagonal mass matrix of the th rfe can be formulated as follows:

In the case of SIGMA electrodes it is relatively easy to formulate the algorithm to calculate elements of matrix using the above formula and angles 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 and the kinetic energy of the whole electrode as where and .

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 are defined in (12). Individual segments of primary element belong to various rfes, whose indexes are determined in (20).

It results from (20) that node belongs to rfe:

Coordinates in the global coordinate system of the electrode   are defined by (17). Having taken into account (23), node coordinates in local system can be defined as follows: where , defined in (28). If rfe has coordinates defined in (22), then coordinates of node in the global system can be calculated from (23):

In order to define the energy of spring deformation of node the node displacements and rotations have to be defined by means of generalized coordinates .

The displacement vector of node belonging to rfe , as in (28), is defined as follows:

After some transformations, when (30), (29), and (17) are taken into account, the following is obtained:

Rotation angles can be calculated from the formula where Finally, the following relation is derived:

Matrix 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 can be written in the following form where is defined in (12) and (13).

Bearing in mind that , and taking into account (36), the above can be written as follows: where .

This leads to a straightforward method of generating stiffness matrix using coordinates as generalized coordinates. The appropriate algorithm is as follows:(1)for , accept: .(2)for ; (i)for compute according to (36)(ii)for ; (a)compute: (b)for ; compute:

It is important to note that matrix 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: where is the acceleration of gravity, , and . Bearing in mind (23) we can calculate derivatives of the potential energy

Since the origin of the local coordinate system of rfe is located at the center of mass of the rfe, that is, , the derivative of potential energy of gravity forces for the th strip of the electrode can be presented in the matrix form where .

Equations of motion of a free electrode can be written in the form and the eigenvalue problem takes the form

Some important differences between the hybrid finite element method and the finite element method applied to modeling of electrodes are listed in Table 1.

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.

The number of elements into which the beams are discretized may be different for both beams. If is the number of electrodes and 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: where , ; , denote top or bottom beam, respectively; and .

The secondary division leads to a system of rigid finite elements with the length equal to

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: where is the mass of the th rfe, is the cross-section area, are inertial moments of the cross-section of the beam, ,  , and . The vector of generalized coordinates of rfe is as follows: where , , are displacements of the center of mass of the rfe and , , are the angles of rotations. The kinetic energy of the th rfe is calculated according to the formula

The energy of spring deformation can be presented in the form where , are the stiffness coefficients calculated according to the formulas given in [9], and is the deformation of spring damping element (sde) , which depends on and .

The potential energy of gravity forces of each rfe can be calculated according to the formula where and is defined in (41).

Kinetic and potential energies of the beams are the sum of energies of all rfes: and the potential energy of spring deformation is the sum of energies of sdes 0 to :

The vector of generalized coordinates of a beam can be presented as

So finally, the following can be written: where is the diagonal mass matrix, is the stiffness matrix with constant coefficients, and is the vector of gravity forces, . Thus, the motion of beams is described by the following number of generalized coordinates:

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 force , and the resulting generalized force follows

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 point by means of spring-damping element (sde) as shown in Figure 8.

We assume the following denotations: denotes the system before deformation and denotes the system rigidly assigned to rfe , where .

Owing to (22) and (49), the motion of rfes is governed by the vectors where , , , , describe the position of the beginning of system in , and , , describe the orientation of system with respect to . After deformation, rfes and move, and thus the position of point in the global coordinate system may vary depending on whether point is considered as a point of rfe or : where is the position vector of the beginning of in the global system , is the vector of coordinates of point in system , is the rotation matrix of direction cosines of system with respect to , when the rotation angles are small, and are coordinates of a point in system rigidly connected with rfe .

Vectors of translational and rotational displacements of sde with respect to the global system can be calculated as follows: where

After transformations the following can be obtained: where It is important to note that the matrices , , , and are constant.

Finally, the energy of spring deformation of the connection can be presented in the following form: where . It can be seen from (64) that the energy of spring deformation of the connection depends on generalized coordinates of rfes and . Thus the following components will appear in the equations of motion where

It should be stressed that the element presented is universal. When stiffness coefficients in matrix 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 Figure 9) and the coefficients of stiffness have been assumed to be N/m for translational displacements and as  N/rad in the case of rotary displacements.

2.4. Equations of Motion

The equations of motion for the system of the beams and electrodes can be written in the following form: where is the number of active electrodes, , , , , , are defined in Section 2.2,   , , are defined in Section 2.1,    denotes the number of the electrode, , and is defined in (58).

It has been assumed that the electrodes are connected, as in Figure 10, 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)).

Finally, the equations of motion after necessary calculations take the following form:

The coupling resulting from the connections between the electrodes is reflected by new stiffness matrices and , while , , and are modified matrices , , and 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: width  m, length  m, and thickness  m, and the following material parameters: Young modulus  N/m2,  kg/m3, and Poisson number . 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 where and .

It can be seen that each of the analyzed frequencies achieves the maximal value for . Thus, for further analysis square elements () 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 , and , elements.

Frequencies of free vibrations obtained from the hybrid model () are compared with results obtained for two models formulated in Abaqus ( and ) in Table 2. First six free vibrations are zero as for a rigid body and thus they are not shown in Figure 11.

Two types of shell elements, S4R of the first order and S8R of the second order, are used for discretisation in the lumped mass formulation () and consistent mass formulation (), respectively. Quantity represents percentage error with respect to the Abaqus solutions: where index and quantity represents percentage error of the lumped mass formulation with respect to the consistent mass formulation in Abaqus solutions:

Results obtained from the hybrid model differ from those obtained for the 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 between the results obtained for and models in Abaqus is larger than 2%. Thus the method proposed enables us to obtain results closer to those from the more exact model . Values presented in Table 2 prove that the method proposed gives reliable results comparable with those obtained from Abaqus.

The influence of numbers and on exactness of the calculations of free vibrations is shown in Figure 12. It can be seen that acceptable exactness is obtained already for and . 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.

Some calculations enabling us to analyse the influence of angles from Figure 2 on frequency of free vibrations of SIGMA plate are also carried out. Figure 13 presents a typical profile of the electrode.

It is assumed that angles , which are multiples of 90° (horizontal and perpendicular segments), do not change. Only values of angles , , , and are changed for , while the width of the strips stays constant.

Figure 14 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 ( and ) influences the frequencies of free vibrations. The more undulated the profile of the electrode (resulting from the change of angles ) 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.

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 modulus  N/m2, and  kg/m3. The course of the force caused by the hammer and inducing the vibrations is presented in Figure 15.

The scheme of the control system (configuration 1) is presented in Figure 16. The system consists of three collecting SIGMA type electrodes (, , and ) suspended on a common beam () and connected at the bottom with a brushing bar () with the anvil () at the end. This system has been used to study the influence of density of discretisation in -axis direction. The respective computation time has been recorded. Table 3 shows the influence of number , into which the plates are divided, on total calculation time.

Figures 17, 18, and 19 show the influence of number , 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 equations of motion have been integrated using the Newmark method with a constant integration step over time interval . Bearing in mind the convergence of the results and a strong dependence of calculation time on number , the results presented below have been obtained for . 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.

Further, we will compare simulation results with those obtained from experimental measurements. The measurements have been carried out at a special test stand (Figure 21) built by a producer of electrostatic precipitators. The system of electrodes was considered during model validation.

The electrodes were 16 meters long and 1.5 millimeters thick. The configuration of checkpoints corresponds (Figure 22) 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 applied to the anvil (Figure 15).

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: where

Correspondence of results has been evaluated by means of factors used in computer fluid mechanics [16, 17].

The Hit Rate. Consider where , are calculated according to (75) (superscript indicates accelerations calculated, while superscript indicates those obtained from experimental measurements), is permissible error, and is the number of points.

Factor 2 (FAC2). Consider

Figure 22 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.

Coefficients and , 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 electrodes and locally, for example, for a single strip of the electrode. Values of coefficients (FAC2) and are presented in Figure 23.

Since maximal values of accelerations are important factors for the producer of electrostatic precipitators, coefficients and are also compared for value defined as follows:

The results are shown in Figure 24.

In both figures hit rate 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 and 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.

List of Symbols

FEM:Finite element method
RFEM:Rigid finite element method
est:Spring-damping element
rfe:Rigid finite element
:Factors
Matrices and vectors are marked by bold letters
:Transposition of matrix or vector
:Coordinates in global coordinate system
:Coordinate in local coordinate system
:Global coordinate system
:Local coordinate system
:Vectors of generalised coordinates with respect to coordinate systems and
:Strain and stress vectors
:Vector of generalised forces
:Stiffness matrices: global and of element
:Mass matrices: global and of element
:Kinetic energy
:Potential energy
:Potential energy of gravity forces
:Potential energy of spring deformation
:Translational displacements for axes, respectively
:Rotational displacements about axes, respectively.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This paper was financially supported by the National Science Centre of Poland under the Grant MAESTRO 2, no. 2012/04/A/ST8/00738, for years 2013–2016.