Sound Propagation Analysis on Sonic Crystal Elastic Structures using the Method of Fundamental Solutions(MFS)

P.G.Santos,J.Carbajo,L.Godinhoand J.Ramis

1　Introduction

The study of periodic structures for sound attenuation,namely sonic crystal type,has been a topic of increasing interest in recent years.Sonic crystals get their name by analogy with ordered structures of semiconductor materials such as silicon crystals,whose feature of allowing certain energy waves to pass through and block others is transposed,in sonic crystals,into the capacity to prevent or limit the propagation of certain sound frequencies.Historically,it is generally considered that the first evidence that it was possible to achieve some effect of acoustic obstruc-tion using structures in periodic arrays was derived fortuitously from a sculptural element,in the gardens of the Fundación Juan March in Madrid[Martínez-Sala,Sancho,Sánchez,Gómez,Llinares and Meseguer(1995)].Since then,different aspects of the behavior of sonic crystals have been studied,some of which from a more theoretical point of view,such as the influence of point defects[Wu,Chen and Liu(2009)]or the formation of waveguides in which the sound propagates with low attenuation[Vasseur,Deymier,Djafari-Rouhani,Pennec and Hladky-Hennion(2008)].In the field of the practical uses of sonic crystals,one which may be regarded perhaps as the most promising is their use for the selective attenuation of sound,for example as traffic noise barriers[Sánchez-Pérez,Rubio,Martínez-Sala,Sánchez-Grandia and Gómez(2002)],since it could become a competitive solution when compared with classic noise barriers[Castiñeira-Ibáñez,Rubio,Romero-García,Sánchez-Pérez and García-Raffi(2012)].

Several methods have been proposed to analyze the sound propagation of sonic crystals,mainly based in 2D approaches.The Plane-Wave Expansion Method[Kushwaha and Halevi(1996);Cao,Hou and Liu(2004);Yuan,Lu and Antoine(2008)]and Wavelet-based Methods[Yan and Wang(2006)]are examples;however those methods could lead,in some situations,to inaccurate results,as they do not account for the transverse waves in the solid elements,that are responsible for the appearance of significant band gaps due to the transverse modes in the scatterers[Sigalas and García(2000)].Even more,in the case of solid scatterers embedded in a fluid medium,those methods,give inaccurate results,as they do not account for the solid- fluid interface interaction[Li,Wang,Zhang and Yu(2013)].Approaches based on the Finite Difference Time Domain Method[Sigalas and García(2000);Cao,Hou and Liu(2004)]have also been used,but they can also lead to inaccurate estimations since the interface interaction in solid- fluid structures is neglected[Li,Wang,Zhang and Yu(2013)].These limitations are currently overpassed by other developed methods as the Multi-Scattering Theory[Mei,Liu and Qiu(2005)],the Dirichlet-to-Neumann Map Method[Wu and Lu(2008),Li,Wang and Zhang(2011)],the Finite Element Method(FEM)[Wang,Wang and Su(2011)]or the Boundary Element Method(BEM)[Li,Wang,Zhang and Yu(2013),Li,Yue-Sheng,Chuanzeng and Gui-Lan(2013)],methods that accounts for the transverse wave propagation in the solid elements as well the solid- fluid interactions.The last two referred methods have some advantages in case of modelling complex geometries of the scatterers.

Other types of numerical strategies may be applicable for the analysis of sonic crystals.One class of numerical methods which may be of interest corresponds to the so-called meshless methods,which require neither domain nor boundary discretization.Techniques such as the MLPG(Meshless-Local-Petrov-Galerkin)[Atluri and Shen(2002)],Kansa’s Method[Kansa(1990a,1990b)]or the Method of Fundamental Solutions(MFS)can be incorporated in this group.The latter(MFS)is a boundary-only technique,in which the solution of a PDE within a given domain is obtained using a linear combination of the effects of virtual sources located outside this domain.Some initial reference papers by Fairweather and Karageorghis(1998),Golberg and Chen(1998),Fairweather,Karageorghis and Martin(2003),or Smyrlis and Karageorghis(2003)describe the MFS and some of its possible applications.More recently,further developments have been published,regarding the analysis of the accuracy and stability of the MFS,such as in the works of Tsai,Lin,Young and Atluri[2006],and also including the application of the MFS for the propagation of electromagnetic waves[Young and Ruan(2005)],for one-dimensional wave equations[Gu,Young and Fan(2009)]or for boundary identification problems[Marin and Karageorghis(2009)],among many others.

