Harnessing inclusions to tune post-buckling deformation and bandgaps of soft porous periodic structures

Journal of Sound and Vibration 459 (2019) 114848

Contents lists available at ScienceDirect

Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi

Harnessing inclusions to tune post-buckling deformation and bandgaps of soft porous periodic structures Jian Li a, b, Yueting Wang c, Weiqiu Chen a, b, d, Yue-Sheng Wang c, e, Ronghao Bao a, * a Key Laboratory of Soft Machines and Smart Devices of Zhejiang Province & Department of Engineering Mechanics, Zhejiang University, Hangzhou, 310027, China b State Key Lab of CAD & CG, Zhejiang University, Hangzhou, 310058, China c Institute of Engineering Mechanics, Beijing Jiaotong University, Beijing, 100044, China d Soft Matter Research Center, Zhejiang University, Hangzhou, 310027, China e School of Mechanical Engineering, Tianjin University, Tianjin, 300350, China

a r t i c l e i n f o

a b s t r a c t

Article history: Received 13 October 2018 Received in revised form 15 June 2019 Accepted 9 July 2019 Available online 10 July 2019 Handling Editor: G. Degrande

The emergence of phononic crystals paves a new way for manipulating elastic waves in structures for their particular bandgap properties. In this paper, the two-dimensional soft porous periodic structures that can be filled with hard inclusions are considered. Both numerical simulations and experiments are conducted to study the effects of inclusions on the buckling modes, post-buckling deformations, and band structures in soft porous periodic structures. It is found that either the number or the arrangement (i.e. filling pattern) of the inclusions has a great influence on the bandgap characteristics. Meanwhile, the material damping affects the wave propagation in soft phononic crystals significantly in the high frequency range. Compared with the unfilled soft porous structure, the sensitivity of the post-buckling deformation to the initial geometrical imperfections can be significantly reduced for the structure filled with inclusions. This means the post-buckling deformation could develop robustly. Further numerical study indicates that the bandgaps can be tuned in a versatile and reversible way when the structure undergoes a large deformation. A more fruitful manner to tune the bandgaps therefore can be achieved by changing the filling pattern of inclusions along with dramatically deforming the structure. The work provides a useful reference for the design of tunable phononic switches and acoustic filters. © 2019 Elsevier Ltd. All rights reserved.

Keywords: Soft porous material Damping Robust post-buckling deformation Inclusion Tunable bandgaps

1. Introduction Recently, the propagation of elastic waves in phononic crystals [1] (periodic structures with bandgap properties) has attracted a great deal of research interest. Phononic crystals exhibit many interesting properties, such as defect states [2], negative refraction [3], and acoustic focusing [4], which can be utilized to design specific acoustic devices, including acoustic filters [5], acoustic waveguides [6], acoustic lenses [7], etc. In addition, the unique bandgap characteristics of phononic

* Corresponding author. E-mail address: [email protected] (R. Bao). https://doi.org/10.1016/j.jsv.2019.114848 0022-460X/© 2019 Elsevier Ltd. All rights reserved.

2

J. Li et al. / Journal of Sound and Vibration 459 (2019) 114848

crystals can be effectively employed to control the propagation of elastic waves in some frequency ranges. Hence, phononic crystals have a wonderful engineering application prospect in vibration, noise reduction, sound insulation [8e10], and so on. In many cases, due to the varying working environments, we hope that the bandgap properties of phononic crystals could be controlled in real time. It is therefore of practical importance to study the tunability of phononic crystals so as to manipulate in real time the position and width of a bandgap for instance. A feasible way to change the bandgaps is to alter the topology and the periodic arrangements (e.g. in square, rectangular, and triangular lattices) of the periodic cells as well as the contrast of material properties between the material phases (e.g. matrix and inclusion). As a special design of phononic crystals, porous periodic structures are easy to be manufactured, and at the same time, the holes can be readily filled with other materials. Therefore, they have been widely considered to clarify the influences of the above-mentioned parameters on the bandgaps. For example, Refs. [11e14] discussed the effects of geometrical configurations on the bandgaps in porous periodic structures with different pore shapes. It is noticed that once a phononic crystal is made, the bandgaps can be tuned, in most cases, only by applying a biasing field so as to locally adjust the geometrical configuration of the structure [15]. In recent years, Bertoldi et al. [16e20] have initiated the exploration of soft porous periodic structures which can undergo reversible elastic deformations under the applied mechanical loads. They have found that the geometrical configuration of a soft porous periodic structure can be dramatically changed by exerting an external load, and consequently, its band structure could evolve in a remarkable way. In particular, once the applied load exceeds a specific critical value, the porous phononic crystal made of soft materials (e.g. rubber or silicone rubber) usually presents a large post-buckling deformation, which in general can effectively change the geometrical configuration as well as the effective stiffness of the structure. Therefore, the bandgaps of soft phononic crystals can be dramatically controlled in a reversible way due to the superelasticity of soft materials [21e23]. Recently, Li et al. [24] proposed a design of soft composite materials consisting of stiff inclusions and voids periodically distributed in a soft matrix. In addition to the appearance of negative Poisson's ratio, they numerically showed that the bandgaps of such soft materials could also be dramatically tuned through deformation. To better assist the design of acoustic devices with tunable dynamic performance, it is of particular significance to develop a deep understanding of the variation with external loads of the bandgaps in a soft phononic crystal. To induce the desirable post-buckling configuration, it is generally necessary to introduce a suitable geometrical imperfection in numerical simulations, so that the results can be compared with the experimental observations. However, if the critical loads of different buckling modes are closely spaced, the induced post-bucking deformation would be very sensitive to the geometrical imperfection [25]. In this case, it will become rather difficult to control the development of post-buckling deformation as one desires. For soft porous periodic structures, filling hard inclusions in the pores periodically could be a feasible way to change the characteristics of the buckling behavior as well as to tune the bandgaps either before or after deformation. In this paper, we will systematically explore the effects of hard inclusions on the post-buckling deformation and the bandgap characteristics of soft porous periodic structures. The aim of the study is actually two-fold: On the one hand, we hope that different arrangements (filling patterns) of hard inclusions could be used to effectively change the bandgap characteristics since this approach seems to be very convenient and also economic; on the other hand, we are interested in attaining robust post-buckling evolution so as to reliably tune the bandgaps via large deformation for soft porous periodic structures with a specially selected filling pattern of hard inclusions. To this end, two-dimensional (2D) soft periodic structures with circular holes arranged in an equilateral triangular pattern are considered here. We will show numerically and experimentally that placing hard inclusions in pores properly can make the post-buckling deformation insensitive to the initial geometrical imperfection. In particular, the effect of damping intrinsic to soft materials on wave propagation will be elucidated. The commercial finite element (FE) software Abaqus is utilized throughout the study to numerically investigate the influences of hard inclusions on the buckling mode, the post-buckling deformation, and the band structure when the soft phononic crystal undergoes a large deformation. Python scripting language is adopted to accomplish the relevant automatic meshing and all the calculation codes in Abaqus [26]. In experiments, frequency response measurements are conducted to compare with the numerical results. The results and conclusions should be very useful, and the present study actually suggests an alternative and efficient way to tune elastic waves in soft porous periodic structures.

