<figure><img src="image_1.png" /></figure>

HAL
open science

# Optimization of spider web-inspired phononic crystals to achieve tailored dispersion for diverse objectives

Vinícius F Dal Poggetto, Federico Bosia, M. Miniaci, Nicola M Pugno

► To cite this version:

Vinícius F Dal Poggetto, Federico Bosia, M. Miniaci, Nicola M Pugno. Optimization of spider web-inspired phononic crystals to achieve tailored dispersion for diverse objectives. Materials & Design, 2021, 209, pp.109980. 10.1016/j.matdes.2021.109980. hal-03440398

HAL Id: hal-03440398

https://hal.science/hal-03440398

Submitted on 22 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Optimization of spider web-inspired phononic crystals to achieve
tailored dispersion for diverse objectives

Vinícius F. Dal Poggetto<fnref n="1" />,*, Federico Bosia<fnref n="2" />, Marco Miniaci<fnref n="3" />, Nicola M. Pugno<fnref n="1" />,<fnref n="4" />,*

Abstract

Spider orb webs are versatile multifunctional structures with optimized mechanical properties for prey capture, but also for transmitting vibrations. The versatility of such a system mainly derives from its variable geometry, which can be effectively used to design phononic crystals, thus inhibiting wave propagation in wide frequency ranges. In this work, the design of spider web-inspired single-phase phononic crystals through selective variation of thread radii and the addition of point masses is proposed, determined through the use of optimization techniques. The obtained results show that spider web geometry displays a rich vibration spectrum, which by varying its geometric characteristics and adding localized masses can be tailored to manipulate wave modes, and the resulting two-dimensional phononic crystals present wide complete band gaps generated by Bragg scattering and local resonances.

**Keywords:** Bioinspired phononic crystals and metamaterials; Spider web-inspired geometry; Optimization; Wave propagation.

**1. Introduction**

Spider webs are the product of evolutionary adaptation, being able to deliver a compromise between demands such as absorbing the impacts of prey while also efficiently transmitting information about the nature and position of vibration sources [1, 2]. Different spider web geometries include the vertical orb web, funnel web, sheet web, and tangle web [3]. Here, a simple plane orb web model is considered, constituted by radial threads connecting the center of the web (hub) to the outer region (frame), and viscid threads (spiral) connecting adjacent radial threads [4, 5]. Radial silk threads are known for their high toughness and tensile strength, with larger Young's modulus and diameters, while viscid threads are known for their large ultimate strain, but have smaller Young's modulus and diameters [6, 7]. Thus, the versatility observed in spider web geometries and the resulting remarkable ability to manipulate waves suggests it can be efficiently used

<footnote>*Corresponding authors</footnote>

<footnote>Email addresses: v.fonsecadalpoggetto@unitn.it (Vinícius F. Dal Poggetto), nicola.pugno@unitn.it (Nicola M. Pugno)</footnote>

<footnote n="1">¹Laboratory of Bio-inspired, Bionic, Nano, Meta Materials & Mechanics, Department of Civil, Environmental and Mechanical Engineering, University of Trento, 38123 Trento, Italy</footnote>

<footnote n="2">²DISAT, Politecnico di Torino, 10129 Torino, Italy</footnote>

<footnote n="3">³CNRS, Centrale Lille, ISEN, Univ. Lille, Univ. Valenciennes, UMR 8520 - IEMN, F-59000 Lille, France</footnote>

<footnote n="4">⁴School of Engineering and Materials Science, Queen Mary University of London, Mile End Road, London E1 4NS, United Kingdom</footnote>
to conceive new artificial materials.

Novel systems designed to exhibit wide attenuation regions have recently been inspired by biological structures [8–13]. In particular, Miniaci et al. [4] have proposed to periodically arrange spider web-inspired structures to achieve phononic crystals (PCs) with efficient wave attenuation performance, through the nucleation of Bragg and local resonance based phononic band gaps (BGs), i.e., frequency regions where wave propagation is inhibited [14–24]. Furthermore, spatially sparse (i.e., lightweight) geometries can be particularly attractive to design systems capable of yielding large attenuation regions [25, 26], which further motivates the investigation of spider web-inspired structures. In particular, the growing use of additive manufacturing and its application in various fields [27] enables the design of structures with complex geometries at a low cost [28], which further encourages the development of single-phase PCs with monolithic structures that can be readily fabricated to yield interesting vibration-control properties [29].

Several types of strategies have been proposed to achieve BGs in two-dimensional structures, such as introducing (i) impedance mismatches through the addition of localized masses [30] or periodic cavities [31], (ii) variations of cross-sections in specific regions [32, 33], or (iii) exploiting the effect of local resonance [34–36]. Spider web structures provide the possibility of adopting all of these strategies due to their versatility, which allows the combination of the impedance modulation realized by graded thread diameters and applied localized masses throughout the web structure, which can be properly determined through the use of optimization techniques. Optimization techniques have been widely used to design PCs to achieve wide BGs, including genetic algorithms [37–39], topology [40–42] and evolutionary optimization [43]. Additionally, when dealing with a given topology, parametric optimization [44] is useful to properly select design variables. This is particularly relevant to spider web-inspired structures.

Although spider web-inspired PC structures have been proposed before [4, 45], their conception was limited to basic configurations, while their full potential for wave manipulation has remained largely unexplored. In this paper, the evaluation of the effects of (i) thread diameter to create impedance mismatches and (ii) the inclusion of point masses to tune the wave filtering characteristics of two-dimensional PCs designed using simple models of spider orb webs and to optimize these parameters to achieve distinctive objectives such as (i) maximum normalized BG widths; (ii) BGs opened at a fundamental frequency and its harmonics; (iii) BGs considering distinct out-of-plane and in-plane polarization modes; (iv) creation of isolated Dirac cones is proposed. This works presents, for the first time, a spider web-inspired single-phase PC which displays great versatility. Its structure can easily be tailored to achieve optimized performance for various different objectives. The presented approach is entirely general and may also be applied to optimize other bio-inspired systems which present similar design freedom. The paper is organized as follows: Section 2 presents the models and methods, Section 3 illustrates the obtained results, and Section 4 presents concluding remarks.

## 2. Models and methods

### 2.1. *Dynamic models*

The development of a two-dimensional PC using a spider orb web-inspired geometry is initially proposed. The relation between the initial spider web-inspired structure, its corresponding periodic hexagonal lattice, and the resulting unit cell are depicted in Figure 1.
<figure><img src="image_3.png" /><figcaption>Figure 1: Spider web-inspired periodic frame: (a) single hexagonal frame structure (radial threads indicated as red lines and viscid threads indicated as blue lines) with distinct composing layers, (b) periodic structure obtained from the repetition of the spider web-inspired frame (colored) in a hexagonal lattice, (c) representative unit cell used to obtain the periodic structure (some frame elements are removed at the edges of the unit cell so they are not duplicated in the corresponding periodic structure).</figcaption></figure>

The initial spider orb web hexagonal geometry with characteristic length L is depicted in Figure 1a, with red lines indicating radial threads and blue lines indicating viscid threads. This structure can be interpreted as a composition of $n_s$ layers, where each layer is constituted by 12 radial and 12 viscid threads, as illustrated for $n_s = 3$. This structure can be repeated to form a hexagonal periodic lattice, as shown in Figure 1b, based on the corresponding unit cell depicted in Figure 1c, where for the i-th layer, the radius of a radial thread is denoted by $r_r^{(i)}$, and of a viscid thread by $r_v^{(i)}$. Thus, for $n_s$ layers, there are $2n_s$ possible different thread diameters. Due to the condition of periodicity, elements connecting corresponding nodes are removed so they are not duplicated [34].

Each structural component is modeled as a 2-node frame element, which superposes the effects of rod, shaft, and beam elements. Considering a frame element in the x-direction, transverse directions are labeled as y and z. The DOFs for the i-th node are represented by the nodal displacement vector $q_i$, and the forces and moments in the same node are represented using vector $f_i$, respectively given by

$$ q_i = \{ u_{xi} \ u_{yi} \ u_{zi} \ \phi_{xi} \ \phi_{yi} \ \phi_{zi} \}^T, \qquad (1a) $$

$$ f_i = \{ F_{xi} \ F_{yi} \ F_{zi} \ M_{xi} \ M_{yi} \ M_{zi} \}^T, \qquad (1b) $$