Recently,the application of the MFS to analyze a sonic crystal composed of circular cylinders was proposed by Martins,Godinho and Picado-Santos(2013).In that work,the authors addressed the application for traffic noise attenuation in a 2D approach,with clear advantages in problem discretization and computational cost when compared to the BEM and FEM approaches.In the formulation,the scatterers were considered as rigid elements,and numerical comparisons with the ones obtained with a BEM formulation(also considering the scatterers as rigid elements)indicated a good level of accuracy of the model.Later,comparisons between numerical results with the ones obtained in an experimental campaign for a reduced scale model(1:5)with different con figurations,also revealed a fine level of accuracy[Martins,Carbajo,Godinho,Mendes and Ramis(2013)].The present paper aims to extend this previously presented methodology to make it suitable for the analysis of non-rigid(elastic)shell cylinders embedded in a fluid medium,and so account for the influence of the scatterer’s stiffness as well as of the fluid-solid interactions on the sound propagation phenomena through the sonic crystal.In this new method,the problem is divided in a series of sub-regions,one of them being the outer region(host fluid),and the remaining regions being defined around each scatterer structure.Since the analytical solutions are known for each of those defined sub-regions(fundamental solutions),it becomes possible to establish a coupled model based on MFS,which accounts for the full interaction between the involved fluids and the solids that compose the scatterer structures,by just establishing the continuity of pressures and displacements along the boundaries connecting the subregions.

Summarily,the method will be presented and validated as follow: first the theoretical formulation in which the numerical analysis methodology is based will be presented;the influence of the model parameters on its accuracy will be evaluated by comparison with FEM in an exemplificative case;a set of results obtained by the proposed model will then be compared againsta previous MFS model developed by Martins,Godinho and Picado-Santos(2013),in which the scatterers are assumed to be rigid.At this point,different combinations of materials and shell thickness of the scatterers will be tested for two types of fluid,in order to vary the structure stiffness as well as the properties’contrast between solids and fluids.The comparison between results from these two models will often show the influence of the scatterer’s stiffness and of the fluid/solid interaction on the sonic crystal behavior.

2　Problem formulation

2.1　Governing equations

The numerical settling for the present work intends to simulate the propagation of sound waves in the presence of a number of elastic scatterers embedded in a fluid media,considering a 2D environment.The scatterers are considered to be cylindrical ring-shaped structures,eventually including a second fluid in their inner part.For that purpose,different governing equations must be considered for either the solid or the fluid parts of the propagation domain,corresponding to the vectorial or scalar wave equations,respectively.

For the case of an isotropic and linear elastic material,with mass density ρS,allowing shear waves to propagate with a velocity βSand compressional waves to propagate with a velocity αS,the vectorial wave equation can be written as:

in which the vector u refers to displacements,ω is the circular frequency and∇=(∂/∂x)ˆ andbeing the unit vectors along thexandydirections.

In the case of a fluid medium,with a mass density ρf,the corresponding scalar wave equation in the frequency domain is the well-known Helmholtz equation.For this type of problems,it can be written as:

in whichpis the pressure andkf= ω/αfis the wave number,with αfbeing the speed of sound in the fluid medium;for this scalar equation,∇2=?∂2/∂x2?+?∂2/∂y2?.For this case,the definition of a fundamental solution corresponding to the effect of a 2D sound source in a fluid medium is possible,and can be easily found in the literature,in the form:

2.2　Analytical solution for a single circular shell scatterer in a fluid medium

Following the work by Godinho,Tadeu and Branco(2004),it is possible to analytically calculate the pressure field within a fluid medium in the presence of an elastic ring shaped scatterer,embedded within this fluid.For this purpose,consider a circular shell solid structure,defined by the internal and external radii,rAandrB,and submerged in a homogenous fluid medium,as illustrated in Fig.1.

This structure is illuminated by a pressure source,placed in the exterior fluid medium,generating sound waves that propagate within the fluid reaching the surface of the embedded structure.Part of the incident energy is then reflected back to the exterior fluid medium,with the rest being transmitted into the solid material,where it propagates as body and guided waves.These waves continue to propagate and may eventually hit the inner surface of the structure,where similar phenomena occur.

The wavefield generated in the outer fluid medium(Fluid 1)depends both on the incident pressure waves and on those coming from the external surface of the shell.The latter propagate away from the cylindrical shell,and can be defined using the following displacement potential when a cylindrical coordinate system is centered on the axis of the circular cylindrical shell:

wherekf1= ω/α1with α1being the pressure wave velocity in the outer fluid,ρf1being the fluid density andrthe distance from the center of the scatterer to the calculation point.Inside the solid material of the shell,two distinct groups of waves exist,corresponding to inward travelling waves,generated at the external surface,and to outward travelling waves,generated at the internal surface of the shell.Each of these groups of waves can be represented using one dilatational and one shear potential:

Figure 1:Circular cylindrical shell structure submerged in a fluid medium

wherekf2= ω/α2with α2being the pressure wave velocity in the inner fluid and ρf2being the fluid density.The term(j=1,6)for each potential of the Eq.5,Eq.6,Eq.7 and Eq.8 are unknown coefficients to be determined by imposing the required boundary conditions.For this case,the necessary boundary conditions are the continuity of normal displacements and stresses,and null tangential stresses on the two solid- fluid interfaces.

Consider,now,that the incident field,in terms of displacement potential,generated by the acoustic source located at(x0,y0)can be de fined at a point(x,y)as:

It should be noted that the pressure field can be obtained from this potential as,leading to Eq.3.In order to establish the appropriate equation system,this incident field must also be expressed in terms of waves centered on the axis of the circular cylindrical shell structure.This can be achieved making use of Graf’s addition theorem,leading to the expression(in cylindrical coordinates):

