Analysis of Acoustic Wave in Homogeneous and Inhomogeneous Media Using Finite Element Method

Even though the propagation of elastic/acoustic waves in inhomogeneous and layered media has been an active research topic for many decades already, new problems and challenges continue to be posed even up to now. In fact, during the last few years, renewed interests have been witnessed by researchers in the various fields of acoustics, such as acoustic mirrors, filters, resonators, waveguides, and other kinds of acoustic devices, in relation to wave propagation in periodic elastic media. In acoustics and applied mechanics, these developments have been triggered by the need for new acoustic devices in order to obtain quality control of elastic/acoustic waves. What sort of material can allow us to have complete control over the elastic/acoustic wave’s propagation? We would like to discuss and answer this question in this chapter. It is well known that the successful applications of photonic band-gap materials have hastened the related researches on phononic band-gap materials. Analysis of Acoustic Wave in Homogeneous and Inhomogeneous Media Using Finite Element Method explores the theoretical road leading to the possible applications of phononic band gaps. It should quickly bring the elastic/acoustic professionals and engineers up to speed in this field of study where elastic/acoustic waves and solid-state physics meet. It will also provide an excellent overview to any course in elastic/acoustic media. Previous research on photonic crystals (Johnson & Joannopoulos, 2001, 2003; Joannopoulos et al., 1995; Leung & Liu, 1990) has sparked rapidly growing interest in the analogous acoustic effects of phononic crystals and periodic elastic structures. The various techniques for band structure calculations were introduced (Hussein, 2009). There are many wellknown methods of calculating the band structures of photonic and phononic crystals in addition to the reduced Bloch mode expansion method: the plane-wave expansion (PWE) method (Huang & Wu, 2005; Kushwaha et al., 1993; Laude et al., 2005; Tanaka & Tamura, 1998; Wu et al., 2004 ; Wu & Huang, 2004), the multiple-scattering theory (MST) (Leung & Qiu, 1993; Kafesaki & Economou, 1999; Psarobas & Stefanou, 2000; Wang et al., 1993), the finite-difference (FD) method (Garica-Pabloset et al., 2000; Sun & Wu, 2005; Yang, 1996), the transfer matrix method (Pendry & MacKinnon, 1992), the meshless method (Jun et al., 2003), the multiple multipole method (Moreno et al., 2002), the wavelet method (Checoury & Lourtioz, 2006; Yan & Wang, 2006), the pseudospectral method (Chiang et al., 2007), the finite element method (FEM) (Axmann & Kuchment, 1999; Dobson, 1999; Huang & Chen,


Theory
In this chapter, based on the theorems of solid-state physics and the finite element method with Bloch calculations, equation of motion of the acoustic modes in two-dimensional inhomogeneous media, phononic band structures, are derived and discussed in detail.In the beginning, the concepts of the real space and k space are introduced while the Brillouin zone is also addressed in the text.Generalized techniques of Bloch calculations in finite element method are used to analyze the acoustic modes in two-dimensional homogeneous and inhomogeneous media, phononic band structures, consisting of materials with general anisotropy.The mixed and transverse polarization modes and quasi-polarization modes are investigated in the text.
where ijk ε is the three-dimensional Levi-Civita completely antisymmetric symbol.The complete set of reciprocal lattice vectors is written as { } bbb , where N 1 , N 2 , and N 3 are integers.Figure 1 shows the primitive unit cell in two-dimensional real space while the Fig. 2 shows the relationship between the real space and k space.A property between the primitive lattice vectors and associated primitive reciprocal lattice vectors is 2 , where ij δ is the kronecker symbol.Note that the associated primitive reciprocal lattice vectors are constructed as k space from the concept of crystal diffraction.We will find that, in following sections, the discrete translational symmetry of a phononic crystal allows us to classify the elastic/acoustic waves with a wave vector k.The propagating modes can be written in "Bloch form," consisting of a plane wave modulated by a function that shares the periodicity of the lattice (Joannopoulos et al., 1995): . ii ee ⋅⋅ == + kr kr kkk P( r ) u( r ) u( r R ) (2) The important feature of the Bloch states is that different values of k do not necessarily lead to different modes.It is clear that a mode with wave vector k and a mode with wave vector k+G are the same mode, where G is a reciprocal lattice vector.The wave vector k serves to specify the phase relationship between the various cells that are described by u.If k is increased by G, then the phase between cells is increased by G⋅R, which we know is 2πn (n= l 1 N 1 +l 2 N 2 + l 3 N 3 is an integer) and not really a phase difference at all.So incrementing k by G results in the same physical mode.This means that we can restrict our attention to a finite zone in reciprocal space in which we cannot get from one part of the volume to another by adding any G.All values of k that lie outside of this zone, by definition, can be reached from within the zone by adding G, and are therefore redundant labels shown in Fig. 3.This zone is the so-called Brillouin zone.π (the lattice constant is a) need to be considered.The Fig. 4 shows the first, second, and third Brillouin zones.For more details, it is best to consult the first few chapters of a solid-state physics text, such as Kittel, 1996, or consult the appendix of popular photonic text like Joannopoulos et al. 1995and Johnson & Joannopoulos, 2001, 2003.

Equation of motion
This section provides a brief introduction of the theory of analyzing acoustic wave propagation in inhomogeneous media like as phononic band structures.The theory in this chapter can also be used to discuss acoustic wave propagation in homogeneous media because a homogeneous medium is symmetric with respect to any periodicity.In an inhomogeneous linear elastic medium with no body force, the equation of motion of the displacement vector (,) t ur can be written as where (,) (,,) zx y z == rx is the position vector, t is the time variable, and () ρ r and ( ) ijmn C r are the position-dependent mass density and elastic stiffness tensor, respectively.The following discussion considers a periodic structure consisting of a two-dimensional periodic array (x-y plane) of material A embedded in a background material B shown in Fig. 5.It is noted that when the properties of materials A and B tend to coincide, the homogeneous case is recovered.To calculate the dispersion diagrams of periodic structures, this study uses COMSOL Multiphysics software to apply the Bloch boundary condition to the unit cell domain in the FEM method.Based on the periodicity of phononic crystals, the displacement and stress components in the periodic structure are expressed as follows: (,) (,) ,  (Tanaka et al., 2000): where R is a lattice translation vector with components of 1 R and 2 R in the x and y directions.The relationships between the original variables (,) .
The Bloch calculations in this study record the variation of the displacements, stress fields, and eigen-frequencies as the wave vector increases.By using the FEM, the unit cell is meshed and divided into finite elements which connect by nodes, and is used to obtain the eigen-solutions and mechanical displacements.The types of finite elements used in this chapter are the default element types, Lagrange-quadratic, in COMSOL Multiphysics.In order to simulate the dispersion diagrams, the wave vectors are condensed inside the first Brillouin zone in the square lattice.According to the above theories, the results of dispersion relations in a band structure along the Γ−Χ−Μ−Γ are characterized and presented in the following sections.Similarly,Fig. 6(b) shows the irreducible part of the Brillouin zone of a rectangular lattice due to the geometric anisotropy, which is a rectangle with vertexes Γ , Χ , Μ , and Y , and the same as discussing the material anisotropy (Wu et al., 2004).The finite element method divides a unit cell with a three-dimensional model into finite elements connected by nodes.The FEM obtains the eigen-solutions and contours of a mode shape.To simulate the dispersion diagrams, the wave vectors are condensed inside the first Brillouin zone in the square and rectangular lattices.Using the theories above, the following section presents the results of dispersion relations in a band structure for the Γ−Χ−Μ−Γ square lattice or isotropic materials, and Y Γ−Χ−Μ− −Γ rectangular lattice or anisotropic materials.Note that the 2D FEM model calculates the dispersion relations of mixed polarization modes, while the 3D FEM model describes the dispersion relations of mixed and transverse polarization modes.

Acoustic wave in homogeneous media
It can be noted that a homogeneous medium is symmetric with respect to any periodicity, and it can be shown that the results for an infinite homogeneous medium can be cast in the form appropriate for a periodic medium.In this section, we introduce the mixed polarization modes and transverse polarization modes in a homogeneous medium.Displacement fields (polarizations) are also investigated and used to distinguish the different modes in the dispersion relations.The aluminum and quartz are adopted for examples and discussed in the section.The wave velocities of different propogating modes are also observed and discussed.

Isotropic medium
In Fig. 5, when the properties of materials A and B tend to coincide, the homogeneous case is recovered.Consider a periodic structure consisting of aluminum (Al) circular cylinders embedded in a background material of Al forming a two-dimensional square lattice with lattice spacing R. It means this is a homogeneous medium in a 3D FEM model.Figure 7 shows the dispersion relations along the boundaries of the irreducible part of the Brillouin zone Γ−Χ−Μ−Γ.The vertical axis is the frequency (Hz) and the horizontal axis is the reduced wave vector * / kk R π = .Here, k is the wave vector along the Brillouin zone.The Young's modulus E, Poisson's ratio ν , and density ρ of the material Al utilized in this example are E=70 GPa, ν =0.33, and ρ =2700 kg/m 3 .
As the elastic waves propagate along the x axis, the nonvanishing displacement fields of the shear horizontal mode (SH), shear vertical mode (SV), and longitudinal mode (L) are u y , u z , and u x respectively.It is noted that wave velocity ,, /2 * , so the slopes of dispersion curves in the Γ−Χ section of Fig. 7 are exactly the straight lines and can be explained as the wave velocities of shear (S) and longitudinal (L) modes.Here, m S,L are the slopes of shear and longitudinal modes in Fig. 7.It is noted that the wave velocities of shear horizontal mode and shear vertical mode are the same in an isotropic material.From the results in Fig. 7, the wave velocities of shear and longitudinal modes are 3119 and 6174 m/s.As we know, the wave velocities of shear and longitudinal modes in an isotropic material can be obtain from

Anisotropic medium
Similarly, the method in this chapter is used to discuss the wave velocities of acoustic modes in an anisotropic material.Consider a periodic structure consisting of quartz circular cylinders embedded in a background material of quartz forming a two-dimensional square lattice with lattice spacing R.This is also a homogeneous medium.The quartz is a piezoelectric and anisotropic material.The density ρ =2651 kg/m 3 .The elastic constants, piezoelectric constants, and relative permittivity of quartz utilized in this example are shown in Tables 1-3.The piezoelectric material, quartz, is a complete structural-electrical material, and thus all piezoelectric material properties were defined and entered into the FEM model.Figure 9 shows the dispersion relations along the boundaries of the irreducible part of the Brillouin zone Y Γ−Χ−Μ− −Γ due to the material anisotropy.In the calculations, the x-y plane is parallel to the (001) plane and the x axis is along the [100] direction of quartz.The vertical axis is the frequency in Hz unit and the horizontal axis is the reduced wave vector.From the discussion, it shows that the method adopted in this chapter can be used to discuss the wave propagations in isotropic and anisotropic media.

Acoustic wave in inhomogeneous media
Previous studies on photonic crystals raise the exciting topic of phononic crystals.This section presents the results of acoustic waves in inhomogeneous media, Al/Ni periodic structures and phononic crystals with reticular geometric structures.It also discusses the tunable band gaps in the acoustic waves of two-dimensional phononic crystals with reticular geometric structures using the 2D and 3D finite element methods.This section Analysis of Acoustic Wave in Homogeneous and Inhomogeneous Media Using Finite Element Method 13 calculates and discusses the band gap variations of the bulk modes due to different sizes of reticular geometric structures.Results show that adjusting the orientation of the reticular geometric structures can increase or decrease the total elastic band gaps for mixed polarization modes.

Periodic structure with two media
It is necessary and worthy to provide evidence supporting the FEM method's (COMSOL Multiphysics) ability to perform Bloch calculations with two media.This chapter compares the dispersion relations of Al/Ni band structure using the PWE method with the results of using the FEM method.Consider a phononic structure consisting of Al circular cylinders embedded in a background material of Ni to form a two-dimensional square lattice with lattice spacing R. Figure 11 shows the dispersion relations along the boundaries of the irreducible part of the Brillouin zone in Fig. 6 The diamond symbols represent the dispersion relations of the transverse polarization modes (shear vertical modes), and the cross symbols represent the mixed polarization modes (shear horizontal mode coupled with longitudinal mode) in the PWE method.The open circles represent the dispersion relations of all modes in the FEM method with a 3D model.The results of the FEM method match well with those of the PWE method.In the similar cases, when the differences of mass densities and elastic constants between the two periodic materials are larger, the convergence of the PWE method is slower and costs more CPU time.As the elastic waves propagate along the x axis, the nonvanishing displacement fields of the shear horizontal mode, shear vertical mode, and longitudinal mode are u y , u z , and u x respectively.For the sequence modes appear, the modes are always the same.When representing the whole wave vector space by the first Brillouin zone alone, they appear as further branches from higher Brillouin zones.In this example, the phase velocities of the SV 0 mode (diamond symbols) are larger than those of the SH 0 mode.The of the Brillouin zone X-M of Fig. 11 represents the dispersion of the bulk waves with propagating direction varied 0 deg~ 45 deg counterclockwise away from the x direction.

Periodic structure with single medium
Figure 12(a) depicts a two-dimensional phononic crystal with the reticular geometric structures of square lattice.These reticular structures are parallel to the z-axis.In a perfect two-dimensional phononic crystal, the periodic structure is constant in the z direction and the size of the structure is infinite in the x and y directions.To analyze the dispersion relations of all bulk acoustic modes in this band structure, the FEM should consider the 3D model in Fig. 12(c).The dimensions of the unit cell in Fig. 12(a) are c=d=0.8Rand R=h=1 in the calculations.Figure 13 shows the dispersion relations of the mixed and transverse polarization modes along the boundaries of the irreducible part of the Brillouin zone in Fig. 6(b) with the scales R=h=1, c=0.8, and a=1.2.The horizontal axis represents the reduced wave number along Y Γ−Χ−Μ− −Γ and the vertical axis represents the frequency (Hz).Note that this band structure shows no full band gap of the mixed and transverse polarization modes.Adopting the 2D FEM model to discuss the mixed polarization modes in this kind of band structure shows that there is only one full frequency band gap in Fig. 13, located at 3311 ~ 3400 Hz. Figure 13 compares the 3D and 2D FEM models.Open circles represent the dispersion relations of mixed polarization modes in the 2D FEM model, while solid circles represent the results of all bulk modes in the 3D FEM model.Figure 14 shows the eigenmode shapes with 4×4 supercell of total displacements for M 1 and M 2 modes indicated Fig. 13.These figures clearly show the phenomena of wave localizations in this reticular geometric structure.Note that the FEM method can easily describe the mode characteristics.In this chapter, M 1 is a shear horizontal mode with mode vibrating displacement along the y direction when the wave propagates along the x direction ( Γ−Χ direction).Also, M 2 is a shear vertical mode with mode vibrating displacement along z direction, and it does not couple with the mixed polarization modes.The following discussion addresses several parameters of the reticular geometric in this chapter.First, the effect of filling fraction is discussed when the parameters c=d varied from 0.1 to 0.9 in Fig. 12(a).Figure 15 shows the distribution of the total band gaps of mixed polarization modes, in which only one total band gap appears at approximately 3560 ~ 3736 Hz in c=d=0.8.The horizontal axis represents the parameter c, and the vertical axis represents frequency (Hz).Figure 15 also shows the 2D diagrams of the reticular geometric structures with c=d=0.1,0.5, and 0.8.Fig. 15.The band gap width with parameters c=d varying from 0.1 to 0.9 when the vertical range is selected from 3500 to 4500 Hz On the other hand, the scale a in Fig. 12(b) varies from 0.1 to 2.0 along the x direction and the width of the unit cell along y direction remains 1.0 in the Bloch calculations.Changing the scale a from 0.1 to 2.0 can tune the full frequency band gaps of mixed polarization modes.Using detailed calculations of dispersion relations of reticular geometric structures with scale a=0.1 to 2.0, Fig. 16 shows the band gap widths with the scale a from 0.1 to 2.0 when the vertical range ranges from 2400 to 5200 Hz.The horizontal axis ranges from 0 to 2.0, and the vertical axis represents frequency (Hz).No full frequency band gap exists when the scale a are 0.1, 0.2, 0.3, 0.4, 0.6, 0.7, 1.5, 1.6, and 1.7.These results clearly show that changing the scale a can increase or decrease the full frequency band gap.
It is noted that the unit cells with a=0.5 and 2.0 are the same in the Bloch calculations.However, the dispersion phenomena is similar except for the scalar of the eigenmode frequencies in the vertical axis of dispersion relations.In both cases, there is only one total band gap of the mixed polarization modes.The location of the band gap ranges from approximately 5009 to 5017.4 Hz with a=0.5, while that for a=2.0 ranges from approximately 2504.5 to 2508.7 Hz.

Conclusion
This chapter examines and discusses the acoustic waves in homogeneous medium and inhomogeneous medium, periodic structures with media and one medium with geometrical periodicity.The wave velocities of shear and longitudinal modes in an isotropic material and those of quasi-SV, quasi-SH, and quasi-L modes in an anisotropic material are obtained using the finite element method.This method also discusses the tunable frequency band gaps of bulk acoustic waves in two-dimensional phononic crystals with reticular geometric structures using the 2D and 3D finite element methods.This study adopts the finite element method to calculate dispersion relations, avoiding the numerical errors, Gibbs phenomenon, from the PWE method.Results show that changing the filling fraction, scale a, and the rotating angles of unit lattices in the reticular geometric structures can increase or decrease the elastic/acoustic band gaps.The effect discussed in this chapter can be utilized to enlarge the phononic band gap frequency and may enable the study of the frequency band gaps of elastic/acoustic modes in special phononic band structures.

a GFig. 3 .Fig. 4 .
Fig. 3.All values of k that lie outside of this zone, by definition, can be reached from within the zone by adding G

Fig. 5 .
Fig. 5. Periodic structures with square lattice.When the properties of materials A and B tend to coincide, the homogeneous case is recovered

Fig. 6 .
Fig. 6.Brillouin regions of the square and rectangular latticesThis chapter considers a periodic homogeneous medium with square lattice and phononic structures with square and rectangular lattices.These lattices consist of periodic structures that form two-dimensional lattices with lattice spacing R (square lattice) and lattice spacing aR (rectangular lattice).The term a is a scale from 0.1 to 2.0 in this chapter.The periodic structures are parallel to the z-axis.Figures6(a) and 6(b) illustrate the Brillouin regions of the square lattice and rectangular lattice, respectively.In the square lattice, Fig. 6(a) shows Fig. 7.The dispersion relations of homogeneous and isotropic material Al along the boundaries of the irreducible part of the Brillouin zone Γ−Χ−Μ−Γ Shown in Γ−Χ section of Fig. 9, the cross symbols represent the quasi shear horizontal (quasi-SH) mode.The square symbols represent the quasi shear vertical (quasi-SV) mode and the open circle symbols represent the quasi longitudinal (quasi-L) mode.The wave velocities of quasi-SH, quasi-SV, and quasi-L modes along x axis are 3306, 5116, and 5741 m/s.Similarly, The wave velocities of quasi-SH, quasi-SV, and quasi-L modes along y axis ( Y Γ− section) are 3922, 4311, and 6009 m/s respectively.Figure 10 also shows the vibration mode shapes of unit cell for quasi-SH, quasi-SV, and quasi-L modes in X point.The arrows shown in Fig. 10 are the polarizations.In this example, the quasi-longitudinal and quasi-transverse waves are almost indistinguishable from the truly longitudinal and truly transverse waves of Fig. 8.

Fig. 9 .
Fig. 9.The dispersion relations of homogeneous material quartz along the boundaries of the irreducible part of the Brillouin zone Y Γ−Χ−Μ− −Γ (a) with filling ratio 0.6.The vertical axis represents the normalized frequency * / t RC ωω = and the horizontal axis represents the reduced wave number * / kk R π = .Here, C t and k are the shear velocity of Ni and the wave vector along the Brillouin zone, respectively.The Young's modulus E, Poisson's ratio ν , and density ρ of the material Ni utilized in this example are E=214 GPa, ν =0.336, and ρ =8905 kg/m 3 .

Fig. 11 .
Fig. 11.Comparison of Bloch calculations between the PWE and FEM methods

Fig. 12 .
Fig. 12.(a) square lattice with lattice spacing R and (b) rectangular lattice with lattice spacing aR along x-axis and R along y-axis, (c) a unit cell with reticular structures in a 3D FEM model The material of the reticular structures in the unit cell in this chapter is aluminum.Figure 12(c) shows a diagram of the unit square lattice in a 3D FEM model.The periodicity of phononic crystals along the z direction is used to calculate the dispersion relations of the mixed and transverse polarization modes.The types of finite elements used for the 2D and 3D cases are the default element types, Lagrange-Quadratic, in COMSOL Multiphysics.Figure13shows the dispersion relations of the mixed and transverse polarization modes along the boundaries of the irreducible part of the Brillouin zone in Fig.6(b) with the scales R=h=1, c=0.8, and a=1.2.The horizontal axis represents the reduced wave number along Y Γ−Χ−Μ− −Γ and the vertical axis represents the frequency (Hz).Note that this band structure shows no full band gap of the mixed and transverse polarization modes.Adopting the 2D FEM model to discuss the mixed polarization modes in this kind of band structure shows that there is only one full frequency band gap in Fig.13, located at 3311 ~ 3400 Hz.Figure13compares the 3D and 2D FEM models.Open circles represent the dispersion

Fig. 13 .Fig. 14 .
Fig. 13.The dispersion relations of the mixed and transverse polarization modes along the boundaries of the irreducible part of the Brillouin zone with the scales R=h=1, c=0.8, and a=1.2

Fig. 16 .
Fig. 16.The band gap widths with the scale a from 0.1 to 2.0 Finally, the rotating angles of reticular geometric structures were changed to analyze the distribution of total band gaps.Figure 17 shows the 2D diagrams of unit rectangular lattices in different rotating angles D=30 deg, 45 deg, 75 deg, and 90 deg.In these cases, the widths of aluminum remain constant, 0.14R, in the reticular geometric structures with different rotating angles in the calculations.Figure 18 shows the band gap widths of rectangular lattices with different rotating angles of reticular geometric structures.Based on the symmetry of the geometry, the different angles in the Bloch calculations were adopted from 15 deg ~ 90 deg.In the calculated results, no band gap is detected from D=5 deg to 65 deg.

Table 3 .
The relative permittivity of quartz