where $u_{xi}$, $u_{yi}$, and $u_{zi}$ represent displacements, $\phi_{xi}$, $\phi_{yi}$, and $\phi_{zi}$ rotations, $F_{xi}$, $F_{yi}$, and $F_{zi}$ forces, and $M_{xi}$, $M_{yi}$, and $M_{zi}$ moments, for the x, y, and z directions, respectively. Quantities $u_x$ and $F_x$ are related using the rod model, $\phi_x$ and $M_x$ using the shaft model, and quantities $u_y$, $\phi_z$, $F_y$, and $M_z$ are related using the beam model, as are also $u_z$, $\phi_y$, $F_z$, and $M_y$. The relations between displacements (rotations) and forces (moments) at nodes are described in the literature for the case of linear elastic behavior using stiffness and mass matrices [46].

## 2.2. Band structure computation

Band diagrams are used to graphically represent dispersion relations, i.e., relations between the propagating waves circular frequency $\omega$ and corresponding wave vectors $\mathbf{k}$. Even though dispersion relations are computed considering infinite periodic structures, they are able to provide important information regarding wave propagation in finite structures [12, 47, 48]. In the case of systems
with two-dimensional periodicity, the wave vector can be restricted to the plane containing the structure and be written in terms of its Cartesian components ($k_x$ and $k_y$) using

$$ \mathbf{k} = k_x \hat{\mathbf{i}} + k_y \hat{\mathbf{j}}. \qquad (2) $$

Due to the system periodicity, the usual process for computing the bands diagram consists in analyzing a characteristic unit cell, representative of the whole structure. In this case, it is possible to restrict the wave vector to the first Brillouin zone (FBZ) [14]. Also, due to symmetry, only a fraction of the FBZ needs to be considered, named the irreducible Brillouin zone (IBZ) [49]. For a two-dimensional hexagonal lattice, the high symmetry points with respective coordinates in the k-space are given by $\Gamma$ (0,0), M ($\pi/L$, $-\pi/L\sqrt{3}$), and K ($4\pi/3L$, 0). These concepts are shown in Figure 2 for the hexagonal geometry.

<figure><img src="image_4.png" /><figcaption>Figure 2: (a) Two-dimensional hexagonal periodic medium with a highlighted unit cell and (b) FBZ with its highlighted IBZ and high symmetry points $\Gamma$ (0,0), M ($\pi/L$, $-\pi/L\sqrt{3}$), and K ($4\pi/3L$, 0).</figcaption></figure>

For the computation of the band diagrams, the finite element (FE) method is used with the application of periodicity to enforce the Bloch-Floquet conditions [50] on the contour of the unit cell. The propagating wave frequencies are computed for fixed directions ($\omega = \omega(\mathbf{k})$), while the evanescent behavior of waves is obtained by computing the complex wave vectors for each frequency of interest ($\mathbf{k} = \mathbf{k}(\omega)$). Further details on the computational implementation of this FE-based band diagram computation technique are given in the supplementary material.

## 2.3. Parameter optimization

### 2.3.1. Parameters and constraints

The number of possible different thread diameters in the proposed unit cell and the addition of localized masses at the junction between radial and viscid threads provide considerable design freedom. A proper selection of design parameters can be obtained using optimization techniques to achieve designs that fulfill specific objectives, such as wide BGs with the lowest possible central frequency [51]. Possible variations of thread diameters in a unit cell and the arrangement of masses are depicted schematically in Figure 3.

Fixed values of $L$ and $n_s$ are used (which are usually limited by physical construction constraints) and keep the values of thread diameters and added masses as possible optimization parameters. Henceforth, the list of the continuous optimization parameters that describe thread radii will be expressed in a vector $\mathbf{r}$ given by

$$ \mathbf{r} = \left\{ r_v^{(1)} \ r_r^{(1)} \ r_v^{(2)} \ r_r^{(2)} \ \dots \ r_v^{(n_s)} \ r_r^{(n_s)} \right\}, \qquad (3) $$
<figure><img src="image_5.png" /><figcaption>Figure 3: Design options used for parameter optimization: (a) variation of radial and viscid thread diameters and (b) introduction of localized masses at the junctions of threads.</figcaption></figure>

where $r_v^{(i)}$ and $r_r^{(i)}$ are the values of viscid and radial thread radii for the i-th layer, respectively, for $i = 1, 2, \dots, n_s$, as illustrated in Figure 1c. For a unit cell with $n_s$ layers, this yields a total of $2n_s$ design parameters. The values of thread radii can be fixed between lower and upper bounds given by $r_{min}$ and $r_{max}$, respectively. This can be summarized as a restriction on $R_r$, stated as

$$R_r: r_v^{(i)}, r_r^{(i)} \in [r_{min}, r_{max}], \forall i = 1, 2, \dots, n_s. \quad (4)$$

Similarly, a mass vector $\mathbf{m}$ can be defined to describe localized masses which can be added to the unit cell. The mass distribution begins from the inner part of the first layer, moving to the center of the orb web, using the continuous variables vector described by

$$\mathbf{m} = \left\{ m^{(1)} \ m^{(2)} \ \dots \ m_v^{(n_s)} \right\}, \quad (5)$$

where $m^{(i)}$ describe the masses added to the inner region of the i-th layer, for $i = 1, 2, \dots, n_s$, yielding a total of $n_s$ design parameters. Instead of setting the values of masses between fixed bounds, the total added mass is limited to $m_{add}$, which can be stated as the restriction $R_m$, written as

$$R_m: \left( 12 \sum_{i=1}^{n_s-1} m^{(i)} \right) + m^{(n_s)} \le m_{add}. \quad (6)$$

### 2.3.2. Objective functions

#### Normalized band gap width

In the following, various objective functions for the dispersion spectrum of the considered structure are presented, showing how each of these alternative configurations can be achieved, and thus demonstrating the versatility of the spider web-inspired geometry.

With the purpose of achieving BGs with maximum widths and minimum central frequencies (typical for structural vibration mitigation applications), the optimization objective is chosen as the maximization of the normalized width of the difference between the n-th and (n + 1)-th bands [52]. Considering a given band diagram, the difference between consecutive propagating frequencies $\omega_n$ and $\omega_{n+1}$ for the IBZ position given by the wave vector $\mathbf{k}$ can be calculated as $\omega_{n+1}(\mathbf{k}) - \omega_n(\mathbf{k})$. Furthermore, if $\mathbf{k}$ is considered over all possible positions along the contour of the IBZ, one can define the difference between the n-th and (n + 1)-th propagating frequencies over this contour, given by $\min_{\mathbf{k}} \omega_{n+1}(\mathbf{k}) - \max_{\mathbf{k}} \omega_n(\mathbf{k})$, which means that, for a positive value of this difference, a BG exists between the n-th and (n + 1)-th bands, which can be further widened. Thus, this
metric can be combined with the corresponding central (mean) frequency to describe the objective function $\varphi_n$, as

$$ \varphi_n^{(s)} = \frac{\min_k \omega_{n+1}(k) - \max_k \omega_n(k)}{(\min_k \omega_{n+1}(k) + \max_k \omega_n(k))/2}, \quad (7) $$

where the superscript s refers to the optimization of a single BG.

The optimization objective described by Eq. (7) is especially useful for the optimization process, since it defines a continuous metric between given bands. Although wave modes may cross and change their ordering in the band diagram, $\varphi_n^{(s)}$ consistently provides the difference between the bands chosen for the optimization, and in the case where a BG is closed, it becomes negative.

## Attenuation of harmonics

Another interesting possibility concerns applications in structural health monitoring problems, where the presence of a defect in a structure may be determined by the detection of its nonlinear response, i.e. the presence of harmonics [53, 54]. Thus, it is important to control both the opening and closing of BGs at a specific frequency (fundamental) and its harmonics, i.e., attenuate a fundamental frequency while preserving its harmonics for detection. With this in mind, a metric can be defined using the minimum attenuation computed at a given fundamental frequency ($\omega_0$) and its harmonics ($n\omega_0$, for a positive integer $n$), also applying weighting factors that account for the greater importance of lower harmonics. This metric can be written using the imaginary part of the wave vector, $\mathrm{Im}(k(\omega))$, which accounts for the wave evanescent behavior for a given circular frequency $\omega$ [55], as

$$ \varphi_{\omega_0}^{(\mathrm{fh})} = \left(1 - e^{-\min|\mathrm{Im}(k(\omega_0))|}\right) + \sum_{n=2}^{N_h} \frac{1}{n} e^{-\min|\mathrm{Im}(k(n\omega_0))|}, \quad (8) $$

