Design of Stretch-Dominated Metamaterials Avoiding Bandgap Resonance

1. Introduction

Owing to the ability to isolate mechanical waves within specific frequency ranges, bandgap mechanical metamaterials are broadly applicable in engineering vibration and sound insulation. Vasconcelos et al. [1] proposed a metamaterial-based interface to attenuate pressure waves induced by the impact hammer acting on the top of an offshore monopile. Zuo et al. [2] applied a star-shaped metamaterial to the shell of an underwater vehicle to reduce engine noise emission. A lattice-structured metamaterial capable of both sound insulation and ventilation was developed by Li et al. [3], and can be used as a roadside noise barrier. Seismic metamaterials composed of periodically arranged building foundations have been validated for shielding seismic waves and train-induced vibrations [4]. Replacing conventional aggregates with local resonant masses [5] enables concrete structural members to function as metamaterials capable of dissipating elastic waves induced by dynamic loads.

Dynamic loads within the bandgap may still induce significant structural responses in the metamaterials, which is known as bandgap resonance [6]. For example, Zhang et al. [7] designed a negative stiffness metamaterial that exhibits a peak in its frequency response curve within the bandgap range. Jiang et al. [8] designed negative stiffness metamaterials whose frequency response curves do not exactly match their bandgaps. The chiral and hexagonal lattices proposed by Li et al. [9,10] also experience large frequency responses within the bandgap ranges. Bandgap resonance can similarly be observed in the frequency response curves of the metamaterials studied in Ref. [11,12,13,14]. The presence of bandgap resonance suggests that certain eigenfrequencies of the metamaterial lie within the bandgap range, thereby compromising its vibration isolation performance. Therefore, it is necessary to prevent this phenomenon in the engineering design of bandgap metamaterials.

Bandgap resonance attracted early attention in the field of solid-state physics. In calculating the specific heat capacity of crystals, periodically arranged atoms are modelled as vibrating one-dimensional atomic chains. Thus, some early studies focused on the correspondence between the eigenfrequency spectra and the dispersion functions of finite atomic chains. Born [15] studied two types of one-dimensional atomic chains and concluded that the spectrum of a finite atomic chain with free ends is consistent with its dispersion function. Wallis [16], however, identified an eigenfrequency within the bandgap while analyzing the spectrum of a finite diatomic chain. Since then, similar phenomena have been widely observed in elastomers [17,18,19,20].

Calculating the dispersion function via Bloch’s theorem assumes that the metamaterials exhibit translational symmetry [21], necessitating that they be either infinitely extended or subjected to periodic boundary conditions [22]. However, finite metamaterials in engineering possess non-periodic boundaries, so their spectra deviate from the dispersion functions. This suggests that the existence of bandgap resonance is determined by the boundary conditions. Bastawrous et al. [6] introduced boundary conditions as perturbations to the dynamic stiffness matrix and analytically derived a closed-form condition for the existence of bandgap resonance in finite diatomic chains. Guo et al. [23] and Sugino et al. [24] solved dynamic differential equations with boundary conditions substituted to obtain the spectra of finite two-phase plates. Jin et al. [25] incorporated boundary conditions in the form of a Fourier series into the dynamic stiffness matrix and solved the spectrum of a two-phase metamaterial plate. By comparing the obtained spectra with the dispersion functions of the metamaterials, the existence of bandgap resonance can be identified.

Analytical calculations of the spectra of metamaterials require explicit expressions for the eigenfrequencies, which poses a substantial mathematical challenge. As a result, existing studies have primarily focused on simple systems, such as two-phase plates or diatomic chains [6,23,24,25]. However, metamaterials designed for specific bandgaps are usually complex, making it difficult to analytically solve their spectra. Therefore, to ensure the vibration isolation performance of a designed metamaterial, bandgap resonance should preferably be avoided without solving the spectrum.

This paper proposes a design method for stretch-dominated metamaterials that achieve expected bandgaps while avoiding bandgap resonance. The paper is organized as follows: Section 2 derives the sufficient condition for avoiding the bandgap resonance in stretch-dominated metamaterials. Section 3 presents the equations for designing the unit cell parameters under symmetry constraint. A design example is illustrated in Section 4, which is used to validate both the sufficient condition and the corresponding design method. Finally, Section 5 draws some conclusions.

2. Sufficient Condition for Avoiding Bandgap Resonance

A two-dimensional stretch-dominated metamaterial is used as an example to demonstrate the derivation of the sufficient condition for avoiding bandgap resonance. As shown in Figure 1, a 3×3 metamaterial, denoted by m, is modelled as a pin-jointed bar structure and subjected to a boundary condition that can be practically realized in engineering. Figure 2a and Figure 3a show the distributions of modal forces and displacements on the boundaries corresponding to any eigenfrequency ω_j of m. Except for the corner nodes, the horizontal components of the node forces are zero on the top and bottom edges, while the vertical components are zero on the left and right edges. Similarly, the vertical components of the node displacements are zero on the top and bottom edges, and the horizontal components are zero on the left and right edges. Thus, it can be derived that