2. Computational model 2.1. Matrix and inclusions A 2D soft periodic structure with circular holes arranged in an equilateral triangular pattern indicated in Fig. 1(a) is considered in this study. The structure is assumed to be long enough in the thickness direction, and an analysis based on the plane-strain model shall be sufficient [16]. In consideration of the computational cost and the deformation characteristics of the structure under an external load, we take in Fig. 1(a) the part enclosed by dashed lines as a representative volume element (RVE), which consists of 2
2 unit cells (one unit cell is indicated by solid lines). It has been verified to be the minimum periodic cell for post-buckling deformation analysis [25]. The RVE has a total of eight circular holes (including those on the boundary in view of the periodic condition), endowed with a sufficient number of arrangements of inclusions to be discussed. Therefore, we merely focus on the arrangements of inclusions with respect to this RVE in the subsequent analysis. Hard cylindrical inclusions or columns made of steel (Young's modulus Es ¼ 194:02 GPa, Poisson's ratio ys ¼ 0:3, and density rs ¼ 7930 kg=m3 ) can be filled in these circular holes. For numerical simulations, a 2D mesh is constructed for inclusions using

J. Li et al. / Journal of Sound and Vibration 459 (2019) 114848

3

Fig. 1. Porous pffiffiffiperiodic structure and filling patterns: (a) The RVE (dashed line) and unit cell (solid line) (the width of the RVE is w ¼ 18.212 mm, the height of the RVE is h ¼ 3w, and the diameter of the circular hole is d ¼ 8 mm); (b) RVEs for eight typical filling patterns; (c) The first Brillouin zone and the irreducible Brillouin zone (gray area).

triangular quadratic elements (element type CPE6 in Abaqus). In accordance with inclusions, the soft host structure will sometimes be referred to as the matrix. To illustrate the arrangements of inclusions, eight typical filling patterns are shown in Fig. 1(b), with the number of inclusions varying from one to eight. More filling patterns can be derived from these typical patterns by adjusting the arrangement of inclusions in each pattern. In the original structure without inclusions (i.e., the matrix only), its geometrical configuration can be changed by adjusting the size (diameter) of the circular hole. Here, the porosity j is used to represent the main geometric feature of the structure, which is defined as

j¼

Acircle ; ARVE

(1)

where Acircle is the area of all the pores in the RVE and ARVE is the total area of the RVE. In this paper, the diameter of the circular hole is set to be d ¼ 8 mm and the porosity is set to be 70%. The soft matrix (the initial density is rm0 ¼ 1300 kg=m3 ) is made of hyperelastic material. Here, the neo-Hookean model, one special kind of the polynomial model (and of the Ogden model as well), is used to describe the constitutive behavior of the soft matrix. The corresponding strain energy density function W is

 K
W ¼ C10 I 1  3 þ m0 ðJ  1Þ2 ; 2

(2)

where C10 ¼ mm0 =2 with mm0 being the initial shear modulus, Km0 is the initial bulk modulus, I 1 is the first deviatoric strain invariant, and J is the local volume change [26]. The initial shear modulus is taken to be mm0 ¼ 0:5 MPa, and hence the material parameter is C10 ¼ 0:25 MPa. For the almost incompressible hyperelastic rubber we used, the ratio of the initial bulk modulus to the initial shear modulus is chosen to be Km0/mm0 ¼ 100, which indicates the initial Poisson's ratio is nm0 ¼ 0:495. The above material parameters are employed in the subsequent numerical simulations so as to compare with the experiments. Furthermore, in the deformation and wave analyses, we assume that the steel columns and the matrix are perfectly connected or bonded, i.e. the interfaces between the matrix and the inclusions are perfect. In addition, for the adopted planestrain model, a fine mesh of 6-node hybrid elements (element type CPE6H in Abaqus) is adopted to ensure the numerical convergence. 2.2. Loading conditions Loading condition is a very important factor that affects the buckling mode and the post-buckling deformation. Here for illustration, both uniaxial and biaxial loadings are considered, which can be uniformly characterized by the deformation gradient F ¼ vx=vX, where X and x are the position vectors of the same material point before and after deformation. For the biaxial compression, we have


 F ¼ ð1 þ εxx Þex 5ex þ 1 þ εyy ey 5ey ;

(3)

where ex and ey are the base vectors in the Cartesian coordinates, and εxx and εyy denote the macroscopic nominal strains in the x and y directions, respectively. Then εxx and εyy can be described as [25]

εxx ¼ l cos q; εyy ¼ l sin q;

(4)

4

J. Li et al. / Journal of Sound and Vibration 459 (2019) 114848

where l (l  0) denotes the load parameter, and the loading path angle q (p  q  3p=2) is related to the ratio between εxx and εyy . For instance, q ¼ 5p=4 corresponds to the equi-biaxial compression, i.e. εxx ¼ εyy . In the case of uniaxial compression, the deformation gradient F may be given by

F ¼ ð1 þ εxx Þex 5ex þ ley 5ey

(5)


 F ¼ lex 5ex þ 1 þ εyy ey 5ey :

