Finite element analysis of bending waves in Mindlin plates with Perfectly Matched Layers

Da Cao 1

Graduate School of Engineering, The University of Tokyo 7-3-1, Hongo, Bunkyo-ku, Tokyo 113-8656, Japan

Naohisa Inoue 2

Faculty of Engineering, Maebashi Institute of Technology 460-1, Kamisadori, Maebashi-shi, Gunma 371-0816, Japan

Tetsuya Sakuma 3

Graduate School of Engineering, The University of Tokyo 7-3-1, Hongo, Bunkyo-ku, Tokyo 113-8656, Japan

ABSTRACT It is important to determine the boundary conditions of walls and floors precisely when simulating the building acoustics. For a certain room, the extension of the spans can be considered as infinite edges. The Perfect Matched Layer (PML) is an artificial absorbing domain for the wave propagations and is widely used in finite element analysis to simulate the acoustical free field conditions right now. In this paper, an e ff ective PML technique for the plate structure will be presented. The PML formulation will be derived based on the Mindlin plate theory and the implementation method will be introduced. This technique will be validated through the numerical experiments. The accuracy and limits of the presented technique will be discussed based on the numerical results compared with the analytical results. The results show that the presented PML technique is e ff ective and accurate to simulate the plate as an infinite large plate. It is expected to implement the technique in the further research of the structure borne sound, such as floor impact sounds.

1. INTRODUCTION

The modelling of boundaries of the structure has a significant e ff ect on the numerical simulation of the floor impact sound and sound insulation of walls. When the simulation target is one specific room in the building, one potential method is to model the floors of neighbor rooms as non-reflective ends (free field) if reflected waves are dominated by the nearest columns and beams of the targeted room.

1 dcao@g.ecc.u-tokyo.ac.jp

2 inoue@maebashi-it.ac.jp

3 sakuma@arch1.t.u-tokyo.ac.jp

a slaty. inter.noise 21-24 AUGUST SCOTTISH EVENT CAMPUS O ¥, ? GLASGOW

Berenger [ 1 ] suggested an artificial absorbing layer, called Perfectly Matched Layer (PML), to absorb plane waves from all directions without reflections. Several research shows that the PML technique can also be applied to elastic waves for both finite di ff erence methods [ 2 ] and finite element methods (FEM) [ 3 ]. However, to the best of authors’ knowledge, no academic study has been published about PML for bending waves in FEM of Mindlin plates. Therefore, in the first part of this paper, the PML formulation for bending waves based on Mindlin plate theory is derived. The implementation method for FEM is also introduced. In the latter part, numerical experiments are done to validate the proposed PML technique. Several numerical factors are also examined to determine the e ffi ciency of this technique.

2. FINITE ELEMENT FORMULATION

The governing equation of Mindlin theory can be derived from the weak form of the generalized elastic wave equation by processing the integral along the thickness of the plate under several kinematic assumptions [ 4 ]. In this section, the weak form of the governing equation of the PML for bending waves in a Mindlin plate is derived by performing the similar algebraic operation to the governing equation of the PML for the elastic field.

2.1. PML Formulation for Generalized Elastic Waves Each component of the displacement vector and the stress tensor are splitted into three axes as follows, where subscripts and the overscript denote original components and splitted components respectively, along three directions 1, 2, 3.

u i = X

k u k i , σ i j X

k σ k i j

Equation of Motion

For an isotropic material with density ρ , the equation of motion is written as

∂σ i j

∂ x j − ρ ∂ 2 u i

X

∂ t 2 = 0 .

j

Splitting the equation of motion into three directions and adding attenuation terms, the splitted equation is written as ∂ 2 u j i ∂ t 2 + τ j ( x j ) ∂ u j i ∂ t = 1

ρ ∂σ i j

∂ x j .

Then, assuming the time harmonic condition which means ∂/∂ t → j ω , the splitted equation can be written in the form as 1 γ j ( x j ) ∂σ ij

