Development of a perfect match system in the improvement of eigenfrequencies of free vibration

ARTICLE IN PRESS

JID: APM

[m3Gsc;March 1, 2017;1:55]

Applied Mathematical Modelling 0 0 0 (2017) 1–26

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Development of a perfect match system in the improvement of eigenfrequencies of free vibration Eric Li a,b,c, Z.C. He a,d,∗ a

State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body, Hunan University, Changsha 410082, PR China Department of Mechanical and Automation Engineering, The Chinese University of Hong Kong, Shatin, New Territories, Hong Kong, China c The State Key Laboratory of Fluid Power and Mechatronic Systems, Zhejiang University, Hangzhou 310027, China d State Key Laboratory of Mechanical System and Vibration, Shanghai Jiao Tong University, Shanghai 200240, PR China b

a r t i c l e

i n f o

Article history: Received 20 January 2016 Revised 3 January 2017 Accepted 9 February 2017 Available online xxx Keywords: Mass redistributed finite element method Finite element method Eigenfrequency Smoothed finite element method

a b s t r a c t It is well known that the standard finite element method (FEM) with overly-stiff effect gives the upper bound solutions of natural frequencies in the free vibration analysis using triangular and tetrahedral elements. In this study, for the first time, this paper aims to improve the prediction of eigenfrequencies through the perfect match between the stiffness and mass matrices. With redistribution of mass in the system, we can tune the balance between stiffness and mass of a discrete model. This can be done by simply shifting the integration points away from the Gaussian locations, while ensuring the mass conservation. A number of numerical examples including 2D and 3D problems have demonstrated that the accuracy of eigenfrequencies is strongly determined by the location of integration points in the mass matrix. With appropriate selection of integration points in the mass matrix, even the exact solution of eigenfrequencies can be obtained in both FEM and smoothed finite element method (SFEM) models. © 2017 Elsevier Inc. All rights reserved.

1. Introduction In the design of dynamic structure, vibration analysis plays a very important role. The prediction of modal frequencies of a structure is very useful to characterize the resonant vibration in machinery and structures [1]. As the analytical solutions for natural frequencies are only available for some simple geometries [2], the finite element method (FEM) has been widely used in the modal analysis [3–5]. It is well known that the displacement-based fully compatible FEM gives the upper bound solutions of eigenfrequencies [6]. This is because the stiffness in the FEM model is stiffer than the exact model [3,4,7]. In order to improve the accuracy of FEM solutions, many efforts have been made. Marangoni et al. [8] developed upper and lower bounds to the natural frequencies of vibration of clamped rectangular orthotropic pates. Wiberg et al. [9] presented the local and global updating techniques to improve FEM solution of the generalized eigenvalue problem in free vibration analysis. An error estimate and an adaptive procedure for generalized linear eigenvalue and eigenvector computation within a framework of the hierarchical finite element method were proposed by Friberg et al. [10]. Mackie [11] modified the standard FEM using the results of dispersion analysis of numerical approximation to improve the accuracy of eignefrequencies of high mode shape. ∗

Corresponding author at: State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body, Hunan University, Changsha 410082, PR China. E-mail address: [email protected] (Z.C. He).

http://dx.doi.org/10.1016/j.apm.2017.02.013 0307-904X/© 2017 Elsevier Inc. All rights reserved.

Please cite this article as: E. Li, Z.C. He, Development of a perfect match system in the improvement of eigenfrequencies of free vibration, Applied Mathematical Modelling (2017), http://dx.doi.org/10.1016/j.apm.2017.02.013

ARTICLE IN PRESS

JID: APM 2

[m3Gsc;March 1, 2017;1:55]

E. Li, Z.C. He / Applied Mathematical Modelling 000 (2017) 1–26

In order to soften the overly-stiff FEM model, Liu and co-workers have developed a G space theory [12] and weakened weak form, leading to the foundation of a family of smoothed finite element method (SFEM) [12–15]. By employing the node-based smoothed technique, node-based smoothed finite element method (NS-FEM) is developed with overly-soft property, which always gives the lower bound solutions of eignefrequencies [16,17]. However, NS-FEM can produce the socalled spurious modes with excessive node-based smoothing technique [16,17]. Due to this reason, an edge-based smoothed finite element method (ES-FEM) is developed in order to achieve a close to exact stiffness [7,18–21]. The ES-FEM gives much more accurate solutions in the prediction of eigenfrequencies compared with overly-soft NS-FEM and overly-stiff FEM [17]. The above efforts to reduce the softening effect of discretized mode are all focused on the modification of stiffness. Recently, Guddati [22,23] proposed a new modified integration rules in the computation of the mass and stiffness for acoustic problems with quadrilateral mesh. Following this, He et al. [24] proposed a mass redistributed finite element method (MR-FEM) using triangular elements in the analysis of acoustic problem in order to balance the discretized model between acoustic mass and stiffness. The implementation of MR-FEM is modified by shifting the integration points away from the usual Gaussian locations in the mass matrix [25]. Motivated by the excellent features of MR-FEM in the acoustic field, for the first time this work analyses the effect of mass redistribution in the prediction of eigenfrequencies for different discretized models including FEM and SFEM. Due to complicated geometry in general engineering problems, the high quality elements such as quadrilateral and brick elements are very difficult (or impossible) to generate. Therefore, the triangular or tetrahedral elements are considered particularly popular in industry. In this work, our study shows that the perfect match of discretized model between elastic mass and stiffness can be achieved with the modification of integration point in the mass matrix, which can improve the computational efficiency of eigenfrequencies in both FEM and SFEM models significantly. In addition, the relationship between the eigenfrequency and integration point t of mass matrix has been derived in this work. More importantly, the numerical results have indicated that the upper and lower bound solutions from NS-FEM and FEM models are furthered improved using coarse mesh. This significant finding is very useful to design the mechanical part or structure to ensure that the resonance leading large oscillation and damage never appears. The paper is organized as follows: Section 2 briefly describes the balanced systems in the discretized model. The formulation of MR-FEM in the computation of eigenfrequencies is shown in Section 3. Numerical examples including 2D and 3D are presented in Section 4 to investigate the effect of integration point of mass matrix on eigenfrequencies using FEM, NS-FEM and ES-FEM with combination of MR-FEM. Finally the conclusions from the numerical results are made in Section 5. 2. The balanced system of discretized model In the prediction of eigenfrequencies of elasticity systems, the global stiffness matrix and mass matrix are system matrices, and the balance between them is critical to yield accurate results. In the standard Galerkin FEM system, the weak form of FEM leads to the differences in “stiffness” between the numerical model and the exact system, and an imbalance between the stiffness matrix and mass matrix emerges. Hence, the eigenfrequencies in the FEM system are always larger than the real ones [6]. Consider the general elasticity problems without boundary conditions, as the stiffness and mass are positive definite matrices using discretized methods, there always exist an eigenvector φ

φ T Kφ = ω 2 φ T Mφ ,

(1)

where φT K φ and φT Mφ are diagonal matrices, and φ is the modal matrix.

⎡ λk1 ⎢ φ T Kφ = ⎢ ⎣

⎤

⎡ λm1 ⎥ ⎢ ⎥, φT Mφ = ⎢ ⎦ ⎣

λk2 ..

.

λkn

λm2 ..

⎤ ⎥ ⎥, ⎦

.

(2)

λmn

where λk1 ∼ λkn are the eigenvalues of stiffness matrix K, λm1 ∼ λmn are the eigenvalues of mass matrix, and the eigenvalues of elasticity are obtained as follows:

⎡

λk1 ⎢ λm1 ⎢ ⎢ T φ K φ ⎢ ω2 I = T =⎢ φ Mφ ⎢ ⎢ ⎣

⎤

λk2 λm2

..

.

⎥ ⎥ ⎥ ⎥ ⎥. ⎥ ⎥ λkn ⎦ λmn

(3)

From Eq. (3), it can be found that the eigenvalues are determined by the stiffness and mass matrices. In order to improve the prediction of eigenfrequencies of free vibrations, one way is to soften the stiffness of the FEM. With a generalized gradient smoothing technique, smoothed finite element method (SFEM) has demonstrated that the stiffness of the discretized model is reduced compared with the model of FEM. In this work, the mass matrix is modified using MR-FEM to improve Please cite this article as: E. Li, Z.C. He, Development of a perfect match system in the improvement of eigenfrequencies of free vibration, Applied Mathematical Modelling (2017), http://dx.doi.org/10.1016/j.apm.2017.02.013

ARTICLE IN PRESS