(6)

Here, Eq. (5) corresponds to the case when the load is applied in the x direction, and l is to be determined from the free boundary condition in the y direction (i.e., the nominal stress in the y direction is zero), while Eq. (6) is for the load in the y direction, and l is to be determined from the free boundary condition in the x direction (i.e., the nominal stress in the x direction is zero). Since the loading conditions in this paper are all compressive in nature, we may write the nominal compression l as l ¼  ε, where ε denotes the nominal strain εxx or εyy in the x or y direction, to facilitate the description and the discussion in the subsequent analysis. 2.3. Deformation and wave analyses In the FE simulations, we first pay attention to the RVE taken from the infinite periodic structure, with periodic boundary conditions and Bloch boundary conditions applied to the boundary of the RVE for deformation analysis and wave analysis, respectively [27]. Based on the general theory of nonlinear elastic deformation [28,29] and the small-on-large theory [30], the problems of linear wave propagation in a deformed body can be solved. In the deformation analysis, two virtual nodes are introduced to model the deformation of the RVE along two periodic directions, while the coupling and constraint equations corresponding to the periodic boundary conditions are established via the motion of these two virtual nodes. The eigenvalue buckling analysis in Abaqus is invoked to obtain the critical loads and the corresponding buckling modes. Then the postbuckling deformation of the structure with the reasonable geometrical imperfections corresponding to the buckling mode can be conducted by the modified pseudo arc-length method (i.e., Riks algorithm) in Abaqus [26]. For Bloch wave analysis, two instances with identical mesh and material properties, representing the real and imaginary parts of the corresponding fields respectively, are taken to deal with the complex-valued problem caused by the Bloch boundary conditions [31]. To acquire the band diagrams, wave vectors are scanned on the boundary of the irreducible Brillouin zone (see Fig. 1(c)). In order to simulate the bandgap behavior of a finite structure, steady-state dynamics analysis in Abaqus may be adopted to obtain the steady-state response of the structure over a frequency range of interest. With the acceleration amplitudes at the output and input sides of the model (Aoutput and Ainput respectively, see Appendix A), we can calculate the transmittance T according to the following formula


  T ¼ 20 log Aoutput Ainput :

(7)

3. Tuning bandgaps by filling inclusions In this section, we are only interested in the undeformed structures and will investigate the effects of the number and arrangement of inclusions on the bandgap characteristics of 2D soft porous structures. 3.1. Numerical analysis As mentioned in Section 2.1, there are many possible arrangements of inclusions in the RVE depicted in Fig. 1(a). Here, for the sake of discussion, we choose four typical filling patterns for which the band diagrams and transmittance spectra are simultaneously given in Fig. 2. Only the results for frequencies lower than 1200 Hz are taken into account. Besides, the 10
1 (or 1
10) enlarged RVEs with periodic conditions applied in a certain direction are considered for the steady-state response analysis (see Appendix A). To make it clear, we use gray and orange bars in Fig. 2 to denote the complete bandgaps and the directional bandgaps, respectively. Firstly, the bandgap characteristics for two filling patterns with two inclusions in a RVE are compared in Fig. 2(a) and (b). Fig. 2(a) is for the case when two hard inclusions are inserted into the holes in a column. It is found from the band diagram that there are abundant complete bandgaps, and the passbands between two adjacent bandgaps (e.g., the first and the second complete bandgaps) are very narrow, which can be further confirmed by the transmittance spectra. For comparison, Fig. 2(b) gives the results for two hard inclusions in a row. It is intriguing that no complete bandgap could be found in Fig. 2(b) in the frequency range lower than 1200 Hz. However, there are plentiful directional bandgaps in the direction of G-Y. It indicates that this structure can be utilized to impede the propagation of elastic waves in one direction while keeping the waves propagable in other directions. The same phenomenon can also be found for the cases when there are three and four inclusions in a RVE as shown in Fig. 2(c) and (d), respectively. The above results indicate that unidirectional suppression of wave propagation can be realized by employing a proper arrangement as well as a certain amount of inclusions.

J. Li et al. / Journal of Sound and Vibration 459 (2019) 114848

5

Fig. 2. The effects of arrangement and number of inclusions on the bandgaps.

It is noteworthy that there are several flat bands in the G-Y direction as shown in Fig. 2(b) and (c), indicating the locally resonant modes. Here, we take the low-order modes for illustration. Fig. 3(a) and (b) show the 4th and 5th modes of the filling pattern with two inclusions, respectively. It is found that the vibrations are localized in the soft matrix. However, when one more inclusion is inserted into the hole, the vibrations become localized around or in the hard inclusions as shown in Fig. 3(c) and (d) for the 4th and 6th modes respectively. In addition, there is no flat band in Fig. 2(d), which implies that the locally resonant mode disappears in the frequency range we consider. Thus, we could use inclusions to tune the resonant mode and change the acoustic properties of the porous structure consequently. 3.2. Experiments To verify the previous numerical results, we fabricated a rubber sample comprising 5
3 RVEs with the same size as the structure shown in Fig. 1(a). The thickness is about 50 mm, which could ensure that the sample deforms approximately in a plane-strain fashion provided the load is appropriately applied. The experimental details are described in Appendix B.

Fig. 3. Contour plots of normalized wave amplitude of different modes at Y point: (a) The 4th mode (at 244.7 Hz) of the RVE with two inclusions; (b) The 5th mode (at 432.9 Hz) of the RVE with two inclusions; (c) The 4th mode (at 254.3 Hz) of the RVE with three inclusions; (d) The 6th mode (at 703.9 Hz) of the RVE with three inclusions.

6

J. Li et al. / Journal of Sound and Vibration 459 (2019) 114848

