*2.2. Material Behaviors*

In this research, we used the polylactic acid (PLA) filament with a diameter of 1.75 mm and glass transition temperature of 65 ◦C. Objects were fabricated using a 3Dgence DOUBLE printer developed by 3Dgence. This is a low-cost desktop 3D printer that extrudes 1.75 mm filaments with a 0.4 mm nozzle. The 3Dgence Slicer software was used to produce G-code files from STL files and command and control parameters of liquefier temperature and printing speed. Beam-like elements were 4D-printed with dimensions of (30 × 1.6 × 1) mm for length, width, and thickness, respectively. Each printing layer was considered to have a 0.1 mm thickness. The print raster was assumed to be along the length direction. The temperature of the liquefier, build tray, and chamber were set 210, 24, and 24 ◦C, respectively. In all the 4D printings, unless otherwise stated, the printing speed was set to 5 mm/s in which the printer does not induce any prestrains in the materials.

Dynamic-mechanical analyzer (DMA, Q800, TA Instruments, New Castle, DE, USA) was first implemented to determine the temperature-dependent mechanical properties of 3D-printed PLAs. DMA test was conducted in an axial tensile way with the frequency of force oscillation 1 Hz and heating rate 5 ◦C/min ranging from 30 to 93 ◦C. The applied dynamic stress was 1.5 times the static stress. Further, the thermomechanical behavior of the 4D-printed sample in terms of storage modulus, *Es*, and phase lag, tan(δ), is shown in Figure 2. It is seen that the storage modulus is reduced drastically as the material is heated to 60 ◦C. It is also observed that the phase lag graph peaks at 65 ◦C, which is assumed to be the glass transition temperature.

Five beam-like elements were 4D-printed at five different speeds (Sp = 5, 10, 20, 40, 70 mm/s). Figure 3 shows the beam configuration after 4D printing. After 4D printing, five samples were heated by dipping into the hot water with a prescribed temperature of 85 ◦C (20 ◦C higher than the transition temperature) and then were cooled down to the room temperature, 24 ◦C. Figure 4 depicts the sample configuration after the heating–cooling process. As can be seen, samples with a straight temporary shape may transform into curved beams. This means that the samples may already be programmed and prestained during the 4D printing process. Figure 4a shows that the sample printed at a low speed of 5 mm/s does not undergo any changes by heating. This implies that this printing speed is not enough to produce any prestrains. However, the configuration presented in Figure 4b indicates that increasing the printing speed to 10 mm/s can cause a slight curvature after heating. Thus, this speed can be considered as a transient speed. The configuration change could be due to an unbalanced FG prestrain regime induced through the thickness direction during the 4D printing process. When dipping the samples into hot water, the unbalanced prestrain, with an increasing trend through the thickness upward from the lower layer to the upper layer, was recovered, leading to a self-bending movement. The mismatch in free-strain recovery inducing curvatures enabled the overall configuration to change toward the top layer. It is worthwhile to mention that an increase in the 4D printing speed increased the bending angle and curvature. The faster the 4D printing, the greater the prestrain and consequently the deformation. Finally, experiments revealed that the FDM 4D printing technology has high potential in fabricating and programming adaptive objects with self-bending features.

**Figure 2.** Dynamic-mechanical analyzer (DMA) test for the 3D-printed polylactic acid (PLA).

**Figure 3.** The beam configuration after 4D printing.

**Figure 4.** The configuration of the samples 4D-printed at different speeds of (**a**) 5, (**b**) 10, (**c**) 20, (**d**) 40, and (**e**) 70 mm/s after the heating–cooling process.

## **3. Theoretical Modeling**

The 4D-printed metastructures with adaptive wave propagation feature can be designed based on the theoretical understanding of the FG programming process and SME of 4D-printed SMP elements. This section is devoted to developing constitutive models and mathematical formulation to predict SME and wave propagation in 4D-printed SMPs.

## *3.1. SMP Model*

In this division, a phenomenological constitutive model presented initially in [17,29] is reformulated to describe shape programming and recovery processes in the 4D-printed structures. Analytical straightforward closed-form solutions are also derived.

