Accelerated Approach for the Band Structures Calculation of Phononic Crystals by Finite Element Method

1 College of Mechanics and Materials, Hohai University, Nanjing 210098, China; yan.zhang@hhu.edu.cn (Y.Z.); lhjiang@hhu.edu.cn (L.-H.J.) 2 College of Civil Engineering, Tongji University, Shanghai 200092, China; changping.mei@163.com 3 College of Harbour, Coastal and Offshore Engineering, Hohai University, Nanjing 210098, China; chenda@hhu.edu.cn * Correspondence: lion_han@hhu.edu.cn; Fax: +86-25-8378-6046


Introduction
The propagation of elastic waves in phononic crystals (PnCs) has received considerable attention in the past decades [1][2][3][4].Most of the studies focus on the band gaps (BGs) in PnCs within which there could be no propagation of elastic waves in these specific frequency ranges.Many technical applications could be made possible owing to the understanding of the BGs property; acoustic filters [1], improvements in the design of transducers [2], noise control [3], vibration absorbers [5], hypersonic device [6], thermal conductance control [7], and optomechanical coupling control [8] are a few of those applications, which could be realized and improved by PnCs.Accurate and efficient modeling and simulation of BGs in PnCs are undoubtedly helpful in designing and manufacturing different devices.
Several methods have been developed in the attempts of calculating the BGs in PnCs.Among those, the transfer matrix (TM) method is the only analytical method, but it is restricted to one-dimensional PnCs [9,10].Other methods, including the plane wave expansion (PWE) method [2,11], finite difference time domain (FDTD) method [12,13], multiple scattering theory (MST) [14,15], finite element method (FEM) [16,17], and lumped-mass (LM) method [18,19].Despite the effort made by these studies, the convergence and computational efficiency are always highly mentioned due to the approximation of finite situation.The sparse discretization provides rough approximation while fine discretization causes better approximation but low computational efficiency.Thus, the models of relevant PnCs are usually somewhat simple for these limitations.Therefore, a method that could provide high precision results and high computational efficiency simultaneously is a better choice to study complex PnCs.
The idea that precise results might be obtained by combining LM and PWE methods is appealing because the results of the two methods converge to the exact solution from the opposite directions [18].However, this idea could unlikely to be implemented because: (1) convenient and efficient computational process cannot be easily realized as their mathematical calibrations are performed by different manners; and (2) it is difficult to control their convergence speeds with the same order.We think FEM is a potential method for solving the above-mentioned problems.The well-known supporting evidence for the finite element theory are: the consistent mass matrix underestimates the mass effects and gives higher Eigen frequencies (upper limit) than the exact ones, while the lumped mass matrix overestimates the mass effects and gives lower Eigen frequencies (lower limit) with similar mathematical calibrations [20].Therefore, it is reasonable to assume that the true solution of the Eigen frequencies should be located between the upper and lower limits after applying the Bloch's theorem.
Several different FEMs have been developed to unveil the relationship between the wave propagation properties and the periodicity of materials [21][22][23][24][25][26].With the help of softwares, such as ATILA [27][28][29], MSC.MARC [30], Comsol [31,32], etc., some convenient approaches with little programming could be achieved.Generally, the common process of BGs calculation in PnCs can be summarized as follows: First, establish the local governing equation in an element with the consistent or lumped mass matrices.Second, assemble the globe governing equation and introduce Bloch's theorem.Third, establish the global governing equation of Eigenvalue problem.Fourth, solve the Eigenvalue.However, the present computational process cannot calculate efficiently.On the one hand, the consistent mass matrix (or the lumped mass matrix) would underestimate (or overestimate) the mass effects.On the other hand, solving a large number of Eigenvalue problems consumes too much time.To overcome the above-mentioned disadvantages, the application of both consistent mass matrix and lumped mass matrix with specific weight function treatment, as well as the optimization of the introduction of Bloch's theorem, will be discussed in this paper.
The present research is organized as follows with the main purpose of establishing an accelerated finite element method for the band structures calculation in PnCs.First, the numerical modeling and conventional FEMs are introduced in detail.Then, the governing equations are formulated for each element through the application of the Bloch's theorem.After the consideration of boundary conditions, the generalized Eigenvalue equation is presented.Next, we give the weight function treatment, determine the weight by the analysis of several instances, and verify the method according to the analytical results.Finally, the computational efficiency of the method is discussed.