∂ x j + ω 2 ρ u j i = 0 , γ j ( x j ) = 1 + τ j ( x j )

j ω .

Finally, merging the equation above by the overscript j , the equation of motion of PML along direction x i is written as X

1 γ j ( x j ) ∂σ i j

∂ x j + ω 2 ρ u i = 0 . (1)

j

Constitutive Equation

For an isotropic material with Lamé’s first constant λ and Lamé’s second constant µ , the constitutive equation

 1 for i = j 0 for i , j

 X

 δ ij + µ ∂ u i

∂ x j + µ ∂ u j

∂ u k ∂ x k

σ ij = λ

∂ x i , δ i j =

k

can also be splitted into each direction as

σ k ij = ∂

∂ x k ( λ u k δ i j + µ u i δ jk + µ u j δ ik ) .

Adding the attenuation term similarly with the equation of motion, the following equation can be obtained, ∂σ k ij ∂ t + τ k ( x k ) σ k ij = ∂ 2

∂ t ∂ x k ( λ u k δ i j + µ u i δ jk + µ u j δ ik ) .

Assuming the time harmonic condition, the equation above is rewritten into the form as

σ k ij = 1 γ k ( x k ) ∂ ∂ x k ( λ u k δ i j + µ u i δ jk + µ u j δ ik ) .

Finally, merging the equation above by the overscript k , the constitutive equation of PML is written as

 δ i j + µ 1 γ j ( x j ) ∂ u i ∂ x j + µ 1 γ i ( x i ) ∂ u j

 X

1 γ k ( x k ) ∂ u k ∂ x k

∂ x i . (2)

σ ij = λ

k

Complex Coordinate System

For each direction, introducing the coordinate transformation as follows

ˆ x i = Z x i

0 γ i ( s ) ds .

The partial derivative of each direction can be written as ∂ ∂ ˆ x i 1 γ i ( x i ) ∂ ∂ x i , (3)

which indicates that the equation of motion and constitutive equation of PML can have the same expressions as in the elastic field by rewriting the equations from the regular Cartesian coordinate system ( x 1 , x 2 . x 3 ) into the complex coordinate system (ˆ x 1 , ˆ x 2 , ˆ x 3 ). Therefore, the following operation will be performed in the complex coordinate system (ˆ x 1 , ˆ x 2 , ˆ x 3 ).

Weak Form Formulation

After multiplying Equation 1 with the weight function δ u i , summing from all the three directions and performing the integral over the entire domain, the following equation can be obtained, Z

 δ u i d ˆ V = 0 .

 X

∂σ i j

X

∂ ˆ x j + ω 2 ρ u i

i

j

Applying the divergence theorem, the weak form of PML for elastic waves can be written as Z

 X

 d ˆ V − Z

i , j δε ij σ ij − ω 2 ρ X

X

i δ u i σ i d ˆ S = 0 , (4)

i δ u i u i

where

∂ ˆ x j + ∂ u j

∂ u i

! , σ i = X

δε ij = 1

j σ i j n j .

∂ ˆ x i

2

2.2. Introduction of Kinematical Assumptions In the following formulation, the neutral plane is defined as the ˆ x 1 ˆ x 2 surface at ˆ x 3 = 0, and ˆ x 3 direction is perpendicular to the plate. The plate thickness is defined as − t / 2 ≤ ˆ x 3 ≤ t / 2 and the plate including PML is a rectangular area. Considering the bending deformation of Mindlin plate, following kinematic assumptions [ 4 ] are introduced into the elastic body:

1. A straight line initially normal to the neutral plane deforms linearly and remains straight.

2. The displacement along ˆ x 3 -direction is independent of ˆ x 3 .

3. The displacement of the neutral plane is of ˆ x 3 -direction exclusively.

4. The normal stress component in ˆ x 3 -direction σ 33 vanishes.