The SMPs consist of two phases, so-called glassy and rubbery phases, which are stable at temperatures above and below the transition temperature, respectively. As SMPs show a mixture of glassy and rubbery phases, volume fractions of the rubbery and glassy phases describe the state of SMPs. The following scalar variables characterize them:

$$
\eta\_{\mathcal{S}} = \frac{V\_{\mathcal{S}}}{V} \qquad \qquad \eta\_{r} = \frac{V\_{r}}{V} \tag{1}
$$

where *Vg* and *Vr* indicate the volume of the glassy and rubbery phases, respectively. The subscripts '*g*' and '*r*' stand for glassy and rubbery phases, respectively, here and henceforth. The change between these phases is considered to be only a function of temperature as a generally well-known assumption. It implies that η*g* and η*r* are just functions of the temperature. Because in every moment the summation of these parameters must satisfy the unity (η*g* + η*r* = 1), the volume fraction of the rubbery phase can be expressed in terms of the glassy one as:

$$
\eta\_r(T) = 1 - \eta\_\mathcal{X}(T) \tag{2}
$$

Considering experimental results from the DMA test, η*g* can explicitly be interpolated by a hyperbolic function as [17,29]:

$$\eta\_{\mathcal{S}} = \frac{\tanh(a\_1 T\_{\mathcal{S}} - a\_2 T) - \tanh(a\_1 T\_{\mathcal{S}} - a\_2 T\_{\mathcal{h}})}{\tanh(a\_1 T\_{\mathcal{S}} - a\_2 T\_{\mathcal{h}}) - \tanh(a\_1 T\_{\mathcal{S}} - a\_2 T\_{\mathcal{l}})} \tag{3}$$

in which *<sup>a</sup>*1and*<sup>a</sup>*2 are determined by adjusting the DMA curve.

The glassy and rubbery phases in SMPs are assumed to be linked in series, i.e., σ = <sup>σ</sup>*g* = σ*r*, where σ denotes the stress. Considering the fact that the 4D-printed objects may experience small strains and moderately large rotations, a small strain regime is assumed. Consequently, an additive strain decomposition is adopted as:

$$
\varepsilon = \eta\_{\mathcal{K}} \varepsilon\_{\mathcal{K}} + (1 - \eta\_{\mathcal{K}}) \varepsilon\_r + \varepsilon\_{in} + \varepsilon\_{th} \tag{4}
$$

where ε denotes the total strain while <sup>ε</sup>*g* and ε*r* designate the strains of the glassy and rubbery phases, respectively. Also, ε*in* is the inelastic strain due to the SMP phase transformation, while ε*th* indicates the thermal strain that can be expressed as ε*th* = *T T*0 α*r* + (<sup>α</sup>*g* − <sup>α</sup>*r*) <sup>η</sup>*g*(*T*) *dT* in which α*r* and <sup>α</sup>*g* denote thermal expansion in rubbery and glassy phases, respectively, and *T*0 is the reference temperature.

Next, ε*in*, associated with the glassy–rubbery phase transformation mechanism, will be formulated. During the cooling step, the rubbery phase changes into the glassy one, and its strain is stored in the SMP material. The strain storage is assumed to be proportional to η*g* based on the rubbery phase strain. In the heating step, the stored strain is assumed to be released gradually, proportional to η*g* concerning the preceding glassy phase and strain storage. Therefore, ε*in* can mathematically be formulated in a rating format as:

$$
\dot{\varepsilon}\_{\dot{m}} = \begin{cases}
\dot{\eta}\_{\mathcal{S}} \, \varepsilon\_r & \dot{T} < 0 \\
\frac{\dot{\eta}\_{\mathcal{S}}}{\eta\_{\mathcal{S}}} \, \varepsilon\_{\dot{m}} & \dot{T} > 0
\end{cases} \tag{5}
$$

By considering ε and *T* as external control variables and <sup>ε</sup>*g*, ε*r*, ε*in*, and η*g* as internal variables, introducing Helmholtz free energy density functions, implementing the second law of thermodynamics in the sense of the Clausius–Duhem inequality, and following a standard argumen<sup>t</sup> [30], the stress state can be obtained as:

$$
\sigma = \sigma\_{\S} = \sigma\_{\Gamma} \tag{6}
$$