$$d_{\mathrm{t}/\mathrm{b}}={[0,0,d_{\mathrm{t}/\mathrm{b}}^{2x},0,\mathrm{...},d_{\mathrm{t}/\mathrm{b}}^{(n_{\mathrm{t}/\mathrm{b}}-1)x},0,0,0]}^{\mathrm{T}}$$

(1)

$$d_{\mathrm{l}/\mathrm{r}}={[0,0,0,d_{\mathrm{l}/\mathrm{r}}^{2y},\mathrm{...},0,d_{\mathrm{l}/\mathrm{r}}^{(n_{\mathrm{l}/\mathrm{r}}-1)y},0,0]}^{\mathrm{T}}$$

(2)

$$f_{\mathrm{t}/\mathrm{b}}={[f_{\mathrm{t}/\mathrm{b}}^{1x},f_{\mathrm{t}/\mathrm{b}}^{1y},0,f_{\mathrm{t}/\mathrm{b}}^{2y},\mathrm{...},0,f_{\mathrm{t}/\mathrm{b}}^{(n_{\mathrm{t}/\mathrm{b}}-1)y},f_{\mathrm{t}/\mathrm{b}}^{(n_{\mathrm{t}/\mathrm{b}})x},f_{\mathrm{t}/\mathrm{b}}^{(n_{\mathrm{t}/\mathrm{b}})y}]}^{\mathrm{T}}$$

(3)

$$f_{\mathrm{l}/\mathrm{r}}={[f_{\mathrm{l}/\mathrm{r}}^{1x},f_{\mathrm{l}/\mathrm{r}}^{1y},f_{\mathrm{l}/\mathrm{r}}^{2x},0,\mathrm{...},f_{\mathrm{l}/\mathrm{r}}^{(n_{\mathrm{l}/\mathrm{r}}-1)x},0,f_{\mathrm{l}/\mathrm{r}}^{(n_{\mathrm{l}/\mathrm{r}})x},f_{\mathrm{l}/\mathrm{r}}^{(n_{\mathrm{l}/\mathrm{r}})y}]}^{\mathrm{T}}$$

(4)

where d_t/b and f_t/b are the displacement and force vectors of the nodes on the top or bottom edge; d_l/r and f_l/r are the displacement and force vectors of the nodes on the left or right edge; $d_{t/b}^{ux}$ and $d_{t/b}^{uy}$ (u=1, …, n_t/b) are the x and y components of the displacement of the u-th node on the top or bottom edge, respectively, $f_{t/b}^{ux}$ and $f_{t/b}^{uy}$ (u=1, …, n_t/b) are the corresponding force components; $d_{l/r}^{ux}$ and $d_{l/r}^{uy}$ (u=1, …, n_l/r) are the x and y components of the displacement of the u-th node on the left or right edge, $f_{l/r}^{ux}$ and $f_{l/r}^{uy}$ (u=1, …, n_l/r) are the corresponding force components; n_t/b is the number of nodes on the top or bottom edge, and n_l/r is that on the left or right edge; the superscript T denotes the transpose of a vector or matrix.

Spatial inversion of m in the x-direction (i.e., mirroring across the y-axis) yields the metamaterial m^x. The distributions of modal forces and displacements on the boundaries are shown in Figure 2b and Figure 3b, and the corresponding frequency remains constant at ω_j. According to Figure 2b and Figure 3b, the boundary node displacement and force vectors of m^x can be written as

$$d_{\mathrm{t}/\mathrm{b}}^{\mathrm{x}}={[0,0,-d_{\mathrm{t}/\mathrm{b}}^{(n_{\mathrm{t}/\mathrm{b}}-1)x},0,\mathrm{...},-d_{\mathrm{t}/\mathrm{b}}^{2x},0,0,0]}^{\mathrm{T}}$$

(5)

$$d_{\mathrm{r}/\mathrm{l}}^{\mathrm{x}}={[0,0,0,d_{\mathrm{l}/\mathrm{r}}^{2y},\mathrm{...},0,d_{\mathrm{l}/\mathrm{r}}^{(n_{\mathrm{l}/\mathrm{r}}-1)y},0,0]}^{\mathrm{T}}$$

(6)