where the superscript fh refers to a first harmonic approach, and $\min|\mathrm{Im}(k(\omega))|$ is the minimum of the imaginary part of all computed wave vectors, which indicates if all waves are attenuated at a given frequency [56], considering the first $N_h$ harmonics. The factor $e^{-\min|\mathrm{Im}(k(\omega))|}$ becomes 1 for $\min|\mathrm{Im}(k(\omega))| = 0$ and 0 for a sufficiently large $\min|\mathrm{Im}(k(\omega))|$, providing a continuous metric for the attenuation considering the chosen harmonics. Thus, this metric yields positive values that increase for a BG opened at the fundamental frequency (1) and for at least one propagating frequency for each of the higher harmonics (1/n). The weighting factor 1/n is somewhat arbitrary and can be generally replaced by any integer-based positive function that reflects the relative importance of each harmonic.

The simplicity in the presented objective function suggests it can be further modified to design structures capable of opening BGs at each harmonic, which can be written as

$$ \varphi_{\omega_0}^{(\mathrm{ah})} = \sum_{n=1}^{N_h} \frac{1}{n} \left( 1 - e^{-\min|\mathrm{Im}(k(n\omega_0))|} \right), \quad (9) $$

where the superscript ah now refers to an all harmonics approach. It is worth noticing that in this case, the opening of multiple BGs is an objective rather diverse from those previously presented, since objective function $\varphi_n^{(s)}$ aims to open and widen a single BG and objective function $\varphi_{\omega_0}^{(fh)}$ opens a single BG and ensures higher harmonics do not lie inside BGs.

Both proposed objective functions ($\varphi_{\omega_0}^{(fh)}$ and $\varphi_{\omega_0}^{(ah)}$) are generic and can be used for optimization procedures applied to structures with frequency-dependent damping, such as viscoelastic structures [10, 57].
Uncoupled mode polarization

The proposed PC is also used to obtain BGs relative to distinct polarization modes, useful for applications in wave controlling and confinement [58, 59]. To this end, a metric is proposed in the form

$$ \varphi^{(\mathrm{p})} = \sum_{n=1}^{N_b} e^{-\frac{1}{N_p} (N_{\mathrm{out}}(\omega_n) N_{\mathrm{in}}(\omega_n))}, \quad (10) $$

where the superscript p refers to a polarization approach, $N_{\mathrm{out}}(\omega_n)$ and $N_{\mathrm{in}}(\omega_n)$ represent the number of out-of-plane and in-plane wave modes contained in the n-th frequency interval centered at $\omega_n$ with width $\Delta\omega$ (i.e., $\omega \in [\omega_n - \Delta\omega/2, \omega_n + \Delta\omega/2]$), respectively. This metric is built in such a way that consecutive intervals do not overlap, $N_p = \sum_{n=1}^{N_b} N_{\mathrm{out}}(\omega_n) + N_{\mathrm{in}}(\omega_n)$ is the total number of computed wave modes, and $N_b = \omega_{\mathrm{max}}/\Delta\omega$ is the number of frequency intervals (bins) for a given maximum frequency $\omega_{\mathrm{max}}$ and width $\Delta\omega$. The sets of corresponding pairs $\omega_n \times N_{\mathrm{out}}(\omega_n)$ and $\omega_n \times N_{\mathrm{in}}(\omega_n)$ have the same meaning as states histograms, where each n-th component represents a frequency bin for propagating frequencies contained in that frequency region. Whenever $N_{\mathrm{out}}(\omega_n)$ or $N_{\mathrm{in}}(\omega_n)$ are zero, the metric yields a 1 for the corresponding bin, decaying to 0 otherwise, while the term $1/N_p$ provides a smoother transition between these limits. The polarization of modes can be computed by comparing out-of-plane and in-plane displacements for each propagating wave mode and categorizing the modes according to their relative magnitudes.

Dirac cones

For the last optimization objective, the use of the variable spider web-inspired PC structure to create isolated Dirac cones [60–62] is proposed. These features may result in elastic systems with localized wave modes that potentially yield topological mode protection and lossless energy transport [63, 64]. Thus, a metric proportional to the differences regarding the propagating frequencies corresponding to the high-symmetry points $\Gamma$, $M$, and $K$ (denoted respectively as $\omega(\Gamma)$, $\omega(M)$, and $\omega(K)$) relative to the consecutive n-th and $(n+1)$-th bands is proposed in the form

$$ \varphi_n^{(\mathrm{D})} = g_n^{(\mathrm{M}\Gamma)} g_n^{(\mathrm{KM})} g_{n,n+1}^{(\Delta\omega)} g_{n+1}^{(\mathrm{MK})} g_{n+1}^{(\Gamma\mathrm{M})}, \quad (11) $$

where the superscript D refers to a Dirac cone approach, and the functions that compose this metric are given by $g_n^{(\mathrm{M}\Gamma)} = g(\omega_n(\mathrm{M}), \omega_n(\Gamma))$, $g_n^{(\mathrm{KM})} = g(\omega_n(\mathrm{K}), \omega_n(\mathrm{M}))$, $g_{n+1}^{(\mathrm{MK})} = g(\omega_{n+1}(\mathrm{M}), \omega_{n+1}(\mathrm{K}))$, $g_{n+1}^{(\Gamma\mathrm{M})} = g(\omega_{n+1}(\Gamma), \omega_{n+1}(\mathrm{M}))$, and $g_{n,n+1}^{(\Delta\omega)} = g(\Delta\omega, \omega_{n+1}(\mathrm{K}) - \omega_n(\mathrm{K}))$. The function $g(a, b) = 1 - e^{-(a/b)^2}$ becomes 0 for $a \ll b$ and 1 for $a \gg b$, thus performing a simple comparison between its arguments. Terms $g_n^{(\mathrm{M}\Gamma)}$, $g_n^{(\mathrm{KM})}$, $g_{n+1}^{(\mathrm{MK})}$, and $g_{n+1}^{(\Gamma\mathrm{M})}$ account for the separation of propagating frequencies of high-symmetry points, while the term $g_{n,n+1}^{(\Delta\omega)}$ closes the distance between consecutive propagating frequencies at the K-point using a threshold value $\Delta\omega$.

2.3.3. Optimization problems

Optimization problems may be formulated by separately optimizing the radii ($2n_s$ variables), the masses ($n_s$ variables), or by simultaneously optimizing both radii and masses ($3n_s$ variables). The proposed objective functions need to be maximized, which may be subject to $\mathbf{R}_r$ and $\mathbf{R}_m$, either simultaneously, yielding the optimization problem

$$ \underset{\mathbf{r}, \mathbf{m}}{\text{maximize}} \quad \varphi, \text{ subject to } \mathbf{R}_r, \mathbf{R}_m, \quad (12) $$
or sequentially, yielding

$$ \underset{\mathbf{r}}{\text{maximize}} \varphi, \text{ subject to } \mathbf{R}_{\mathbf{r}}, \qquad (13a) $$

$$ \underset{\mathbf{m}}{\text{maximize}} \varphi, \text{ subject to } \mathbf{R}_{\mathbf{m}}, \qquad (13b) $$

where the optimization variables **r** and **m** are given, respectively by Eqs. (3) and (5), restrictions $\mathbf{R}_{\mathbf{r}}$ and $\mathbf{R}_{\mathbf{m}}$ are given, respectively by Eqs. (4) and (6); the generic objective function, $\varphi$, may refer to $\varphi_n^{(s)}$ (Eq. (7)), $\varphi_{\omega_0}^{(fh)}$ (Eq. (8)), $\varphi_{\omega_0}^{(ah)}$ (Eq. (9)), $\varphi^{(p)}$ (Eq. (10)), or $\varphi_n^{(D)}$ (Eq. (11)), depending on which optimization objective is considered.

The resulting constrained nonlinear optimization problems can be solved using a sequential quadratic programming (SQP) algorithm, which yields excellent efficiency and percentage of successful solutions when compared to several other methods [65]. The SQP algorithm uses a quadratic approximation of a Lagrangian function to formulate a quadratic programming subproblem [66]. Candidate solutions are iterated using approximations for the Hessian matrix of the Lagrangian function, conveniently updated through quasi-Newton methods [67]. This type of method is also particularly useful for constrained parameters, and is currently available in commercial software such as Matlab [68]. The interested reader is referred to the general overviews presented in [69, 70].

# 3. Results

## 3.1. Initial considerations

To compute the band diagrams of the structure shown in Figure 1c, the material properties of Digital ABS Plus (a commonly used material for 3D printing) are considered for both the radial and viscid threads (which simplifies the fabrication process): Young's modulus $E = 2800$ MPa, Poisson's ratio $\nu = 0.35$, and mass density $\rho = 1175$ kg/m$^3$. Before applying the optimization process for distinct objectives, the band diagrams yielded by a PC with all threads having the same diameters are first investigated.

### 3.1.1. Initial band diagrams