This equation is consistent with the taken assumption that the stress in each phase is equal in the series model. The glassy- and rubbery-phase stresses are also derived as:

$$
\sigma\_{\mathcal{S}} = \mathbb{C}\_{\mathcal{S}} \varepsilon\_{\mathcal{S}} \; , \; \sigma\_{r} = \mathbb{C}\_{r} \varepsilon\_{r} \tag{7}
$$

in which *C* signifies the elasticity sti ffness matrix defined as:

$$\mathbf{C} = \frac{E}{(1+\nu)(1-2\nu)} \begin{bmatrix} 1-\nu & \nu & \nu & 0 & 0 & 0\\ \nu & 1-\nu & \nu & 0 & 0 & 0\\ \nu & \nu & 1-\nu & 0 & 0 & 0\\ 0 & 0 & 0 & \frac{(1-2\nu)}{2} & 0 & 0\\ 0 & 0 & 0 & 0 & \frac{(1-2\nu)}{2} & 0\\ 0 & 0 & 0 & 0 & 0 & \frac{(1-2\nu)}{2} \end{bmatrix} \tag{8}$$

where *E* and *v* denote Young's modulus and Poisson's ratio, respectively. The stress–strain relationship can be derived by substituting Equation (7) into Equation (4) as:

$$
\sigma = \left( S\_r + \eta\_{\mathcal{K}} (S\_{\mathcal{K}} - S\_r) \right)^{-1} (\varepsilon - \varepsilon\_{\rm in} - \varepsilon\_{\rm th}) \tag{9}
$$

where *S* indicates the compliance matrix defined as *C*−1.

From a numerical viewpoint, the nonlinear SMP behaviors can be treated in an explicit time-discrete stress/strain-temperature-driven framework. The time domain [0,*t*] is discretized into increments, and the evolution equation is solved over the local band [*t n*,*t <sup>n</sup>*+<sup>1</sup>]. The superscript *n* + 1 shows the current step, while the superscript *n* denotes the previous step. The inelastic strain can be computed by applying the linearized implicit backward-Euler integration method to the flow rule (5) as:

$$\varepsilon\_{in}^{n+1} = \begin{cases} \varepsilon\_{in}^{n} + \Delta \eta\_{\mathcal{S}}^{n+1} \varepsilon\_{r}^{n+1} & \dot{T} < 0 \\\ \varepsilon\_{in}^{n} + \frac{\Delta \eta\_{\mathcal{S}}^{n+1}}{\eta\_{\mathcal{S}}^{n+1}} \varepsilon\_{in}^{n+1} & \dot{T} > 0 \end{cases} \tag{10}$$

where

$$
\Delta \eta\_{\mathcal{S}}^{n+1} = \eta\_{\mathcal{S}}^{n+1} - \eta\_{\mathcal{S}}^{n} \tag{11}
$$

By using Equations (7) and (9) to substitute the rubbery strain and the stress, and performing some mathematical simplifications, ε*in* defined in (10) can be updated for cooling and heating steps in stress and strain control ways as:

**Cooling:**

$$\text{stress control mode} \rightarrow \varepsilon\_{in}^{n+1} = \varepsilon\_{in}^{n} + \Delta \eta\_{\mathcal{S}}^{n+1} S\_r^{n+1} \sigma^{n+1} \tag{12}$$

$$\text{strain control mode} \rightarrow \boldsymbol{\varepsilon}\_{\text{in}}^{n+1} = \left(\boldsymbol{I} + \Delta \boldsymbol{\eta}\_{\text{g}}^{n+1} \mathbf{S}\_{r}^{n+1} \mathbf{C}\_{\varepsilon}^{n+1}\right)^{-1} \left(\boldsymbol{\varepsilon}\_{\text{in}}^{n} + \Delta \boldsymbol{\eta}\_{\text{g}}^{n+1} \mathbf{S}\_{r}^{n+1} \mathbf{C}\_{\varepsilon}^{n+1} \left(\boldsymbol{\varepsilon}^{n+1} - \boldsymbol{\varepsilon}\_{\text{th}}^{n+1}\right)\right) \tag{13}$$

**Heating:**