From Assumption 1-3, the three-dimensional displacement vector in the elastic body can be written as follows, u α (ˆ x 1 , ˆ x 2 , ˆ x 3 ) = ψ α (ˆ x 1 , ˆ x 2 )ˆ x 3 , u 3 (ˆ x 1 , ˆ x 2 , ˆ x 3 ) = w (ˆ x 1 , ˆ x 2 )

where the subscript α takes 1 and 2, and β is defined as β ≡ mod( α, 2) + 1. ψ α is the rotational angle around x α which is the axis of the cross section perpendicular to the x β axis of the infinitesimal element, and w is the out-of-plane displacement of the neutral plane. Each component of the strain tensor is written as follows,

ε αα = ∂ψ β

∂ ˆ x α + ∂ψ β

∂ψ α

ψ β + ∂ w

! ˆ x 3 , ε α 3 = 1

! ˆ x 3

∂ ˆ x α ˆ x 3 , ε αβ = 1

∂ ˆ x β

∂ ˆ x α

2

2

However, ε 33 will be zero according to the definition of linear strain, which means Equation 2 will be in conflict with Assumption 4 if Poisson’s ratio is not zero. To satisfy Assumption 4, ε 33 in plate theory is defined as ε 33 = − ν 1 − ν ( ε 11 + ε 22 ) ,

where ν denotes to Poisson’s ratio. Inserting the kinematic assumptions into Equation 4 , the weak form of bending wave equation for PML in Mindlin plate can be obtained by performing the integral over the plate thickness. Firstly, the sti ff ness matrix can be obtained by performing the integral over the surface domain as follows,

T 

∂δψ 2

∂ψ 2 ∂ ˆ x 1 ∂ψ 1 ∂ ˆ x 2 ∂ψ 1 ∂ ˆ x 1 + ∂ψ 2







 d ˆ S



B B ν 0

Z t / 2

∂ ˆ x 1 ∂δψ 1

i , j δε ij σ ij dˆ x 3 d ˆ S = Z

Z

X

B ν B 0

∂ ˆ x 2 ∂δψ 1

− t / 2

A

A

∂ ˆ x 1 + ∂δψ 2

0 0 B 1 − ν

(5)

∂ ˆ x 2

∂ ˆ x 2

2

T  ψ 2 + ∂ w

 δψ 2 + ∂δ w



 d ˆ S ,

Z

∂ ˆ x 1 δψ 1 + ∂δ w

∂ ˆ x 1 ψ 1 + ∂ w

A kt µ

∂ ˆ x 2

∂ ˆ x 2

where B = µ t 3

12(1 − ν ) and k denotes to shear correction factor. Similarly, the mass matrix can be obtained as follows, Z

T  ψ y ψ ′ x

Z t / 2

 δψ y δψ ′ x



 + Z

i ρδ u i u i dˆ x 3 d ˆ S = Z

X

A ρ L δ ww d ˆ S , (6)

A ρ J

− t / 2

A

where ρ J = ρ t 3 / 12 and ρ L = ρ t . Finally, d ˆ S can be written as follows according to Equation 3 ,

d ˆ S = dˆ x 1 dˆ x 2 = γ 1 ( x 1 ) γ 2 ( x 2 )d x 1 d x 2 .

During the implementation, γ i ( x i ) should be included in the surface integral. Although the analytical process is di ffi cult, the surface integral can be performed by Gauss-Legendre quadrature.

3. NUMERICAL EXPERIMENTS

3.1. Numerical Model The numerical model used to validate the proposed PML theory is a square plate whose length of sides and thickness are 2 m and 0 . 02 m respectively. The plate has a density ρ of 2500 kg / m 3 , Poisson’s ratio of 0 . 22 and Young’s modulus of 78 GPa. The mesh of the plate has 40 × 40 elements which means all of the elements are squares whose lengths of sides are 50 mm. The numerical computation is performed at frequencies of 100, 500 and 1000 Hz. The function γ α can be defined artificially. Here, Equation 7 is used to define γ α according to Bermúdez’s definition [ 5 ].