We first measured the frequency responses of the sample with two rigid inclusions in a column in the RVE, for which the band diagram was already given in Fig. 2(a). The theoretical predictions (based on the RVE) and the experimental observations (based on the finite sample) are compared in Fig. 4 for waves propagating both in the G-X and G-Y directions. It is seen that the two results match quite well in terms of the bandgap characteristics. However, the experimental curves exhibit obvious attenuation at the high frequencies outside the bandgap. It is worth mentioning that most soft materials like rubber are rate dependent and whose damping in general cannot be neglected [32]. Therefore, it seems necessary to include the damping effect in the frequency domain analysis. To this end, we took the 5
3 enlarged RVEs into account (see Appendix A). Through dynamic mechanical analysis (DMA) test on the same specimen with the fabricated sample, we obtained the loss factor as h ¼ 0:061 at room temperature (25 + C) (see Appendix B). In the simulations incorporating damping effect, the same phenomenon was observed numerically, that is, there is significant attenuation within the passbands at high frequencies. By comparing the transmittance spectra for h ¼ 0 (no damping) and h ¼ 0:061, we also find that the resonant peaks in bandgaps, which correspond to natural frequencies of the finite structure [33], are weakened or even vanish due to the damping effect. To further illustrate the unidirectional suppression of elastic waves in the structure shown in Fig. 2(b), the frequency response measurements were then conducted for waves propagating both in the G-X and G-Y directions. Here, only the results for frequencies below 800 Hz are considered and plotted in Fig. 5. From the experimental curves and numerical transmittance spectra with damping effect, we find that the two results are in good agreement and there is significant attenuation for wave propagation in the G-Y direction at high frequencies. Besides, the transmittance in the G-X direction is obviously larger than that in the G-Y direction in most frequency ranges, which might be feasible to design a switch for elastic waves by rotating the structure by 90+. The results without damping effect are displayed in Fig. 5(d), which may be compared with Fig. 5(c) to clarify the effect of damping. It is noteworthy that it is very difficult to identify the bandgaps from the transmittance spectrum in the G-Y direction as indicated in Fig. 5(d). That is because, on the one hand, the resonant peaks in

Fig. 4. Comparison between the calculated band diagrams and the transmittance spectra in the directions of G-X (a) and G-Y (b).

Fig. 5. Experimental transmission coefficients in the G-X and G-Y directions: (a) Band diagram of the structure shown in Fig. 2(b); (b) Experimental transmittance curves; (c) Numerical transmittance curves with damping effect (the loss factor h ¼ 0:061); (d) Numerical transmittance curves without damping effect (h ¼ 0).

J. Li et al. / Journal of Sound and Vibration 459 (2019) 114848

7

the transmittance spectrum, which correspond to the natural frequencies of the finite structure, will blur the boundaries of the bandgaps. On the other hand, the excitation source adopted in the computational and experimental models (see the similar cases in Figs. A.1(b) and A.1(c) and the experimental set-up in Fig. B.3) cannot guarantee the generation of pure plane waves, which will also reduce the possibility of detecting the directional bandgaps accurately. Nevertheless, it is still a common practice to experimentally characterize the wave propagation behavior in phononic crystals [25,34]. Further comparing the numerical results (transmittance spectra in the G-Y direction) in Fig. 5(c) and (d), we can find that the peaks in the bandgaps will be weakened sharply when the damping effect is involved. It is therefore intriguing to design practical acoustic devices when material damping presents. 4. Tuning bandgaps by large deformation 4.1. Effects of inclusions on buckling modes In this section, we are going to show how the distribution of inclusions affects the buckling modes, and to demonstrate the effect of large deformation caused by external loads on the bandgaps as well. Actually, Shan et al. have investigated the buckling modes of the structure in Fig. 1(a) under biaxial compression and found that the change of the loading path angle q could induce the transition of buckling modes due to the closeness of the first two critical buckling loads. They have also attained a chiral mode by a linear combination of the first two buckling modes [25]. Inserting some hard inclusions into the pores of the soft porous structure may restrain the transition of buckling modes, making the post-buckling deformation insensitive to the initial geometrical imperfections and enabling the chiral mode easy to be derived. Since inserting more steel inclusions into the holes may not be conducive to the deformation of the structure, we here take one specially selected filling pattern shown in Fig. 6(a) as an example to perform the buckling analysis. Firstly, the case of uniaxial compression in the x or y direction is considered, for which the first two buckling modes and the corresponding critical loads are given in Fig. 6. It is seen that the first-order critical loads in both cases are about 0.056, while the difference between the first two critical loads in each case is quite obvious. Hence, in the subsequent post-buckling analysis, initial geometrical imperfection corresponding to the first-order buckling mode could be introduced to induce the post-buckling deformation in Abaqus. Next, we consider the first two buckling modes of the structure under a biaxial compression. The corresponding critical loads are calculated when the loading path angle q changes. Fig. 7 depicts the variations of the critical load lcr with the loading path angle q for the RVE shown in Fig. 6(a). For the porous structure without inclusions, Shan et al. [25] have concluded that the first two critical loads are very close to each other and the transition of buckling modes appears when q ¼ 5p=4 where the two modes by chance have the same critical load. Therefore, the buckling mode is quite sensitive to the initial geometrical imperfection, and hence it is difficult to induce the desired post-buckling deformation either numerically or experimentally. However, as shown above, placing the hard inclusions in the pores will make the critical loads of the first two buckling modes differ from each other more obviously when the structure is subjected to a biaxial compression. Actually, there is no transition of buckling modes that will occur when the loading path angle q varies. Interestingly, the buckling modes shown in Fig. 6(b) and (d) and the A-mode in Fig. 7 are almost of the same chiral configuration. These buckling modes can all be denoted as A-mode to facilitate the discussion. We now qualitatively analyze the effect of inclusions on the formation of A-mode. In general, the formation of buckling modes of a porous periodic structure may be understood by the folding mechanism developed in Ref. [25]. Fig. 8(a) plots the expanded structure corresponding to the RVE in Fig. 6(a). For the low-order buckling modes, the part enclosed by solid lines in Fig. 8(a) can be simplified as a rigid triangle, while the inclusion and the part surrounding it (dashed line) can be simplified as a rigid hexagram. Hence, the original periodic structure can be simplified to the network as shown in Fig. 8(b), which is built with rigid triangles and hexagrams. The reconfiguration of this network can only be achieved by rotating the triangles and/or the hexagrams. In other