$$
\kappa\_{in}^{n+1} = \frac{\eta\_{\mathcal{S}}^{n+1}}{\eta\_{\mathcal{S}}^n} \kappa\_{in}^n \tag{14}
$$

Finally, by considering updated inelastic strain, the stress–strain constitutive Equation (9) can be discretized and unified for heating and cooling processes as:

$$
\sigma^{n+1} = \mathbb{C}\_D^{n+1} (\varepsilon^{n+1} - \varepsilon \varepsilon\_{in}^{n+1} - \varepsilon\_{th}^{n+1}) \tag{15}
$$

where the so-called unified stiffness matrix *CD* and parameter ς for cooling and heating processes are defined as:

$$\begin{cases} \begin{array}{c} \mathsf{C}\_{D}^{n+1} = \left(I + \Lambda \eta\_{\mathcal{G}}^{n+1} S\_{r}^{n+1} \mathbf{C}\_{\varepsilon}^{n+1}\right)^{-1} \mathbf{C}\_{\varepsilon}^{n+1}, \zeta = 1 & \dot{T} < 0\\ \updownarrow \mathcal{C}\_{D}^{n+1} = \mathsf{C}\_{\varepsilon \prime} \ \zeta = \frac{\eta\_{\mathcal{G}}^{n+1}}{\eta\_{\mathcal{G}}^{n}} & \dot{T} > 0 \end{array} \tag{16}$$

#### *3.2. Wave Propagation Model*

The wave propagation analysis of the architected periodic structures can be carried out by using the Bloch's theorem for local resonance. Based on this theory, the displacement of each node in a chosen unit cell in the region of a periodic structure depends only on the displacement field of the equal node in the reference unit cell ( → *<sup>U</sup>*Ref<sup>→</sup>*r*). The formulation of this theorem is implemented to solve the equations of motion in the periodic boundary conditions. This concept is stated as [31]:

$$\overrightarrow{LI}\left(\overrightarrow{r} + \overrightarrow{R}, t\right) = \overrightarrow{\mathcal{U}}\_{\text{Ref}}\left(\overrightarrow{r}\right) \exp\left(i\left[\overrightarrow{\kappa}.\left(\overrightarrow{r} + \overrightarrow{R}\right) - \omega t\right]\right) \tag{17}$$

where →*r* and →*R* are position and lattice vectors, respectively. Also, the Bloch wave vector →κ in the 2D periodicity is considered as →κ = <sup>κ</sup>*x*, <sup>κ</sup>*y*, where κ*x* and <sup>κ</sup>*y* denote the phase constants which are measures of the phase variations over one unit cell in two directions of the periodicity (X–Y). Also, ω and *t* refer to frequency and time, respectively. As illustrated in Figure 5, the periodicity in two directions is defined by the direct vectors →*a x* and →*a y*. Hence, the lattice vector is described as:

$$
\overrightarrow{\vec{R}} = n\_x \overrightarrow{\vec{a}}\_x + n\_y \overrightarrow{\vec{a}}\_y \quad , \quad n\_{\ge \prime} n\_y = 0, \pm 1, \pm 2, \cdots \tag{18}
$$

**Figure 5.** Direct vectors in a periodic model.

Since the Bloch wave vector varies in the Brillouin zone, the procedures to calculate reciprocal vectors and the first Brillouin zone are explained. Equation (19) formulates the reciprocal vectors → *b x* and → *b y* as [32]:

$$\begin{cases} \begin{array}{l} \overrightarrow{b}\_{X} = \frac{2\pi}{\left| \overrightarrow{a}\_{x} \times \overrightarrow{a}\_{y} \right|} \end{array} \begin{bmatrix} a\_{y} \chi \\ -a\_{yX} \end{bmatrix} \\\ \overrightarrow{b}\_{Y} = \frac{2\pi}{\left| \overrightarrow{a}\_{x} \times \overrightarrow{a}\_{y} \right|} \begin{bmatrix} -a\_{x} \chi \\ a\_{yX} \end{bmatrix} \end{cases} \tag{19}$$