wherer0is the distance from the source to the axis of the circular cylindrical shell;εnis 1 ifn=0 and 2 in the remaining cases.

Using the previously presented equation,a linear equation system,with 6 equations and 6 unknowns,can be assembled for each value ofn.After the solution of the equation system is computed,the unknown values(j=1,6)can be used to determine the final wave fields.For the outer fluid,the pressure field at a point(x,y)can be written as:

withr0being the distance from the source to the field point,andNbeing the number of terms required to attain convergence

2.3　Hybrid solution for multiple circular shell scatterers

To extend the applicability of the previously defined analytical solution to the case of multiple scatterers embedded in a fluid medium,the meshless Method of Fundamental Solutions will be used,based in the work of Godinho,Amado-Mendes and Pereira(2011).For that purpose,consider an in finite fluid medium and the presence of an arbitrary number of circular ring shaped structures made of elastic materials,and filled with a fluid material,as illustrated in Fig.2.

Consider that,in the presence ofNRscatterers,the problem is divided inNR+1 sub-regions,one of them being the outer region(host fluid medium),and each of the remainingNRsubregions corresponding to a limited region which includes one of the scatterers.This con figuration is also represented in Fig.2.

For each of the regions described above,the corresponding fundamental solutions are known,and can be written,in terms of acoustic pressure,by means of Eq.3 and Eq.11.Using those solutions,it becomes possible to establish a coupled model accounting for the full interaction between the involved fluids and the solids that compose the scatterers,by just establishing the continuity of pressures and particle velocities along the boundaries connecting the sub-regions that include scatterers

Figure 2:Schematic representation of the problem.

with the outer sub-region,corresponding to the host fluid.If the MFS is used,the acoustic field in the outer sub-region,containing the host in finite fluid,can be defined by considering a number of virtual sources,placed within each of the remaining sub-regions,and linearly combining their effects as:

whererepresents a field point of coordinateslocated in the host(external)fluid,is the position of the real source exciting the system,andis the position of each of theNVSjvirtual sources placed within sub-regionj.Additionally,are the fundamental solutions for the host fluid at a pointoriginated by a source positioned atand,respectivelly,which is indeed the same presented in Eq.3.The coefficientsaj,lcorrespond to the amplitudes of each of the virtual sources used to simulate the pressure field in the host medium,and are “a-priori”unknown.

If a point located within the outer fluid of thekthinner sub-region is considered,the pressure field can be written as:

To determine the unknown coefficientsaj,landbk,ldefined above,it becomes necessary to establish a system of linear equations,enforcing the continuity of pressures and normal particle velocities along each of theNRboundaries separating the outer sub-region from each internal sub-region.Assuming that the boundary conditions are enforced atNVSkcollocation points along thekthinterface(as illustrated in Fig.2),the continuity equations at themthcollocation pointxc,kmof that boundary can be written in terms of acoustic pressure as:

and of normal particle velocity as:

Writing these two equations for each collocation point,aNxNsystem of linear equations,can be built.Once this system of equations is solved,one may obtain the pressure at any point in the host fluid by applying Eq.12.

Considering,as an example,the case of a fluid with four embedded scatterers,each of them within a sub-region with an interface defined withMk(k=1,4)collocation points and incudingNVSkvirtual sources,the corresponding equation system

would assume the form:

in which each of the sub-matrices of the system matrix related to the host fluid can be de fined as:

with:

The sub-matrices related to each scatterer can be written as:

The unknown vector in Eq.16 is given by:

and the right-hand-side term is defined as:

An important point that should be highlighted concerning this formulation is that the coupling between sub-regions is enforced in fluid- fluid interfaces,at some distance from the interfaces between each scatterer and the host fluid.Thus,the coupling conditions are imposed in a region with smooth variations of the pressure,which greatly improves the global performance of the strategy.Additionally,since the interface between sub-regions is virtual,it can assume a smooth shape,such as that of a circle,which has been demonstrated in previous works to lead to very accurate results[Godinho,Tadeu,and Mendes(2007)].Since the fundamental solutions are computed analytically for both the host fluid and for the sub-regions containing the scatterers,high accuracy may be obtained in the calculations.

2.4　Time domain responses

Using the above defined method,all the responses are initially computed in the frequency domain,and additional steps are required in order to calculate time domain signals.For the purpose of computing time signals,in all presented results it is assumed that the source emits a Ricker pulse,which is defined in the frequency domain as:

where Ω = ωt0/2,Ais the amplitude,tsis the time when the maximum occurs and πt0is the dominant wavelet period.The choice of the Ricker pulse to simulate the temporal variation of the acoustic source is related to the fact that it decays rapidly both in time and frequency domains,reducing computational effort and allowing easier interpretation of the computed time series and synthetic waveforms After calculating the frequency domain responses for a full range of discrete frequencies,these responses must be scaled using the frequency spectrum defined in Eq.26,and after this an inverse Fast Fourier Transform(iFFT)must be applied to the result.In this process,it becomes necessary to make use of complex frequencies,defined as ωc= ω −iξ,with ξ =0.7∆ω where ∆ω represents the frequency increments,to reduce the aliasing effect that occurs due the frequency discretization.In the time domain,this shift is later taken into account by applying an exponential windoweξtto the response[Kausel and Roësset(1992)]Using this process,it is possible to adequately compute time signals within a time window defined byT=2π/∆ω.