The Model
We use two-dimensional PnCs to illustrate the method.Two-dimensional PnCs have several forms, including infinitely and finitely long rods embedded in a matrix with different elastic properties and lattices.The cross-section shape of the scatterer could be arbitrarily selected.We consider a system consisting of infinitely long cylinders embedded in a homogeneous material.A cross section of such a system is shown in Figure 1.The lattice constant is a.We consider the z-axis to be parallel to the axis of the cylinders and the wave propagation in the x-y plane.Because of the periodicity, a cell contains all the properties of PnCs and we use the quadrilateral isoperimetric element with four nodes to mesh it.
If no damping and forcing terms exist, the governing equation for the motion of any point p in the homogeneous elastic medium is given by [33]: where ρ is the density; u p is the displacement vector; λ is the Lamé constant; µ is the shear modulus; and ∇ is the nabla differential operator.Considering the wave propagation in the x-y plane, the wave equation for the z component is decoupled from the equations for the x and y component (hereinafter referred to as the Z mode and XY mode, respectively).ρ ..
where u xy is the displacement in the x-y plane and u z is the displacement along the z-axis.Considering the wave propagation in the x-y plane, the wave equation for the z component is decoupled from the equations for the x and y component (hereinafter referred to as the Z mode and XY mode, respectively).
where uxy is the displacement in the x-y plane and uz is the displacement along the z-axis.

The Conventional FEM
According to the finite element theory [34], at any point in each element, the displacements can be expressed by the interpolation of nodal displacements.
  where ux, uy and uz are the displacements in x, y and z directions at any point in an element; u , e iy u , e iz u are the displacement components of the node i in x, y and z directions, respectively; and i = 1,2,3,4.
The equilibrium equation in an element is given as: e e e e    m u k u 0 (6) where m e and k e are the element mass matrix and element stiffness matrix and u e is the nodal displacement vector.m e can be a consistent or lumped mass matrix.The lumped mass matrix is generally adopted for software.The method of mass lumping or diagonalization is introduced in many books, and more details can be seen in Chapter 17 of reference [34].
Finally, using the well-known assembly principle, the equilibrium equation of the entire structure could be given.

The Conventional FEM
According to the finite element theory [34], at any point in each element, the displacements can be expressed by the interpolation of nodal displacements.
where u x , u y and u z are the displacements in x, y and z directions at any point in an element; where m e and k e are the element mass matrix and element stiffness matrix and u e is the nodal displacement vector.m e can be a consistent or lumped mass matrix.The lumped mass matrix is generally adopted for software.The method of mass lumping or diagonalization is introduced in many books, and more details can be seen in Chapter 17 of reference [34].
Finally, using the well-known assembly principle, the equilibrium equation of the entire structure could be given.

Application of the Bloch's Theorem in an Element
Generally, the Bloch's theorem is applied to the global matrices for the whole model.For the method that uses just one kind of mass matrix, the treatment is suitable.However, one can see that we will use two kinds of mass matrices below; therefor, the global mass matrices should be generated twice.Thus, the repetition of introducing the Bloch's theorem to the global matrices would lead to cumbersome calculation.To facilitate the calculation, we apply the Bloch's theorem in each element, which requires less computational operations, as well as physical memory.
We just give the derivation in XY mode and a similar process in Z mode can be easily achieved.Because of the periodic array of materials in PnCs, the system satisfies the Bloch's theorem [35].In addition, the free-vibration motion is simple harmonic [36].Thus, the displacements of any point p in the cell are given as: / -expi `kx x `ky y ´ωt ˘ (7) where v p is the amplitude vector, it has the same translational periodicity as the PnC; k = (k x ,k y ) is the wave vector; r = (x,y) is the coordinate; ω is the circular frequency; t is time.
Thus, in XY mode, the displacements of the four nodes in Equation ( 4) are given as: where , v e ix and v e iy are the amplitude components of node i in x and y directions, respectively.Substituting Equation (8) into Equation ( 4) gives the displacement of any point in an element as: where Substituting Equation ( 9) into Equation ( 6), the equilibrium equation in an element is given as: Evidently, the matrices k e T and m e T depend on the selected coordinate system in Equation ( 10).The following treatment can solve the problem.