where (*axX*, *axY*) and *ayX*, *ayY*are components of direct vectors →*a x* and →*a y* along the X- and Y-directions, respectively. The first Brillouin zone (FBZ) and irreducible Brillouin zone (IBZ) are depicted in Figure 6. The IBZ is the smallest part of the FBZ representing all the high-symmetry points in which the wave propagation is analyzed [33].

**Figure 6.** First Brillouin zone (FBZ, dashed square) and irreducible Brillouin zone (IBZ, colored triangle).

In this research, numerical simulations have been produced by employing an FE software. Each unit cell is discretized into finite elements. The dispersion relation is an eigenvalue problem that is solved By COMSOL Multiphysics 4.3. The problem has three unknown parameters (<sup>κ</sup>*x*, <sup>κ</sup>*y*, and ω) in which the wave vectors are examined in the IBZ (G-X-M-G), and the problem is determined to calculate the eigenfrequencies. Convergence studies of the mesh size have also been performed to achieve the converged results accurately to three significant digits.

## *3.3. FE Solution*

#### 3.3.1. In-House FE Code

In order to replicate thermomechanical behaviors of FG 4D-printed elements, a Ritz-based FE solution is developed in MATLAB software. A 3D 20-node quadratic serendipity hexahedron element is considered in this problem. It has 20 nodes so that 8 corner nodes are augmented with 12 side nodes located at the side center. The element also has 3 degrees of freedom per node (*ui*, i = 1, 2, 3). More details on the FE solution and numerical programming can be found in [17].

#### 3.3.2. COMSOL Multiphysics FE Modeling

In this section, thermomechanical behaviors of 4D-printed samples are analyzed by implementing a simple method in COMSOL Multiphysics software. For this purpose, by using the DMA test data, the temperature-dependent Young's modulus is implemented in the COMSOL Multiphysics. Table 1 shows the dependency of Young's modulus on the temperature.

**Table 1.** The temperature-dependent Young's modulus from the DMA test.


As described already, during 4D printing, a prestrain that varies through the thickness is induced in the object producing an FG structure. For modeling 4D-printed structures in COMSOL Multiphysics, the object can be divided into multiple sections with variable thermal expansion. In this study, the 4D-printed beam-like structures are divided into six sections. Figure 7 illustrates a discretized form of a printed beam with different thermal expansion coefficients.

**Figure 7.** Discretized 4D-printed sample.

The thermal expansion of each layer is chosen to replicate a configuration similar to the experiments as depicted in Figure 4 for a specific printing speed. Table 2 indicates the thermal expansion coefficient assumed for each layer.

**Table 2.** The thermal expansion coefficient of each layer for different printing speeds.


A relationship between thermal expansion coefficients and printing speed can be formulated as:

$$a\_i = C\_1 S^3 p + C\_2 S^2 p + C\_3 S p + C\_4 \tag{20} \le S p \le 70 \tag{20}$$

where C1, C2, and C3 are constants defined for each layer in Table 3.


**Table 3.** The constant of thermal expansion interpolation function for each printing layer.

Figure 8 shows the deformed configuration obtained from the FE COMSOL Multiphysics simulation of self-bending 4D-printed beams after the heating–cooling process. In order to characterize the configuration of the printed samples after the heating–cooling process, three geometric parameters are considered. For this purpose, we use parameters *R*1,*R*2, and *R*3 which describe the outer length, opening, and depth of mid surface, respectively, as shown in Figure 8c. In order to determine the

accuracy and efficiency of the simple method implemented in COMSOL Multiphysics, the geometric features obtained from the experiments, FE COMSOL Multiphysics, and in-house FE code are compared in Table 4. It can be concluded that the simulation results of FE COMSOL Multiphysics are in good agreemen<sup>t</sup> with those measured from experiments and calculated by the in-house FE solution. This way, the reliability and accuracy of setting variable thermal expansion coefficients in the COMSOL Multiphysics in replicating the self-bending feature observed in the 4D-printed samples is validated.

**Figure 8.** Finite element (FE) COMSOL Multiphysics simulation of the samples 4D-printed with different speeds of (**a**) 10, (**b**) 20, (**c**) 40, and (**d**) 70 mm/s after the heating–cooling process.


**Table 4.** The geometric parameters of the beam-like 4D structures after the heating–cooling process.

#### *3.4. Periodic Structural Design*