JID: APM

[m3Gsc;March 1, 2017;1:55]

E. Li, Z.C. He / Applied Mathematical Modelling 000 (2017) 1–26

3

Fig. 1. Illustration of two types of S-FEM.

y

C

c

a

b

A

x

(a) Parent mesh

B (b) Regular mesh

Fig. 2. Transformation of triangular mesh.

its eigenvalues to match properly with the current stiffness matrix. The theoretical analysis is studied to verify how the mass balances FEM and SFEM models with adjustment of integration points in the mass system. Two typical SFEM models are discussed in this study: Edge-based smoothed finite element method (ES-FEM) [7,19] and node-based smoothed finite element method (NS-FEM) [16,26]. The formulations of ES-FEM and NS-FEM are quite similar to the standard FEM. By introducing different smoothing domain based on edges of elements [7,27,28] or node of elements [16,29–31] as shown in Fig. 1, the compatible strain B matrix is replaced by the smoothed strain B. Hence, the stiffness in Eq. (3) using S-FEM can be formulated as follows:



KES−FEM =

T



 KNS−FEM =

∇ BES D∇ BES d The ES − FEM stiffness matrix, T



∇ BNS D∇ BNS d The NS − FEM stiffness matrix.

(4) (5)

In general, the softening effect of discretized system is proportional to the smoothing domain. With excessive node-based smoothing domain, NS-FEM behaves overly-soft, leading to the lower bound solution of eigenfrequency. While FEM is found overly-stiff, it provides upper solutions of eigenfrequency using consistent mass. In this work, we introduced the flexible integration point in the mass to improve the predication of eigenfrequency in the ES-FEM, NS-FEM and FEM models. 3. Theoretical analysis of eigenfrequency using MR-FEM 3.1. 2D eigenfrequency using MR-FEM In the 2D solid mechanics using linear triangular element, it is very general and customary to map the finite elements from a parent triangular mesh to the regular triangular mesh by area coordinates form, as shown in Fig. 2 and thus the integration of stiffness and mass matrices can be transformed into the natural coordinate system as follows:



Ke =

+1 −1



1 −ξ

−1

BT DB det (J )dξ dη, Me = ρ



+1 −1



1 −ξ

−1

NT N det (J )dξ dη,

(6)

Please cite this article as: E. Li, Z.C. He, Development of a perfect match system in the improvement of eigenfrequencies of free vibration, Applied Mathematical Modelling (2017), http://dx.doi.org/10.1016/j.apm.2017.02.013

ARTICLE IN PRESS

JID: APM 4

[m3Gsc;March 1, 2017;1:55]

E. Li, Z.C. He / Applied Mathematical Modelling 000 (2017) 1–26

Fig. 3. Area coordinates for a linear triangular element.

where ξ , η are local coordinates, D is the matrix of materials constant, J is Jacobian. Ni (x) are the linearly shape functions, which can be expressed in the area coordinates form, as shown in Fig. 3. By connecting point P with the three vertices a, b and c, three sub-areas A1, A2, and A3 corresponding to the triangles Pab, Pbc, and Pac respectively can be determined for any given point p within the element. The shape functions for each node are then defined as

N1 =

A1 = 1 − ξ − η; A

N2 =

A2 = ξ; A

N3 =

A3 = η. A

(7)

The above integrals are traditionally evaluated using Gauss quadrature



+1



−1

+1 −1

φ dξ dη =

ng

wi φ (ξi , ηi ),

(8)

i=1

where ng is the number of Gauss points, wi are the weights. Then the generalized stiffness and mass matrices can be obtained by the integration rule.

Ke = wi (∇ N(ξi , ηi ) ) D(∇ N(ξi , ηi ) ) det (J ) T

Me =

ρ

3



wki N

ξik , ηik

T

N

 ξik , ηik det (J ).

(9)

i=1

It can be seen from Eq. (9) that the gradient of shape function of triangular mesh is constant within the elements and thus the integration of stiffness matrix in Eq. (9) equals to the constant multiplied by the Jacobian, and no numerical integration is needed.

⎡

1 ⎢ 2 (1 + v )

⎢ ⎢ 0 ⎢ ⎢ ⎢ −1 ⎢ ⎢ 2 (1 + v ) e K = E⎢ ⎢ −1 ⎢ ⎢ 2 (1 + v ) ⎢ ⎢ ⎢ 0 ⎢ ⎣ 1 2 (1 + v )

0 −1 −1 + v2

v −1 + v2 1 −1 + v2

v −1 + v2 0

−1 2 (1 + v )

v −1 + v2 v−3





2 −1 + v2 −1 2(−1 + v ) 1 −1 + v2 −1 2 (1 + v )

−1 2 (1 + v ) 1 −1 + v2 −1 2(−1 + v ) v−3



2 −1 + v2



v

−1 + v2 −1 2 (1 + v )

−1 ⎤ 2 (1 + v ) ⎥

0

v −1 + v2 1 −1 + v2

v −1 + v2 −1 −1 + v2

⎥ ⎥ ⎥ ⎥ −1 ⎥ ⎥ 2 (1 + v ) ⎥ ⎥, −1 ⎥ ⎥ 2 (1 + v ) ⎥ ⎥ ⎥ ⎥ 0 ⎥ ⎦ 0

(10)

1 2 (1 + v )

0

where E and v are the Young’s modulus and Poisson’s ratio respectively. The eigenvalues of stiffness matrix Ke using single triangular element can be obtained as follows:

λk1 = λk2 = λk3 = 0; λk4 =

2E 1+v

λk5 = E

√ −2 + 1 + 3v2 −1 + v2

λk6 = E

√ −2 − 1 + 3v2 . −1 + v2

(11)

On the other hand, the consistent mass matrix is evaluated by three following Gauss points:

Point

α :ξi = 1/6, ηi = 1/6; Point β : ξi = 2/3, ηi = 1/6; Point γ : ξi = 1/6, ηi = 2/3.

(12)

The lumped mass matrix is achieved by the following three integration points.

Point α :ξi = 0, ηi = 0; Point

β : ξi = 1, ηi = 0; Point γ : ξi = 0, ηi = 1.

(13)

Please cite this article as: E. Li, Z.C. He, Development of a perfect match system in the improvement of eigenfrequencies of free vibration, Applied Mathematical Modelling (2017), http://dx.doi.org/10.1016/j.apm.2017.02.013

ARTICLE IN PRESS

JID: APM

[m3Gsc;March 1, 2017;1:55]

E. Li, Z.C. He / Applied Mathematical Modelling 000 (2017) 1–26

5

Fig. 4. Adjustable integration points for the mass matrix of linear triangular element.

As the mass matrix is influenced by the location of integration points, the following new locations of integration points for the general mass matrix are proposed:

Point α :

ξi =

1−t , 2

ηi =

1−t ; 2

Point β :

ξi = t, ηi =

1−t ; 2

Point γ :

ξi =

1−t , 2

ηi = t ,

(14)

where t ∈ [0, 1], and weight wi = 1/3 so that 3i=1 wi = 1. Note that when t is varied from 0 to 1, the integration point will move from a triangular vertex to midpoint of the edge opposite to the vertex, as shown in Fig. 4. Using these three sampling points, the mass matrix can be expressed in a parameter t as follows:

⎡

a11 ⎢0 ⎢a31 1 Me = ρ Ae ⎢ ⎢0 3 ⎣ a51 0

0 a22 0 a42 0 a62

a13 0 a33 0 a53 0

0 a24 0 a44 0 a64

a15 0 a35 0 a55 0

a11 = a22 = a33 = a44 = a55 = a66 = t 2 +

⎤

0 a26 ⎥ 0 ⎥ ⎥, a46 ⎥ ⎦ 0 a66

(15)

(1 − t )2 2

a13 = a15 = a24 = a26 = a31 = a35 = a42 = a46 = a51 = a53 = a62 = a64 = t (1 − t ) +

(1 − t )2 4

.

(16)

It can be easily confirmed that the summation of all matrix elements of Me equals to the area of the element regardless of value of t: 36 1 e e ρ A Mi = 3 i=1

2ρ Ae

 

.

(17)

total area of the element

It has been proven [27] that the total mass will always be conserved, if (1) the shape functions satisfy the partitions of unity property; and (2) the sum of the weights at all the sampling points is equal to 1. Although the formulation of the mass matrix will be influenced by the location of the integration points, the mass conservation is always guaranteed, as long as these two conditions are satisfied. This gives us the flexibility to alter the location of the sampling points in the mass matrix to improve the prediction of eigenfrequency of discretized system. After obtaining the general mass matrix for 2D elasticity problems using triangular mesh, the eigenvalues of mass system in the discretized system can be derived as follows:

1 6

λm1 = λm2 = h2 ρ , λm3 = λm4 = λm5 = λm6 =

 h2 ρ 2 9t − 6t + 1 , 24

(18)

where h is the size of element. Based on Eq. (11), it is noted that the eigenvalues for the stiffness are constant. Thus, the eigenvalues for the whole system are only controlled by the mass system. Apparently, Eq. (18) indicates that the location of integration point t controls eigenvalues of mass system by the term (9t2 − 6t + 1). Fig. 5 plots that the term (9t2 − 6t + 1)varies with the integration point t. As outlined in Fig. 5, it can be found that when t equals to 1/3, the largest eigenvalue is obtained. As t > 1/3, the eigenvalue is inversely proportional to the location of integration point t. While 0 < t < 1/3, with the increase of t, the eigenvalue becomes larger and larger. 3.2. 3D eigenfrequency using MR-FEM In the formulation of 3D MR-FEM using linear tetrahedral element, the mapping from a parent tetrahedral mesh to the regular tetrahedral mesh is based on volume coordinates. Thereby, the integration of stiffness and mass matrix is Please cite this article as: E. Li, Z.C. He, Development of a perfect match system in the improvement of eigenfrequencies of free vibration, Applied Mathematical Modelling (2017), http://dx.doi.org/10.1016/j.apm.2017.02.013

ARTICLE IN PRESS

JID: APM 6

[m3Gsc;March 1, 2017;1:55]

E. Li, Z.C. He / Applied Mathematical Modelling 000 (2017) 1–26

Fig. 5. Effect of integration points on the eigenvalue using single triangle element.

Fig. 6. Transformation of tetrahedral mesh.

transformed into the natural coordinate system as follows:

 Ke =

+1 −1



+1 −1



+1

−1



BT DB det (J )dξ dηdζ , Me = ρ



+1 −1



+1

−1



+1 −1

NT N det (J )dξ dηdζ ,

(19)

where Ni (x) are the linear shape functions, which can also be expressed in the volume coordinates form, as shown in Fig. 6. By connecting point P with the four vertices a, b, c and d, four sub-areas VPabc , VPabd , VPbcd and VPacd can be determined for any given point P within the element. Thus, the shape functions for each node are then defined as follows:

N1 =

VPbcd = 1 − ξ − η − ζ; Vabcd

N2 =

VPacd = ξ; Vabcd

N3 =

VPabd = η; Vabcd

N4 =

VPabc = ζ. Vabcd

(20)

Please cite this article as: E. Li, Z.C. He, Development of a perfect match system in the improvement of eigenfrequencies of free vibration, Applied Mathematical Modelling (2017), http://dx.doi.org/10.1016/j.apm.2017.02.013

ARTICLE IN PRESS

JID: APM

[m3Gsc;March 1, 2017;1:55]

E. Li, Z.C. He / Applied Mathematical Modelling 000 (2017) 1–26

7

With single tetrahedral element, the stiffness is expressed as follows:

⎡a

11

⎢0 ⎢0 ⎢ ⎢0 ⎢0 ⎢ ⎢a61 K = E⎢ ⎢a71 ⎢ ⎢0 ⎢a ⎢ 91 ⎢0 ⎣ 0 0

0 a22 0 0 0 0 0 a82 a92 0 0 a122

0 0 a33 a43 0 0 a73 a83 a93 0 a113 0

0 0 a34 a44 0 0 a74 a84 a94 0 a114 0

0 0 0 0 a55 0 a75 a85 0 a105 0 0

a16 0 0 0 0 a66 a76 0 a96 0 0 0

a17 0 a37 a47 a57 a67 a77 a87 a97 a107 a117 0

0 a28 a38 a48 a58 0 a78 a88 a98 a108 a118 a128

a19 a29 a39 a49 0 a69 a79 a89 a99 0 a119 a129

0 0 0 0 a510 0 a710 a810 0 a1010 0 0

0 0 a311 a411 0 0 a711 a811 a911 0 a1111 0

⎤

0 a212 ⎥ 0 ⎥ ⎥ 0 ⎥ ⎥ 0 ⎥ 0 ⎥ ⎥, 0 ⎥ ⎥ a812 ⎥ a912 ⎥ ⎥ 0 ⎥ ⎦ 0 a1212

(21)

where

1 2(1+v ) = a76 = a69 = a96

a11 = a22 = a55 = a66 = a1010 = a1212 = a16 = a61 = a212 = a122 = a510 = a105 = a17 = a17 = a19 = a19 = a82 = a28 = a29 = a92 = a57 = a75 = a58 = a85 = a67 −1 = a710 = a107 = a810 = a108 = a912 = a129 = 2 (1 + v ) (−1 + v) a33 = a44 = a11 = , (1+v )( − 1 + 2v ) −v a34 = a43 = a311 = a113 = a34 = a411 = a114 = (1+v )( − 1 + 2v )

v (1+v )( − 1 + 2v ) −(−1 + v ) (−2 + 3v) = a88 = a99 = , a = a93 = a47 = a74 = a811 = a118 = (1+v )( − 1 + 2v ) 39 (1+v )( − 1 + 2v )

a37 = a73 = a38 = a83 = a48 = a84 = a49 = a94 = a711 = a117 = a911 = a119 = a77

a78 = a87 = a79 = a97 = a89 = a98 =

−1 . 2(1 + v )(−1 + 2v )

(22)

The eigenvalues for the stiffness are:

1

5E

λk1 = λk2 = λk3 = λk4 = λk5 = λk6 = 0, λk7 = λk8 = , λ = λk10 = 1 + v k9 2 (1 + v )



 √ √ 2E −5 + 4v + 9 − 24v + 48v2 −2E 5 − 4v + 9 − 24v + 48v2 λk11 = , λk12 = . −1 + v + 2v2 −1 + v + 2v2

(23)

It is clearly observed that the eigenvalues of stiffness derived in this work are constant within the elements from Eq. (23). While different locations of integration points will lead to different types of mass matrix, the consistence mass matrix is evaluated by using the following four integration points with weight wi = 1/4:

Point α : ξi = 0.1381966, ηi = 0.1381966, ζi = 0.1381966 Point β : ξi = 0.1381966, ηi = 0.5854102, ζi = 0.1381966 Point γ : ξi = 0.5854102, ηi = 0.1381966, ζi = 0.1381966 Point o: ξi = 0.1381966, ηi = 0.1381966, ζi = 0.5854102,

(24)

and the lump mass matrix is evaluated by using the following four integration points with weight wi = 1/4:

Point α : ξi = 0, ηi = 0, ζi = 0, Point β : ξi = 0, ηi = 1, ζi = 0 Point γ : ξi = 0, ηi = 1, ζi = 0, Point o : ξi = 0, ηi = 0, ζi = 1.

(25)

As the mass matrix varies with the locations of the integration points, the general new locations of integration points for the mass matrix are proposed as shown in Fig. 7:

Point α :

1−t 1−t 1−t 1−t 1−t , ηi = , ζi = , Point β : ξi = , ηi = t, ζi = 3 3 3 3 3 1−t 1−t 1−t 1−t = t, ηi = , ζi = , Point o :ξi = , ηi = , ζi = t, 3 3 3 3

ξi =

Point γ : ξi

(26)

Please cite this article as: E. Li, Z.C. He, Development of a perfect match system in the improvement of eigenfrequencies of free vibration, Applied Mathematical Modelling (2017), http://dx.doi.org/10.1016/j.apm.2017.02.013

ARTICLE IN PRESS

JID: APM 8

[m3Gsc;March 1, 2017;1:55]

E. Li, Z.C. He / Applied Mathematical Modelling 000 (2017) 1–26

Fig. 7. Adjustable integration points for the mass matrix of linear tetrahedral element.

wheret ∈ [0, 1], and weight wi = 1/4 so that as follows:

⎡b

11

⎢ 0 ⎢ 0 ⎢ ⎢ b41 ⎢ 0 ⎢ ⎢ 0 1 Me = ρ V e ⎢ ⎢ b71 4 ⎢ ⎢ 0 ⎢ 0 ⎢ ⎢b101 ⎣ 0 0

0 b22 0 0 b52 0 0 b82 0 0 b112 0

0 0 b33 0 0 b63 0 0 b93 0 0 b123

3

b14 0 0 b44 0 0 b74 0 0 b104 0 0