$$f_{\mathrm{t}/\mathrm{b}}^{\mathrm{x}}={[-f_{\mathrm{t}/\mathrm{b}}^{(n_{\mathrm{t}/\mathrm{b}})x},f_{\mathrm{t}/\mathrm{b}}^{(n_{\mathrm{t}/\mathrm{b}})y},0,f_{\mathrm{t}/\mathrm{b}}^{(n_{\mathrm{t}/\mathrm{b}}-1)y},\mathrm{...},0,f_{\mathrm{t}/\mathrm{b}}^{2y},-f_{\mathrm{t}/\mathrm{b}}^{1x},f_{\mathrm{t}/\mathrm{b}}^{1y}]}^{\mathrm{T}}$$

(7)

$$f_{\mathrm{r}/\mathrm{l}}^{\mathrm{x}}={[-f_{\mathrm{l}/\mathrm{r}}^{1x},f_{\mathrm{l}/\mathrm{r}}^{1y},-f_{\mathrm{l}/\mathrm{r}}^{2x},0,\mathrm{...},-f_{\mathrm{l}/\mathrm{r}}^{(n_{\mathrm{l}/\mathrm{r}}-1)x},0,-f_{\mathrm{l}/\mathrm{r}}^{(n_{\mathrm{l}/\mathrm{r}})x},f_{\mathrm{l}/\mathrm{r}}^{(n_{\mathrm{l}/\mathrm{r}})y}]}^{\mathrm{T}}$$

(8)

Spatial inversion of m in the y-direction (i.e., mirroring across the x-axis) yields the metamaterial m^y. The distributions of modal forces and displacements on the boundaries are shown inFigure 2c and Figure 3c, and the corresponding frequency remains constant at ω_j. According to Figure 2c and Figure 3c, the boundary node displacement and force vectors of m^y can be written as

$$d_{\mathrm{b}/\mathrm{t}}^{\mathrm{y}}={[0,0,d_{\mathrm{t}/\mathrm{b}}^{2x},0,\mathrm{...},d_{\mathrm{t}/\mathrm{b}}^{(n_{\mathrm{t}/\mathrm{b}}-1)x},0,0,0]}^{\mathrm{T}}$$

(9)

$$d_{\mathrm{l}/\mathrm{r}}^{\mathrm{y}}={[0,0,0,-d_{\mathrm{l}/\mathrm{r}}^{(n_{\mathrm{l}/\mathrm{r}}-1)y},\mathrm{...},0,-d_{\mathrm{l}/\mathrm{r}}^{2y},0,0]}^{\mathrm{T}}$$

(10)

$$f_{\mathrm{b}/\mathrm{t}}^{\mathrm{y}}={[f_{\mathrm{t}/\mathrm{b}}^{1x},-f_{\mathrm{t}/\mathrm{b}}^{1y},0,-f_{\mathrm{t}/\mathrm{b}}^{2y},\mathrm{...},0,-f_{\mathrm{t}/\mathrm{b}}^{(n_{\mathrm{t}/\mathrm{b}}-1)y},f_{\mathrm{t}/\mathrm{b}}^{(n_{\mathrm{t}/\mathrm{b}})x},-f_{\mathrm{t}/\mathrm{b}}^{(n_{\mathrm{t}/\mathrm{b}})y}]}^{\mathrm{T}}$$

(11)

$$f_{\mathrm{l}/\mathrm{r}}^{\mathrm{y}}={[f_{\mathrm{l}/\mathrm{r}}^{(n_{\mathrm{l}/\mathrm{r}})x},-f_{\mathrm{l}/\mathrm{r}}^{(n_{\mathrm{l}/\mathrm{r}})y},f_{\mathrm{l}/\mathrm{r}}^{(n_{\mathrm{l}/\mathrm{r}}-1)x},0,\mathrm{...},f_{\mathrm{l}/\mathrm{r}}^{2x},0,f_{\mathrm{l}/\mathrm{r}}^{1x},-f_{\mathrm{l}/\mathrm{r}}^{1y}]}^{\mathrm{T}}$$

(12)

Spatial inversion of m^x in the y-direction yields the metamaterial m^xy. The distributions of modal forces and displacements on the boundaries are shown in Figure 2d and Figure 3d, and the corresponding frequency remains constant at ω_j. According to Figure 2d and Figure 3d, the boundary node displacement and force vectors of m^xy can be written as

$$d_{\mathrm{b}/\mathrm{t}}^{\mathrm{xy}}={[0,0,-d_{\mathrm{t}/\mathrm{b}}^{(n_{\mathrm{t}/\mathrm{b}}-1)x},0,\mathrm{...},-d_{\mathrm{t}/\mathrm{b}}^{2x},0,0,0]}^{\mathrm{T}}$$

(13)

$$d_{\mathrm{r}/\mathrm{l}}^{\mathrm{xy}}={[0,0,0,-d_{\mathrm{l}/\mathrm{r}}^{(n_{\mathrm{l}/\mathrm{r}}-1)y},\mathrm{...},0,-d_{\mathrm{l}/\mathrm{r}}^{2y},0,0]}^{\mathrm{T}}$$