The band diagram (real and imaginary parts of the wave vector) of a hexagonal lattice with length $L = 50$ mm, number of layers $n_s = 3$, and elements with radius $r = 1.5$ mm, using the material properties of Digital ABS Plus are initially computed, which yields a unit cell with a mass of $4.84 \times 10^{-3}$ kg. The resulting band diagram, unit cell, and corresponding wave mode displacements are shown in Figure 4 with the real part of wave vectors marked in blue; since the imaginary part of wave vectors serves as an indicative of the type of mechanism associated with BG formation, these are also shown in red, indicated as $\text{Im}(k_{\Gamma M}a/2\pi)$ and $\text{Im}(k_{\Gamma K}a/2\pi)$ for the $\Gamma M$ and $\Gamma K$ directions at the left and right portions of the band diagram, respectively. Also, the normalized frequency $\bar{\omega} = \omega a/2\pi c_L$ is used, where $a = L\sqrt{3}$ is the hexagonal lattice characteristic length and $c_L = \sqrt{E/\rho}$ is the speed of longitudinal waves in the material.

The band diagram in Figure 4a shows no BGs. Wave mode shapes are marked using letters in Figure 4b, where a rigid-body-like translational behavior can be noticed in wave modes V (out-of-plane displacements) and H1/H2 (in-plane displacements); wave mode pairs B1/B2, BB1/BB2, and BB3/BB4 represent first-, second-, and third-order-like bending modes, respectively; and wave mode T is torsional. Although, as observed, the lack of impedance mismatches indicates that no BGs are formed in this initial structure, these initial results allow further comprehension on the
<figure><img src="image_6.png" /><figcaption>Figure 4: Initial results for the hexagonal lattice structure: (a) band diagram with the real and imaginary parts of the wave vector marked in blue and red, respectively, showing no BGs; (b) mode displacements profiles (colored according to the magnitude of displacement and normalized according to the maximum of displacement in each case): V is an out-of-plane displacement mode, B1 and B2 represent bending modes, BB1 and BB2 represent second-order-like bending modes, BB3 and BB4 represent third-order-like bending modes, H1 and H2 are in-plane displacement modes, and T is a torsional mode.</figcaption></figure>

influence that the parameter optimization has on altering the band diagram and wave modes, as presented in the subsequent sections.

## 3.2. Optimization of normalized band gap width

### 3.2.1. Optimization of thread radii

Here, the results of the optimization process considering the diameters of threads (Figure 3a) with the objective of opening and widening various BGs (Eq. (7)) through the selection of proper values of $n$ using the optimization problem stated in Eq. (13a) are presented. The same initial structure as in the previous section is chosen, with a lattice length of $L = 50$ mm and $n_s = 3$ layers, which yields 6 design variables.

In view of future experiments, the bounds for the optimization variable $r$ are chosen compliant to typical manufacturing restrictions such as the resolution of the additive process. Here, the bounds $r_{min} = 0.5$ mm and $r_{max} = 3.0$ mm are set. Also, due to the large number of converging threads at the unit cell center (see Figure 1a), an increase in the thread diameters in this vicinity might be unfeasible for the manufacturing process of the PC. Thus, to facilitate the connection between these elements, the innermost thread diameters are removed from the list of parameters and fix it at the lower limit (i.e., $r_r^{(n_s)} = r_{min}$).
The first bands that could be successfully separated using the optimization process were the 6th and 7th bands. The optimized parameter vector that widens the corresponding normalized BG width is obtained as

$$ \mathbf{r}_{\text{opt}}^{(6,7)} = \{ 0.500 \ 0.594 \ 3.000 \ 3.000 \ 3.000 \ 0.500 \} \text{mm}, \quad (14) $$

which yields the corresponding band diagram, unit cell, and wave modes shown in Figure 5. Henceforth, the obtained parameters representing the unit cell with the thickness of elements proportional to the optimized diameters are depicted, also using colors to indicate the added masses at nodes, ranging from white (no added mass) to black (all added mass).

<figure><img src="image_7.png" /><figcaption>Figure 5: Results for optimized structure for maximum normalized BG width between 6th and 7th bands: (a) band diagram with BGs between 2296 Hz – 2788 Hz, 2998 Hz – 3025 Hz, 3464 Hz – 3742 Hz, 4433 Hz – 4640 Hz, and 6633 Hz – 7000 Hz and its corresponding unit cell obtained using optimization objective $\varphi_6^{(s)}$ (Eq. (7)). (b) Wave modes are separated between low-frequency (B1, B2, V, H1, H2, and T, dominated by rigid-body-like inclusion behavior) and mid- to high-frequency wave modes (L1, L2, L3, and L4) with localized motion at the outer threads.</figcaption></figure>

The band diagram computed using optimized geometric parameters for widening the normalized BG width between the 6th and 7th bands (Figure 5a) has the first BG located between 2296 Hz – 2788 Hz (19% normalized width) which separates wave modes (presented in Figure 5b) associated with the motion of a rigid-body-like inclusion (wave modes B1, B2, V, H1, H2, and T) from wave modes with higher frequencies and localized motion (L1, L2, L3, and L4). The second-order bending modes (BB1 and BB2, Figure 4b) are not present.

As side effects, additional BGs are opened at 2998 Hz – 3205 Hz, 3464 Hz – 3742 Hz, 4433 Hz – 4640 Hz, and 6633 Hz – 7000 Hz frequency ranges, which are associated with the localized wave
modes. Thus, although above 2788 Hz no wide BGs are noticed, one can expect limited energy propagation, given the low group velocity presented by these wave modes.

The imaginary parts of wave vectors have a fairly symmetrical distribution around the central frequency of each BG, indicating they are formed through Bragg scattering, although localized modes are visible in the diagram (e.g., a locally resonant BG appears around 7000 Hz). This indicates that the lattice length can be scaled to match the desired wavelengths associated with a given vibration attenuation application. Also, since Bragg scattering is associated with the lattice length, the opening of BGs at frequencies corresponding to multiple wavelengths is usually expected. The occurrence of several Bragg scattering BGs for wavelengths in the same order of magnitude indicates another advantage of the proposed structure.

Regarding the optimal radii distribution $\mathbf{r}_{\text{opt}}^{(6,7)}$, it is interesting to notice that the first parameter is equal to the lower bound (0.5 mm), and the second parameter is close to this lower bound, while all other optimized parameters (viscid and radial threads starting from the second layer, with exception of the innermost radial threads) are equal to the upper bound (3.0 mm). This indicates that the optimal configuration to achieve the maximum normalized BG width consists in maximizing the stiffness and mass of the inner regions, while minimizing the stiffness of the outer regions, which leads to a spring-mass-like system with the outer region acting as soft springs and the inner region acting as a rigid and heavy core. This result agrees with what has been previously presented in literature [25, 26] and validates the proposed optimization scheme.

With this configuration, the unit cell has a total mass of 9.65 g, which represents a 99% weight increase in the mass of the unit cell when compared to the initial configuration. This further reinforces the spring-mass behavior for the low-frequency wave modes of the system, i.e., low stiffness with large mass values.

It is also interesting that the second parameter is not exactly equal to the lower bound. To understand why, it is necessary to analyze the band diagram considering the radii vector $\mathbf{r} = \{ 0.5 \ 0.5 \ 3.0 \ 3.0 \ 3.0 \ 0.5 \} \text{ mm}$ (not shown here, for the sake of brevity). For these parameters, the first BG has a slightly narrower normalized width (15%), which suggests that a smoother radii transition between the outer and inner regions of the PC is beneficial for BG formation, as opposed to excessively large differences in thread radii.

The described optimization procedure can also be used to manipulate other characteristics of the band diagram, for example to obtain a BG associated with the separation of the 8th and 9th bands. In this case, the optimization process yields the optimized radii vector given by

$$ \mathbf{r}_{\text{opt}}^{(8,9)} = \{ 0.680 \ 0.500 \ 3.000 \ 3.000 \ 3.000 \ 0.500 \} \text{ mm}, \quad (15) $$

showing a smooth variation between the optimized radii (as opposed to setting the first 2 variables at the lower bound), reinforcing the previous observation that a smoother radii transition is beneficial for BG formation, which would not be possible if an optimization method using a discrete set of variables were chosen. The corresponding computed band diagram, unit cell, and some wave modes are depicted in Figure 6.