3　Model verification

In order to analyze the proposed MFS model efficiency and accuracy,the results for the Insertion Loss obtained for a specific example are compared with the ones provided by the Finite Element Method(FEM),as no analytical solution is known for this kind of problem.The proposed example consists in a set of four scatterers embedded in a fluid medium,water with a density of 1000 kg/m3and allowing sound waves to propagate at 1500 m/s The scatterers are assumed to be made of PolyVinyl Chloride(PVC)with a density of 1400 Kg/m3a Young’s modulus of 3.0146 GPa and a Poisson’s ration of 0.40622.The scatterers,consisting of a ring shaped structure,have an external radius of 0.1 m and a thickness of 0.025 m,with their centres equally spaced of 0.4 m between them in a rectangular lattice con figuration.An acoustic source placed 3 m away,aligned with the middle of the scatterers illuminates the system and the response is evaluated as the average value obtained over a square grid of 25 receivers,positioned behind the scatterers as depicted in Fig.3a).

The response was evaluated in the frequency range from 10 Hz to 5000 Hz with a frequency step of 10 Hz.In the FEM model a very re fined mesh was used in this study in order to allow the analysis of higher frequencies(smaller wavelengths)of 5000 Hz as depicted in Fig.3b).

Figure 3:a)Con figuration of the problem;b)FEM mesh around the scatterers.

In order to evaluate the influence of the proposed model parameters in the accuracy of the results,namely the radius of the interface delimiting the sub-regions(R),the distance of the virtual sources to the interface(D)and the number of collocation points(n),a set of combinations of these parameters were defined:Rvaried between 1.1 and 1.2 times the external scatterers radius;the ratio betweenDandRvaried from 0.10 to 0.50 with increments of 0.10,andnassumed values of 10,15,20 and 30.Fig.4 presents the results of the Insertion Loss obtained by the reference FEM formulation and the MFS with the variations on those model parameters.

Analysing Fig.4 for all the cases,it becomes clear that increasing the number of collocation points improves the accuracy of the results,as it allows a better definition of the virtual interface,with more points in which continuity is enforced.For the considered maximum frequency(5000Hz)and fluid medium sound speed(1500m/s),if we assume that at least 5 to 6 collocation points per wavelength should be used(as in[Godinho,Amado-Mendes and Pereira(2011)]),it is concluded that at least 12 to 15 points are necessary to obtain more reasonable results,particularly when R=1.2.This can contribute to the fact that the results are not satisfactory forn=10,independently of the interface radius and of theR/Dratio(Fig.4a,b);this inaccuracy is also particularly detected by observing a false resonance peak that appears approximately at 4500 Hz in all the cases in which 10 collocation points are used(Fig.4a,b).In addition,for the smaller numbers of points,the distance of the sources to the boundary becomes small when compared to the spacing between sources,a situation which is also not advisable in the MFS,leading to an inaccurate result.

For the case of 15 collocation points and for both interface radius,the results seems to stabilize forD/Requal or above 0.40(Fig.4c,d).Similarly,forn=20 the results stabilize forD/Requal or above 0.30(Fig.4e,f),while forn=30 this occurs forD/Requal or above 0.20(Fig.4g,h).In short,the results show that the number of colocation points is closely correlated with theD/Rratio in what concerns the accuracy of the results.However,the interface position seems to have no relevant influence on the accuracy of the results in the presented cases,assuming that it is not too far away from the scatterers(which in some cases could lead to very closely spaced collocation points and virtual sources corresponding to different scatterers),or far too close to the scatterer(which implies stronger spatial variations of the pressure,thus creating problems in correctly simulating this field).

Figure 4:Insertion Loss obtained as a function of the D/R ratio,the number of collocation points(n)and interface radius(R):a)n=10,R=1.1x0.1;b)n=10,R=1.2x0.1;c)n=15,R=1.1x0.1;d)n=15,R=1.2x0.1;e)n=20,R=1.1x0.1;f)n=20,R=1.2x0.1;g)n=30,R=1.1x0.1;h)n=30,R=1.2x0.1.

Observing the results provided by the two models(MFS and FEM),it becomes evident that when an adequate choice of the MFS model parameters is done,an excellent agreement exists.However,very small differences are still visible originated by the fact that the FEM model truncates the discretization at some distance from the scatterers,and then simulates the in finite fluid by means of a PML.This approximation introduces some errors,and can be responsible for those small differences.It should also be noted that the MFS model in this example only requires a total of 80 collocation points(considering 20 colocation points per scatterer withD/R=0.30),and thus a matrix with just 160x160 elements had to be formed.In comparison,the number of nodes in the FEM model was of several thousands,and thus the computational performance is considerably better in the MFS.

4　Numerical simulations