In this section, two periodic architected metastructures with adaptive dynamical characteristics are conceptually proposed. These metastructures are made of passive mainframes printed at a low speed so that no prestrain is induce, and active beam-like members printed at high speed with induced prestrains and self-folding features. They have the potential to be 4D-printed by setting different printing speeds for two nozzles. The PLA material as characterized in Sections 2.2 and 3.3 are considered for 4D printing. The passive main frame of the metastructure is printed at a low speed such as 5 mm/s, while active elements with self-bending features are fabricated with three different 4D printing speeds. Different arrangements have been designed, and their dynamic performance has been examined numerically. Two metastructures have shown a high dynamic performance that will further be examined numerically in the next section. The first adaptive structure consists of active elements connected in parallel and diagonal to the frame, while the second adaptive structure consists of active elements connected in parallel with the frame. For convenience, they are called diagonal and parallel metastructures, respectively. Figure 9 shows the designed structures in which yellow and blue colors signify active and passive elements, respectively.

**Figure 9.** Periodic metastructures with active and passive components: (**a**) diagonal structure; (**b**) parallel structure (the red dashed oval shows the fixed-fixed beam used for the frequency normalization).

#### **4. Results and Discussion**

Periodic active structures are associated with granting propagating and bandgap ranges. In the propagating frequency ranges, the elastic wave propagation is done in all directions, but the bandgap represents a frequency range where the elastic wave propagation is stopped. In this section, COMSOL-based numerical results are presented, revealing how to design adaptive periodic structures with the ability to optimize the dynamical functionality without embedding any additional resonating components. Since it is difficult to actuate all mode shapes in experimental studies, experimental works have been considered in-plane or out-of-plane mode shapes only. However, in the present numerical study, a 3D dynamic case is investigated by FE COMSOL Multiphysics, and all modes of vibrations (i.e., bending, torsion, and elongation) can be measured without any limitation.

In order to verify the FE simulation, the bandgap of a triangular structure is simulated, as shown in Figure 10, and compared with the results of [25]. It is seen that there is an excellent agreemen<sup>t</sup> between the dispersion curves from the FE simulation and [25].

Eigenfrequencies of adaptive periodic structures, as shown in Figure 9, are computed by imposing periodic boundary conditions in different elastic wave vectors which are estimated based on IBZ detailed in the previous section. These eigenfrequencies are finally normalized against Ω = ω /<sup>ω</sup>0, in which ω0 = 22.4 *EI*/(*mL*40) is the first natural frequency of a fixed-fixed single beam, as shown in Figure 9a. The dispersion diagram and some out-of-plane mode shapes of the diagonal metastructure

are illustrated in Figure 11, whereas the mode shapes are depicted in high symmetry points (G, X, or M). As can be seen, there is a wide bandgap (gray square) in the range of Ω = 1.902 − 2.043 which is 4.68% of the frequency ranging between 0 and 3. As shown in Figure 11, the cause of bandgap in this range is the resonance of active parts. Further, there is a flat eigenfrequency before the bandgap which helps to assure its local resonance nature.

**Figure 10.** Verification of the band structure of the triangular topology: (**a**) the current study; (**b**) Ref. [25].

**Figure 11.** Band structure and mode shapes of the diagonal metastructure.

The configuration of the diagonal structure after the heating–cooling process for three printing speeds of active elements is illustrated in Figure 12. As it is expected, the elements 4D-printed faster produce more curvature after thermal activation. It is worth mentioning that the heating could also be another variable to control the curvature. i.e., by partially heating the metastructure with active elements printed at 70 mm/s speed, configurations become like those of thoroughly heated structures with active elements printed at 20 and 40 mm/s. Dispersion curves of diagonal structure with active elements of different 4D printing speeds after the heating–cooling process are depicted in Figures 13–15. Comparing the results presented in Figures 11 and 13 reveal that the bandgap area is adapted from Ω = 1.902 − 2.043 to Ω = 1.751 − 1.812, and the expanse of the bandgap area decreases from 4.68% to 2.01% as well. This phenomenon shows the significant effect of the self-bending feature on the locally resonant filtration. As can be seen, the functional range changes by tuning the 4D printing speed of active elements.