In this case, the band diagram (Figure 6a) shows that the first BG, formed between the 8th and 9th bands, appears between 1980 Hz – 2570 Hz, which indicates a better result (26% normalized width) than the one obtained using the previous configuration, and smaller BGs appear at higher frequencies (4116 Hz – 4202 Hz and 4205 Hz – 4768 Hz). Also, in this configuration the unit cell has a mass of 9.68 g (a 100% increase when compared to the initial design). Thus, a wider
<figure><img src="image_8.png" /><figcaption>Figure 6: Results for optimized structure for maximum normalized BG width between 8th and 9th bands: (a) band diagram with BGs between 1980 Hz – 2570 Hz, 4116 Hz – 4202 Hz, 4205 Hz – 4768 Hz, and 7180 Hz – 7327 Hz and its corresponding unit cell obtained using optimization objective $\phi_8^{(s)}$ (Eq. (7)). (b) Localized wave modes (L1, L2, L3, and L4) control BG limits and a vertical localized wave mode L5 with zero group velocity can be noticed.</figcaption></figure>

normalized BG width is obtained using practically the same amount of added mass, indicating a more efficient use of mass addition.

The opening mechanisms of the computed BGs can be explained analyzing the wave modes depicted in Figure 6b. The slight increase in the outermost viscid threads radius yields several localized wave modes. In the case of L1, a negligible group velocity is noticed, and its lower frequency controls the first BG opening, which closes at the localized wave mode depicted as L2. A similar mechanism occurs between L3 and L4, where a wave mode with zero group velocity (L5) can be noticed, which is explained due to the localized vibration of the innermost threads. In fact, the presence of mode L5 is barely noticeable when analyzing the imaginary part of the wavenumbers, which further reinforces the conclusion that in practice, a single BG is achieved between L3 and L4. It should also be noticed that the mechanisms associated with the opening of BGs are based on localized wave modes not obviously foreseen using simple design guidelines, which further supports the use of an optimization process for a more efficient use of design variables.

### 3.2.2. Optimization of added masses

The optimization of the distribution of added masses (Figure 3b), considering the same objective function as the previous section (Eq. (7)) is now investigated. Initially, only mass addition is discussed (without variations in thread radii), as stated by Eq. (13b). For an unbiased comparison, the initial structure ($r = 1.5$ mm) is considered, and the total added mass is limited to the difference between the initial design and the maximum mass addition obtained in the previous
cases, i.e., $m_{add} = 4.84$ g.

The obtained results show that when such restrictions are imposed, no complete (i.e, for all types of waves) or full (i.e., for all wave vector values) BGs can be opened considering up to the 10th band. These results also suggest that the simple addition of masses, i.e., without the manipulation of local stiffness, is not effective in creating complete BGs.

To further investigate the effects of mass addition, the necessary allocated mass is computed increasing $m_{add}$ until a first complete and full BG is obtained using two distinct configurations: (i) at the PC center and (ii) with equally distributed masses located at the radial-viscid junctions closest to the PC center. Results are shown in Figure 7.

<figure><img src="image_9.png" /><figcaption>Figure 7: Band diagrams of structures with added masses: (a) concentrated at the PC center (BG around 630 Hz) and (b) distributed around the first junctions around PC center (BG around 4410 Hz).</figcaption></figure>

The band diagram shown in Figure 7a was obtained for $m = \{ 0.000 \ 0.000 \ 305 \}$ g and indicates that the addition of mass at the PC center is able to separate the vertical and horizontal rigid-body-like wave modes from the first bending modes. However, an excessively large mass addition (nearly $63 \times m_{add}$) is necessary to nucleate this BG. Meanwhile, in Figure 7b, the band diagram obtained using masses allocated around the PC center (radial-viscid junctions closest to the PC center) is shown, obtained for $m = \{ 0.000 \ 0.687 \ 0.000 \}$ g, indicating that localized wave modes can be more efficiently used to nucleate BGs at higher frequencies and require less mass addition (in this case, $1.7 \times m_{add}$).

To further investigate the effects of mass addition, the allocation of masses using a previously optimized structure presented in Section 3.2.1 is now considered. The same restrictions on the addition of masses are imposed, i.e., $m_{add} = 4.84$ g, and consider the optimized radii given by $\mathbf{r}_{opt}^{(6, 7)}$.

Since this optimization process is independent of the previous one, different BGs may be widened. The BGs between two different sets of bands are then chosen: (i) the 6th and 7th and (ii) the 7th and 8th bands. The obtained optimized mass vectors are respectively given by

$$ \mathbf{m}_{opt}^{(6, 7)} = \{ 0.000 \ 0.403 \ 0.000 \} \text{ g}, \qquad (16) $$
which represents masses allocated at the junction between radial and viscid threads closest to the PC center, and

$$ \mathbf{m}_{\text{opt}}^{(7,8)} = \{ 0.000 \ 0.000 \ 4.84 \} \text{ g,} \qquad (17) $$

which represents all of the mass allocated at the PC center. The corresponding band diagrams and unit cells are presented in Figure 8.

<figure><img src="image_10.png" /><figcaption>Figure 8: Band diagrams computed (considering a previously optimized threads configuration) using a mass distribution for maximizing the normalized BG width between (a) the 6th and 7th bands (first BG between 1868 Hz – 2759 Hz) and (b) the 7th and 8th bands (first BG between 1705 Hz – 2788 Hz) and their corresponding unit cells obtained using optimization objectives $\phi_n^{(s)}$ (Eq. (7)) for $n = 6$ and $n = 7$, respectively.</figcaption></figure>

Considering the band diagram shown in Figure 8a, obtained using the mass vector given by Eq. (16), BGs are noticed between 1868 Hz – 2759 Hz, 2985 Hz – 3098 Hz, 3353 Hz – 3684 Hz, 4425 Hz – 4532 Hz, 4565 Hz – 4640 Hz, and 6633 Hz – 7000 Hz. The first BG presents a 39% normalized width, which represents a significant improvement over its previous version, i.e., without added masses. Interestingly, even though the optimization constraints are based on an inequality (i.e., not all masses must be allocated), these results indicate that the allocation of total masses is capable of yielding improved results.

Moreover, the band diagram shown in Figure 8b, obtained using the mass vector given by Eq. (17), displays BGs between 1705 Hz – 2788 Hz, 2998 Hz – 3205 Hz, 4429 Hz – 4640 Hz, and 6633 Hz – 7000 Hz, which indicate a first BG with a considerable improvement in its normalized width (48%). However, this result also indicates a shortcoming of the optimization objective, since a better result was obtained considering different bands to be separated in consecutive optimization processes, i.e., 8th and 9th bands for the mass addition and 6th and 7th bands for the thread diameters processes, instead of using the same pair of bands in both cases. The interplay between
optimization objectives relative to different bands to be separated can present intricate solutions, in which case, the procedure shown here is meant to be an initial approach in demonstrating the versatility of the proposed structure. Further investigations are beyond the scope of this paper and left for future work.

### 3.3. Optimization of fundamental frequency and harmonics

The results of the optimization process considering the diameters of threads and distribution of added masses with the objective of attenuating a fundamental frequency and preserving or attenuating its higher harmonics using the problem stated in Eq. (12) and the objective functions given in Eqs. (8) and (9) are now presented. The same initial geometry, material properties, and optimization variables constraints as presented in the last section are considered. Also, the considered harmonics are restricted to the maximum frequency of 7000 Hz. Two distinct fundamental harmonics are considered: (i) 2000 Hz and (ii) 2500 Hz, which were chosen to demonstrate the versatility of the structure, since no harmonics are common to both fundamental frequencies in the analyzed frequency range. The lowest fundamental frequency was chosen to reflect the lowest limit of previously obtained BGs.

Even though the BG opening mechanisms have already been explained in the last section, the optimized parameters results are presented for additional insight on the obtained band diagrams, comparing the results that yield a single BG at the fundamental frequency to those that also attenuate higher harmonics. For case (i), the optimized thread radii and masses to open a BG only at the fundamental frequency are obtained as

$$ \mathbf{r}_{\text{opt, f}}^{(f_0=2000 \text{ Hz})} = \{ 0.565 \ 0.602 \ 0.500 \ 0.684 \ 2.720 \ 0.500 \} \text{mm}, \quad (18a) $$

$$ \mathbf{m}_{\text{opt, f}}^{(f_0=2000 \text{ Hz})} = \{ 0.403 \ 0.000 \ 0.000 \} \text{g}, \quad (18b) $$

while in the case of opening BGs at each harmonic, the optimization yields

$$ \mathbf{r}_{\text{opt, a}}^{(f_0=2000 \text{ Hz})} = \{ 0.900 \ 0.744 \ 2.994 \ 0.501 \ 0.559 \ 0.500 \} \text{mm}, \quad (19a) $$

$$ \mathbf{m}_{\text{opt, a}}^{(f_0=2000 \text{ Hz})} = \{ 0.000 \ 0.402 \ 0.000 \} \text{g}, \quad (19b) $$