(14)

$$f_{\mathrm{b}/\mathrm{t}}^{\mathrm{xy}}={[-f_{\mathrm{t}/\mathrm{b}}^{(n_{\mathrm{t}/\mathrm{b}})x},-f_{\mathrm{t}/\mathrm{b}}^{(n_{\mathrm{t}/\mathrm{b}})y},0,-f_{\mathrm{t}/\mathrm{b}}^{(n_{\mathrm{t}/\mathrm{b}}-1)y},\mathrm{...},0,-f_{\mathrm{t}/\mathrm{b}}^{2y},-f_{\mathrm{t}/\mathrm{b}}^{1x},-f_{\mathrm{t}/\mathrm{b}}^{1y}]}^{\mathrm{T}}$$

(15)

$$f_{\mathrm{r}/\mathrm{l}}^{\mathrm{xy}}={[-f_{\mathrm{l}/\mathrm{r}}^{(n_{\mathrm{l}/\mathrm{r}})x},-f_{\mathrm{l}/\mathrm{r}}^{(n_{\mathrm{l}/\mathrm{r}})y},-f_{\mathrm{l}/\mathrm{r}}^{(n_{\mathrm{l}/\mathrm{r}}-1)x},0,\mathrm{...},-f_{\mathrm{l}/\mathrm{r}}^{2x},0,-f_{\mathrm{l}/\mathrm{r}}^{1x},-f_{\mathrm{l}/\mathrm{r}}^{1y}]}^{\mathrm{T}}$$

(16)

As can be seen from Figure 2 and Eqs. (1), (2), (5), (6), (9), (10), (13) and (14), the boundary node displacement vectors satisfy

$$d_{\mathrm{r}}=d_{\mathrm{l}}^{\mathrm{x}}\text{，}{d}_{\mathrm{r}}^{\mathrm{y}}=d_{\mathrm{l}}^{\mathrm{xy}}\text{，}{d}_{\mathrm{b}}=d_{\mathrm{t}}^{\mathrm{y}}\text{，}{d}_{\mathrm{b}}^{\mathrm{x}}=d_{\mathrm{t}}^{\mathrm{xy}}$$

(17)

Moreover, Figure 3 and Eqs. (3), (4), (7), (8), (11), (12), (15) and (16) indicate that the boundary node force vectors satisfy the following relationships except for corner nodes:

$$f_{\mathrm{r}}+f_{\mathrm{l}}^{\mathrm{x}}=0\text{，}{f}_{\mathrm{r}}^{\mathrm{y}}+f_{\mathrm{l}}^{\mathrm{xy}}=0\text{，}{f}_{\mathrm{b}}+f_{\mathrm{t}}^{\mathrm{y}}=0\text{，}{f}_{\mathrm{b}}^{\mathrm{x}}+f_{\mathrm{t}}^{\mathrm{xy}}=0$$

(18)

Let nodes 1, 2, 3 and 4 denote the lower right corner node of m, the lower left corner node of m^x, the upper right corner node of m^y, and the upper left corner node of m^xy, respectively. Then, their node forces satisfy

$$f_1+f_2+f_3+f_4=\left\{\begin{array}{c}{f}_b^{x{n}_{\mathrm{b}}}\\ f_b^{y{n}_b}\end{array}\right\}+\left\{\begin{array}{c}-f_b^{x{n}_{\mathrm{b}}}\\ f_b^{y{n}_b}\end{array}\right\}+\left\{\begin{array}{c}{f}_b^{x{n}_{\mathrm{b}}}\\ -f_b^{y{n}_b}\end{array}\right\}+\left\{\begin{array}{c}-f_b^{x{n}_{\mathrm{b}}}\\ -f_b^{y{n}_b}\end{array}\right\}=0$$

(19)

Eqs. (17)– (19) indicate that the corresponding nodes on the bottom edge of m and the top edge of m^y have identical displacements, and their node forces are in equilibrium. The same relationship holds between the right edge of m and the left edge of m^x, the top edge of m^xy and the bottom edge of m^x, and the right edge of m^y and the left edge of m^xy. This suggests that m, m^x, m^y, and m^xy can be merged into a new structure M. Figure 4a and Figure 4b show the distributions of its modal forces and displacements on the boundaries, and the corresponding frequency is still ω_j.

The boundary node displacement and force vectors of M can be written as

$$d_{\mathrm{t}}^{\mathrm{M}}={[0,0,d_{\mathrm{t}}^{2x},0,\mathrm{...},d_{\mathrm{t}}^{(n_{\mathrm{t}}-1)x},0,0,0,-d_{\mathrm{t}}^{(n_{\mathrm{t}}-1)x},0,\mathrm{...},-d_{\mathrm{t}}^{2x},0,0,0]}^{\mathrm{T}}$$