Fig. 6. The first two buckling modes and the corresponding critical loads in the case of uniaxial compression: (a) RVE with two inclusions (b) lxcr;1 ¼ 0:056; (c) lxcr;2 ¼ 0:165; (d) lycr;1 ¼ 0:056; (e) lycr;2 ¼ 0:093. lxcr;i denotes the critical load, where x ¼ x; y means the loading direction and i ¼ 1; 2 signifies the order of the buckling mode.

8

J. Li et al. / Journal of Sound and Vibration 459 (2019) 114848

Fig. 7. The relation between the buckling critical load lcr and the loading path angle q.

Fig. 8. Formation of the chiral buckling mode (A-mode): (a) Expanded structure; (b) Equivalent structure; (c) The evolution of the circular hole; (d) Deformed structure and the A-mode.

words, the deformation of the structure is actually converted into the relative rotations between the triangles and hexagrams. On the other hand, the circular hole could be equivalent to a hexagonal hole with six hinged joints, but due to the rigid inclusions, two of them should be changed to the rigid ones. Thus, the degrees of freedom of the motion of the hexagonal hole are reduced by two, becoming the same as a parallelogram with four hinged joints, as shown in Fig. 8(c). Therefore, the structure has only one degree of freedom in view of its deformation, and only the chiral pattern shown in Fig. 8(d), which corresponds to the A-mode, can be attained by rotating the triangles. In this sense, the chiral modes are prone to occur among the low-order buckling modes of the porous structure filled with hard inclusions.

4.2. Experiments To illustrate that more stable post-buckling deformation will occur when properly inserting hard inclusions into the rubber sample, both uniaxial and equi-biaxial compression tests were conducted. Firstly, the porous structure without hard inclusions was considered. We took two samples from the same batch to do the uniaxial compression experiments. Although these two samples are almost the same in terms of material properties and configuration, more or less different geometrical imperfections exist in them due to the inevitable randomness. Fig. 9 shows the deformations of these two rubber samples under uniaxial compression when the load parameter is l ¼ 0:20. It is found that the post-buckling deformation of the sample in Fig. 9(a) perfectly follows the first buckling mode (see mode I in Fig. 9) [25]. However, some parts of the deformation of the second sample in Fig. 9(b) are obviously triggered by another buckling mode (see mode II in Fig. 9), which leads to the nonuniform deformation of the whole structure. Next, we inserted the hard inclusions into the holes of the samples and considered the finite structures comprising the RVEs indicated in Fig. 6(a). Both numerical and experimental results are shown and compared in Fig. 10, which indicates an excellent agreement. Notably, the first buckling modes under the uniaxial and biaxial compressive loadings are almost the

J. Li et al. / Journal of Sound and Vibration 459 (2019) 114848

9

Fig. 9. Comparison between two samples of the porous structure with different geometrical imperfections under the same loading.

Fig. 10. Numerical and experimental results for the filling pattern under uniaxial and equi-biaxial compressions: (a) Uniaxial compression in the x direction with l ¼ 0:15; (b) Uniaxial compression in the y direction with l ¼ 0:15; (c) Equi-biaxial compression with l ¼ 0:15.

same. And indeed, the post-buckling deformations shown in Fig. 10 are all induced by the chiral mode or A-mode, which indicates the deformation of the porous structure filled with hard inclusions is also insensitive to the loading. Hence, the experimental observation here confirms the analysis based on Fig. 8.

4.3. Tunability of bandgaps through post-buckling deformation In the previous discussion, we have attained the robust post-buckling deformation by selecting the filling pattern shown in Fig. 6(a). When a soft periodic structure undergoes post-buckling deformation, its geometrical configuration as well as the effective material properties will be altered, which in turn may change its bandgap characteristics. Here we pay more attention to the tunability of bandgaps through this robust post-buckling deformation. Since there is no shear deformation with respect to the RVE during the compression, the lattice remains as a rectangular lattice after deformation. Thus, for the calculation of the band structure, the wave vectors are still selected to be on the boundary G-X-M-Y-G of the irreducible Brillouin zone, but the size of lattice changes. Fig. 11 shows the band diagrams of the structure which undergoes the large deformations depicted in Fig. 10. When the structure is load-free (i.e., undeformed), it possesses an ultra-broad complete bandgap in the frequency range of 314e1143 Hz. In each deformed configuration, the original complete bandgap is narrowed to the range of about 450e965 Hz and new bandgaps occur. Notably, the band diagrams in Fig. 11(b) and (c) have the similar bandgap characteristics due to the same load magnitude and the robustness of deformation. As mentioned before, the damping in soft matrix has great effect on wave propagation in the high frequency region. Here, we introduce the loss factor obtained from DMA test into our computational model and calculate the transmittance spectra in

10

J. Li et al. / Journal of Sound and Vibration 459 (2019) 114848

Fig. 11. Transmittance spectrum, band diagram, and the corresponding deformation of the structure under applied mechanical load: (a) Undeformed structure; (b) Uniaxial compression in the x direction with l ¼ 0:15; (c) Uniaxial compression in the y direction with l ¼ 0:15; (d) Equi-biaxial compression with l ¼ 0:15.

Fig. 12. Variation of the bandgap characteristics with the applied load in the y direction: (a) The relation between the distribution of bandgaps and the applied load; (b) The relation between the relative bandgap width fr and the applied load l.

the G-X direction with the 10
1 enlarged RVEs (see Appendix A). As shown in Fig. 11, the damping has negligible effect on the transmission coefficients in the bandgap as well as in the low frequency region. It is noteworthy that in Fig. 11(b) and (c), several new bandgaps are generated, which contributes to the existence of narrow passbands in the high frequency range, where obvious attenuation is observed from the transmission coefficients. This means it is hard to identify the narrow

J. Li et al. / Journal of Sound and Vibration 459 (2019) 114848

11