In order to observe the effect of modelling the elastic behaviour of the scatterers on the sound propagation through a sonic crystal using the proposed MFS model,results for a set of different examples were compared with those obtained when considering perfectly rigid scatterers[Martins,Godinho and Picado-Santos(2013)].In the presented examples the two basic layouts consisted of a set of shell cylinders arranged in rectangular or triangular lattice con figurations as depicted in Fig 5.

In both layouts the shell cylinders had their geometric centres equally spaced at 0.4 m with a fixed radius of 0.1 m.The acoustic point source was placed 3 m apart from the scatterers,and varied its vertical position considering two situations:centred with the periodic structure,at y=1.8 m(middle),and close to the top,at y=3.24 m(top).In both cases the response was evaluated at a square grid of 25 points equally spaced of 0.25 m,and positioned 0.5 m away from the array of scatterers.

Figure 5:Tested layout con figurations:a)rectangular lattice;b)triangular lattice.

The scatterers were placed within a fluid medium and filled with the same fluid.A set of variations were established:the type of fluid medium(air and water);the scatterer’s material(PVC,Concrete and Steel);the scatterer’s thickness(5,25 and 95 mm,which correspond to 5%,25%and 95%of the radius),as well as the source position(middle and top).The fluid medium properties(density,ρf,and sound wave velocity,αf)are presented in Tab.1,and the scatterer’s material properties(density,compressional wave velocity,α and shear wave velocity,β)are presented in Tab.2.In the numerical models,the model parameters were set toR=1.1 times the external radius of the cylinders(to de fine the virtual fluid- fluid interface)andD/R=0.40(to de fine the position of the virtual sources),while 30 collocation points were used in the cases involving water and 60 points in the cases in which the fluid is air(indeed,fewer points are required when the fluid is water due to the larger wavelengths involved in the analysis).In the analysis of the Insertion Loss,both damping in the solid and sound absorption in the host fluid were set to zero.

Table 1:Fluid medium properties.

Table 2:Scatterer’s material properties.

In what follows,two sub-sections are presented,the first corresponding to a frequency domain analysis of the behaviour of the sonic crystal structures in terms of their Insertion Loss,and the second to the analysis of time-domain responses,trying to better illustrate the propagation of the sound waves in the presence of those structures.

4.1　Frequency domain Insertion Loss

The results for the acoustic attenuation provided by the different structures are presented in Fig.6 to Figure 8,where ‘Rigid’refers to the results obtained with the model in[Martins,Godinho and Picado-Santos(2013)].

As can be seen,and as a general statement,for a fluid with the properties of air,the Insertion Loss results are almost independent of the scatterers stiffness as the results provided by the two models match almost perfectly.It can be concluded that in this case,as expected,the strong contrast of properties between the fluid and the solid elements leads to a quasi-rigid behaviour of all the tested scatterers,and thus the assumption of rigid boundaries should be a sufficiently good approximation.On the other hand,in the case of a fluid with the properties of water,the results provided by the new numerical model are quite different from those of the rigid model,and the differences become even more pronounced as the cylinders’stiffness is progressively reduced.

Scrutinizing the presented results in more detail,it can be seen that,for rigid scatterers in a rectangular lattice(Fig.6),a band of large attenuation(usually designated as band gap)is clearly identifiable around 400 Hz when the fluid is air,and around 1800 Hz when it is water;this shift for higher frequencies in the case of water is clearly related with the higher sound speed in that fluid,which originates longer wavelengths and translates the band gap phenomenon to those higher frequencies.A similar behaviour is identifiable in Fig.7,for the case of a triangular lattice,although the band gaps occur at slightly higher frequencies.When elastic scatterers are analysed,the same band gaps occur when the fluid medium has the properties of air,and considerable attenuation can be observed;however when the fluid medium is water,and in particular for the lower thickness(5 mm),the provided attenuation is drastically reduced,and the Insertion Loss reaches almost zero;exception is made to specific frequencies that correspond to resonance peaks of the scatterers,where considerable attenuation peaks are visible.

Focusing in the case in which the host fluid is water,it is worth analysing the significant differences that occur between the different elastic materials.In the case of scatterers made of PVC,it can be seen that the Insertion Loss curve differs dramatically from the one obtained with the model in[Martins,Godinho and Picado-Santos(2013)];for that material,the sound reduction effect is only visible when the thicker scatterers are modelled,and even in that case the computed Insertion Loss has a distinctly different character when compared to the rigid structure.For the other materials,concrete and steel,it is observed that as the stiffness increases,the Insertion Loss curve computed by the proposed model tends to approach the one of the rigid model;again,exception is made to some specific frequencies,where relatively high values of attenuation occur,corresponding to resonance peaks of the scatterers.In spite of the marked approximation of the computed curves to the rigid case,it can be noted that,even for steel scatterers,the Insertion Loss never reaches those levels,and thus modelling the real properties of the elastic materials is mandatory in order to provide physically meaningful results when the host fluid is water.

Figure 6:Insertion Loss for the rectangular lattice when the source is at middle and the fluid medium is:air-a1)PVC;a2)Concrete;a3)Steel;water-b1)PVC;b2)Concrete;b3)Steel.