which indicate that, for multiple BGs, the thick viscid layer is shifted from the center to the periphery of the unit cell, while allocated masses are shifted in the opposite direction. The corresponding band diagrams and unit cells are shown in Figure 9, where the corresponding fundamental frequencies and their harmonics are marked using orange dashed lines.

Results presented in Figure 9 indicate the successful opening of BGs only for a fundamental frequency (Figure 9a) or also attenuating higher harmonics (Figure 9b). In both cases, a large number of zero group velocity waves are located below the first BG, explained due to the small value of radii in the first viscid layer. It is worth noticing that in Figure 9a, at least one propagating mode (i.e., zero attenuation) is obtained for each of the higher harmonics, thus ensuring detectable wave modes. Also, the analysis of the imaginary part of the obtained wave vectors in Figure 9b indicates decreasing attenuation values for the harmonics as frequency increases.

For case (ii), the optimized thread radii and masses to open a BG only at the fundamental frequency are given by

$$ \mathbf{r}_{\text{opt, f}}^{(f_0=2500 \text{ Hz})} = \{ 0.657 \ 0.544 \ 2.317 \ 1.877 \ 2.109 \ 0.500 \} \text{mm}, \quad (20a) $$

$$ \mathbf{m}_{\text{opt, f}}^{(f_0=2500 \text{ Hz})} = \{ 0.340 \ 0.011 \ 0.542 \} \text{g}, \quad (20b) $$
<figure><img src="image_11.png" /><figcaption>Figure 9: Band diagrams and their corresponding unit cells obtained using the optimization objective that (a) opens a BG only at the fundamental frequency ($\varphi_{\omega_0}^{(\text{fh})}$, Eq. (8)) and (b) attenuates all harmonics ($\varphi_{\omega_0}^{(\text{ah})}$, Eq. (9)) for $f_0 = 2000$ Hz. Harmonics are indicated using dashed lines (- -).</figcaption></figure>

while in the case of opening a BGs at each harmonic, the optimization yields

$$ \mathbf{r}_{\text{opt, a}}^{(f_0=2500 \text{ Hz})} = \{ 0.676 \ 0.636 \ 2.996 \ 1.813 \ 2.698 \ 0.500 \} \text{mm}, \quad (21a) $$

$$ \mathbf{m}_{\text{opt, a}}^{(f_0=2500 \text{ Hz})} = \{ 0.004 \ 0.000 \ 4.330 \} \text{g}, \quad (21b) $$

which indicate somehow similar distributions, with a more noticeable difference in the shifting of added mass between the periphery and the center of the unit cell. The corresponding band diagrams and unit cells are shown in Figure 10.

The band diagrams presented in Figures 10a and 10b show that some BGs are opened at frequencies unrelated to the desired harmonics, which do not impede the proposed objective from being achieved. Also, in Figure 10a, the second harmonic is close to a BG unrelated to the optimization objective, which is shifted up in Figure 10b. This slight change in the dispersion relation is enough to achieve the objective of attenuating all harmonics, which helps to understand why both objectives could be achieved with a small change in the parameters.

### 3.4. Optimization of in-plane and out-of-plane mode polarization

The results of the optimization process considering the diameters of threads and distribution of added masses to obtain BGs with distinct (out-of-plane and in-plane) polarization modes using the problem stated in Eq. (12) and the objective function given in Eq. (10) are now presented. The same previous initial geometry, material properties, and optimization constraints are considered.
<figure><img src="image_12.png" /><figcaption>Figure 10: Band diagrams and their corresponding unit cells obtained using the optimization objective that (a) opens a BG only at the fundamental frequency ($\varphi_{\omega_0}^{(fh)}$, Eq. (8)) and (b) attenuates all harmonics ($\varphi_{\omega_0}^{(ah)}$, Eq. (9)) for $f_0 = 2500$ Hz. Harmonics are indicated using dashed lines (- -).</figcaption></figure>

The optimization process yields the radii and mass vectors given by

$$ \mathbf{r}_{\text{opt}}^{(p)} = \{ 0.500 \ 0.500 \ 0.502 \ 3.000 \ 3.000 \ 0.500 \} \text{ mm}, \qquad (22a) $$

$$ \mathbf{m}_{\text{opt}}^{(p)} = \{ 0.403 \ 0.000 \ 0.000 \} \text{ g}, \qquad (22b) $$

indicating a single layer of thick viscid and radial threads, with masses distributed around its outer nodes. Results for the optimized unit cell, initial and final band diagrams, and density of states (DOS, shown in arbitrary units) histograms are shown in Figure 11, calculated using 20 equally spaced frequency bins between 0 and 7000 Hz.

The band diagram and DOS histogram (Figures 11a and 11b) show that initially, every frequency bin contains wave modes from both types of polarization (out-of-plane and in-plane). The final results indicate how the optimization procedure allows to open wide BGs relative to both types of polarization (Figure 11c), with non-overlapping polarization modes between 350 Hz – 1750 Hz, 2100 Hz – 4550 Hz and 6300 Hz – 7000 Hz (Figure 11d). The final band diagram also shows several almost zero group velocity out-of-plane modes starting below 2000 Hz and going up to approximately 4500 Hz. This suggests that the thin outer threads play a major role in determining this type of behavior. Curiously, the case $N_p = 1$ was also tested in the metric proposed in Eq. (10) and obtained very similar results, even though in the latter case, each component of the metric is practically restricted to binary values.
<figure><img src="image_13.png" /><figcaption>Figure 11: Band diagrams and histograms accounting for (●) out-of-plane and (●) in-plane polarization modes. The initial band diagram (a) and its corresponding propagating frequency histograms (b) show each of the 20 frequency bins with wave modes from each type of polarization. After the optimization, the final band diagram and unit cell (c) – obtained using the optimization objective $\phi^{(P)}$ (Eq. (10)) – with its corresponding histogram (d) show 13 frequency bins with a single type of polarization and a few with non polarized modes.</figcaption></figure>

### 3.5. Dirac cones

Finally, the results of the optimization process considering the diameters of threads and distribution of added masses to create Dirac cones using the problem stated in Eq. (12) and the objective function given in Eq. (11) are presented, with $\Delta\omega = 0.2\pi$ rad/s. The same previous initial geometry, material properties, and optimization constraints are considered. Using the proposed optimization process, the radii and mass vectors that create a Dirac cone ($D_1$) between the 10th and 11th bands are obtained as

$$ \mathbf{r}_{\text{opt}}^{(D_1)} = \{ 0.944 \ 0.712 \ 0.500 \ 0.741 \ 2.391 \ 0.500 \} \text{ mm,} \qquad (23a) $$

$$ \mathbf{m}_{\text{opt}}^{(D_1)} = \{ 0.000 \ 0.348 \ 0.007 \} \text{ g.} \qquad (23b) $$

It is also possible to experiment on creating another Dirac cone ($D_2$) between the 17th and
18th bands, yielding the radii and mass vectors

$$ \mathbf{r}_{\text{opt}}^{(D_2)} = \{ 1.252 \ 0.500 \ 2.950 \ 0.863 \ 2.937 \ 0.500 \} \text{ mm,} \qquad (24a) $$

$$ \mathbf{m}_{\text{opt}}^{(D_2)} = \{ 0.000 \ 0.000 \ 0.080 \} \text{ g.} \qquad (24b) $$

The corresponding band diagrams, unit cells, and dispersion surfaces indicating the 6 cones corresponding to the vertices of the FBZ (high-symmetry K point) representing the Dirac cones $D_1$ and $D_2$ are shown in Figure 12.

<figure><img src="image_14.png" /><figcaption>Figure 12: Band diagrams, corresponding unit cells, and dispersion surfaces computed for the generated Dirac cones obtained using the optimization objectives $\phi_n^{(D)}$ (Eq. (11)) for $n = 10$ and $n = 17$, respectively: $D_1$ (2839 Hz) with its (a) band diagram, unit cell, (b) and dispersion surface and $D_2$ (4298 Hz) with its (c) band diagram, unit cell, and (d) dispersion surface.</figcaption></figure>

When comparing band diagrams shown in Figures 12a and 12c, it is worth noticing that the Dirac cones $D_1$ and $D_2$ seem to be associated with the same localized wave modes (see Figure 6 for the description of wave modes), but at different frequencies (2839 Hz and 4298 Hz, respectively).
The resemblance between these wave modes is further reinforced by the imaginary part of the wave vectors, which have analogous partial BGs in the $\Gamma\text{M}$ direction, and the dispersion surfaces (Figures 12b and 12d, respectively). Also, the difference in the band ordering ($D_1$ is formed between the 10th and 11th bands, while $D_2$ is formed between the 17th and 18th bands) is due to the addition of several localized modes in the low frequency range below $D_2$. Furthermore, the optimized parameters values (Eqs. (23) and (24)) show that, in this case, mass addition is less necessary to create Dirac cones for increasing values of frequency (4.19 g for $D_1$ and 0.080 g for $D_2$, respectively).