passbands in the high frequency region from the transmittance curves. As a consequence, the tunability of bandgaps at high frequencies will be limited in consideration of the damping in actual soft materials (e.g. rubber). We further consider the relation between the distribution of bandgaps and the load magnitude. Here we take the structure exhibited in Fig. 6(a) under a uniaxial compression for instance. As shown in Fig. 12(a), there is a bandgap between the 6th and 7th bands with low midgap position and wide gap width at the initiation of the loading. Additionally, with the increase of the load, the frequency of the lower boundary of this bandgap increases while that of the upper boundary decreases, except for some particular and small periods of the load magnitude. This result indicates that we can tune the width of this bandgap sharply through large deformation. Moreover, we calculate the ratio of the gap width fw to the midgap frequency fc, namely, the relative bandgap width fr , which is an important parameter in the design of acoustic devices [35]. The relation between the relative bandgap width fr and the applied load l is shown in Fig. 12(b), which is obtained according to the results in Fig. 12(a). It is inconceivable that the relative bandgap width fr of the bandgap between the 6th and 7th bands is more than 100% when the applied load l is less than 8.5%. Although it generally decreases with the increase of the applied load, it is always larger than those of the other bandgaps as shown in Fig. 12(b).

5. Conclusions In this paper, we performed numerical simulations as well as experiments to study the effects of inclusions on the buckling modes, post-buckling deformations, and band structures in soft porous periodic structures. We also numerically investigated the tunability of bandgaps of the structures undergoing large deformation. We now arrive at the following conclusions: 1. The arrangement and the amount of inclusions have remarkable effects on the bandgap characteristics of the porous structure. A reasonable arrangement of the inclusions can lead to ultra-broad and ultra-low bandgaps. In addition, resonant modes may exist in a porous structure with inclusions, which can be altered by properly changing the filling pattern. Hence, the property of unidirectional suppression of elastic waves can be tuned as well. 2. The hard inclusions placed in the pores will affect the buckling modes of the structure and the chiral mode can be triggered more easily. The inclusions will restrain the transition of buckling modes and make the post-buckling deformation insensitive to the initial geometrical imperfections. Therefore, the post-buckling deformation could develop in a robust and controllable way. 3. The elastic bandgaps in the porous structure can be tuned in a very fruitful manner by harnessing inclusions. In general, the bandgap characteristics are more sensitive to the applied load, e.g., some new bandgaps will be opened when the load varies. These observations could be a useful reference for the design of phononic switches for elastic waves. 4. The damping in soft materials has remarkable suppressive effect on elastic waves propagating in passbands at high frequencies while it has negligible influence on waves in the low frequency region. Thus, for the actual soft materials, the frequency range where the tunability of bandgaps could be effectively exploited is restricted. Nevertheless, making use of damping reasonably will pave a new way for designing acoustic devices.

Acknowledgments The work was supported by the National Natural Science Foundation of China [grant numbers 11532001, 11621062 and 11872329]. Partial supports from the Fundamental Research Funds for the Central Universities [grant number 2016XZZX00105] and the Shenzhen Scientific and Technological Fund for R & D [grant number JCYJ20170816172316775] are also acknowledged.

Appendix A. Simulation details Fig. A.1 shows five cases for carrying out the frequency domain analysis in Abaqus, including three structures (Figs. A.1(a), A.1(d) and A.1(e)) with periodic boundary conditions in a certain direction and two finite structures (Figs. A.1(b) and A.1(c)). In Fig. A.1, a sinusoidal mechanical load is applied at the left (or bottom) sides highlighted by red lines and the average acceleration amplitude (Ainput ) at each side will be calculated as the input signal. Besides, the right (or top) sides highlighted by red lines will be seen as the receiving ends and the average acceleration amplitude (Aoutput ) at each side will be calculated as the output signal. In most cases, FE models consisting of 10
1 (or 1
10) RVEs illustrated in Figs. A.1(a), A.1(d) and A.1(e) are considered in our calculations to compare with the band diagram. However, it is difficult to apply periodic boundary conditions in experiments since the sample is finite in both the horizontal and vertical directions. Therefore, the models shown in Figs. A.1(b) and A.1(c) are adopted for the steady-state response analysis to compare with the experimental results.

12

J. Li et al. / Journal of Sound and Vibration 459 (2019) 114848

Fig. A.1. FE models for frequency domain analysis: (a) Model comprising 10
1 RVEs (periodic boundary conditions are applied to the horizontal (top and bottom) edges) for calculating the transmittance in the G-X direction; (b) Model comprising 5
3 RVEs mimicking the experimental set-up shown in Fig. B.3(a); (c) Model comprising 5
3 RVEs mimicking the experimental set-up shown in Fig. B.3(b); (d) Deformed model comprising 10
1 RVEs (periodic boundary conditions are applied to the top and bottom edges) for calculating the transmittance in the G-X direction; (e) Model comprising 1
10 RVEs (periodic boundary conditions are applied to the left and right edges) for calculating the transmittance in the G-Y direction.

Appendix B. Experimental details B.1. DMA test In order to explore the effect of material damping on wave propagation in the composite structure, dynamic mechanical analysis (DMA) test was conducted to quantify the damping in rubber which is the soft material used in experiments. The storage modulus, loss modulus and the loss tangent (tan d) are shown in Fig. B.1 as functions of temperature. Here, we employ the loss tangent at room temperature (25 + C) as the loss factor h [36,37] (which is about 0.061) in the numerical simulations when damping is taken into consideration.

Fig. B.1. DMA test of the rubber sample.

B.2. Frequency response measurements To experimentally investigate the propagation behavior of elastic waves in the soft periodic structures, tests on two undeformed samples were conducted. The devices (see Fig. B.2) used for tests were placed on the optical table (POT18-12, LSXPT). Four situations of frequency response measurements as shown in Fig. B.3 were considered to compare with the numerical results in Fig. 2(a) and (b). Firstly, to reduce the effects of gravity on the in-plane deformation and the in-plane wave propagation, each sample was laid on the metal block shown in Fig. B.2. There were four nuts under the four

J. Li et al. / Journal of Sound and Vibration 459 (2019) 114848

13