(20)

$$d_{\mathrm{b}}^{\mathrm{M}}={[0,0,d_{\mathrm{t}}^{2x},0,\mathrm{...},d_{\mathrm{t}}^{(n_{\mathrm{t}}-1)x},0,0,0,-d_{\mathrm{t}}^{(n_{\mathrm{t}}-1)x},0,\mathrm{...},-d_{\mathrm{t}}^{2x},0,0,0]}^{\mathrm{T}}$$

(21)

$$d_{\mathrm{l}}^{\mathrm{M}}={[0,0,0,d_{\mathrm{l}}^{2y},\mathrm{...},0,d_{\mathrm{l}}^{(n_{\mathrm{l}}-1)y},0,0,0,-d_{\mathrm{l}}^{(n_{\mathrm{l}}-1)y},\mathrm{...},0,-d_{\mathrm{l}}^{2y},0,0]}^{\mathrm{T}}$$

(22)

$$d_{\mathrm{r}}^{\mathrm{M}}={[0,0,0,d_{\mathrm{l}}^{2y},\mathrm{...},0,d_{\mathrm{l}}^{(n_{\mathrm{l}}-1)y},0,0,0,-d_{\mathrm{l}}^{(n_{\mathrm{l}}-1)y},\mathrm{...},0,-d_{\mathrm{l}}^{2y},0,0]}^{\mathrm{T}}$$

(23)

$$f_{\mathrm{t}}^{\mathrm{M}}={[f_{\mathrm{t}}^{1x},f_{\mathrm{t}}^{1y},0,f_{\mathrm{t}}^{2y},\mathrm{...},,0,f_{\mathrm{t}}^{(n_{\mathrm{t}})y},\mathrm{...},0,f_{\mathrm{t}}^{2y},f_{\mathrm{t}}^{1x},f_{\mathrm{t}}^{1y}]}^{\mathrm{T}}$$

(24)

$$f_{\mathrm{b}}^{\mathrm{M}}=-{[-f_{\mathrm{t}}^{1x},f_{\mathrm{t}}^{1y},0,f_{\mathrm{t}}^{2y},\mathrm{...},,0,f_{\mathrm{t}}^{(n_{\mathrm{t}})y},\mathrm{...},0,f_{\mathrm{t}}^{2y},-f_{\mathrm{t}}^{1x},f_{\mathrm{t}}^{1y}]}^{\mathrm{T}}$$

(25)

$$f_{\mathrm{l}}^{\mathrm{M}}={[f_{\mathrm{l}}^{1x},f_{\mathrm{l}}^{1y},f_{\mathrm{l}}^{2x},0,\mathrm{...},f_{\mathrm{l}}^{(n_{\mathrm{l}/\mathrm{r}})x},0,\mathrm{...},f_{\mathrm{l}}^{2x},0,f_{\mathrm{l}}^{1x},f_{\mathrm{l}}^{1y}]}^{\mathrm{T}}$$

(26)

$$f_{\mathrm{r}}^{\mathrm{M}}=-{[-f_{\mathrm{l}}^{1x},f_{\mathrm{l}}^{1y},0,f_{\mathrm{l}}^{2y},\mathrm{...},0,f_{\mathrm{l}}^{(n_{\mathrm{l}/\mathrm{r}})y},\mathrm{...},0,f_{\mathrm{l}}^{2y},-f_{\mathrm{l}}^{1x},f_{\mathrm{l}}^{1y}]}^{\mathrm{T}}$$

(27)

According to Eqs. (20)– (27), the node displacements and forces on the boundaries of M satisfy the periodic boundary condition [26]. However, M may not be a translationally symmetric metamaterial because m is not necessarily identical to m^x, m^y and m^xy. Conversely, if m remains unchanged after the spatial inversions in both the x- and y-directions, M is a metamaterial satisfying the periodic boundary condition. This requires that the unit cell of m exhibits spatial inversion symmetry along x- and y-axes. In this case, the spectrum of M must coincide with the dispersion function obtained via Bloch’s theorem, indicating that none of the eigenfrequencies of M, including ω_j, lie within the bandgap. Since any eigenfrequency ω_j of m lies outside the bandgap, the metamaterial m can avoid bandgap resonance.

The above derivation can be generalized to three-dimensional metamaterials and some other boundary conditions. For two- and three-dimensional cases, the applicable boundary conditions are listed in Table 1 and Table 2, respectively.

In Table 1 and Table 2, x, y and z represent nodes on the boundaries normal to the x-, y- and z-directions, respectively; c₁, c₂ and c₃ represent constraints applied in the x-, y- and z- directions, respectively.

In summary, the sufficient condition for avoiding bandgap resonance can be express as