Figure 7:Insertion Loss for the triangular lattice when the source is at middle and the fluid medium is:air-a1)PVC;a2)Concrete;a3)Steel;water-b1)PVC;b2)Concrete;b3)Steel.

Figure 8:Insertion Loss when the source is at top and the fluid medium is air(a)or water(b),and when the elastic material is PVC,for a source positioned close to the top of the structure.In a1)and b1),a rectangular lattice is considered,while in a2)and b2)a triangular lattice is modeled.

Fig.8 illustrates the Insertion Loss results computed when the acoustic source is positioned closer to the top of the sonic crystal structures made on PVC,for both propagation media and for the two lattice con figurations.Observing the behaviour registered when rigid structures are considered,it can be observed that small changes occur in the position of the first band gap and that the corresponding attenuation becomes smaller,for both the rectangular and triangular lattices.A second gap seems to be formed in this case,with very significant attenuations being registered specially for the triangular lattice.It should be noted that,in this case,the dominant incidence of sound waves is oblique to the structure,and thus different interference phenomena are generated.As for the different materials that may compose the structures,the comments given above for the first source position are still valid.

4.2　Time domain analysis

In order to better analyse the obtained results,time domain responses were calculated for the triangular lattice con figuration with 25 mm thickness,considering the host fluid medium with the properties of water and the source centred with the periodic structure,for the three types of material.In the presented cases,the pressure field in the host fluid was evaluated over a regular grid,with a spacing of 0.1 m,spreading over a rectangular area between(x=-3 m,y=-2 m)and(x=3 m,y=5.6 m).

Fig.9 presents a sequence of snapshots taken at different time instants,for the case of steel scatterers,and assuming that the central frequency of the pulse emitted by the source is 2000 Hz.

To allow a better identification of the effect of the structure,the chosen frequency is located within the band gap identified in the previous section.

For an initial time instant of t=2.35 ms,the sound wave has propagated away from the source,and is reaching the scatterers,but still evidencing a perfectly circular wavefront.For t=3.13 ms,the wavefront has already entered the structure,and the scattering effect becomes visible,with a fraction of the incident energy being reflected back to the source side of the structure.When t=3.90 ms,the scattering effect is now very clearly,as well as a reduction of the wave amplitude in the region behind the scatterers.Indeed,this reduction evidences the band gap effect of the sonic crystal,which occurs precisely around the dominant frequency of the Ricker pulse.The presence of multiple scattering from the structures is further revealed for the snapshot taken at t=4.69 ms,with multiple wavefronts radiating from the scatterers.Part of this effect may also be related to resonances occurring in the elastic structures,which were also identified in Fig.7.

To further scrutinize the dynamic behaviour of the structure,the time signals computed at three receivers,all of them at y=2.7 m and with x=-3 m(R1),x=-1 m(R2)and x=2 m(R3),are illustrated in Fig.10;in the same figure,the corresponding frequency spectra are also presented.

Figure 9:Snapshots displaying the pressure wave field calculated using the proposed model when the fluid medium is water and considering steel scatterers,for a central frequency of 2000Hz at:a)t=2.35 ms;b)t=3.13 ms;c)t=3.90 ms;d)t=4.69 ms.

Clearly,at R1 there is very little effect of the periodic structure,and only some minor disturbance is visible starting at t=4 ms as a sequence of oscillations with low amplitude.Similarly,the FFT of this signals exhibits a shape very similar to the Ricker wavelet,with a small disturbance between 1500 Hz and 2500 Hz.In fact,at this receiver the pulse coming from source strongly dominates the response,and is responsible for a very large fraction of the energy reaching the receiver.At R2 and R3,the effect of the structure is much stronger,and affects both time and frequency responses.Indeed,although at R2 the incident pulse is still dominant,the sequence of pulses arriving later has now a comparatively more significant amplitude.In the frequency spectrum,the effect of these oscillations is also very clear,affecting the range around the central frequency of the Ricker pulse.At R3,this effect is now extremely marked,and only a trace of the incident pulse can be identified.The filtering effect of the sonic crystal can be clearly identified also in the frequency spectrum,in which the frequency range of the band gap is greatly attenuated(see also Fig.7 for a clearer identification of the band gap).The final part of the time signal at R3 also reveals further oscillations,which are displayed in more detail in Fig.11.

Figure 10:Responses at receivers R1(a),R2(b)and R3(c).Left column exhibit time signals(a1,b1,c1),and right column the corresponding FFT(a2,b2,c2).

The analysis of this detailed view and of the corresponding FFT response reveal that these oscillations are strongly related with peaks and dips observable also in Fig.7,the first(and dominant)being related to the resonance of the elastic scatterers.Snapshots of the pressure wave field at the time instant t=3.90 ms were also taken for the cases in which the scatterers are composed of the various elastic materials,and are shown in Fig.12,again for a central frequency of 2000 Hz.

Figure 11:Detail of the time signal at the receiver R3 with the resonance effect from the scatterers(a).In(b)the FFT of the signal between 6.5 ms and 10 ms is shown.