### 4. Concluding remarks

In summary, several two-dimensional PCs inspired by distinctive geometric characteristics of spider web structures were designed and optimized, incorporating additional features for improved dynamic behavior for applications. Starting from simple spider orb web models, the effects of impedance modulation using variable diameters of radial and viscid threads and the addition of point masses at the junctions of threads to manipulate the dispersion characteristics were investigated. The proposed designs consider a single material, which allows realizations of single-phase PCs amenable to additive manufacturing. In the case of hexagonal periodic lattices with all elements having the same diameters, low- to mid-frequency wave modes are dominated by rigid-body-like (translational, bending, and torsional) wave modes and mid- to high-frequencies usually show localized wave modes. The design freedom in the proposed lattice yields a great number of possible combinations, and optimization techniques were employed to achieve objectives such as (i) wide normalized BG widths, (ii) opening a BG at a fundamental frequency while preserving or attenuating its higher harmonics, (iii) BGs with respect to distinct polarization modes, and (iv) generating Dirac cones.

The first of these objectives is relatively standard and has been used to validate our optimization process. The optimization of thread diameters confirms expected results, i.e., low-frequency BGs can be obtained using thick central threads with a soft boundary configuration, resembling a heavy central resonator. Less obvious configurations can also be achieved and yield wider BGs through the manipulation of localized wave modes. The addition of masses alone is not efficient in opening BGs, but can enhance BGs that were previously opened using optimized thread diameters.

When optimizing the structure to achieve BGs at a chosen fundamental frequency, similar configurations may yield distinct objectives, such as preserving or attenuating higher harmonics. The shifting of thick viscid layers in one direction while shifting masses in the opposite direction may open several BGs at higher harmonics where none existed before. In the cases where band diagrams are not very different and harmonics are already close to similar BGs between different objectives, mass reallocation seems more necessary to achieve each objective.

It was also demonstrated that the highly tailorable impedance of the system can be used to achieve configurations that yield BGs with respect to distinct polarization modes. A proposed metric was tested and yielded optimization parameters indicating that wide frequency ranges with separate polarization modes, especially due to highly localized wave modes at the outer threads. This aspect will be explored more in depth in future work.

Finally, it was shown that the presented PC structure may present configurations that yield Dirac cones, revealed using a newly proposed metric, also suggesting that mass addition is not as necessary as in the case of BG nucleation.
In conclusion, the presented optimization process is entirely general and can be used to manip-
ulate several dispersion characteristics according to the requirements of the desired application but
requires a structure with tailorable characteristics. The proposed spider web-inspired PC provides
a versatile structure suitable for optimization procedures, given the richness of available vibration
modes and their tunability through the variation of geometric parameters and localized masses.

Acknowledgments

VFDP, FB, MM, and NMP are supported by the EU H2020 FET Open “Boheme” grant No. 863179.

References

[1] W. Eberhard, Spider webs: Behavior, Function, and Evolution, University of Chicago Press, 2020.

[2] G. Greco, V. Mastellari, C. Holland, N. M. Pugno, Comparing modern and classical perspectives on spider silks and webs, Perspectives on Science 29 (2) (2021) 133–156.

[3] E. J. Roberson, M. J. Chips, W. P. Carson, T. P. Rooney, Deer herbivory reduces web-building spider abundance by simplifying forest vegetation structure, PeerJ 4 (2016) e2538.

[4] M. Miniaci, A. Krushynska, A. B. Movchan, F. Bosia, N. M. Pugno, Spider web-inspired acoustic metamaterials, Applied Physics Letters 109 (7) (2016) 071905.

[5] S. W. Cranford, A. Tarakanova, N. M. Pugno, M. J. Buehler, Nonlinear material behaviour of spider silk yields robust webs, Nature 482 (7383) (2012) 72–76.

[6] B. Mortimer, A. Soler, C. R. Siviour, R. Zaera, F. Vollrath, Tuning the instrument: sonic properties in the spider's web, Journal of The Royal Society Interface 13 (122) (2016) 20160341.

[7] G. Greco, M. F. Pantano, B. Mazzolai, N. M. Pugno, Imaging and mechanical characterization of different junctions in spider orb webs, Scientific Reports 9 (2019) 5776.

[8] P. Zhang, A. C. To, Broadband wave filtering of bioinspired hierarchical phononic crystal, Applied Physics Letters 102 (12) (2013) 121910.

[9] Y. Chen, L. Wang, Multiband wave filtering and waveguiding in bio-inspired hierarchical composites, Extreme Mechanics Letters 5 (2015) 18–24.

[10] M. Miniaci, A. Krushynska, A. S. Gliozzi, N. Kherraz, F. Bosia, N. M. Pugno, Design and fabrication of bioinspired hierarchical dissipative elastic metamaterials, Physical Review Applied 10 (2) (2018) 024012.

[11] Z. Shen, T. R. Neil, D. Robert, B. W. Drinkwater, M. W. Holderied, Biomechanics of a moth scale at ultrasonic frequencies, Proceedings of the National Academy of Sciences 115 (48) (2018) 12200–12205.
[12] M. Rupin, G. Lerosey, J. de Rosny, F. Lemoult, Mimicking the cochlea with an active acoustic metamaterial, New Journal of Physics 21 (2019) 093012.

[13] Y. Chen, L. Wang, Bio-inspired heterogeneous composites for broadband vibration mitigation, Scientific Reports 5 (2016) 17865.

[14] L. Brillouin, Wave Propagation in Periodic Structures: Electric Filters and Crystal Latices, Chemische Reihe. Lehrbucher und Monographien aus dem Gebiete der exakten Wissenschaften, Dover Publications, 1953.

[15] M. Sigalas, E. N. Economou, Band structure of elastic waves in two dimensional systems, Solid State Communications 86 (3) (1993) 141-143.

[16] R. Martínez-Sala, J. Sancho, J. V. Sánchez, V. Gómez, J. Llinares, F. Meseguer, Sound attenuation by sculpture, Nature 378 (6554) (1995) 241.

[17] Z. Liu, X. Zhang, Y. Mao, Y. Y. Zhu, Z. Yang, C. T. Chan, P. Sheng, Locally resonant sonic materials, Science 289 (5485) (2000) 1734-1736.

[18] A. Khelif, B. Aoubiza, S. Mohammadi, A. Adibi, V. Laude, Complete band gaps in two-dimensional phononic crystal slabs, Physical Review E 74 (4) (2006) 046610.

[19] U. Lee, Spectral Element Method in Structural Dynamics, John Wiley & Sons, 2009.

[20] R. V. Craster, S. Guenneau, Acoustic Metamaterials: Negative Refraction, Imaging, Lensing and Cloaking, Springer Series in Materials Science, Springer Netherlands, 2012.

[21] J. F. Doyle, Wave Propagation in Structures: Spectral Analysis Using Fast Discrete Fourier Transforms, Springer Science & Business Media, 2012.

[22] P. A. Deymier, Acoustic Metamaterials and Phononic Crystals, Springer Series in Solid-State Sciences, Springer Berlin Heidelberg, 2013.

[23] V. Laude, Phononic Crystals: Artificial Crystals for Sonic, Acoustic, and Elastic Waves, De Gruyter Studies in Mathematical Physics, De Gruyter, 2015.

[24] G. Ma, P. Sheng, Acoustic metamaterials: From local resonances to broad horizons, Science Advances 2 (2) (2016) e1501595.

[25] L. D'Alessandro, E. Belloni, R. Ardito, A. Corigliano, F. Braghin, Modeling and experimental verification of an ultra-wide bandgap in 3D phononic crystal, Applied Physics Letters 109 (22) (2016) 221907.

[26] L. D'Alessandro, B. Bahr, L. Daniel, D. Weinstein, R. Ardito, Shape optimization of solid-air porous phononic crystal slabs with widest full 3D bandgap for in-plane acoustic waves, Journal of Computational Physics 344 (2017) 465-484.

[27] L. E. Murr, Frontiers of 3D printing/additive manufacturing: from human organs to aircraft fabrication, Journal of Materials Science & Technology 32 (10) (2016) 987-995.
[28] R. Hague, I. Campbell, P. Dickens, Implications on design of rapid manufacturing, Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science 217 (1) (2003) 25–30.