Theorem 1.A d-dimensional metamaterial can avoid bandgap resonance if both of the following conditions are satisfied: 1. The unit cell exhibits spatial inversion symmetry along all d coordinate axes; 2. One of the boundary conditions listed inTable 1orTable 2is imposed.

3. Perturbation of Bandgap

The following generalized characteristic equation [26,27,28,29] can be used to solve the dispersion function of the metamaterial consisting of unit cells with spatial inversion symmetry:

$$\widehat{D}(k)\widehat{d}(k)=0$$

(28)

where $\widehat{D}(k)=T_{\mathrm{d}}{(k)}^{\mathrm{H}}[K-\omega {(k)}^{2}M]T_{\mathrm{d}}(k)$; $K$ and $M$ are the d·n×d·n stiffness and mass matrices of the unit cell, respectively, and their explicit forms can be found in Ref. [29]; n is the number of nodes in the unit cell; the superscript H denotes the Hermite transpose; $T_{\mathrm{d}}(k)$ is the d·n×d·r reduction matrix introducing Bloch’s theorem [29], and r is the number of inequivalent nodes [29]; ω(k) is the circular frequency, and k is the d×1 wave vector [30]. Once k is given, the generalized eigenvalues ω_j(k)² (j=1, 2, …, d·r) and their corresponding eigenvectors ${\widehat{d}}_j(k)$ (d·r×1) can be solved via Eq. (28). Thus, the dispersion function of the metamaterial, i.e., the k‒ω relation, is established with a total of d·r branches of k‒ω_j.

According to matrix perturbation theory [28,29], the perturbation of branch j caused by the increments of unit cell parameters, including node coordinates and element cross-sectional areas, can be calculated by

$$\Delta {\omega }_j(k)=O_j(k)\Delta \epsilon$$

(29)

where $\Delta \epsilon ={\{A_1,A_2,\mathrm{...},A_b,x_1^1,\mathrm{...},x_1^d,\mathrm{...},x_n^1,\mathrm{...},x_n^d\}}^{\mathrm{T}}$ is the (b+d·n)×1 incremental unit cell parameter vector; b is the number of elements in the unit cell; A_i (i=1,…, b) is the cross-sectional area of the i-th element; $x_p^k$( p=1, …, n, k=1, .., d) is the k-th coordinate of the p-th node. O_j(k) is a 1×(b+d·n) row vector whose v-th component is