corners of the sample to reduce the contact area between the sample and the support as small as possible. In order to stimulate the sample, a signal generator (DG4162, RIGOL) was used to generate the Gaussian white noise which was further amplified by the power amplifier (2718, Brüel & Kjær) to control the shaker (4809, Brüel & Kjær). To increase the contact area between the shaker and the sample so as to generate a uniform in-plane wave as far as possible, a PLA (Polylactic acid) block was placed between the shaker and the sample. Additionally, the contact interface was made bonded using glue. Besides, an accelerometer (4507 B, Brüel & Kjær) was stuck onto the top of the PLA block to record the input signals. On the receiving end of the sample, a thin PTFE (Polytetrafluorethylene) block was stuck onto the sample. The PTFE block and another accelerometer (4518-003, Brüel & Kjær) which was used to record the output signals were also glued together. Then the input and output signals were acquired through a data acquisition module (3050-A-060, Brüel & Kjær) and processed by the computer. The transmission coefficient T is finally calculated according to Eq. (7).

Fig. B.2. Devices for wave propagation test.

14

J. Li et al. / Journal of Sound and Vibration 459 (2019) 114848

Fig. B.3. Four situations of frequency response measurements: (a) Measurements in the G-X direction corresponding to the RVE in Fig. 2(a); (b) Measurements in the G-Y direction corresponding to the RVE in Fig. 2(a); (c) Measurements in the G-X direction corresponding to the RVE in Fig. 2(b); (d) Measurements in the G-Y direction corresponding to the RVE in Fig. 2(b).

B.3. Uniaxial and equi-biaxial compression In order to carry out the uniaxial compression test, a device was designed as shown in Fig. B.4(a). It consists of a rail with two movable blocks and a digital caliper composing a digital display (see the inset in Fig. B.4(d)) for the measurement of the relative displacement between the two blocks. The resolution of the digital caliper is 0.01 mm. The figure displayed on the digital display gives the displacement of one block, and the relative displacement is the double of this figure due to the symmetry. For equi-biaxial compression, another type of movable blocks was designed as shown in Figs. B.4(b), B.4(c) and B.4(d), which consists of a pair of vertical blocks and a pair of horizontal blocks. The vertical and the horizontal blocks were designed such that the same nominal strain is applied to the sample simultaneously through them. Actually, there will be inevitable fabrication error associated with the movable blocks, making it difficult to achieve the exact equi-biaxial compression state. Nonetheless, we find the fabrication error is very small, and its effect on the experimental results is negligible. During the tests, lubricant was used on the contact surfaces between the movable blocks and the rubber sample to reduce the friction and the boundary effects as less as possible.

J. Li et al. / Journal of Sound and Vibration 459 (2019) 114848

15

Fig. B.4. Devices for uniaxial and equi-biaxial compression tests: (a) Device for uniaxial compression test; (b) Device for equi-biaxial compression test; (c) Axonometric drawing of the device for equi-biaxial compression; (d) Process of the equi-biaxial compression test.

References [1] M.S. Kushwaha, P. Halevi, L. Dobrzynski, B. Djafari-Rouhani, Acoustic band structure of periodic elastic composites, Phys. Rev. Lett. 71 (1993) 2022e2025. https://doi.org/10.1103/PhysRevLett.71.2022. [2] M.M. Sigalas, Defect states of acoustic waves in a two-dimensional lattice of solid cylinders, J. Appl. Phys. 84 (1998) 3026e3030. https://doi.org/10. 1063/1.368456. [3] X. Zhang, Z. Liu, Negative refraction of acoustic waves in two-dimensional phononic crystals, Appl. Phys. Lett. 85 (2004) 341e343. https://doi.org/10. 1063/1.1772854. [4] B.C. Gupta, Z. Ye, Theoretical analysis of the focusing of acoustic waves by two-dimensional sonic crystals, Phys. Rev. E. 67 (2003) 36603. https://doi. org/10.1103/PhysRevE.67.036603. nchez-Pe rez, L.M. Garcia-Raffi, Tunable wideband bandstop acoustic filter based on two-dimensional multiphysical phe[5] V. Romero-García, J.V. Sa nomena periodic systems, J. Appl. Phys. 110 (2011) 14904. https://doi.org/10.1063/1.3599886. [6] J.O. Vasseur, M. Beaugeois, B. Djafari-Rouhani, Y. Pennec, P.A. Deymier, Transmission of acoustic waves through waveguide structures in twodimensional phononic crystals, Phys. Status Solidi. 1 (2004) 2720e2724. https://doi.org/10.1002/pssc.200405287. [7] Y. Kanno, K. Tsuruta, K. Fujimori, H. Fukano, S. Nogi, Phononic-crystal acoustic lens by design for energy-transmission devices, Electron. Commun. Jpn. 97 (2013) 22e27. https://doi.org/10.1002/ecj.11491. [8] F. Casadei, L. Dozio, M. Ruzzene, K.A. Cunefare, Periodic shunted arrays for the control of noise radiation in an enclosure, J. Sound Vib. 329 (2010) 3632e3646. https://doi.org/10.1016/j.jsv.2010.04.003. [9] N.S. Gao, H. Hou, Y.H. Mu, Low frequency acoustic properties of bilayer membrane acoustic metamaterial with magnetic oscillator, Theor. Appl. Mech. Lett. 7 (2017) 252e257. https://doi.org/10.1016/j.taml.2017.06.001. [10] N.S. Gao, H. Hou, Sound absorption characteristic of micro-helix metamaterial by 3D printing, Theor. Appl. Mech. Lett. 8 (2018) 63e67. https://doi.org/ 10.1016/j.taml.2018.02.001. [11] S. Mohammadi, A.A. Eftekhar, A. Khelif, H. Moubchir, R. Westafer, W.D. Hunt, A. Adibi, Complete phononic bandgaps and bandgap maps in twodimensional silicon phononic crystal plates, Electron. Lett. 43 (2007) 898e899. https://doi.org/10.1049/el:20071159. [12] C. Goffaux, J.P. Vigneron, Theoretical study of a tunable phononic band gap system, Phys. Rev. B. 64 (2001) 75118. https://doi.org/10.1103/PhysRevB.64. 075118. [13] Y. Liu, J.Y. Su, Y.L. Xu, X.C. Zhang, The influence of pore shapes on the band structures in phononic crystals with periodic distributed void pores, Ultrasonics 49 (2009) 276e280. https://doi.org/10.1016/j.ultras.2008.09.008. [14] Y. Chen, T. Li, F. Scarpa, L. Wang, Lattice metamaterials with mechanically tunable Poisson's ratio for vibration control, Phys. Rev. Appl. 7 (2017) 24012. https://doi.org/10.1103/PhysRevApplied.7.024012. [15] Y. Huang, C.L. Zhang, W.Q. Chen, Tuning band structures of two-dimensional phononic crystals with biasing fields, J. Appl. Mech. 81 (2014) 91008. https://doi.org/10.1115/1.4027915. [16] K. Bertoldi, M.C. Boyce, Mechanically triggered transformations of phononic band gaps in periodic elastomeric structures, Phys. Rev. B. 77 (2008) 52105. https://doi.org/10.1103/PhysRevB.77.052105. [17] K. Bertoldi, M.C. Boyce, Wave propagation and instabilities in monolithic and periodically structured elastomeric materials undergoing large deformations, Phys. Rev. B. 78 (2008) 1e16. https://doi.org/10.1103/PhysRevB.78.184107. [18] K. Bertoldi, M.C. Boyce, S. Deschanel, S.M. Prange, T. Mullin, Mechanics of deformation-triggered pattern transformations and superelastic behavior in periodic elastomeric structures, J. Mech. Phys. Solids 56 (2008) 2642e2668. https://doi.org/10.1016/j.jmps.2008.03.006. [19] K. Bertoldi, Stability of periodic porous structures, in: D. Bigoni (Ed.), Extrem. Deform. Struct., Springer Vienna, 2015, pp. 157e177. https://doi.org/10. 1007/978-3-7091-1877-1_4. [20] J. Shim, P. Wang, K. Bertoldi, Harnessing instability-induced pattern transformation to design tunable phononic crystals, Int. J. Solids Struct. 58 (2015) 52e61. https://doi.org/10.1016/j.ijsolstr.2014.12.018.