[29] Z. Jia, Y. Chen, H. Yang, L. Wang, Designing phononic crystals with wide and robust band gaps, Physical Review Applied 9 (2018) 044021.

[30] C. Claeys, K. Vergote, P. Sas, W. Desmet, On the potential of tuned resonators to obtain low-frequency vibrational stop bands in periodic panels, Journal of Sound and Vibration 332 (6) (2013) 1418–1436.

[31] L. Tang, L. Cheng, Periodic plates with tunneled acoustic-black-holes for directional band gap generation, Mechanical Systems and Signal Processing 133 (2019) 106257.

[32] A. Bibi, H. Liu, J. L. Xue, Y. X. Fan, Z. Y. Tao, Manipulation of the first stop band in periodically corrugated elastic layers via different profiles, Wave Motion 88 (2019) 205–213.

[33] V. F. Dal Poggetto, J. R. F. Arruda, Widening wave band gaps of periodic plates via shape optimization using spatial Fourier coefficients, Mechanical Systems and Signal Processing 147 (2021) 107098.

[34] C. Claeys, E. Deckers, B. Pluymers, W. Desmet, A lightweight vibro-acoustic metamaterial demonstrator: Numerical and experimental investigation, Mechanical Systems and Signal Processing 70–71 (2016) 853–880.

[35] E. J. P. Miranda Jr., E. D. Nobrega, A. H. R. Ferreira, J. M. C. Dos Santos, Flexural wave band gaps in a multi-resonator elastic metamaterial plate using Kirchhoff-Love theory, Mechanical Systems and Signal Processing 116 (2019) 480–504.

[36] P. Gao, A. Climente, J. Sánchez-Dehesa, L. Wu, Single-phase metamaterial plates for broadband vibration suppression at low frequencies, Journal of Sound and Vibration 444 (2019) 108–126.

[37] V. Romero-García, J. V. Sánchez-Pérez, L. M. García-Raffi, J. M. Herrero, S. García-Nieto, X. Blasco, Hole distribution in phononic crystals: Design and optimization, The Journal of the Acoustical Society of America 125 (6) (2009) 3774–3783.

[38] Z. F. Liu, B. Wu, C. F. He, Band-gap optimization of two-dimensional phononic crystals based on genetic algorithm and FPWE, Waves in Random and Complex Media 24 (3) (2014) 286–305.

[39] X. K. Han, Z. Zhang, Topological optimization of phononic crystal thin plate by a genetic algorithm, Scientific Reports 9 (2019) 8331.

[40] E. Andreassen, J. S. Jensen, Topology optimization of periodic microstructures for enhanced dynamic properties of viscoelastic composite materials, Structural and Multidisciplinary Optimization 49 (2014) 695–705.
[41] L. Xie, B. Xia, J. Liu, G. Huang, J. Lei, An improved fast plane wave expansion method for topology optimization of phononic crystals, International Journal of Mechanical Sciences 120 (2017) 171–181.

[42] L. Xie, B. Xia, G. Huang, J. Lei, J. Liu, Topology optimization of phononic crystals with uncertainties, Structural and Multidisciplinary Optimization 56 (2017) 1319–1339.

[43] M. I. Hussein, K. Hamza, G. M. Hulbert, R. A. Scott, K. Saitou, Multiobjective evolutionary optimization of periodic layered materials for desired wave dispersion characteristics, Structural and Multidisciplinary Optimization 31 (2006) 60–75.

[44] V. F. Dal Poggetto, A. L. Serpa, J. R. F. Arruda, Optimization of local resonators for the reduction of lateral vibrations of a skyscraper, Journal of Sound and Vibration 446 (2019) 57–72.

[45] A. O. Krushynska, F. Bosia, M. Miniaci, N. M. Pugno, Spider web-structured labyrinthine acoustic metamaterials for low-frequency sound control, New Journal of Physics 19 (10) (2017) 105001.

[46] R. D. Cook, D. S. Malkus, M. E. Plesha, R. J. Witt, Concepts and Applications of Finite Element Analysis, 4th Edition, Wiley, 2001.

[47] C. C. Claeys, P. Sas, W. Desmet, On the acoustic radiation efficiency of local resonance based stop band materials, Journal of Sound and Vibration 333 (14) (2014) 3203–3213.

[48] M. Miniaci, M. Mazzotti, M. Radzieński, N. Kherraz, P. Kudela, W. Ostachowicz, B. Morvan, F. Bosia, N. M. Pugno, Experimental observation of a large low-frequency band gap in a polymer waveguide, Frontiers in Materials 5 (2018) 8.

[49] F. Maurin, C. Claeys, E. Deckers, W. Desmet, Probability that a band-gap extremum is located on the irreducible Brillouin-zone contour for the 17 different plane crystallographic lattices, International Journal of Solids and Structures 135 (2018) 26–36.

[50] F. Bloch, Über die Quantenmechanik der Elektronen in Kristallgittern, Zeitschrift für Physik 52 (1929) 555–600.

[51] G. Yi, B. D. Youn, A comprehensive survey on topology optimization of phononic crystals, Structural and Multidisciplinary Optimization 54 (5) (2016) 1315–1344.

[52] Z. Zhang, G. M. Ji, X. K. Han, Optimization design of a novel zigzag lattice phononic crystal with holes, International Journal of Modern Physics B 33 (13) (2019) 1950124.

[53] M. Miniaci, A. S. Gliozzi, B. Morvan, A. Krushynska, F. Bosia, M. Scalerandi, N. M. Pugno, Proof of concept for an ultrasensitive technique to detect and localize sources of elastic non-linearity using phononic crystals, Physical Review Letters 118 (2017) 214301.

[54] A. S. Gliozzi, M. Miniaci, A. O. Krushynska, B. Morvan, M. Scalerandi, N. M. Pugno, F. Bosia, Proof of concept of a frequency-preserving and time-invariant metamaterial-based nonlinear acoustic diode, Scientific Reports 9 (2019) 9560.
[55] V. Laude, Y. Achaoui, S. Benchabane, A. Khelif, Evanescent Bloch waves and the complex band structure of phononic crystals, Physical Review B 80 (2009) 092301.

[56] A. O. Krushynska, M. Miniaci, F. Bosia, N. M. Pugno, Coupling local resonance with bragg band gaps in single-phase mechanical metamaterials, Extreme Mechanics Letters 12 (2017) 30–36.

[57] A. O. Krushynska, V. G. Kouznetsova, M. G. D. Geers, Visco-elastic effects on wave dispersion in three-phase acoustic metamaterials, Journal of the Mechanics and Physics of Solids 96 (2016) 29–47.

[58] Y. Lai, Y. Wu, P. Sheng, Z.-Q. Zhang, Hybrid elastic solids, Nature Materials 10 (2011) 620–624.

[59] G. Ma, C. Fu, G. Wang, P. Del Hougne, J. Christensen, Y. Lai, P. Sheng, Polarization bandgaps and fluid-like elasticity in fully solid elastic metamaterials, Nature Communications 7 (2016) 13536.

[60] A. K. Geim, K. S. Novoselov, The rise of graphene, in: Nanoscience and Technology, World Scientific, 2010, pp. 11–19.

[61] T. O. Wehling, A. M. Black-Schaffer, A. V. Balatsky, Dirac materials, Advances in Physics 63 (1) (2014) 1–76.

[62] J. Wang, S. Deng, Z. Liu, Z. Liu, The rare two-dimensional materials with Dirac cones, National Science Review 2 (1) (2015) 22–39.

[63] M. Miniaci, R. K. Pal, B. Morvan, M. Ruzzene, Experimental observation of topologically protected helical edge modes in patterned elastic plates, Physical Review X 8 (2018) 031074.

[64] M. Miniaci, R. K. Pal, R. Manna, M. Ruzzene, Valley-based splitting of topologically protected helical waves in elastic plates, Physical Review B 100 (2019) 024304.

[65] K. Schittkowski, NLPQL: A FORTRAN subroutine solving constrained nonlinear programming problems, Annals of Operations Research 5 (1986) 485–500.

[66] J. Nocedal, S. J. Wright, Numerical Optimization, Springer Science & Business Media, 2006.

[67] R. Fletcher, Practical Methods of Optimization, John Wiley & Sons, 2013.

[68] The MathWorks, Inc., MATLAB Optimization Toolbox, Natick, MA, US (2020). URL https://www.mathworks.com/help/optim/

[69] P. E. Gill, W. Murray, M. H. Wright, Practical Optimization, Academic Press New York, 1981.

[70] M. J. D. Powell, Variable metric methods for constrained optimization, in: Mathematical Programming the State of the Art, Springer, 1983, pp. 288–311.