$${(O_j(k))}_v=\left\{\begin{array}{c}{\displaystyle \frac{1}{2{\omega }_j(k)}}{\widehat{d}}_{\overline{j}}{(k)}^{\mathrm{H}}{\displaystyle \frac{\partial \widehat{D}(k)}{\partial A_v}}{\widehat{d}}_{\overline{j}}(k),v\le b\\ {\displaystyle \frac{1}{2{\omega }_j(k)}}{\widehat{d}}_{\overline{j}}{(k)}^{\mathrm{H}}{\displaystyle \frac{\partial \widehat{D}(k)}{\partial x_p^k}}{\widehat{d}}_{\overline{j}}(k),v\text{>}b\end{array}\right.$$

(30)

where v=b+d·(p-1)+k for the second term; ${\widehat{d}}_{\overline{j}}(k)$ is the eigenvector normalized by ${\widehat{d}}_{\overline{j}}{(k)}^{\mathrm{H}}{T}_{\mathrm{d}}{(k)}^{\mathrm{H}}M{T}_{\mathrm{d}}(k){\widehat{d}}_{\overline{j}}(k)=1$; ${\displaystyle \frac{\partial \widehat{D}(k)}{\partial A_v}}=T_{\mathrm{d}}{(k)}^{\mathrm{H}}({\displaystyle \frac{\partial K}{\partial A_v}}-{\omega }_j{(k)}^2{\displaystyle \frac{\partial M}{\partial A_v}})T_{\mathrm{d}}(k)$, ${\displaystyle \frac{\partial \widehat{D}(k)}{\partial x_p^k}}=T_{\mathrm{d}}{(k)}^{\mathrm{H}}({\displaystyle \frac{\partial K}{\partial x_p^k}}-{\omega }_j{(k)}^2{\displaystyle \frac{\partial M}{\partial x_p^k}})T_{\mathrm{d}}(k)$. The explicit forms of ${\displaystyle \frac{\partial K}{\partial A_v}}$,${\displaystyle \frac{\partial K}{\partial x_p^k}}$, ${\displaystyle \frac{\partial M}{\partial A_v}}$ and ${\displaystyle \frac{\partial M}{\partial x_p^k}}$ are given in Ref. [29].

An illustrative bandgap between the lower branch ω_l(k) and the upper branch ω_u(k) is shown in Figure 5. The range of the bandgap is [ω_l(k_l), ω_u(k_u)], where k_l and k_u are the wave vectors corresponding to the maximum value of ω_l(k) and the minimum value of ω_u(k), respectively. Substituting k_l and k_u into Eq. (29), the perturbation of ω_l(k_l) and ω_u(k_u) can be expressed as

$$\Delta \omega =\left\{\begin{array}{c}\Delta {\omega }_{\mathrm{u}}(k_{\mathrm{u}})\\ \Delta {\omega }_{\mathrm{l}}(k_{\mathrm{l}})\end{array}\right\}=\left\{\begin{array}{c}{O}_{\mathrm{u}}(k_{\mathrm{u}})\\ O_{\mathrm{l}}(k_{\mathrm{l}})\end{array}\right\}\Delta \epsilon =O\Delta \epsilon$$

(31)

The difference between the expected bandgap [ω_l, ω_u] shown in Figure 5 and the current bandgap is Δω_t={ω_u-ω_u(k_u), ω_l-ω_l(k_l)}^T. By substituting Δω_t into Eq. (31), Δε for adjusting unit cell parameters and tuning bandgap can be solved. However, Δε must satisfy spatial inversion symmetry, or it will break the symmetry of the unit cell. Figure 6 shows a unit cell with 8 nodes and 12 elements. An increment vector adjusting the coordinates of nodes 5, 6, 7 and 8 with $\Delta x_5^1=1$, $\Delta x_6^1=-1$, $\Delta x_7^1=-1$ and $\Delta x_8^1=1$ satisfies spatial inversion symmetry, and is denoted as Δε_g. Similarly, ΔA₅=1, ΔA₆=1, ΔA₇=1 and ΔA₈=1 can be assembled into another Δε_g satisfying spatial inversion symmetry. If there are at most τ linearly independent Δε_g, then any incremental unit cell parameter vector satisfying spatial inversion symmetry can be expressed as

$$\Delta {\epsilon }_{\mathrm{sym}}=\{\Delta {\epsilon }_{g1},\Delta {\epsilon }_{g1},\mathrm{...}\Delta {\epsilon }_{g\tau }\}{\{{\kappa }_1,{\kappa }_2,\mathrm{...},{\kappa }_{\tau }\}}^{\mathrm{T}}={{\rm E}}_g\Delta \kappa$$

(32)

where Δκ is an arbitrary τ×1 coefficient vector.

Substituting Eq. (32) and Δω_t into Eq. (31), it can be obtained that

$$\Delta {\omega }_{\mathrm{t}}=O{{\rm E}}_g\Delta \kappa$$

(33)

The least squares solution [31] of Eq. (33) is

$$\Delta \kappa ={(O{{\rm E}}_g)}^{\mathrm{T}}{[(O{{\rm E}}_g){(O{{\rm E}}_g)}^{\mathrm{T}}]}^{-1}\Delta {\omega }_{\mathrm{t}}$$

(34)

By substituting Eq. (34) into Eq. (32), the resulting Δε_sym can be used to adjust the unit cell parameters and tune the bandgap under the constraint of spatial inversion symmetry.

The unit cell parameters should be modified in small steps due to the strong nonlinearity in their relationship with the dispersion function. Thus, Δε_sym is scaled by

$$\Delta {\epsilon }_{\mathrm{s}}=s{\displaystyle \frac{\Delta {\epsilon }_{\mathrm{sym}}}{\left|\right|\Delta {\epsilon }_{\mathrm{sym}}\left|\right|}}$$

(35)

where s is a specified scaling factor; ||•|| denotes the 2-norm of the vector. The expected bandgap can be obtained by iteratively adjusting the unit cell parameters with Δε_s until ||Δω_t|| falls below a specified tolerance c.

4. Example

Figure 7a shows a square unit cell with a side length of 10 mm that exhibits spatial inversion symmetry along the x- and y-axes. The origin of the coordinate system is at the center of the unit cell with the axes shown in Figure 7a. The node coordinates are listed in Table 3. All the elements of the unit cell have a circular cross-section with a radius of 0.3 mm (area A=0.283 mm²) and are made of plastic with a Young’s modulus of E=1.2 GPa and a density of ρ=0.9 g/cm³.

For this unit cell, the dispersion function obtained using Eq. (28) consists of 52 branches, whose significant overlap makes it difficult to depict their surfaces on a plane. To illustrate the band structure, a frequency distribution function for branch j can be constructed as follows:

$$D_j(\omega )={\displaystyle \underset{\Pi }{\iint }{\left.\mathrm{d}A\right|}_{{\omega }_j(k)=\omega }}(j=1,2,\mathrm{...},d\cdot r)$$

(36)

where Π denotes the first Brillouin zone [30]; dA=dk_xdk_y is the area of the infinitesimal element of Π; and k_x and k_y are the x and y components of k, respectively. According to Eq. (36), D_j(ω) represents the total area of the infinitesimal elements for which ω_j=ω. If D_j(ω)≠0, there must exist a wave vector k such that ω_j(k)=ω. Therefore, ω_j(k) lies within the frequency domain with nonzero D_j(ω). The frequency distribution of the dispersion function can be visualized by plotting D_j(ω) for all branches, with ω as the horizontal axis and area as the vertical axis. Figure 7b shows 27 lower-frequency branches distributed between 0 and 3100 Hz.

This example aims to achieve a bandgap of 600–1000 Hz by designing the coordinates of nodes 5–29 and the cross-sectional areas of all elements. According to Figure 7b, the expected bandgap is most likely generated between branches 3 and 4. The unit cell parameters are iteratively adjusted using Eqs. (31)– (35) with a tolerance of c=2. The resulting unit cell is shown in Fib. 8a, whose element cross-sectional areas and node coordinates are listed in Table 4 and Table 5, respectively. The frequency distribution of the updated dispersion function is shown in Figure 8b, where a bandgap of 600.0–998.6 Hz can be found.

A metamaterial consisting of 4×4 designed unit cells is shown in Figure 9a. The metamaterial is denoted as PBC, FBC or SBC when the periodic, free, or specific boundary condition is imposed. The specific boundary condition is illustrated in Figure 9b and corresponds to case 3 in Table 1. In addition, the symmetry of the designed unit cell can be slightly broken by adjusting the coordinate of node 15 to (-4.03, 1.00), resulting in a bandgap of 601.0–904.7 Hz. The 4×4 metamaterial, consisting of this asymmetric unit cell and subjected to the specific boundary condition, is denoted as SBCA.

The eigenfrequencies and frequency-response curves of PBC, FBC, SBC, and SBCA are shown in Figure 10, Figure 11, Figure 12 and Figure 13. The frequency-response function is defined as

$$y(\omega )=20{\log }_{10}({\displaystyle \frac{\left|\right|d_{\mathrm{out}}(\omega )\left|\right|}{\left|\right|d_{\mathrm{in}}(\omega )\left|\right|}})$$

(37)

where ||d_in(ω)|| is the displacement amplitude at input node 1, and ||d_out(ω)|| is that at output node 2 or 3. These nodes are labeled in Figure 9a. With a modal damping of 0.02, ||d_in(ω)|| and ||d_out(ω)|| are calculated using the modal superposition method (MS) [32] and the finite element software ANSYS, respectively. The MS is implemented in MATLAB, and the procedure includes the following steps: (1) Set up the dynamic equation for the metamaterial consisting of 4×4 unit cells. (2) Calculate the eigenfrequencies and eigenmodes of the metamaterial. (3) Calculate the displacements at the input and output nodes, i.e., d_in and d_out, under a simple harmonic load using the displacement equation [32] of the modal superposition method. The metamaterials are also modelled in ANSYS. All the elements are simulated by Link180. A harmonic response analysis is conducted to obtain the displacements d_in and d_out, with the solution method [33] set to FULL.

As illustrated in Figure 10, the eigenfrequencies and frequency-response curves of PBC are consistent with the bandgap, confirming that the metamaterial with expected bandgap has been successfully designed. However, when the metamaterial is subjected to the free boundary condition, eigenfrequencies emerge within the bandgap, as illustrated in Figure 11. Figure 12 shows that the bandgap resonances in SBC are effectively suppressed by imposing the specific boundary condition, thereby validating the second condition of Theorem 1. For the asymmetric metamaterial, even when the specific boundary condition is imposed, bandgap resonances still emerge in the frequency-response curves of Figure 13. This validates the first condition of Theorem 1.

5. Conclusions

This paper proposes a design method for stretch-dominated metamaterials that achieve expected bandgaps while avoiding bandgap resonance. The conclusions are summarized as follows:

(1) A d-dimensional metamaterial can avoid bandgap resonance if both of the following conditions are satisfied: 1. The unit cell exhibits spatial inversion symmetry along all d coordinate axes; 2. One of the boundary conditions listed in Table 1 and Table 2is imposed.

(2) The perturbation expression of the bandgap is derived. By introducing symmetry constraint, the incremental unit cell parameter vector can be solved via this perturbation expression. This vector is then used to adjust the node coordinates and element cross-sectional areas to tune the bandgap.

(3) A two-dimensional metamaterial example is used to validate the proposed sufficient condition and design method. The designed metamaterial, which satisfies spatial inversion symmetry and the boundary condition in Table 1, achieves the expected bandgap and avoids bandgap resonance, thus validating the proposed design method. However, bandgap resonances emerge once either of these two conditions is not satisfied. This validates the sufficient condition for avoiding bandgap resonance.