These plots help con firming the effect of the different elastic materials in the scattering patterns,revealing that very little energy is reflected back when the PVC scatterers are used,and the scattering becomes more significant as stiffer materials are considered.Indeed,for the case of PVC scatterers,almost no distortion or loss of amplitude is visible in the propagating wavefront,and only traces of scattered energy are observable on the source side of the structure.

5　Conclusions

In this paper,a model based on the Method of Fundamental Solutions that models the sound propagation through a sonic crystal structure is presented.The model makes use of known analytical solutions for the case of a shell solid cylinder embedded in a fluid medium,performing a coupling between these solutions and those for an in finite fluid along a virtual interface located at some distance from the solid,creating a fluid- fluid virtual interface.With this procedure,it is possible to avoid the discretization of the solid material and of the solid- fluid interfaces,reducing computational efforts compared to other approaches such as BEM or FEM,but maintaining a fine level of accuracy.

A set of numerical simulations comparing the proposed model and a previous MFS formulation where the structure behaviour is assumed to be rigid are presented.The simulations included variations on the basic structure geometry as well as the fluid and solid properties.In the simulated cases it was found that when the fluid medium has the properties of air,a considerable band of attenuation could be observed;however when the fluid medium is water,and specially for the lower stiffness scatterers,almost zero level of attenuation is observed,with exception of the frequencies that corresponds to resonance peaks of the scatterers.It was also concluded that in case of high contrast between solid and fluid,the sound propagation through the sonic crystal is not widely affected by the structure stiffness;however in case of low contrast between fluid and solids,strong contrasts between results occur.Indeed,as the stiffness decreases,the more the results tends to diverge,as the Insertion Loss curve of the proposed model starts differing more from the one obtained with the rigid structure model.

In conclusion,the computed results are indicative that the effect of the elastic material of the scatterers cannot be neglected whenever the contrast between the properties of both media is not large.If,however,a very strong contrast occurs(as for example when the fluid is air),the elastic cylinders behave as rigid structures,and almost no differences could be identified.

Acknowledgement:The authors would like to acknowledge the financial support by FCT(Fundação para a Ciência e a Tecnologia,Portugal)and COMPETE,through research project PTDC/ECM-COM/1438/2012.

Atluri,S.N.;Shen,S.(2002):The meshless local Petrov-Galerkin(MLPG)method:a simple&less-costly alternative to the finite element and boundary element methods.CMES:Computer Modeling in Engineering&Sciences,vol.3,no.1,pp.11-51.

Cao,Y.J.;Hou,Z.L.;Liu,Y.Y.(2004):Convergence problem of plane-wave expansion method for phononic crystals.Phys.Lett.A,vol.327,pp.247–253.

Cao,Y.J.;Hou,Z.L.;Liu,Y.Y.(2004):Finite difference time domain method for band-structure calculations of two-dimensional phononic crystals.Solid State Commun,vol.132,pp.539-543.

Castiñeira-Ibáñez,S.;Rubio,C.;Romero-García,V.;Sánchez-Pérez,J.V.;

García-Raffi,L.M.(2012):Design,Manufacture and Characterization of an Acoustic Barrier Made of Multi-Phenomena Cylindrical Scatterers Arranged in a Fractal-Based Geometry.Archives of Acoustics,vol.37,no.4,pp.455–462.

Fairweather,G.;Karageorghis,A.(1998):The method of fundamental solutions for elliptic boundary value problems.Advances in Computational Mathematics,vol.9,pp.69–95.

Fairweather,G.;Karageorghis,A.;Martin,P.A.(2003):The method of fundamental solutions for scattering and radiation problems.Engineering Analysis with Boundary Elements,vol.27 pp.759–69.

Godinho,L.;Amado-Mendes,P.;Pereira,A.(2011):A Hybrid Analytical-Numerical Model Based on the Method of Fundamental Solutions for the Analysis of Sound Scattering by Buried Shell Structures.Mathematical Problems in Engineering,vol.2011,pp.710623

Godinho,L.;Tadeu,A.;Branco,F.(2004):Dynamic analysis of submerged fluidfilled pipelines subjected to a point pressure load.Journal of Sound and Vibration,vol.271,no.1-2,pp.257–277.

Godinho,L.;Tadeu,A.;Mendes,P.(2007):Wave propagation around thin structures using the MFS.Computers,Materials and Continua(CMC),vol.5,no.2,pp.117–127.

Golberg,M.;Chen,C.S.(1998):The method of fundamental solutions for potential,Helmholtz and diffusion problems.Boundary Integral Methods:Numerical and Mathematical Aspects.Comput.Mech.,vol.1,pp.103–176.

Gu,M.H.;Young,D.L.;Fan,C.M.(2009):The method of fundamental solutions for one-dimensional wave equations.Computers,Materials&Continua(CMC),vol.11,no.3,pp.185-208.

Kansa,E.(1990a):Multiquadrics-A scattered data approximation scheme with applications to computational fluid-dynamics–I:surface approximations and partial derivative estimates.Computers&Mathematics with Applications,vol.19,pp.127-145.