, we get
Tk e Tv e " ω 2 Tm e Tv e (11) or k e v e " ω 2 m e v e ( 12) Crystals 2016, 6, 11 where k e and m e are the new element stiffness and mass matrices, respectively.They are both Hermitian positive definite matrices.

The Boundary Treatment
The formation and settlement of the Eigenvalue problem also depend on the treatment of boundary conditions.A sparse mesh of the cell in Figure 1b is illustrated in Figure 2. A, B, C and D are the four corner points of the cell.P is a mesh node located on the left boundary, and Q is the mesh node on the right boundary, which is one period away from P. Mesh nodes M and N have a similar situation as Q and P. f the basis vectors of the cell are a1 and a2, due to the periodicity of v in Equation ( 8), we g a a e m and n are both integers.hus, for the points shown in Figure 2, we have the following rules.
e the subscript denotes the corresponding mesh node.hus, the degree of freedom (DOF) information of P is the same as that of Q.Similarly, information of M is the same as that of N. A, B, C and D all have the same DOF informa the global matrix is assembling, for example, the DOF numbers of Q should be set as tho at is the boundary treatment method.In order to apply the method, the mesh must main dicity on the boundary.Uniform distribution of the same number of nodes on the opp is the easiest way.
he Generalized Eigenvalue Problem onsidering the same assembly principle of global matrix as the conventional FEM and -mentioned boundary treatment rules, the generalized Eigenvalue problem could be give If the basis vectors of the cell are a 1 and a 2 , due to the periodicity of v in Equation ( 8), we get v prq " v pr `ma 1 `na 2 q (13) where m and n are both integers.Thus, for the points shown in Figure 2, we have the following rules.
where the subscript denotes the corresponding mesh node.Thus, the degree of freedom (DOF) information of P is the same as that of Q.Similarly, The DOF information of M is the same as that of N. A, B, C and D all have the same DOF information.When the global matrix is assembling, for example, the DOF numbers of Q should be set as those of P. That is the boundary treatment method.In order to apply the method, the mesh must maintain periodicity on the boundary.Uniform distribution of the same number of nodes on the opposite sides is the easiest way.

The Generalized Eigenvalue Problem
Considering the same assembly principle of global matrix as the conventional FEM and the above-mentioned boundary treatment rules, the generalized Eigenvalue problem could be given as: where K and M are the new global stiffness matrix and global mass matrix, respectively, and they are both Hermitian positive definite matrices.v is the global amplitude vector.
Because of the property of the Hermitian positive definite matrix [37], Equation ( 15) has only real Eigenvalues, and could be transformed into a real symmetric generalized Eigenvalue problem.Let where the subscripts R and I denote the real and imaginary parts of a complex number, respectively.Substituting Equation ( 16) into Equation (15) gives: Solving the Eigenvalue problem, the frequency dispersion relation of wave propagation could be given.Considering the Bloch's theorem and symmetrical feature of cell, the wave vectors, captured along the edge of the irreducible Brillouin zone, could be enough to determine the frequency ranges of BGs.

The Weight Function Treatment
The element mass matrix in the above process could be the consistent mass matrix or the lumped mass matrix, corresponding to the CM-FEM and LM-FEM, respectively.According to the finite element theory, the consistent mass matrix gives higher Eigen frequencies than the exact ones, while the lumped mass matrix causes lower Eigen frequencies.
For a given wave vector, one can get two circular frequencies ω a and ω b corresponding to the generalized Eigenvalue problems using the consistent and lumped mass matrices, respectively.We prefer to use the natural frequencies, which are given as f a = ω a /2π and f b = ω b /2π.It is assumed that when the mesh approaches the limiting case, f a Ñf and f b Ñf, where f is the exact solution.
Due to the exact solution must be located between f a and f b , it is reasonable to give the exact solution f as: f " q f a `p1 ´qq f b (18) where q is the weight, q P r0, 1s.
Choosing the suitable weight, Equation (18) could give the higher accuracy solution than f a or f b , even in the case of rather sparse mesh.