i=1

wi = 1. Using these four sampling points, the mass matrix can be expressed

0 b25 0 0 b55 0 0 b85 0 0 b115 0

0 0 b36 0 0 b66 0 0 b96 0 0 b126

b17 0 0 b47 0 0 b77 0 0 b107 0 0

0 b28 0 0 b58 0 0 b88 0 0 b118 0

0 0 b39 0 0 b69 0 0 b99 0 0 b129

b110 0 0 b410 0 0 b710 0 0 b1010 0 0

0

0 0

b211 0 0 b511 0 0 b811 0 0 b1111 0

⎤ ⎥

b312 ⎥ ⎥ 0 ⎥ ⎥ 0 ⎥ b612 ⎥ ⎥, 0 ⎥ ⎥ 0 ⎥ b912 ⎥ ⎥ 0 ⎥ ⎦ 0 b1212

(27)

where

b11 = b22 = b33 = b44 = b55 = b66 = b77 = b88 = b99 = b1010 = b1111 = b1212 = t 2 +

(1 − t )2 3

.

2

) The remaining parameters in Eq. (27) equal to 2t (1−t + 2(19−t ) 3 e The summation of all matrix elements of M equals to the volume of the element Ve regardless of value of t: 144





 ρV e 3t 2 + (1 − t )2 + 6t (1 − t ) + 2(1 − t )2

i=1

=

144



 ρV e 4t 2 + 1 − 2t + 6t − 6t 2 + 2 − 4t + 2t 2

i=1

=

3ρ V e

.

 

(28)

total mass of the element

The conservation of mass matrix at the element level whose individual entries can “tuned” using the “knob” t. After obtaining the general mass matrix for 3D elasticity problems using tetrahedral mesh, the eigenvalues of the mass system can be derived as follows:

λm1 = λm2 = λm3 =

1 3 h 24

λm4 = λm5 = λm6 = λm7 = λm8 = λm9 = λm10 = λm11 = λm12 =

1 3 h (16t 2 − 8t + 1 ). 216

(29)

It can be seen from Eq. (29) that the location of integration points controls the eigenvalues of free vibration problems by the term (16t2 − 8t + 1). Fig. 8 plots that the term (16t2 − 8t + 1) varies with the location of integration point t. It can be found that when t equals to 1/4, the largest eigenvalue is obtained. When 0 < t < 1/4, with the increase of t, the eigenvalues will become bigger and bigger. While t is greater than 1/4, the eigenvalues are proportional to t value. Since the new location of integration points can modify the mass matrix to alter its eigenvalues, it gives us the flexibility to balance the stiffness and mass matrices in the discretized model. This offers an excellent approach for accurate solutions Please cite this article as: E. Li, Z.C. He, Development of a perfect match system in the improvement of eigenfrequencies of free vibration, Applied Mathematical Modelling (2017), http://dx.doi.org/10.1016/j.apm.2017.02.013

JID: APM

ARTICLE IN PRESS E. Li, Z.C. He / Applied Mathematical Modelling 000 (2017) 1–26

[m3Gsc;March 1, 2017;1:55] 9

Fig. 8. Effect of integration points on the eigenvalue using single tetrahedral element.

Fig. 9. 2D cantilever beam.

using tetrahedral elements to predict the solutions of eigenfrequency. In the following section, the performance of 3D MRFEM to evaluate the eigenfrequency is illustrated in detail. 4. Numerical examples In this work, the triangular or tetrahedral elements are focused due to the easy creation for complicated boundary. The excellent performance of MR-FEM in the analysis of eigenfrequency of free vibration can be easily extend to quadrilateral and brick elements. 4.1. 2D cantilever beam This section provides a quantitative analysis in the prediction of eigenfrequencies for free vibration problems using different types of numerical methods with combination of MR-FEM. The first 2D cantilever beam with length l = 1 m and width d = 0.2 m is studied, which is subjected to fixed boundary condition at left hand side as outlined in Fig. 9. The material properties are: Young’s Modulus E = 100 MPa, Poisson’s ratio v = 0.3 and the density ρ = 103 kg/m3 . The reference solution is obtained using FEM with 80 0 0 quadrilateral elements. 4.1.1. Effect of integration point on the eigenfrequencies using triangular element As studied in Section 3, the theoretical analysis has indicated that the integration point in the mass matrix has a great effect on the prediction of eigenvalue with a single triangular element. It is expected that the numerical solutions of eigenfrequencies can be improved using coarse mesh in the FEM and SFEM models with combination of MR-FEM. As shown in Figs. 10–12, it is noticed that the eigenfrequencies obtained from FEM, NS-FEM and ES-FEM modes (105 nodes) vary with the integration point in the mass matrix for different mode shapes. Due to the length limit of the paper, only the results for the first four mode shapes are presented. Apparently, the mass matrix using integration point t = 1/3 gives the maximum eigenfrequency in FEM, NS-FEM and ES-FEM models, which is consistent with the analytical solution using single triangular Please cite this article as: E. Li, Z.C. He, Development of a perfect match system in the improvement of eigenfrequencies of free vibration, Applied Mathematical Modelling (2017), http://dx.doi.org/10.1016/j.apm.2017.02.013

ARTICLE IN PRESS

JID: APM 10

[m3Gsc;March 1, 2017;1:55]

E. Li, Z.C. He / Applied Mathematical Modelling 000 (2017) 1–26

Fig. 10. Effect of integration point in the mass matrix on natural frequency for different mode shapes using FEM with combination of MR-FEM. Table 1 Numerical results of eigenfrequency using FEM. Mode shape

Reference

FEM (consistent)

FEM (t = 1/3)

FEM (lumped)

1 2 3 4 5 6 7 8 9 10

9.9291 53.639 79.223 128.10 213.37 236.95 304.16 391.88 395.28 481.53

10.95 58.71 79.37 139.96 233.48 238.17 334.55 396.34 436.20 522.59

10.95 58.91 79.38 140.91 235.98 238.68 340.07 398.27 445.90 533.76

10.93 58.12 79.32 137.24 226.17 236.90 319.34 390.62 410.18 490.86

elements. As NS-FEM is overly-soft model, it always predicts the smallest eigenfrequency compared with ES-FEM and FEM models using the same integration point. On the contrary, the overly-stiff FEM model always gives the maximum eigenfrequency. This demonstrates again that the eigenfrequency is proportional to the stiffening effect of discretized systems. As t equals to 1 corresponding to the lumped mass matrix, it is found that the lowest solutions of eigenfrequency are achieved among FEM, NS-FEM and ES-FEM models respectively. From Figs. 10–12, it is easily understood that the softening effect is inversely proportional to integration point t as t ranges between 0 and 1/3. When t > 1/3, the softening effect is becoming stronger and stronger in all numerical discretized modes including FEM, NS-FEM and ES-FEM models. For the first time, it is found that the integration point in the mass matrix acts as a knob to control the softening and stiffening effects in the prediction of eigenfrequencies. Tables 1–3 list the first ten natural eigenfrequencies obtained from FEM, NS-FEM and ES-FEM using lumped and consistent mass matrices with 105 nodes. As FEM, NS-FEM and ES-FEM always achieve the largest eigenfrequencies with integration point t = 1/3 in the mass matrix, this special integration point with unique property is studied in the following sections Please cite this article as: E. Li, Z.C. He, Development of a perfect match system in the improvement of eigenfrequencies of free vibration, Applied Mathematical Modelling (2017), http://dx.doi.org/10.1016/j.apm.2017.02.013

ARTICLE IN PRESS

JID: APM

[m3Gsc;March 1, 2017;1:55]

E. Li, Z.C. He / Applied Mathematical Modelling 000 (2017) 1–26

11

Fig. 11. Effect of integration point in the mass matrix on natural frequency for different mode shapes using ES-FEM with combination of MR-FEM. Table 2 Numerical results of eigenfrequency using NS-FEM. Mode shape

Reference

NS-FEM (consistent)

NS-FEM (t = 1/3)

NS-FEM (lumped)

1 2 3 4 5 6 7 8 9 10

9.9291 53.639 79.223 128.10 213.37 236.95 304.16 391.88 395.28 481.53

9.4114 50.7519 79.1609 120.4988 174.0303 198.5493 235.8044 278.2750 296.3904 311.0289

9.4171 50.9343 79.1762 121.4037 201.1306 236.2252 284.2326 298.6494 365.4603 388.6063

9.3943 50.2151 79.1153 100.8029 117.9031 171.0610 177.4087 188.1259 191.0808 206.0504

to examine the accuracy and convergence rate for all FEM, NS-FEM and ES-FEM models in detail. The reference solutions are also listed in the table. As indicated in Table 1, the upper bound solutions of eigenfrequencies are always obtained in the FEM model. Apparently, the lumped mass matrix gives the best solutions in FEM model. The solutions from NS-FEM tabulated in Table 2 indicate that NS-FEM always gives the lower bound solutions regardless of integration point in the mass matrix. However, the accuracy of NS-FEM model using integration point t = 1/3 is much better than the consistent and lumped mass matrices. In the Table 3, it is found that the ES-FEM always provides the upper bound solutions in the first 3 mode shapes. However, ES-FEM using lumped mass matrix gives the lower bound solutions after mode shape 3, which implies that the softening effect using lumped mass matrix dominates in the larger mode shape. The ES-FEM with consistent mass matrix and mass matrix using integration point t = 1/3 still behaves the stiffening effect, which predicts the upper Please cite this article as: E. Li, Z.C. He, Development of a perfect match system in the improvement of eigenfrequencies of free vibration, Applied Mathematical Modelling (2017), http://dx.doi.org/10.1016/j.apm.2017.02.013

ARTICLE IN PRESS

JID: APM 12

[m3Gsc;March 1, 2017;1:55]

E. Li, Z.C. He / Applied Mathematical Modelling 000 (2017) 1–26

Fig. 12. Effect of integration point in the mass matrix on natural frequency for different mode shapes using NS-FEM with combination of MR-FEM. Table 3 Numerical results of eigenfrequency using ES-FEM. Mode shape

Reference

ES-FEM (consistent)

ES-FEM (t = 1/3)

ES-FEM (lumped)

1 2 3 4 5 6 7 8 9 10

9.9291 53.639 79.223 128.10 213.37 236.95 304.16 391.88 395.28 481.53

10.03 54.27 79.28 129.90 216.91 237.18 310.15 392.50 404.43 493.49

10.03153 54.46516 79.29773 130.8209 219.5078 237.5931 315.7484 394.4862 414.5215 507.5337

10.01 53.71 79.24 127.23 209.59 235.96 294.83 377.44 386.64 454.40

bound solutions for all mode shapes. With a proper integration point in the large mode shapes, ES-FEM can even bound the exact value from upper and lower sides. Using these bound properties of ES-FEM with MR-FEM, one can now effectively certify a numerical solution and conduct elegant adaptive analyses for solutions of desired accuracy in the analysis of eigenfrequencies of free vibration problems. 4.1.2. Convergence rate of 2D MR-FEM This section further discusses the convergence properties of different numerical methods using MR-FEM. Fig. 13 shows the convergence rate of eigenfrequencies using FEM for different mode shapes. Without loss of generalization, mode shapes for 7, 10, 16, and 20 are randomly selected to analyse the properties of MR-FEM with combination of FEM. It is clearly shown that FEM achieves the upper bounds for all eigenfrequencies, which is consistent with overly-stiff FEM property. Please cite this article as: E. Li, Z.C. He, Development of a perfect match system in the improvement of eigenfrequencies of free vibration, Applied Mathematical Modelling (2017), http://dx.doi.org/10.1016/j.apm.2017.02.013

JID: APM

ARTICLE IN PRESS

[m3Gsc;March 1, 2017;1:55]

E. Li, Z.C. He / Applied Mathematical Modelling 000 (2017) 1–26

13

Fig. 13. Convergence rate of natural frequency for different mode shapes in FEM with combination of MR-FEM.

With increasing number of nodes, the solutions from FEM converge the exact solution from above. From Fig. 13, it is observed that FEM using lumped mass matrix gives the best solution. The consistent mass matrix and mass matrix using integration point t = 1/3 increase the stiffening effect of FEM model further, which gives the poor solutions in the prediction of eigenfrequencies. In particular, the mass matrix using integration point t = 1/3 provides the stiffest model resulting the worst solutions. Fig. 14 presents the convergence rate of eigenfrequency for mode shapes 7, 10, 16, and 20 using NS-FEM. As NS-FEM is overly-soft model [14], the lower bounds of eigenfrequencies obtained using triangular elements with various integration points are observed. In addition, it is found that the lumped mass matrix is able to soften the discretized model further, leading to the lowest eigenfrequencies for all mode shapes. The integration point t = 1/3 gives the largest eigenfrequency in NS-FEM, which provides the best solutions in NS-FEM for mode shape. That implies mass matrix using integration point t = 1/3 reduces the softening effect of NS-FEM, making NS-FEM mode more close to exact system. The softening effect of NSFEM using consistent mass matrix is between the lumped mass matrix and mass matrix using integration point t = 1/3. This 2D example shows clearly that NS-FEM can provide more close solutions from lower sides with adjustment of integration point in the mass matrix. The convergence rate of eigenfrequencies using ES-FEM for mode shapes 7, 10, 16 and 20 is outlined in Fig. 15. From Fig. 15, it is observed that ES-FEM using lumped mass matrix converges the reference solutions from below. On the other hand, ES-FEM using consistent mass matrix and mass matrix with integration point t = 1/3 approach the reference solutions from above. In terms of accuracy, ES-FEM using the consistent mass matrix gives the best solutions. In the previous study, ES-FEM is considered as a softened model which always gives the upper bound solutions in the computation of eigenfrequency [15]. However, the mass redistributed method has provided an alternative way to soften ES-FEM model. It is interesting to note that ES-FEM provides the bound solutions in the prediction of larger eigenfrequencies with MR-FEM. This 2D numerical example clearly indicates that the softening and stiffening effects in the prediction of eigenfrequencies can be controlled by the integration point in the mass matrix. The quantitative analysis of integration point in the mass Please cite this article as: E. Li, Z.C. He, Development of a perfect match system in the improvement of eigenfrequencies of free vibration, Applied Mathematical Modelling (2017), http://dx.doi.org/10.1016/j.apm.2017.02.013

JID: APM 14

ARTICLE IN PRESS

[m3Gsc;March 1, 2017;1:55]

E. Li, Z.C. He / Applied Mathematical Modelling 000 (2017) 1–26

Fig. 14. Convergence rate of natural frequency for different mode shapes in NS-FEM with combination of MR-FEM.

matrix provides a solution to improve the accuracy of eigenfrequency for free vibration problems. Using MR-FEM, we can bind the exact solution in a more effective way.

4.2. Analysis of 2D bracket using MR-FEM The second 2D bracket model as shown in Fig. 16 is further analysed in this section. The bottom is subjected to fixed boundary condition. The material properties are the same as Section 4.1. FEM model using 4117 quadrilateral elements is employed to provide the reference solutions. It is expected that the lower bound solutions are always achieved in the NS-FEM model as shown in Fig. 17. Using integration point t = 1/3 in the mass matrix, NS-FEM can give much better lower bound solutions. The lumped and consistent mass matrices corresponding t = 1 and t = 2/3 give very poor solutions in the NS-FEM model. Using these lower bound properties of NS-FEM and upper bound properties of FEM, we can always effectively certify a numerical solution. In addition, it is noticed that ES-FEM can also bind the solutions from above and lower using consistent and lumped mass matrices respectively. This exciting finding is very useful to conduct elegant adaptive analysis for solutions of desired accuracy.

4.3. Performance of 3D MR-FEM in eigenfrequencies prediction Lured by the excellent features of mass redistributed method in the 2D model, the 3D mass redistributed method is further studied in this section to analyse the eigenfrequency. As shown in Fig. 18, a three dimensional cantilever beam with dimension 1 × 0.2 × 0.2m subjected to fixed boundary condition is presented. The material properties are: Young’s Modulus E = 1 GPa, Poisson’s ratio v = 0.4 and the density ρ = 103 kg/m3 . The reference solution using 40,0 0 0 brick elements are obtained from FEM model for the purpose of comparison. Please cite this article as: E. Li, Z.C. He, Development of a perfect match system in the improvement of eigenfrequencies of free vibration, Applied Mathematical Modelling (2017), http://dx.doi.org/10.1016/j.apm.2017.02.013

JID: APM

ARTICLE IN PRESS

[m3Gsc;March 1, 2017;1:55]

E. Li, Z.C. He / Applied Mathematical Modelling 000 (2017) 1–26

15

Fig. 15. Convergence rate of natural frequency for different mode shapes in ES-FEM with combination of MR-FEM.

Fig. 16. 2D bracket model.

4.3.1. Effect of integration point on the eigenfrequencies using tetrahedral element In this 3D free vibration problem, it is still seen that the integration point of mass matrix controls the accuracy of eigenfrequency of 3D FEM, NS-FEM and ES-FEM models with 106 nodes as depicted in Figs. 19–21 respectively. Due to the length limit of the paper, only the eigenfrequencies for mode shapes 1, 2, 3 and 4 are presented. From Figs. 19–21, it is shown that the softening effect in FEM, NS-FEM and ES-FEM models is proportional to the integration point t as t is greater than 1/4, which leads the eigenfrequency to decrease with increasing integration point t. On the contrary, the stiffening Please cite this article as: E. Li, Z.C. He, Development of a perfect match system in the improvement of eigenfrequencies of free vibration, Applied Mathematical Modelling (2017), http://dx.doi.org/10.1016/j.apm.2017.02.013

JID: APM 16

ARTICLE IN PRESS

[m3Gsc;March 1, 2017;1:55]

E. Li, Z.C. He / Applied Mathematical Modelling 000 (2017) 1–26

Fig. 17. Effect of integration point t on the eigenfrequency from different numerical methods.

Fig. 18. Three dimensional cantilever beam.

effect is becoming stronger and stronger in 3D FEM, NS-FEM and ES-FEM models when the integration point t in the mass matrix is greater than 0 but less than 1/4. Apparently, the integration point t = 1/4 in the mass matrix results the stiffest model, which is consistent with the analytical solutions using singe tetrahedral element as mentioned in Section 3. The overly-stiff FEM always gives the largest eigenfrequencies among all discretized model with the same integration point in the mass matrix. On the other hand, the smallest eigenfrequencies is always achieved by NS-FEM model with the same integration point in the mass matrix due to overly-soft effect. This 3D example has confirmed that the integration point t = 1/4 results the stiffest model including FEM, NS-FEM and ES-FEM models. Please cite this article as: E. Li, Z.C. He, Development of a perfect match system in the improvement of eigenfrequencies of free vibration, Applied Mathematical Modelling (2017), http://dx.doi.org/10.1016/j.apm.2017.02.013

ARTICLE IN PRESS

JID: APM

[m3Gsc;March 1, 2017;1:55]

E. Li, Z.C. He / Applied Mathematical Modelling 000 (2017) 1–26

17

Fig. 19. Effect of integration point in the mass matrix on natural frequency for different mode shapes using FEM with combination of MR-FEM. Table 4 Numerical results of eigenfrequency using FEM. Mode shape

FEM

FEM (consistent)

FEM (t = 1/4)

FEM (lumped)

1 2 3 4 5 6 7 8 9 10

31.738 31.738 137.24 169.68 169.68 252.47 401.97 401.97 411.60 665.18

43.39 44.14 218.63 223.24 224.90 256.82 512.11 514.52 655.13 770.88

43.45 44.21 224.62 225.40 234.01 256.98 519.34 521.79 700.45 774.85

43.15 43.90 181.18 217.00 218.00 256.22 485.71 488.19 534.83 753.72

The first ten natural eigenfrequencies computed from FEM, NS-FEM and ES-FEM with 106 nodes are tabulated in Tables 4–6 respectively. For the purpose of comparison, the reference solutions are also given in the Table. As indicated in Table 4, the upper bound solutions of FEM model are obtained regardless of integration point in the mass matrix. Again, it is found that the lumped mass matrix in FEM model gives more accurate solutions compared with consistent mass matrix and mass matrix using integration point t = 1/4. The stiffest FEM corresponding to integration point t = 1/4 results the poorest solutions. From Table 5, it is expected that 3D NS-FEM still gives the lower bound solutions for all mode shapes regardless of integration point. The accuracy of NS-FEM has been improved significantly by using integration point t = 1/4 in the mass matrix, which implies that the softening effect in NS-FEM has been decreased. In the NS-FEM model, the lumped mass matrix giving the most softened model provides the worst solutions compared with consistent mass matrix and mass matrix using integration point t = 1/4. In the Table 6, it is found the ES-FEM always provides the upper bound solutions in the first Please cite this article as: E. Li, Z.C. He, Development of a perfect match system in the improvement of eigenfrequencies of free vibration, Applied Mathematical Modelling (2017), http://dx.doi.org/10.1016/j.apm.2017.02.013

ARTICLE IN PRESS

JID: APM 18

[m3Gsc;March 1, 2017;1:55]

E. Li, Z.C. He / Applied Mathematical Modelling 000 (2017) 1–26

1

3 Fig. 20. Effect of integration point in the mass matrix on natural frequency for different mode shapes using NS-FEM with combination of MR-FEM. Table 5 Numerical results of eigenfrequency using NS-FEM. Mode shape

FEM

NS-FEM (consistent)

NS-FEM (t = 1/4)

NS-FEM (lumped)

1 2 3 4 5 6 7 8 9 10

31.738 31.738 137.24 169.68 169.68 252.47 401.97 401.97 411.60 665.18

26.68 27.31 125.13 144.26 146.42 251.95 339.79 345.16 366.95 550.21

26.72 27.35 133.10 145.57 147.73 252.10 346.37 351.58 393.03 570.57

26.51 27.14 103.35 139.35 141.49 251.35 298.05 316.43 322.35 381.97

two mode shape. After mode shape two, the lower bound solutions of ES-FEM using consistent mass matrix are observed. The consistent mass matrix and mass matrix with integration point t = 1/4 always gives stiffened effect in the ES-FEM model in the prediction of eigenfrequencies, leading the solutions larger than exact one. This 3D example has confirmed again that the integration point in the mass matrix has altered the stiffening or softening effect in the numerical discretized model. 4.3.2. Convergence rate of 3D MR-FEM The convergence properties of different numerical methods using 3D MR-FEM is analysed in this section. First, the convergence rate of eigenfrequencies using FEM for different mode shape is studied. Without loss of generalization, different mode shapes for 5, 10, 15, and 20 are randomly selected. As shown in Fig. 22, the upper bound solutions of FEM are always Please cite this article as: E. Li, Z.C. He, Development of a perfect match system in the improvement of eigenfrequencies of free vibration, Applied Mathematical Modelling (2017), http://dx.doi.org/10.1016/j.apm.2017.02.013

ARTICLE IN PRESS

JID: APM

[m3Gsc;March 1, 2017;1:55]

E. Li, Z.C. He / Applied Mathematical Modelling 000 (2017) 1–26

19

Fig. 21. Effect of integration point in the mass matrix on natural frequency for different mode shapes using ES-FEM with combination of MR-FEM.

Table 6 Numerical results of eigenfrequency using ES-FEM. Mode shape

Reference

ES-FEM (consistent)

ES-FEM (t = 1/4)

ES-FEM (lumped)

1 2 3 4 5 6 7 8 9 10

31.738 31.738 137.24 169.68 169.68 252.47 401.97 401.97 411.60 665.18

33.23 33.31 155.57 178.06 179.23 254.22 425.08 426.03 465.13 708.37

33.28 33.36 165.49 179.56 180.72 254.37 432.07 432.97 497.71 728.29

33.03 33.11 128.51 172.39 173.59 253.65 378.50 400.51 401.34 619.44

observed. The lumped mass matrix performs better than consistent mass matrix and mass matrix using integration point t = 1/4. On the other hand, the stiffening effect becomes stronger in the mass matrix using integration point t = 1/4, which gives the numerical solutions less accurate than the consistent and lumped mass matrices. The convergence rate of eigenfrequency for 5th, 10th, 15th and 20th mode shapes using NS-FEM is illustrated in Fig. 23. It is clearly indicated that 3D NS-FEM using tetrahedral element gives the lower bound solutions due to its overly-soft effect again. The mass matrix with integration point t = 1/4 stands out clearly compared with lumped and consistent mass matrices. This is because mass matrix using integration t = 1/4 makes the NS-FEM less softened. Thus, the numerical solutions from mass matrix with integration point t = 1/4 in NS-FEM model are close to the reference one. From Fig. 23, it is also noticed that the NS-FEM model with lumped mass matrix provides the worst solutions due to overly-softened effect. Please cite this article as: E. Li, Z.C. He, Development of a perfect match system in the improvement of eigenfrequencies of free vibration, Applied Mathematical Modelling (2017), http://dx.doi.org/10.1016/j.apm.2017.02.013

JID: APM 20

ARTICLE IN PRESS

[m3Gsc;March 1, 2017;1:55]

E. Li, Z.C. He / Applied Mathematical Modelling 000 (2017) 1–26

Fig. 22. Convergence rate of natural frequency for different mode shapes in FEM with combination of MR-FEM.

Fig. 24 shows the convergence rate of eigenfrequency using ES-FEM for 5th, 10th, 15th and 20th mode shapes using various integration points. It is seen that for larger mode shapes, ES-FEM can approach the exact solution from above and below with different integration points. With refinement of mesh, both lumped and consistent mass matrix can converge the exact solution from below and above respectively, which provides an alternative approach to bind the exact solution using coarse mesh. The lumped mass matrix in the larger mode shape results the ES-FEM model softer, while the consistent mass matrix provides the stiffening effect of ES-FEM model. The stiffest ES-FEM model corresponding to integration point t = 1/4 gives the worst solution. This 3D example illustrates clearly one very important fact that we now can bind the exact eigenfrequencies from both sides in the analysis of free vibration problems in the large mode shape using ES-FEM. It is expected that the alteration of softened and stiffened discretized model can be achieved by adjustment of integration point in the 3D mass matrix based on current study. However, in the previous study, the softening effect depends on the number of node to create the smoothing domain [13]. The more elements, the larger smoothing effects are [13]. For the first time, it is found that mass redistribution in the mass matrix has a softening or stiffening effect on the discretized system in the evaluation of eigenfrequency. This important finding opens a new window for the development of a new class of numerical methods via manipulating the integration point in the mass matrix to achieve a desirable accuracy in the study of free vibration problems. 4.4. Analysis of 3D bracket using MR-FEM The fourth 3D bracket model as shown in Fig. 25 is further studied in this section. The bottom is subjected to fixed boundary condition. The material properties are taken as: Young’s Modulus E = 71 GPa, Poisson’s ratio v = 0.3, density ρ = 2700 kg/m3 . The reference solutions are obtained from FEM using FEM with 38,862 brick elements. Please cite this article as: E. Li, Z.C. He, Development of a perfect match system in the improvement of eigenfrequencies of free vibration, Applied Mathematical Modelling (2017), http://dx.doi.org/10.1016/j.apm.2017.02.013

JID: APM

ARTICLE IN PRESS

[m3Gsc;March 1, 2017;1:55]

E. Li, Z.C. He / Applied Mathematical Modelling 000 (2017) 1–26

21

Fig. 23. Convergence rate of natural frequency for different mode shapes in NS-FEM with combination of MR-FEM.

The effect of integration point of mass matrix on the prediction of eigenfrequencies for mode shapes 1, 4, 7 and 10 using different numerical methods is plotted in Fig. 26 (The mode shape number is randomly selected in all numerical examples). Again, it is observed that the integration points of mass matrix play an important role in the computation of eigenfrequencies. The ES-FEM with lumped mass matrix corresponding to t = 1 stands out clearly. As the mode shape increases, the lumped mass matrix makes ES-FEM softer, which gives the lower bound solutions. The NS-FEM with mass matrix using integration point t = 1/4 gives more close solutions to the reference ones compared with other integration points in the mass matrix. Again, the upper bound solutions of FEM are always observed regardless of integration point.

4.5. Discussion In the numerical examples described above, the locations of integration points in the mass matrix are strictly refined inside the element of triangular and tetrahedral mesh. In other words, the adjustable parameter t controlling the locations of integration points in the mass matrix is: 0 ≤ t ≤ 1. Thus, a general question naturally arises: can we develop a method to improve the accuracy of the eigenfrequency by further increasing the integration point t? If this can be done, it might provide a new approach to obtain more accurate and even nearly exact solutions. In the following section, the integration point t > 1 as shown in Fig. 27 is formulated in the MR-FEM. For the sake of simplicity, the numerical results (105 nodes) of eigenfrequencies for 4th and 10th model shapes (randomly selected) presented in example 2.1 are shown in Fig. 28 as 0 ≤ t ≤ 2. It is noticed that t = 0.9 and 1.6 in the ES-FEM and FEM respectively give the exact solution in the computation of eigenfrequency of 4th mode shape. Additionally, the exact solutions are obtained in the FEM and ES-FEM models as t = 0.8 and 1.07 are applied in the prediction of eigenfrequency of Please cite this article as: E. Li, Z.C. He, Development of a perfect match system in the improvement of eigenfrequencies of free vibration, Applied Mathematical Modelling (2017), http://dx.doi.org/10.1016/j.apm.2017.02.013

JID: APM 22

ARTICLE IN PRESS

[m3Gsc;March 1, 2017;1:55]

E. Li, Z.C. He / Applied Mathematical Modelling 000 (2017) 1–26

Fig. 24. Convergence rate of natural frequency for different mode shapes in ES-FEM with combination of MR-FEM.

0.1 0.04

0.05

0.16

Fig. 25. 3D bracket model.

Please cite this article as: E. Li, Z.C. He, Development of a perfect match system in the improvement of eigenfrequencies of free vibration, Applied Mathematical Modelling (2017), http://dx.doi.org/10.1016/j.apm.2017.02.013

ARTICLE IN PRESS

JID: APM

E. Li, Z.C. He / Applied Mathematical Modelling 000 (2017) 1–26

[m3Gsc;March 1, 2017;1:55] 23

Fig. 26. Effect of integration point in the different numerical methods.

10th mode shape. This interesting phenomenon indicates that we can even compute the exact solution using coarse mesh with adjustment of integration point in the mass matrix. Next, the effect of mass redistribution using MR-FEM on the eigenmodes is also studied. The modal assurance criterion (MAC), which is used to evaluate the consistency degree of the eigenmodes between the numerical results and the reference, is defined as follows:

Maci,i =

  φ T φi 2 i 

, φiT φi φiT φi



(30)

where φi is the ith model shape obtained using numerical methods, and φi is the ith reference mode shape obtained using very fine mesh. φiT and φiT are the transposition of φi and φi , respectively. The beam is also discretized into 156 nodes and 250 triangular elements. Fig. 29(a) and (b) plot the effect of integration point using MR-FEM on modal assurance criterion. It can be seen from the figure that eigenmodes of NS-FEM using MR-FEM is very sensitive to the integration point t of mass redistribution. The mass matrix corresponding to integration point t = 1/3 gives the best solutions in the NS-FEM model. In addition, it is noticed that the effect of integration point of mass matrix is smaller in both FEM and ES-FEM models for lower mode shape, and the numerical results in both FEM and ES-FEM models are almost identical for 5th mode shape as shown in Fig. 29(a). However, the variations of eigenmodes in both FEM and ES-FEM modes for high mode shape (10th mode shape) are also sensitive to the integration point t with mass redistribution as outlined in Fig. 29(b). More importantly, it is worth noting that integration point t = 1.1 (outside of the element) in the FEM model gives the best MAC among all the integration point t. On the other hand, it is observed that the numerical results in the ES-FEM model are always better than those from FEM and NS-FEM models as 0 ≤t ≤ 1 for 10th mode shape. Please cite this article as: E. Li, Z.C. He, Development of a perfect match system in the improvement of eigenfrequencies of free vibration, Applied Mathematical Modelling (2017), http://dx.doi.org/10.1016/j.apm.2017.02.013

JID: APM 24

ARTICLE IN PRESS

[m3Gsc;March 1, 2017;1:55]

E. Li, Z.C. He / Applied Mathematical Modelling 000 (2017) 1–26

Fig. 27. Illustration of integration t outside of the element.

Fig. 28. Effect of integration point on the eigenfrequencies as 0 ≤ t ≤ 2.

Please cite this article as: E. Li, Z.C. He, Development of a perfect match system in the improvement of eigenfrequencies of free vibration, Applied Mathematical Modelling (2017), http://dx.doi.org/10.1016/j.apm.2017.02.013

JID: APM

ARTICLE IN PRESS E. Li, Z.C. He / Applied Mathematical Modelling 000 (2017) 1–26

[m3Gsc;March 1, 2017;1:55] 25

Fig. 29. Effect of integration point on modal assurance criterion with mass redistribution.

However, determining how to select the correct integration point in the mass matrix is a challenge task as t > 1. This could be beyond our study which needs a rigorous proof, and we hope that these interesting results can draw some attention from mathematical communality to help us accomplish a class of more efficient computational methods in the analysis of free vibration problems. 5. Conclusion In this study, the mass redistributed finite element method (MR-FEM) is formulated to balance the stiffness and mass matrix in order to improve the accuracy of the eigenfrequency in the 2D and 3D free vibration problems. With a proper selection of integration points in the mass matrix, the accuracy for standard finite element method (FEM) and smoothed finite element methods (SFEM) including NS-FEM and ES-FEM models can be improved. Based on the numerical results, the following remarks can be made: (1) With adjustment of integration point in the mass matrix, the softening and stiffening effects can be altered in both FEM and SFEM models. (2) The integration points t = 1/3 and t = 1/4 in the mass matrix for 2D and 3D problems give the best accuracy in NS-FEM model. (3) The lumped mass matrix corresponding to t = 1 always predicts the best eigenfrequencies in 2D and 3D FEM models as the integration point is refined within the element (0 < t < 1). (4) In the prediction of eigenfrequency for large mode shape, the ES-FEM can bound the solutions by consistent and lumped mass matrix as 0 < t < 1. Please cite this article as: E. Li, Z.C. He, Development of a perfect match system in the improvement of eigenfrequencies of free vibration, Applied Mathematical Modelling (2017), http://dx.doi.org/10.1016/j.apm.2017.02.013

JID: APM 26

ARTICLE IN PRESS

[m3Gsc;March 1, 2017;1:55]

E. Li, Z.C. He / Applied Mathematical Modelling 000 (2017) 1–26

(5) The numerical procedure of MR-FEM without increasing the computational cost is very easy and straightforward. (6) The exact solutions using coarse mesh can be obtained in FEM and ES-FEM models with a proper selection of integration point in the mass matrix as the location of integration point is not strictly inside the element. Although these preliminary results have lacked some rigorous mathematical proof, it has revealed an important fact that the alteration of integration point in the mass matrix plays an important role to match the stiffness matrix. That interesting finding could be extended to some complicated cases such as coupling problems and structures to improve the prediction of eigenfrequency. Acknowledgments The project is supported by the Science Fund of State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body Nos. 31315002 and 31615002. The authors also wish to thank Research Project of State Key Laboratory of Mechanical Systems and Vibration MSV 201613 and MSV201711, and Open Foundation of the State Key Laboratory of Fluid Power and Mechatronic Systems GZK F-201601. Reference [1] Z.C. He, G.Y. Li, Z.H. Zhong, A.G. Cheng, G.Y. Zhang, E Li, An improved modal analysis for three-dimensional problems using face-based smoothed finite element method, Acta Mech. Solida Sin. 26 (2013) 140–150. [2] I. Senjanovic´ , M. Tomic´ , N. Vladimir, N Hadžic´ , An approximate analytical procedure for natural vibration analysis of free rectangular plates, Thin-Walled Struct. 95 (2015) 101–114. [3] O.C. Zienkiewicz, R.L Taylor, The Finite Element Method, fifth ed, Butterworth-Heinemann, Oxford, Boston, 20 0 0. [4] G.R. Liu, S.S Quek, The Finite Element Method: A Practical Course, Butterworth-Heinemann, Oxford, 2003. [5] A.B. Sabir, M.S. Djoudi, A Sfendji, Natural frequencies of arch bridges with slab decks by the finite element method, Thin-Walled Struct. 18 (1994) 31–45. [6] W.L. Wood, Note on finite element predictions of natural frequencies, Commun. Numer. Methods Eng. 10 (1994) 731–734. [7] G.R. Liu, T. Nguyen-Thoi, K.Y. Lam, An edge-based smoothed finite element method (ES-FEM) for static, free and forced vibration analyses of solids, J. Sound Vib. 320 (2009) 1100–1130. [8] R.D. Marangoni, L.M. Cook, N Basavanhally, Upper and lower bounds to the natural frequencies of vibration of clamped rectangular orthotropic plates, Int. J. Solids Struct. 14 (1978) 611–623. [9] N.-E. Wiberg, R. Bausys, P Hager, Improved eigenfrequencies and eigenmodes in free vibration analysis, Comput. Struct. 73 (1999) 79–89. [10] O. Friberg, P. Moller, D. Makovicka, N.E Wiberg, An adaptive procedure for eigenvalue problems using the hierarchical finite-element method, Int. J. Numer. Methods Eng. 24 (1987) 319–335. [11] R.I. Mackie, Improving finite-element predictions of modes of vibration, Int. J. Numer. Methods Eng. 33 (1992) 333–344. [12] G.R. Liu, On G space theory, Int. J. Comput. Methods 6 (2009) 257–289. [13] G.R. Liu, G.Y. Zhang, A normed G space and weakened weak (W2) formulation of a cell-based smoothed point interpolation method, Int. J. Comput. Methods 06 (2009) 147–179. [14] G.R. Liu, A G space theory and a weakened weak (W2) form for a unified formulation of compatible and incompatible methods: Part II applications to solid mechanics problems, Int. J. Numer. Methods Eng. 81 (2010) 1127–1156. [15] G.R. Liu, A generalized gradient smoothing technique and the smoothed bilinear form for Galerkin formulation of a wide class of computational methods, Int. J. Comput. Methods 5 (2008) 199–236. [16] G.R. Liu, T. Nguyen-Thoi, H. Nguyen-Xuan, K.Y Lam, A node-based smoothed finite element method (NS-FEM) for upper bound solutions to solid mechanics problems, Comput. Struct. 87 (2009) 14–26. [17] Z.Q. Zhang, G.R. Liu, Upper and lower bounds for natural frequencies: a property of the smoothed finite element methods, Int. J. Numer. Methods Eng. 84 (2010) 149–178. [18] E. Li, G.R. Liu, V. Tan, Simulation of hyperthermia treatment using the edge-based smoothed finite-element method, Numer. Heat Transfer Part A-Appl. 57 (2010) 822–847. [19] Z.C. He, G.Y. Li, Z.H. Zhong, A.G. Cheng, G.Y. Zhang, E. Li, et al., An ES-FEM for accurate analysis of 3D mid-frequency acoustics using tetrahedron mesh, Comput. Struct. 106 (2012) 125–134. [20] Z.C. He, G.Y. Li, Z.H. Zhong, A.G. Cheng, G.Y. Zhang, G.R. Liu, et al., An edge-based smoothed tetrahedron finite element method (ES-T-FEM) for 3D static and dynamic problems, Comput. Mech. 52 (2013) 221–236. [21] E. Li, Z.P. Zhang, C.C. Chang, G.R. Liu, Q Li, Numerical homogenization for incompressible materials using selective smoothed finite element method, Compos. Struct. 123 (2015) 216–232. [22] M.N. Guddati, B. Yue, Modified integration rules for reducing dispersion error in finite element methods, Comput. Methods Appl. Mech. Eng. 193 (2004) 275–287. [23] B. Yue, M.N. Guddati, Dispersion-reducing finite elements for transient acoustics, J. Acoust. Soc. Am. 118 (2005) 2132–2141. [24] Z. He, G. Li, G.Y. Zhang, G.R. Liu, Y. Gu, E Li, Acoustic analysis using a mass-redistributed smoothed finite element method with quadrilateral mesh, Eng. Comput. 32 (2015) 2292–2317. [25] E. Li, Z.C. He, Y. Jiang, B. Li, 3D mass-redistributed finite element method in structural–acoustic interaction problems, Acta Mech. 227 (2016) 857–879. [26] H. Nguyen-Xuan, L.V. Tran, C.H. Thai, T Nguyen-Thoi, Analysis of functionally graded plates by an efficient finite element method with node-based strain smoothing, Thin-Walled Struct. 54 (2012) 1–18. [27] G.R. Liu, On partitions of unity property of nodal shape functions: rigid-body-movement reproduction and mass conservation, Int. J. Comput. Methods 13 (2016) 1640 0 03. [28] E. Li, Z.C. He, X. Xu, An edge-based smoothed tetrahedron finite element method (ES-T-FEM) for thermomechanical problems, Int. J. Heat Mass Transf. 66 (2013) 723–732. [29] E. Li, G.R. Liu, V. Tan, Z.C. He, Modeling and simulation of bioheat transfer in the human eye using the 3D alpha finite element method (aFEM), Int. J. Numer. Meth. Biomed. 26 (2010) 955–976. [30] E. Li, Z.C. He, G. Wang, An exact solution to compute the band gap in phononic crystals, Comput. Mater. Sci. 122 (2016) 72–85. [31] E. Li, G.R. Liu, V. Tan, Z.C. He, An efficient algorithm for phase change problem in tumor treatment using alpha FEM, Int. J. Therm. Sci. 49 (2010) 1954–1967.

Please cite this article as: E. Li, Z.C. He, Development of a perfect match system in the improvement of eigenfrequencies of free vibration, Applied Mathematical Modelling (2017), http://dx.doi.org/10.1016/j.apm.2017.02.013