γ α ( x α ) = 1 + ( jk b L ( x α )) , (7)

where k b denotes to the bending wavenumber and L ( x α ) denotes to the distance from point x α to the end of PML. The numerical results are compared with the analytical solutions of the point excited infinite Kirchho ff plate and the solutions [ 6 ] are written as follows

w a ( r , t ) = w 0 h H 2 0 ( kr ) − H 2 0 ( − ikr ) i · e i ω t , w 0 = F − 8 ik 2 D , (8)

where r denotes to the length from the excitation point, k denotes to the wavenumber, H 2 0 denotes to the Hankel function of the second kind, F denotes to the excitation force and D denotes to the bending sti ff ness. To avoid the shear locking in finite element analysis, the Mixed Interpolation of Tensorial Components (MITC) method [ 7 ] is implemented in the following computations.

3.2. Validation of Directivity First of all, as a validation of the proposed PML theory and implementation method, a qualitative observation is made on the displacement distribution. Figure 1 shows the numerical results of displacement fields. It can be clearly seen that spherical wavefronts propagate from excitation points at all frequencies. The distribution is slightly distorted as the bending wave spreads away from the vibration point at low frequencies 100 Hz, which is considered to be due to the numerical reflection from PML.

45°

27°

-8.3 14 -5 0 5 10 displacement Z [10 -9 m]

-4.9 7.3 -2 0 2 4 6 displacement Z [10 -8 m]

-4.3 7.2 -2 0 2 4 displacement Z [10 -9 m]

Figure 1: Numerical results of displacement fields with PML at (a) 100 Hz, (b) 500 Hz, (c) 1000 Hz.

In addition, Figure 2 shows the comparison between numerical results and analytical solutions in three directions 0 ◦ , 27 ◦ and 45 ◦ . Since the numerical results are in good accordance with analytical

solutions in all directions, it can be confirmed that the reflected wave from PML is considerably small. From the observations above, it is validated that the proposed PML theory and implementation method can be used to simulate non-reflective ends for Mindlin plates.

10 -8

10 -8

10 -9

1.5

8

8

Analytical

0 ° 27 ° 45 °

6

6

1

4

4

Displacement [m]

Displacement [m]

Displacement [m]

0.5

2

2

0

0

0

-2

-2

-0.5

-4

-4

-6

-1

-6

0 0.5 1 1.5 Distance [m]

0 0.5 1 1.5 Distance [m]

0 0.5 1 1.5 Distance [m]

Figure 2: Numerical results of displacement fields with PML at (a) 100 Hz, (b) 500 Hz, (c) 1000 Hz.

3.3. Evaluation of Accuracy Then, as a quantitative measure, the relative error between the numerical results and analytical solutions are computed in the L 2 norm by

v tR R | w n − w a | 2 d x d y R R | w a | 2 d x d y , (9)

Relative Error =

where w n and w a are numerical results and analytical solutions respectively. Table 1 shows the relative errors of numerical results under di ff erent element numbers N for PML domain. The PML width d and the numbers of sample points in Gauss-Legendre quadrature are set to 0 . 25 m and 4 points respectively. It can be seen that the relative errors of numerical results become smaller as the element numbers increases, which indicates that the accuracy of the proposed PML implementation largely depends on the density of the mesh for PML domain.

Table 1: Relative erro rs of displacement fields with di ff erent element num bers for PML domain

Element numbers 100 Hz 500 Hz 1000 Hz

N = 2 54.5% 160% 38.8%

N = 4 51.5% 24.2% 12.5%

N = 8 17.5% 7.48% 14.2%

N = 16 8.34% 4.74% 14.0%

Then, the influence of the PML width d on the accuracy is studied by setting the element number N = 16. Since the element number is fixed here, which means the degree of freedom is also fixed, the