Determination of the Weight
We use three PnCs examples, including square lattice of cylindrical inclusions in a matrix, square lattice of square prism inclusions in a matrix and triangular lattice of cylindrical inclusions in a matrix, which are shown in Figures 3 and 4. The lattice constant are all 0.02 m, the diameter of cylindrical inclusion are all 0.012 m, and the side length of the cross section of square prism inclusion is also 0.012 m.We consider the cases with lead inclusions and vulcanized rubber matrix.The material parameters [38] of lead and rubber employed in the calculations are λ Lead = 4.23 ˆ10 10 Pa, µ Lead = 1.49ˆ10 10 Pa, ρ Lead = 11,600 kg/m 3 and λ Rubber = 5.44 ˆ10 6 Pa, µ Rubber = 3.4 ˆ10 5 Pa, ρ Rubber = 1300 kg/m 3 .In a recent paper [39], the elastic modulus of silicon rubber was studied and the P-wave modulus has been given as a value of the order of GPa.However, the vulcanized rubber has never been done the similar test, thus we still use the material parameters in [38].The choose of material parameters is also based on the consideration of possibly high contrast between scatterers and matrix, in order to give a severe test condition.In a recent paper [39], the elastic modulus of silicon rubber was studied and the P-wave modulus has been given as a value of the order of GPa.However, the vulcanized rubber has never been done the similar test, thus we still use the material parameters in [38].The choose of material parameters is also based on the consideration of possibly high contrast between scatterers and matrix, in order to give a severe test condition.We calculate the band structures of the three cases by CM-FEM and LM-FEM, respectively, in order to give the results under the consistent and lumped mass matrices.The frequency dispersion relations of wave propagation in XY and Z modes are shown in Figures 5-7.The results in Figures 5-7, respectively, correspond to the meshes that contains 393 nodes, 379 nodes and 367 nodes.In these cases with sparse mesh, one can notice that although an agreement is found among the CM-FEM and LM-FEM results, but the difference is becoming larger as the increase of frequency.One can even see the apparent differences of the results locate in the fourth, fifth and sixth curves.Thus, FEM with only one kind of mass matrix cannot give high accuracy results in the relatively high frequency range, especially for the sparse mesh.Therefore the band structure calculation usually requires rather fine mesh, in order to obtain precise enough results, while the computational efficiency cannot be considered simultaneously.We calculate the band structures of the three cases by CM-FEM and LM-FEM, respectively, in order to give the results under the consistent and lumped mass matrices.The frequency dispersion relations of wave propagation in XY and Z modes are shown in Figures 5-7.The results in Figures 5-7 respectively, correspond to the meshes that contains 393 nodes, 379 nodes and 367 nodes.In these cases with sparse mesh, one can notice that although an agreement is found among the CM-FEM and LM-FEM results, but the difference is becoming larger as the increase of frequency.One can even see the apparent differences of the results locate in the fourth, fifth and sixth curves.Thus, FEM with only one kind of mass matrix cannot give high accuracy results in the relatively high frequency range, especially for the sparse mesh.Therefore the band structure calculation usually requires rather fine mesh, in order to obtain precise enough results, while the computational efficiency cannot be considered simultaneously.We calculate the band structures of the three cases by CM-FEM and LM-FEM, respectively, in order to give the results under the consistent and lumped mass matrices.The frequency dispersion relations of wave propagation in XY and Z modes are shown in Figures 5-7.The results in Figures 5-7, respectively, correspond to the meshes that contains 393 nodes, 379 nodes and 367 nodes.In these cases with sparse mesh, one can notice that although an agreement is found among the CM-FEM and LM-FEM results, but the difference is becoming larger as the increase of frequency.One can even see the apparent differences of the results locate in the fourth, fifth and sixth curves.Thus, FEM with only one kind of mass matrix cannot give high accuracy results in the relatively high frequency range, especially for the sparse mesh.Therefore the band structure calculation usually requires rather fine mesh, in order to obtain precise enough results, while the computational efficiency cannot be considered simultaneously.For the determination of the weight, we need the convergence behavior of results.We still use the above-mentioned PnCs examples but consider different mesh levels.The convergence of the results given by CM-FEM and LM-FEM are compared in Figures 8-10.The behaviors of the nine results of the Z modes at the M point (points M1, M2, M3, M4, M5 and M6 in Figures 5b and 6b) and X point (points X1, X2 and X3 in Figure 7b) are shown as functions of the number of nodes employed in a cell with CM-FEM and LM-FEM.The dash lines illustrate the exact results that CM-FEM and LM-FEM results finally both approach.The convergence tendencies verify that the CM-FEM results indeed converge to the exact solution from the direction opposite to the LM-FEM results.Furthermore, clearly, the convergence speed of both CM-FEM and LM-FEM results are very close.Therefore q = 0.5 is a reasonable weight.So we give an easy estimate of exact solution as Figures 8-10 also show the results given by Equation (19).Evidently, the new results have much faster convergence speed than that of the other two methods.
Actually, in the conventional FEM, using the proper element type, such as the quadrilateral element, the consistent mass formulation and the lumped mass formulation have the same order of convergence [40].From the numerical results, we can see that this feature still exists in the FEM for the band structure calculation of PnCs.The same order of convergence of the consistent and lumped mass formulations gives the equal weight.4b for a triangular lattice of lead columns with a cylindrical cross section embedded in vulcanized rubber.The lattice constant is 0.02 m and the diameter of the cylindrical inclusion is 0.012 m.
For the determination of the weight, we need the convergence behavior of results.We still use the above-mentioned PnCs examples but consider different mesh levels.The convergence of the results given by CM-FEM and LM-FEM are compared in Figures 8-10.The behaviors of the nine results of the Z modes at the M point (points M 1 , M 2 , M 3 , M 4 , M 5 and M 6 in Figures 5b and 6b) and X point (points X 1 , X 2 and X 3 in Figure 7b) are shown as functions of the number of nodes employed in a cell with CM-FEM and LM-FEM.The dash lines illustrate the exact results that CM-FEM and LM-FEM results finally both approach.The convergence tendencies verify that the CM-FEM results indeed converge to the exact solution from the direction opposite to the LM-FEM results.Furthermore, clearly, the convergence speed of both CM-FEM and LM-FEM results are very close.Therefore q = 0.5 is a reasonable weight.So we give an easy estimate of exact solution as Figures 8-10 also show the results given by Equation (19).Evidently, the new results have much faster convergence speed than that of the other two methods.
Actually, in the conventional FEM, using the proper element type, such as the quadrilateral element, the consistent mass formulation and the lumped mass formulation have the same order of convergence [40].From the numerical results, we can see that this feature still exists in the FEM for the Crystals 2016, 6, 11 9 of 14 band structure calculation of PnCs.The same order of convergence of the consistent and lumped mass formulations gives the equal weight.
Figures 8-10 also show the results given by Equation (19).Evidently, the new results have much faster convergence speed than that of the other two methods.
Actually, in the conventional FEM, using the proper element type, such as the quadrilateral element, the consistent mass formulation and the lumped mass formulation have the same order of convergence [40].From the numerical results, we can see that this feature still exists in the FEM for the band structure calculation of PnCs.The same order of convergence of the consistent and lumped mass formulations gives the equal weight.