Kansa,E.(1990b):Multiquadrics-A scattered data approximation scheme with applications to computational fluid-dynamics—II solutions to parabolic,hyperbolic and elliptic partial differential equations.Computers&Mathematics with Applications,vol.19,pp.147-161.

Kausel,E.;Roësset,J.M.(1992):Frequency domain analysis of undamped systems.Journal of Engineering Mechanics,vol.118,no.4,pp.724–734.

Kushwaha,M.S.;Halevi,P.(1996):Giant acoustic stop bands in two-dimensional periodic arrays of liquid cylinders.Appl.Phys.Lett.,vol.69,pp.31-33.

Li,F.L.;Wang,Y.S.;Zhang,C.(2011):Bandgap calculation of two-dimensional mixed solid– fluid phononic crystals by Dirichlet-to-Neumann maps.Phys.Scr.,vol.84,pp.055402.

Li,F.L.;Wang,Y.S.;Zhang,C.;Yu,G-L.(2013):Bandgap calculations of two-dimensional solid– fluid phononic crystals with the boundary element method.Wave Motion,vol.50,pp.525–541.

Li,F-L.;Yue-Sheng,W.;Chuanzeng,Z.;Gui-Lan,Y.(2013):Boundary element method for band gap calculations of two-dimensional solid phononic crystals.Engineering Analysis with Boundary Elements,vol.37,pp.225–235.

Marin,L.;Karageorghis,A.(2009):Regularized MFS-based boundary identification in two-dimensional Helmholtz-type equations.Computers,Materials&Continua(CMC),10:259–293

Martínez-Sala,R.;Sancho,J.;Sánchez,J.V.;Gómez,V.;Llinares,J.;Meseguer,F.(1995):Sound attenuation by sculpture.Nature,vol.378,pp.241.

Martins,M.;Carbajo,J.;Godinho,L.;Mendes,P.;Ramis,J.(2013):Insertion Loss Provided by a Periodic Structure–Numerical and Experimental Evaluation.Proceedings of TecniAcustica Valladolid(2013).

Martins,M.;Godinho,L.;Picado-Santos,L.(2013):Numerical Evaluation of Sound Attenuation Provided by Periodic Structures.Archives of Acoustics,vol.38,no.4,pp.503–516.

Mei,J.;Liu,Z.Y.;Qiu,C.Y.(2005):Multiple-scattering theory for out-ofplane propagation of elastic waves in two-dimensional phononic crystals.J.Phys.:Condens.Matter,vol.17,pp.3735-3757.

Sánchez-Pérez,J.V.;Rubio,C.;Martínez-Sala,R.;Sánchez-Grandia,R.;

Gómez,V.(2002):Acoustic barriers based on periodic arrays of scatterers.Appl.Phys.Lett.,vol.81,pp.5240.

Sigalas,M.M.;García,N.(2000):Importance of coupling between longitudinal and transverse components for the creation of acoustic band gaps:the aluminum in mercury case,Appl.Phys.Lett.,vol.76,pp.2307-2309.

Smyrlis,Y.S.;Karageorghis,A.(2003):Some aspects of the method of fundamental solutions for certain biharmonic problems.CMES:Computer Modeling in Engineering&Sciences,vol.4,no.5,pp.535-550.

Tsai,C.C.;Lin,Y.C.;Young,D.L.;Aturi,S.N.(2006):Investigations on the accuracy and condition number for the method of fundamental solutions.CMES:Computer Modeling in Engineering&Sciences,vol.16,no.2,pp.103-114.

Vasseur,J.O.;Deymier,P.A.;Djafari-Rouhani,B.;Pennec,Y.;Hladky-

Hennion,A.C.(2008):Absolute forbidden bands and waveguiding in two-dimensional phononic crystal plates.Phys.Rev.B,vol.77,pp.085415.

Wang,Y.F.;Wang,Y.S.;Su,X.X.(2011):Large band gaps of two-dimensional phononic crystals withcross-likeholes.J.Appl.Phys.,vol.110,pp.113-520.

Wu,L.Y.;Chen,L.W.;Liu,C.M.(2009):Acoustic pressure in cavity of variously sized two-dimensional sonic crystal with various filling fraction.Phys.Lett.A,vol.373,pp.1189–1195.

Wu,Y.M.;Lu,Y.Y.(2008):Dirichlet-to-Neumann map method for analyzing interpenetrating cylinder arrays in a triangular lattice.J.Opt.Soc.Am.B,vol.25,pp.1466–1473.

Yan,Z.Z.;Wang,Y.S.(2006):Wavelet-based method for calculating elastic band gaps of two-dimensional phononic crystals,Phys.Rev.B,vol.74,pp.224-303.

Young,D.L.;Ruan,J.W.(2005):Methods of fundamental solutions for scattering problems of electromagnetic waves.CMES:Computer Modeling in Engineering&Sciences,vol.7,no.1,pp.223-232.

Yuan,J.H.;Lu,Y.Y.;Antoine,X.(2008):Modeling photonic crystals by boundary integral equations and Dirichlet-to-Neumann maps.J.Comput.Phys.,vol.227,pp.4617-4629.