16

J. Li et al. / Journal of Sound and Vibration 459 (2019) 114848

[21] Y.L. Huang, N. Gao, W.Q. Chen, R.H. Bao, Extension/compression-controlled complete band gaps in 2D chiral square-lattice-like structures, Acta Mech. Solida Sin. 31 (2018) 51e65. https://doi.org/10.1007/s10338-018-0004-z. [22] B. Wu, W.J. Zhou, R.H. Bao, W.Q. Chen, Tuning elastic waves in soft phononic crystal cylinders via large deformation and electromechanical coupling, J. Appl. Mech. 85 (2018) 31004. https://doi.org/10.1115/1.4038770. [23] N. Gao, Y.L. Huang, R.H. Bao, W.Q. Chen, Robustly tuning bandgaps in two-dimensional soft phononic crystals with criss-crossed elliptical holes, Acta Mech. Solida Sin. 31 (2018) 573e588. https://doi.org/10.1007/s10338-018-0044-4. [24] J. Li, V. Slesarenko, S. Rudykh, Auxetic multiphase soft composite material design through instabilities with application for acoustic metamaterials, Soft Matter 14 (2018) 6171e6180. https://doi.org/10.1039/C8SM00874D. [25] S. Shan, S.H. Kang, P. Wang, C. Qu, S. Shian, E.R. Chen, K. Bertoldi, Harnessing multiple folding mechanisms in soft periodic structures for tunable control of elastic waves, Adv. Funct. Mater. 24 (2014) 4935e4942. https://doi.org/10.1002/adfm.201400665. mes Simulia Corp., Proidence, 2014. [26] Abaqus Analysis User's Guide, Dassault Syste [27] Y.L. Huang, N. Gao, R.H. Bao, W.Q. Chen, A novel algorithm for post-buckling analysis of periodic structures and its application, Chin. J. Comput. Mech. 33 (2016) 509e515 (in Chinese), https://doi.org/10.7511/jslx201604014. [28] G.A. Holzapfel, Nonlinear Solid Mechanics: A Continuum Approach for Engineering Science, Wiley, Chichester, 2000. [29] R.W. Ogden, Non-linear Elastic Deformations, Dover Publications, New York, 1997. [30] J.D. Achenbach, Wave Propagation in Elastic Solids, North-Holland, Amsterdam, 1973. [31] M. Åberg, P. Gudmundson, The usage of standard finite element codes for computation of dispersion relations in materials with periodic microstructure, J. Acoust. Soc. Am. 102 (1997) 2007e2013. https://doi.org/10.1121/1.419652. [32] J.R. Raney, N. Nadkarni, C. Daraio, D.M. Kochmann, J.A. Lewis, K. Bertoldi, Stable propagation of mechanical signals in soft media using stored elastic energy, Proc. Natl. Acad. Sci. 113 (2016) 9722e9727. https://doi.org/10.1073/pnas.1604838113. [33] Y.J. Chen, Y. Huang, C.F. Lü, W.Q. Chen, A two-way unidirectional narrow-band acoustic filter realized by a graded phononic crystal, J. Appl. Mech. 84 (2017) 091003. https://doi.org/10.1115/1.4037148. [34] F. Javid, P. Wang, A. Shanian, K. Bertoldi, Architected materials with ultra-low porosity for vibration control, Adv. Mater. (2016) 5943e5948. https://doi. org/10.1002/adma.201600052. [35] Z. Jia, Y. Chen, H. Yang, L. Wang, Designing phononic crystals with wide and robust band gaps, Phys. Rev. Appl. 9 (2018) 44021. https://doi.org/10.1103/ PhysRevApplied.9.044021. [36] M.V. Barnhart, X. Xu, Y. Chen, S. Zhang, J. Song, G. Huang, Experimental demonstration of a dissipative multi-resonator metamaterial for broadband elastic wave attenuation, J. Sound Vib. 438 (2019) 1e12. https://doi.org/10.1016/j.jsv.2018.08.035. [37] Y. Chen, F. Qian, L. Zuo, F. Scarpa, L. Wang, Broadband and multiband vibration mitigation in lattice metamaterials with sinusoidally-shaped ligaments, Extrem. Mech. Lett. 17 (2017) 24e32. https://doi.org/10.1016/j.eml.2017.09.012.