Analytical Solution of Z Mode Wave Propagation
Although the dash lines in Figures 8-10 illustrate the exact results according to the convergence tendencies of CM-FEM and LM-FEM, one should note that the analytical solution of the band structure of PnCs that consists of different materials cannot be obtained, but that of the unique material can be obtained.Thus, the analytical results of a unique material can be the references.Here, we just consider the Z mode propagation of elastic waves in a unique material.The unique material is obtained from replacing different materials in PnC by one material.Then,

Analytical Solution of Z Mode Wave Propagation
Although the dash lines in Figures 8-10 illustrate the exact results according to the convergence tendencies of CM-FEM and LM-FEM, one should note that the analytical solution of the band structure of PnCs that consists of different materials cannot be obtained, but that of the unique material can be obtained.Thus, the analytical results of a unique material can be the references.Here, we just consider the Z mode propagation of elastic waves in a unique material.The unique material is obtained from replacing different materials in PnC by one material.Then,

Analytical Solution of Z Mode Wave Propagation
Although the dash lines in Figures 8-10 illustrate the exact results according to the convergence tendencies of CM-FEM and LM-FEM, one should note that the analytical solution of the band structure of PnCs that consists of different materials cannot be obtained, but that of the unique material can be obtained.Thus, the analytical results of a unique material can be the references.Here, we just consider the Z mode propagation of elastic waves in a unique material.The unique material is obtained from replacing different materials in PnC by one material.Then, we can surely consider the unique material also has the same periodicity and reciprocal lattice vectors as the original PnC.
According to the governing equation of motion using the harmonic solution u z pr, tq " v z prq exp p´iωtq (21) the Eigenvalue problem is given as: According to the Bloch's theorem, where w(r) has the same periodicity as the original PC.Clearly, exp(iGr) can be the solution of w(r), where G is the reciprocal lattice vector.Finally, the frequency dispersion relation is given as: We consider the same square lattice of cylindrical inclusions in the matrix, which is shown in Figure 3a.The matrix and inclusions are both set as lead.The computation using the presented method is based on the mesh of a cell in Figure 3a, which contains 934 nodes.The comparison of the frequency dispersion relations obtained by Equation (24) and the presented method is shown in Figure 11.Evidently, the results given by the presented method have a good agreement with the analytical results.The high frequency results still have high accuracy.The maximum relative errors of the two values corresponding to M point in Figure 11 where w(r) has the same periodicity as the original PC.Clearly, exp(iGr) can be the solution of w(r), where G is the reciprocal lattice vector.Finally, the frequency dispersion relation is given as: We consider the same square lattice of cylindrical inclusions in the matrix, which is shown in Figure 3a.The matrix and inclusions are both set as lead.The computation using the presented method is based on the mesh of a cell in Figure 3a, which contains 934 nodes.The comparison of the frequency dispersion relations obtained by Equation (24) and the presented method is shown in Figure 11.Evidently, the results given by the presented method have a good agreement with the analytical results.The high frequency results still have high accuracy.The maximum relative errors of the two values corresponding to M point in Figure 11 are just 0.04% and 0.19%.