**Figure 12.** The configuration of adaptive periodic diagonal metastructure after heating–cooling process for three different printing speeds: (**a**) 20, (**b**) 40, and (**c**) 70 mm/s.

**Figure 13.** Band structure and mode shapes of the diagonal structure with active elements printed at 20 mm/s after the heating–cooling process.

**Figure 14.** The counterpart of Figure 12 for 40 mm/s.

**Figure 15.** The counterpart of Figure 12 for 70 mm/s.

Figures 13 and 14 reveal that increasing the 4D printing speed from 20 to 40 mm/s does not affect the bandgap range much, changes the stop-band area from 2.01% to 2.07%, and moves the range of bandgap to Ω = 1.801 − 1.863. However, Figures 12c and 15 show an interesting point: after heating–cooling of the diagonal metastructure with the printing speed of 70 mm/s, the whole stop-band area is vanished, and this model allows all the frequencies in the range of 0–3 times of reference frequency to pass. Also, the mode shapes of the structure in some of the high symmetry points and before the bandgap area are depicted in Figure 13. The results presented in Figures 11–15 imply that varying 4D printing speed or heating temperature in the diagonal structure changes the dispersion behaviors significantly and can be manipulated to find an appropriate locally resonant vibration filter. The phenomenon of bandgap switch caused by changing the natural frequency of the active part in the periodic structure is diminished by local resonance changing frequency in different printing speeds.

The parallel metastructure, as shown in Figure 9b displays different dispersion behaviors than the diagonal metastructure (Figure 9a). The dispersion curve for this model is illustrated in Figure 16, where there is no bandgap area, meaning that this structure allows all elastic waves in the frequency range of 0 to 3 to pass. Like diagonal metastructures, the active elements of the parallel structure are also manufactured with three different printing speeds. The configuration of parallel metastructures after the heating–cooling process is depicted in Figure 17, where parts a–c represent the self-bending metastructures 4D-printed at the speeds of 20, 40, and 70 mm/s, respectively. Furthermore, the band structure and mode shapes of parallel structures with self-bending elements of different printing speeds 20, 40, and 70 mm/s after the heating–cooling process are depicted in Figures 18–20, respectively.

As it can be seen in Figures 16 and 18, Figures 19 and 20, the dynamic behaviors of these metastructures are remarkable such that the bandgap area has an increasing–decreasing trend as the structure is heated, revealing self-bending features. Figure 18 shows that the actuated metastructure with self-bending elements 4D-printed at 20 mm/s has a narrow bandgap area in the range of Ω = 1.945 − 1.989 with the amount of 1.48%. However, by using active elements with a higher printing speed of 40 mm/s, the dynamic behaviors change. It is found that the amount of bandgap area increases to 12.32%, and the system exhibits stop-bands in multiple frequencies (Figure 19). These ranges are read as Ω = 2.172 − 2.231, Ω = 2.371 − 2.505, Ω = 2.421 − 2.430, Ω = 2.441 − 2.569, and Ω = 2.594 − 2.765. This implies that this design has a better performance than the others. It can be concluded that this type of 4D-printed architected metastructure has excellent potential in adapting its locally resonant filters from 0 to a significant value such as 12.32%. These bandgaps are generated by Bragg scattering within the medium. In this type of adaptive periodic structure, there is no locally resonant bandgap and the bandgaps are the Bragg type.

**Figure 16.** Band structure and mode shapes of the parallel metastructure.

**Figure 17.** The configuration of adaptive periodic parallel metastructure after heating–cooling process for three different printing speeds: (**a**) 20, (**b**) 40, and (**c**) 70 mm/s.

**Figure 18.** Band structure and mode shapes of the parallel metastructure with active elements printed at 20 mm/s after the heating–cooling process.

**Figure 19.** The counterpart of Figure 17 for 40 mm/s.

**Figure 20.** The counterpart of Figure 17 for 70 mm/s.

Finally, the numerical results presented in Figure 20 reveal that by using active elements with a 4D printing speed of 70 mm/s, the bandgap area vanishes. This means that this metastructure, such as the structure before heating–cooling, as shown in Figure 16, propagates all the locally resonant vibration in all directions.