computational cost remains the same for di ff erent PML width. The sample points in Gauss-Legendre quadrature are set to 4 points. Table 2 shows the results. Although the relative errors tend to decrease as the PML is widened, some instability can also be also seen at 500 Hz. It is inferred that the accuracy may depend on the dimensions and elastic properties of the numerical model.

Table 2: Relative errors of displacement fields with di ff erent wid th of PML domain

PML width 100 Hz 500 Hz 1000 Hz

d = 0 . 125 m 8.66% 7.57% 13.4%

d = 0 . 25 m 8.34% 4.74% 14.0%

d = 0 . 5 m 6.91% 5.86% 14.0%

d = 1 m 2.94% 4.35% 13.5%

Since the function γ is included in the computation of element matrices, the numerical integral of Gauss-Legendre quadrature is an approximation. The influence of numbers of sample points in Gauss-Legendre quadrature is studied with fixed element numbers N = 16 and fixed PML width d = 1 m. Table 3 shows the results.

Table 3: Relative e rrors of displacement fields with di ff erent numbe rs of sample points

Sample points 100 Hz 500 Hz 1000 Hz

4 2.940% 4.352% 13.47%

9 2.911% 4.347% 13.49%

16 2.907% 4.342% 13.50%

(a) Physical: unstructured, PML: structured (b) Physical: structured, PML: unstructured

-8.9 14 -5 0 5 10 displacement Z [10 -9 m]

Figure 3: Numerical results of displacement fields with distorted quadrilateral meshes.

Theoretically, MITC elements can also be implemented with distorted quadrilateral meshes. The influence of the distorted elements is studied at the frequency of 500 Hz. Figure 3 (a) shows the numerical results when the physical plate is divided into distorted quadrilateral mesh and the PML

are divided into regular rectangle mesh. The relative error in this case is 5 . 91%, which is slightly larger than the case where there is no distortion in meshes. Figure 3 (b) shows the numerical results when the physical plate is divided into regular rectangle mesh and the PML are divided into distorted quadrilateral mesh. The relative error in this case increases to 9 . 10%, which is due to numerical reflections from PML.

4. CONCLUSION

In this paper, PML for bending waves in Mindlin plate theory is formulated and its implementation method in finite element analysis is introduced. The proposed PML theory and implementation method is validated through numerical experiments. It is found that PML width and element numbers have the most influence on the accuracies of the numerical results. To have deep understanding on the accuracies of the proposed PML theory, more numerical experiments taking elastic properties into consideration are needed. In the future, the proposed PML theory and implementation methods will be applied to the numerical simulation of floor impact sounds.

ACKNOWLEDGE

This work has been supported by JST SPRING, Grant Number JPMJSP2108.

REFERENCES

[1] J. Berenger. A perfectly matched layer for the absorption of electromagnetic waves. Journal of Computational Physics , 114(2):185–200, 1994. [2] W.C. Chew and Q.H. Liu. Perfectly matched layers for elastodynamics: A new absorbing boundary condition. Journal of Computational Acoustics , 04(04):341–359, 1996. [3] I. Harari and U. Albocher. Studies of fe / pml for exterior problems of time-harmonic elastic waves. Computer Methods in Applied Mechanics and Engineering , 195(29):3854–3879, 2006. [4] L. Nagler. Simulation of Sound Transmission through Poroeleastic Plate-like Structures . Technischen Universität Graz, 2011. [5] A. Bermúdez, L. Hervella-Nieto, A. Prieto, and R. Rodrı´guez. An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems. Journal of Computational Physics , 223(2):469–488, 2007. [6] A. Nilsson and B. Liu. Vibro-Acoustics, Volume 1 . Springer Berlin, Heidelberg, 2015. [7] K. Bathe and E. Dvorkin. A four-node plate bending element based on mindlin / reissner plate theory and a mixed interpolation. International Journal for Numerical Methods in Engineering , 21(2):367–383, 1985.