Computational Efficiency
Besides, we discuss the computational efficiency.The presented method needs to solve two generalized Eigenvalue problems for one wave vector, and nearly consumes double time.However, the convergence speed of the presented method is much faster than those of CM-FEM or LM-FEM,

Computational Efficiency
Besides, we discuss the computational efficiency.The presented method needs to solve two generalized Eigenvalue problems for one wave vector, and nearly consumes double time.However, the convergence speed of the presented method is much faster than those of CM-FEM or LM-FEM, so a sparser mesh could give even the higher precision results.From Figures 8-10 one can clearly see that the results of the presented method by employing no more than 400 nodes are fairly close to the final convergence results, and they have even higher precision than the CM-FEM or LM-FEM results given under the mesh with more than 900 nodes.We give the comparison of the time of calculating all Eigenvalues in both XY and Z modes in Tables 1-3.The generalized Eigenvalue problems are solved with the help of international mathematics and statistics library (IMSL).By using the presented method, it achieves higher accuracy results and nearly 90% reduction in computational time, compared with CM-FEM and LM-FEM.Evidently, the presented method has apparent advantages.

Extension To Three Dimensional Problem
Finally, we discuss the band structures of sonic wave propagation in water.The governing equation of sonic wave has a similar expression to the governing equation of Z mode wave propagation above.We first build a cubic lattice cell with a sphere inclusion embedded in matrix.Then, the materials of matrix and inclusion are both set as water.A sparse mesh with tetrahedron element and 266 nodes is shown in Figure 12.As discussed above, this case has analytical solution ω 2 " c 2 0 pG `kq 2 , while the calculations of the present method, CM-FEM and LM-FEM are performed simultaneously.The results are plotted in Figure 13.Clearly, the present method has better agreement with analytical solutions in wide frequency range.One interesting phenomena is observed that weight function q is kept constant (0.5) with tetrahedron element.With the same mesh, the maximum relative error by the present method can be neglected while the maximum relative error by CM-FEM or LM-FEM exceeds 20%.
  0 c    G k , while the calculations of the present method, CM-FEM and LM-FEM are performed simultaneously.The results are plotted in Figure 13.Clearly, the present method has better agreement with analytical solutions in wide frequency range.One interesting phenomena is observed that weight function q is kept constant (0.5) with tetrahedron element.With the same mesh, the maximum relative error by the present method can be neglected while the maximum relative error by CM-FEM or LM-FEM exceeds 20%.13.Clearly, the present method has better agreement with analytical solutions in wide frequency range.One interesting phenomena is observed that weight function q is kept constant (0.5) with tetrahedron element.With the same mesh, the maximum relative error by the present method can be neglected while the maximum relative error by CM-FEM or LM-FEM exceeds 20%.

Conclusions
The accelerated FEM we introduced in the band structures calculation is confirmed to have a higher convergence speed and computational efficiency compared with other FEMs for similar situations.The computational time can be reduced using the relative sparse mesh because of the fast convergence of the presented method.Besides, the presented method is shown to be suitable for cases with very high contrast in the material parameters between scatterers and matrix, such as lead and vulcanized rubber.
The accelerated FEM is mainly based on giving suitable weights to the results under the consistent and lumped mass matrices.Using the simple element type, such as the quadrilateral isoperimetric element with four nodes utilized in this study, the consistent and lumped mass formulations could have the same order of convergence.So we choose the weight as 0.5 and the quadrilateral isoperimetric element is recommended in the band structures calculation of 2D PnCs.In the case of the above-mentioned 3D problem, the weight of 0.5 is also selected.Furthermore, the introduction of the Bloch's theorem to the local governing equation improves the computational efficiency as well.
In this paper, the accelerated method is mainly illustrated and discussed for two-dimensional solid/solid PnCs.Like other FEMs, this method is not limited to these cases.It can also be used for one-and three-dimensional PnCs, as well as the drilled solid PnCs, fluid/solid and fluid/fluid PnCs.Especially, the three-dimensional PnCs with complex geometry and material properties, as well as the cases of the supercell which usually lead to the complex meshing with heavy computational job [41].These problems could be solved due to the present method with a sparse mesh.Considering the CM-FEM/LM-FEM are widely applied in software, while the presented average treatment and weighing values fixed as 0.5 can be easily implemented, the presented method provides a simple, efficient and workable choice for BGs computation in PnCs.To perform the presented method, the periodical mesh condition at the facing boundaries is obligated.

Figure 1 .
Figure 1.(a) A cross section of the PnC consisting of infinitely long cylinders embedded in a homogeneous material; (b) a cell of the PnC; and (c) the meshed cell.

Figure 1 .
Figure 1.(a) A cross section of the PnC consisting of infinitely long cylinders embedded in a homogeneous material; (b) a cell of the PnC; and (c) the meshed cell.

2 .
m are the new element stiffness and mass matrices, respectively.They are itian positive definite matrices.heBoundary Treatment he formation and settlement of the Eigenvalue problem also depend on the treatme dary conditions.A sparse mesh of the cell in Figure1bis illustrated in Figure2.A, B, C a e four corner points of the cell.P is a mesh node located on the left boundary, and Q i node on the right boundary, which is one period away from P. Mesh nodes M and N ha r situation as Q and P. The sparse mesh of a cell.A, B, C and D are four corner points.P is a mesh node located on he left boundary, and Q is the mesh node on the right boundary which is one period away from P. esh nodes M and N have a similar situation.

Figure 2 .
Figure 2. The sparse mesh of a cell.A, B, C and D are four corner points.P is a mesh node located on the left boundary, and Q is the mesh node on the right boundary which is one period away from P. Mesh nodes M and N have a similar situation.

Figure 3 .
Figure 3. (a) The mesh of a cell with a square lattice of cylindrical inclusions in a matrix.(b) The mesh of a cell with a square lattice of square prism inclusions in a matrix.The corresponding irreducible Brillouin zone is shown as the shadow region in (c).

Figure 3 .Figure 4 .
Figure 3. (a) The mesh of a cell with a square lattice of cylindrical inclusions in a matrix.(b) The mesh of a cell with a square lattice of square prism inclusions in a matrix.The corresponding irreducible Brillouin zone is shown as the shadow region in (c).Crystals 2016, 6, 11 7 of 14

Figure 5 . 6 Figure 4 .
Figure 5. Phononic band structures in (a) XY mode and (b) Z mode along the edge of the irreducible Brillouin zone shown in Figure3cfor a square lattice of lead columns with a cylindrical cross section embedded in vulcanized rubber.The lattice constant is 0.02 m and the diameter of the cylindrical inclusion is 0.012 m.

Figure 4 .
Figure 4. (a) The mesh of a cell with a triangular lattice of cylindrical inclusions in a matrix.The corresponding irreducible Brillouin zone is shown as the shadow region in (b).

Figure 5 . 6 Figure 5 .Figure 5 .Figure 6 .
Figure 5. Phononic band structures in (a) XY mode and (b) Z mode along the edge of the irreducible Brillouin zone shown in Figure3cfor a square lattice of lead columns with a cylindrical cross section embedded in vulcanized rubber.The lattice constant is 0.02 m and the diameter of the cylindrical inclusion is 0.012 m.

Figure 6 .Figure 7 .
Figure 6.Phononic band structures in (a) XY mode and (b) Z mode along the edge of the irreducible Brillouin zone shown in Figure 3c for a square lattice of lead columns with a square cross section embedded in vulcanized rubber.The lattice constant is 0.02 m and the side length of the cross section of inclusion is 0.012 m.Crystals 2016, 6, 11 8 of 14

Figure 7 .
Figure 7. Phononic band structures in (a) XY mode and (b) Z mode along the edge of the irreducible Brillouin zone shown in Figure4bfor a triangular lattice of lead columns with a cylindrical cross section embedded in vulcanized rubber.The lattice constant is 0.02 m and the diameter of the cylindrical inclusion is 0.012 m.

Figure 8 .
Figure 8. Convergence results corresponding to the points M1, M2 and M3 in Figure 5b.Figure 8. Convergence of the three results corresponding to the points M 1 , M 2 and M 3 in Figure 5b.

Figure 8 . 14 Figure 9 .
Figure 8. Convergence results corresponding to the points M1, M2 and M3 in Figure 5b.Figure 8. Convergence of the three results corresponding to the points M 1 , M 2 and M 3 in Figure 5b.Crystals 2016, 6, 11 9 of 14

Figure 10 .
Figure 10.Convergence of the three results corresponding to the points X1, X2 and X3 in Figure 7b.

Figure 9 . 14 Figure 9 .
Figure 9. Convergence of the three results corresponding to the points M 4 , M 5 and M 6 in Figure 6b.

Figure 10 .
Figure 10.Convergence of the three results corresponding to the points X1, X2 and X3 in Figure 7b.

Figure 10 .
Figure10.Convergence of the three results corresponding to the points X 1 , X 2 and X 3 in Figure7b.

Figure 11 .
Figure 11.The comparison of the frequency dispersion relations of the Z mode propagation of elastic waves in lead.The wave vectors are captured along the edge of the irreducible Brillouin zone shown in Figure 3c.

Figure 11 .
Figure 11.The comparison of the frequency dispersion relations of the Z mode propagation of elastic waves in lead.The wave vectors are captured along the edge of the irreducible Brillouin zone shown in Figure 3c.

Figure 12 .
Figure 12.The (a) appearance and (b) inner part of the mesh of the cubic lattice cell with a sphere inclusion embedded in matrix.

Figure 13 .
Figure 13.The comparison of the frequency dispersion relations of the sonic wave propagation in water.

Figure 12 .
Figure 12.The (a) appearance and (b) inner part of the mesh of the cubic lattice cell with a sphere inclusion embedded in matrix.

Figure 12 .
Figure 12.The (a) appearance and (b) inner part of the mesh of the cubic lattice cell with a sphere inclusion embedded in matrix.

Figure 13 .
Figure 13.The comparison of the frequency dispersion relations of the sonic wave propagation in water.Figure 13.The comparison of the frequency dispersion relations of the sonic wave propagation in water.

Figure 13 .
Figure 13.The comparison of the frequency dispersion relations of the sonic wave propagation in water.Figure 13.The comparison of the frequency dispersion relations of the sonic wave propagation in water.

Table 1 .
Comparison of the computational time in the case of square lattice of cylindrical inclusions in the matrix (s).

Table 2 .
Comparison of the computational time in the case of square lattice of square prism inclusions in the matrix (s).

Table 3 .
Comparison of the computational time in the case of triangular lattice of cylindrical inclusions in the matrix (s).