Transient analysis of nonlinear locally resonant metamaterials via computational homogenization

Abstract

In this paper, the transient computational homogenization scheme is extended to allow for nonlinear elastodynamic phenomena. The framework is used to analyze wave propagation in a locally resonant metamaterial containing hyperelastic rubber-coated inclusions. The ability to properly simulate realistic nonlinearities in elasto-acoustic metamaterials constitutes a step forward in metamaterial design as, so far, the literature has focused only on academic nonlinear material models and simple lattice structures. The accuracy and efficiency of the framework are assessed by comparing the results with direct numerical simulations for transient dynamic analysis. It is found that the band gap features are adequately captured. The ability of the framework to perform accurate nonlinear transient dynamic analyses of finite-size structures is also demonstrated, along with the significant computational time savings achieved.

1. Introduction

Metamaterials are engineered materials designed to exhibit extraordinary properties. Such materials have been attracting attention in many fields, from optics to acoustics and mechanics [1]. Metamaterials properties originate from their artificially designed microstructures, which may yield unconventional properties, such as negative Poisson’s ratio [2, 3], negative mass density [4], negative refractive index [5–7], negative bulk modulus [8], among others. This results in a broad range of potential applications, most of them involving wave manipulation or attenuation, as used in, e.g., filtering devices devices [9], superlenses [10, 11], shielding structures [12], and cloaking devices [13, 14]. In this paper, the focus is on metamaterials designed to manipulate elasto-acoustic wave propagation.

The most striking feature of such metamaterials is the emergence of band gaps, i.e. frequency bands in which waves cannot propagate or are strongly attenuated. The local resonance phenomenon is one of the mechanisms underlying these forbidden frequency zones. The local resonance band gaps may occur in subwavelength regimes and are typically induced by a negative effective dynamic mass [15]. A convincing experimental demonstration of this phenomenon, within the framework of metamaterials, was designed by Liu et al. [4]. Their design consisted of heavy rubber-coated lead spherical inclusions stacked in a simple cubic lattice embedded in an epoxy matrix. The soft rubber ensures that the heavy core can move independently of the matrix material, acting as a local resonator through the matrix material, waves propagate and may interact with the inclusions. Close to the local resonance frequency, the strong interaction between the matrix and the inclusions and the periodic arrangement of inclusions give rise to wave reflection and attenuation, i.e. triggering the so-called local resonance band gap.

The subwavelength feature of locally resonant elasto-acoustic metamaterials makes them suitable for homogenization. Homogenization theories are powerful tools for the computation of effective properties of composites by averaging their micro-scale behavior [16–18]. In problems where the underlying microstructure should be modeled in detail, computational homogenization [19] can be used. In general, the key feature of computational homogenization schemes is their generality and accuracy compared with other homogenization schemes, while still being efficient compared with direct numerical simulations (i.e. the numerical solution, in space and time, for the dynamics of the full system, including all degrees of freedom in the system, without performing any homogenization or model reduction). Both analytical and computational homogenization frameworks have been recently extended to account for complex transient interactions and micro-inertia effects. However, most proposed frameworks are based on frequency-domain formulations [20–22], which are only convenient for linear elastic materials and for the straightforward computation of band diagrams in periodic materials, and thus band gap identification. For more general cases, involving material nonlinearity and aperiodic time-dependent boundary conditions, it is often necessary to solve the multi-scale problem directly in space and time. A step forward in this direction was made by Pham et al. [23], who extended the classical computational homogenization scheme [19, 24], also called FE², in specific situations, to incorporate microdynamics. Both macro- and micro-scales are governed by the full balance of linear momentum; solving the multi-scale problem results in solving an initial boundary value problem at both scales [23, 25]. Although the new micro–macro kinematics and dual relations for the effective stress tensor and momentum have been derived in a general sense, their numerical implementation has still been restricted to the special case of a linear elasto-acoustic locally resonant metamaterial. More recently, Sridhar et al. [26] have extended that work by using a static–dynamic superposition under the linearity assumption, leading to an enriched reduced-order model, thus enhancing the performance of the framework in incorporating linear microdynamics in the subwavelength regime of elastic metamaterials.

Recently, the role of nonlinearity on the metamaterial response has raised attention as it would enable the design of tunable structures [1]. Nonlinearities may induce amplitude-dependent dispersion and group velocities, and also the propagation of emergent waveforms, such as solitons [27–29]. However, only a few papers in the literature have addressed the effect of nonlinearity on locally resonant metamaterials [30–32]. Furthermore, the material nonlinearity considered so far was usually restricted to academic material models, which do not correspond to any realistic nonlinear material.

This paper aims to demonstrate that the dynamic computational homogenization scheme extended from Pham et al. [23] is applicable to metamaterial microstructures involving nonlinear constituents. For this purpose, a semi-1D locally resonant metamaterial of Liu-type is considered, whereby the computational homogenization scheme is derived in a general sense. The coating layer, responsible for the interaction between the matrix material and the heavy inclusion in Liu’s material, is modeled as a neo-Hookean hyperelastic material, commonly used for rubber materials. The presented implementation of the computational homogenization framework involves a finite-element description for the spatial discretization and the Newmark method for time integration. Unlike the linear implementation provided by Pham et al. [23], to incorporate material nonlinearities, the macroscopic tangents are derived and computed at each integration point. Since the approach makes use of a finite-element discretization, complex geometries can be modeled, thus providing a step forward in the analysis of nonlinear locally resonant elasto-acoustic metamaterials, which have been restricted so far to lattice systems.

This paper is organized as follows. In Section 2, the dynamic computational homogenization framework is briefly presented. In Section 3, an extended nonlinear numerical implementation of the dynamic computational homogenization framework is developed. Next, numerical simulations are carried out to compare the computational homogenization framework results with direct numerical simulations in terms of accuracy and efficiency, as presented in Section 4. Two types of numerical analysis are performed: a free wave propagation analysis and a transient structural dynamic analysis. Finally, in Section 5, conclusions are drawn.

In the following, scalars, vectors, and second-order tensors are denoted by a (or A), $\overrightarrow{a}$ , and $\mathbf{A}$ , respectively. A general m -th order tensor is denoted ${}^m\mathbf{A}$ . Standard tensor operations are denoted as follows: double contraction, $\mathbf{A}:\mathbf{B}=A_{\mathrm{pq}}{B}_{\mathrm{qp}}$ ; dyadic product, $\overrightarrow{a}\otimes \overrightarrow{b}=a_p{b}_q{\overrightarrow{e}}_p\otimes {\overrightarrow{e}}_q$ ; dot product, $\mathbf{A}\cdot \overrightarrow{b}=A_{\mathrm{pq}}{b}_q{\overrightarrow{e}}_p$ , where ${\overrightarrow{e}}_k$ , and $k=1,2,3$ are the Cartesian basis vectors. A column matrix is denoted $\left(\underset{̃}{\bullet }\right)$ and matrices are denoted $\left(\underset{\_}{\bullet }\right)$ . A right superscript is used to index quantities belonging to a group and to denote sub-matrices of a matrix, for instance, for $\underset{\_}{a}$ and $\underset{̃}{b}$ , by ${\underset{\_}{a}}^{\left(\mathrm{mn}\right)}$ and ${\underset{̃}{b}}^{\left(m\right)}$ , respectively. The transpose operation ${\left(\bullet \right)}^{\mathrm{T}}$ of a second-order tensor is defined as $\mathbf{A}=A_{\mathrm{pq}}{\overrightarrow{e}}_p\otimes {\overrightarrow{e}}_q$ , ${\mathbf{A}}^{\mathrm{T}}=A_{\mathrm{qp}}{\overrightarrow{e}}_p\otimes {\overrightarrow{e}}_q$ . The first and second time derivatives are denoted $\left(\stackrel{.}{\bullet }\right)$ and $\left(\stackrel{··}{\bullet }\right)$ , respectively.

2. Dynamic computational homogenization

The full-scale problem is decomposed into two scales: the macroscopic (or coarse) scale, and the microscopic (or fine) scale, which may present complex topology and nonlinear, evolving material behavior. These scales are coupled by scale transition relations, which usually involve the application of a localization operator to set the macro-to-micro kinematic relations, and the Hill–Mandel condition, which sets the micro-to-macro homogenization relation through energy consistency across the scales. The dynamic multi-scale model is set up by extending the classical first-order computational homogenization framework to the transient dynamic case. In the following, the subscript “M” refers to a macroscopic quantity, whereas the subscript “m” refers to a microscopic quantity.

2.1. Separation of scales

Consider a general heterogeneous material with n microscopic constituents. The classical computational homogenization requires the typical size of each microscopic constituent $l_{\mathrm{m},k}$ to be much smaller than the shortest characteristic wavelength ${\mathrm{\Lambda }}_k$ in this microscopic constituent k for a given applied excitation, i.e.

$$l_{\mathrm{m},k}\ll {\mathrm{\Lambda }}_k,\forall k=1,\dots ,n.$$

(1)

This equation is generally known as the long-wavelength approximation. Under this assumption, the microscopic representative volume element (RVE) can be considered to behave quasistatically. Therefore, micro-inertial effects are negligible. Since a locally resonant elasto-acoustic metamaterial operating in a subwavelength regime is designed in such a way that it obeys the long-wavelength approximation, a more relaxed assumption has been proposed by Pham et al. [23]. Let ${\mathrm{\Lambda }}_{\mathrm{mat},i}$ and ${\mathrm{\Lambda }}_{\mathrm{het},j}$ be the shortest characteristic wavelength in the i -th and j -th constituents of the matrix (host structure) and heterogeneities (resonators), respectively. The relaxed separation of scales is then expressed as

$$\{\begin{array}{cc}\mathrm{Matrix}\left(\mathrm{long}\mathrm{wavelength}\mathrm{approximation}\right):l_{\mathrm{m},i}\ll {\mathrm{\Lambda }}_{\mathrm{mat},i},\hfill & \forall i=1,\dots ,n_{\mathrm{mat}},\hfill \\ \mathrm{Heterogeneities}:l_{\mathrm{m},j}\lesssim {\mathrm{\Lambda }}_{\mathrm{het},j},\hfill & \forall j=1,\dots ,n_{\mathrm{het}}.\hfill \end{array}$$

(2)

Here, $n_{\mathrm{mat}}$ and $n_{\mathrm{het}}$ are the number of microscopic matrix and heterogeneity constituents, respectively. Using this relaxed separation of scales, the matrix is still assumed to react instantaneously; however, wavelengths corresponding to the heterogeneities can now be of the same order as the typical size of these resonators, taking into account micro-inertia effects.

2.2. Macroscopic problem

At the macroscopic scale, in the absence of body forces, the system is governed by the momentum balance equation:

$${\overrightarrow{\nabla }}_{0\mathrm{M}}\cdot {\mathbf{P}}_{\mathrm{M}}^{\mathrm{T}}-{\stackrel{.}{\overrightarrow{p}}}_{0\mathrm{M}}=\overrightarrow{0},$$

(3)

where ${\overrightarrow{\nabla }}_{0\mathrm{M}}$ is the gradient operator with respect to the reference configuration, ${\mathbf{P}}_{\mathrm{M}}$ is the first Piola–Kirchhoff stress tensor, and ${\overrightarrow{p}}_{0\mathrm{M}}$ is the momentum. In computational homogenization, the constitutive behaviour is not known a priori, and it is obtained from the microscopic response at each macroscopic material point.

2.3. Microscopic problem

The locally resonant elasto-acoustic metamaterial considered in this work is a periodic structure. The choice of microscopic RVE is therefore straightforward and is identical to the unit cell (see Figure 1).

Figure 1.

2D representative volume element used for dynamic computational homogenization.

The displacements at the microscopic scale can be related to the macroscopic kinematics (i.e., position or displacement, velocity, and acceleration fields) by means of a Taylor expansion in terms of the position vector around a micro-scale reference position vector. Restricting the approximation to the first-order derivative, the resulting kinematic relation that should hold at each time $t\in \left[0,t_{\mathrm{end}}\right]$ , with $t_{\mathrm{end}}$ the end time of the simulation, yields

$${\overrightarrow{u}}_{\mathrm{m}}\left({\overrightarrow{X}}_{0\mathrm{m}},t\right)={\overrightarrow{u}}_{\mathrm{M}}\left(t\right)+\left({\overrightarrow{X}}_{0\mathrm{m}}-{\overrightarrow{X}}_{0\mathrm{m}}^{\mathrm{r}}\right)\cdot \left({\mathbf{F}}_{\mathrm{M}}^{\mathrm{T}}\left(t\right)-\mathbf{I}\right)+{\overrightarrow{w}}_{\mathrm{m}}\left({\overrightarrow{X}}_{0\mathrm{m}},t\right),\forall t\in \left[0,t_{\mathrm{end}}\right].$$

(4)

Here, ${\overrightarrow{u}}_{\mathrm{m}}$ and ${\overrightarrow{u}}_{\mathrm{M}}$ are the displacements at the microscopic and macroscopic scales, respectively. ${\overrightarrow{X}}_{0\mathrm{m}}$ is the microscopic position vector in the undeformed configuration and ${\overrightarrow{X}}_{0\mathrm{m}}^{\mathrm{r}}$ a micro-scale reference position vector, taken here to be the geometric center of the RVE; ${\mathbf{F}}_{\mathrm{M}}={\left({\overrightarrow{\nabla }}_{0\mathrm{M}}\otimes {\overrightarrow{x}}_{\mathrm{M}}\right)}^{\mathrm{T}}$ is the macroscopic deformation gradient tensor, with ${\overrightarrow{x}}_{\mathrm{M}}$ the macroscopic position vector in the current configuration, and ${\overrightarrow{w}}_{\mathrm{m}}$ are the microscopic fluctuations. These fluctuations provide the necessary degrees of freedom to account for the locally resonant behaviour. At the RVE level, the system is governed by the momentum balance equation, which, in the absence of body forces, reads

$${\overrightarrow{\nabla }}_{0\mathrm{m}}\cdot {\mathbf{P}}_{\mathrm{m}}^{\mathrm{T}}-{\stackrel{.}{\overrightarrow{p}}}_{0\mathrm{m}}=\overrightarrow{0}.$$

(5)

At the microscopic scale, the constitutive behaviour of each material constituent $\alpha$ (i.e. matrix, coating, or inclusion) is assumed to be known and described by constitutive laws, which may be time- or history-dependent:

$${\mathbf{P}}_{\mathrm{m}}^{\alpha }\left(t\right)={\mathcal{F}}^{\alpha }\left\{{\mathbf{F}}_{\mathrm{m}}^{\alpha }\left(\tau \right),{\stackrel{.}{\mathbf{F}}}_{\mathrm{m}}^{\alpha }\left(\tau \right),\tau \in \left[0,t\right]\right\},$$

(6)

$${\stackrel{.}{\overrightarrow{p}}}_{0\mathrm{m}}^{\alpha }\left(t\right)={\rho }_{0\mathrm{m}}^{\alpha }{\stackrel{··}{\overrightarrow{u}}}_{\mathrm{m}}^{\alpha }\left(t\right).$$

(7)

Note that the rate of change of the microscopic linear momentum is directly related to the acceleration ${\stackrel{··}{\overrightarrow{u}}}_{\mathrm{m}}^{\alpha }$ by the density ${\rho }_{0\mathrm{m}}^{\alpha }$ . Initial and boundary conditions are needed to solve the microscopic problem. These are derived by introducing the scale transition relations, as given in the following section.

2.4. Scale transition relations

2.4.1. Macro-to-micro: kinematics

Kinematic averaging is used to establish the downscaling relation from the macroscopic deformation gradient tensor ${\mathbf{F}}_{\mathrm{M}}$ , as

$${\mathbf{F}}_{\mathrm{M}}=\frac{1}{V_{0\mathrm{m}}}{\int }_{V_{0\mathrm{m}}}{\mathbf{F}}_{\mathrm{m}}\mathrm{d}{V}_{0\mathrm{m}}.$$

(8)

Here, $V_{0\mathrm{m}}$ is the microscopic volume in the reference configuration.

The microscopic deformation gradient tensor ${\mathbf{F}}_{\mathrm{m}}$ can be determined from the kinematics relation by taking the gradient of equation (4) with respect to the undeformed microscopic configuration, which yields

$${\mathbf{F}}_{\mathrm{m}}^{\mathrm{T}}={\mathbf{F}}_{\mathrm{M}}^{\mathrm{T}}+{\overrightarrow{\nabla }}_{0\mathrm{m}}\otimes {\overrightarrow{w}}_{\mathrm{m}}.$$

(9)

Substitution of this relation in equation (8) and application of the divergence theorem results in the following constraint equation on the microfluctuation field:

$${\int }_{{\mathrm{\Gamma }}_{0\mathrm{m}}}{\overrightarrow{n}}_{0\mathrm{m}}\otimes {\overrightarrow{w}}_{\mathrm{m}}\mathrm{d}{\mathrm{\Gamma }}_{0\mathrm{m}}=\mathbf{0},$$

(10)

where ${\mathrm{\Gamma }}_{0\mathrm{m}}$ denotes the undeformed microscopic boundary with outward normal ${\overrightarrow{n}}_{0\mathrm{m}}$ . Among different options, periodic boundary conditions are the most natural choice for periodic microstructures and have been shown to perform adequately, even for random microstructures in static homogenization [33]. Periodic boundary conditions require the microscopic fluctuations ${\overrightarrow{w}}_{\mathrm{m}}$ to be equal at the corresponding points at opposite boundaries with opposite normals (see Figure 1). These relations can be expressed as

$${\overrightarrow{w}}_{\mathrm{m}}\left({\overrightarrow{X}}_{0\mathrm{m}}^{+}\right)={\overrightarrow{w}}_{\mathrm{m}}\left({\overrightarrow{X}}_{0\mathrm{m}}^{-}\right),\forall {\overrightarrow{X}}_{0\mathrm{m}}^{+}\mathrm{and}{\overrightarrow{X}}_{0\mathrm{m}}^{-}\mathrm{such}\mathrm{that}{\overrightarrow{n}}_{0\mathrm{m}}^{+}=-{\overrightarrow{n}}_{0\mathrm{m}}^{-}.$$

(11)

Note, that for transient problems, care should be taken when applying periodic boundary conditions. For the considered RVE, the locally resonant behaviour is localized in the center and fully confined in the quasistatic matrix, which is in accordance with the relaxed scale separation condition (equation (2)). Hence, the use of periodic boundary conditions is allowed for the considered locally resonant elasto-acoustic metamaterial RVE.

2.4.2. Micro-to-macro: Hill–Mandel condition

The macroscopic stress, momentum and tangent stiffness can be determined from the microscopic scale based on the Hill–Mandel condition [34]. This requires the volume average of the microscopic virtual work in the RVE to equal the virtual work per unit of volume at the macro-scale. In a transient framework, both work of deformation and work of displacement should be taken into account, leading to

$${\mathbf{P}}_{\mathrm{M}}:\delta {\mathbf{F}}_{\mathrm{M}}^{\mathrm{T}}+{\stackrel{.}{\overrightarrow{p}}}_{0\mathrm{M}}\cdot \delta {\overrightarrow{u}}_{\mathrm{M}}=\frac{1}{V_{0\mathrm{m}}}{\int }_{V_{0\mathrm{m}}}{\mathbf{P}}_{\mathrm{m}}:\delta {\mathbf{F}}_{\mathrm{m}}^{\mathrm{T}}\mathrm{d}{V}_{0\mathrm{m}}+\frac{1}{V_{0\mathrm{m}}}{\int }_{V_{0\mathrm{m}}}{\stackrel{.}{\overrightarrow{p}}}_{0\mathrm{m}}\cdot \delta {\overrightarrow{u}}_{\mathrm{m}}\mathrm{d}{V}_{0\mathrm{m}}.$$

(12)

By making use of the virtual variations of the downscaling relations, equations (4) and (9), and after applying the chain rule, the right hand side of equation (12) can be written as

$$\begin{array}{c}\frac{1}{V_{0\mathrm{m}}}\left[\right({\int }_{V_{0\mathrm{m}}}{\mathbf{P}}_{\mathrm{m}}\mathrm{d}{V}_{0\mathrm{m}}):\delta {\mathbf{F}}_{\mathrm{M}}^{\mathrm{T}}+{\int }_{V_{0\mathrm{m}}}{\overrightarrow{\nabla }}_{0\mathrm{m}}\cdot ({\mathbf{P}}_{\mathrm{m}}^{\mathrm{T}}\cdot \delta {\overrightarrow{w}}_{\mathrm{m}})\mathrm{d}{V}_{0\mathrm{m}}+({\int }_{V_{0\mathrm{m}}}{\stackrel{.}{\overrightarrow{p}}}_{0\mathrm{m}}\mathrm{d}{V}_{0\mathrm{m}})\cdot \delta {\overrightarrow{u}}_{\mathrm{M}}\\ +\left({\int }_{V_{0\mathrm{m}}}{\stackrel{.}{\overrightarrow{p}}}_{0\mathrm{m}}\otimes \left({\overrightarrow{X}}_{0\mathrm{m}}-{\overrightarrow{X}}_{0\mathrm{m}}^{\mathrm{r}}\right)\mathrm{d}{V}_{0\mathrm{m}}\right):\delta {\mathbf{F}}_{\mathrm{M}}^{\mathrm{T}}+{\int }_{V_{0m}}\left({\stackrel{.}{\overrightarrow{p}}}_{0\mathrm{m}}-{\overrightarrow{\nabla }}_{0\mathrm{m}}\cdot {\mathbf{P}}_{\mathrm{m}}^{\mathrm{T}}\right)\cdot \delta {\overrightarrow{w}}_{\mathrm{m}}\mathrm{d}{V}_{0\mathrm{m}}].\end{array}$$

(13)

In this equation, the last term vanishes, owing to the momentum balance at the microscopic scale (cf. equation (5)). The second term can be rewritten using Gauss’s theorem, after which, the tractions on the boundary ${\overrightarrow{t}}_{0\mathrm{m}}={\overrightarrow{n}}_{0\mathrm{m}}\cdot {\mathbf{P}}_{\mathrm{m}}^{\mathrm{T}}$ can be substituted. A direct consequence of the periodic boundary conditions is the appearance of anti-periodic tractions [33]. By making use of equation (11), it is trivial to show that this term also vanishes (for detailed developments, see Kouznetsova [33]). The Hill–Mandel condition then reads

$$\begin{array}{c}{\mathbf{P}}_{\mathrm{M}}:\delta {\mathbf{F}}_{\mathrm{M}}^{\mathrm{T}}+{\stackrel{.}{\overrightarrow{p}}}_{0\mathrm{M}}\cdot \delta {\overrightarrow{u}}_{\mathrm{M}}\\ =\frac{1}{V_{0\mathrm{m}}}\left[\left({\int }_{V_{0\mathrm{m}}}\left({\mathbf{P}}_{\mathrm{m}}+{\stackrel{.}{\overrightarrow{p}}}_{0\mathrm{m}}\otimes \left({\overrightarrow{X}}_{0\mathrm{m}}-{\overrightarrow{X}}_{0\mathrm{m}}^{\mathrm{r}}\right)\right)\mathrm{d}{V}_{0\mathrm{m}}\right):\delta {\mathbf{F}}_{\mathrm{M}}^{\mathrm{T}}+\left({\int }_{V_{0\mathrm{m}}}{\stackrel{.}{\overrightarrow{p}}}_{0\mathrm{m}}\mathrm{d}{V}_{0\mathrm{m}}\cdot \delta {\overrightarrow{u}}_{\mathrm{M}}\right)\right].\end{array}$$

(14)

This condition must hold for any virtual deformation and displacement ( $\delta {\mathbf{F}}_{\mathrm{M}}^{\mathrm{T}}$ , $\delta {\overrightarrow{u}}_{\mathrm{M}}$ ). Hence,

$${\mathbf{P}}_{\mathrm{M}}=\frac{1}{V_{0\mathrm{m}}}{\int }_{V_{0\mathrm{m}}}\left({\mathbf{P}}_{\mathrm{m}}+{\stackrel{.}{\overrightarrow{p}}}_{0\mathrm{m}}\otimes \left({\overrightarrow{X}}_{0\mathrm{m}}-{\overrightarrow{X}}_{0\mathrm{m}}^{\mathrm{r}}\right)\right)\mathrm{d}{V}_{0\mathrm{m}},$$

(15)

$${\stackrel{.}{\overrightarrow{p}}}_{0\mathrm{M}}=\frac{1}{V_{0\mathrm{m}}}{\int }_{V_{0\mathrm{m}}}{\stackrel{.}{\overrightarrow{p}}}_{0\mathrm{m}}\mathrm{d}{V}_{0\mathrm{m}}.$$

(16)

These upscaling relations describe the dependency of the macroscopic stress and momentum on volume-averaged microscopic quantities. From a computational efficiency point of view, it is more convenient to rewrite the volume integrals as boundary integrals. This can be done by rewriting the microscopic stress using the microscopic equilibrium equation (equation (5)), and the identity tensor $\mathbf{I}={\overrightarrow{\nabla }}_{0\mathrm{m}}\otimes \left({\overrightarrow{X}}_{0\mathrm{m}}-{\overrightarrow{X}}_{0\mathrm{m}}^{\mathrm{r}}\right)$ . Application of the chain rule then leads to

$${\mathbf{P}}_{\mathrm{m}}={\overrightarrow{\nabla }}_{0\mathrm{m}}\cdot \left({\mathbf{P}}_{\mathrm{m}}^{\mathrm{T}}\otimes \left({\overrightarrow{X}}_{0\mathrm{m}}-{\overrightarrow{X}}_{0\mathrm{m}}^{\mathrm{r}}\right)\right)-{\stackrel{.}{\overrightarrow{p}}}_{0\mathrm{m}}\otimes \left({\overrightarrow{X}}_{0\mathrm{m}}-{\overrightarrow{X}}_{0\mathrm{m}}^{\mathrm{r}}\right).$$

(17)

Substituting equation (17) into equation (15), and applying Gauss’s theorem gives

$${\mathbf{P}}_{\mathrm{M}}^{\mathrm{T}}=\frac{1}{V_{0\mathrm{m}}}{\int }_{{\mathrm{\Gamma }}_{0\mathrm{m}}}\left({\overrightarrow{X}}_{0\mathrm{m}}-{\overrightarrow{X}}_{0\mathrm{m}}^{\mathrm{r}}\right)\otimes {\overrightarrow{t}}_{0\mathrm{m}}\mathrm{d}{\mathrm{\Gamma }}_{0\mathrm{m}}.$$

(18)

Similarly, by making use of equation (5) and the divergence theorem, the macroscopic momentum rate can be written as

$${\stackrel{.}{\overrightarrow{p}}}_{0\mathrm{M}}=\frac{1}{V_{0\mathrm{m}}}{\int }_{{\mathrm{\Gamma }}_{0\mathrm{m}}}{\overrightarrow{t}}_{0\mathrm{m}}\mathrm{d}{\mathrm{\Gamma }}_{0\mathrm{m}}.$$

(19)

The macroscopic stress and momentum can now be determined via integration over the RVE boundary only. The numerical implementation of this general dynamic homogenization framework is given in the next section.

3. Numerical implementation

3.1. Macroscopic system of equations

Adopting a standard finite-element discretization scheme, the macroscopic balance of momentum, equation (3), can be expressed in terms of its residual as follows:

$${\underset{̃}{\overrightarrow{r}}}_{\mathrm{M}}={\underset{̃}{\overrightarrow{f}}}_{\mathrm{M},\mathrm{inertia}}+{\underset{̃}{\overrightarrow{f}}}_{\mathrm{M},\mathrm{int}}-{\underset{̃}{\overrightarrow{f}}}_{\mathrm{M},\mathrm{ext}}=\underset{̃}{\overrightarrow{0}},$$

(20)

which vanishes if the balance of forces is satisfied; ${\underset{̃}{\overrightarrow{f}}}_{\mathrm{M},\mathrm{inertia}}$ , ${\underset{̃}{\overrightarrow{f}}}_{\mathrm{M},\mathrm{int}}$ , and ${\underset{̃}{\overrightarrow{f}}}_{\mathrm{M},\mathrm{ext}}$ are the macroscopic nodal inertia forces, internal forces, and external forces, respectively, and can be expressed as

$${\underset{̃}{\overrightarrow{f}}}_{\mathrm{M},\mathrm{inertia}}={\int }_{V_{0\mathrm{M}}^h}\underset{̃}{N}\left({\overrightarrow{X}}_{0\mathrm{M}}\right){\stackrel{.}{\overrightarrow{p}}}_{0\mathrm{M}}\mathrm{d}{V}_{0\mathrm{M}},$$

(21)

$${\underset{̃}{\overrightarrow{f}}}_{\mathrm{M},\mathrm{int}}={\int }_{V_{0\mathrm{M}}^h}{\overrightarrow{\nabla }}_{0\mathrm{M}}\underset{̃}{N}\left({\overrightarrow{X}}_{0\mathrm{M}}\right)\cdot {\mathbf{P}}_{\mathrm{M}}^{\mathrm{T}}\mathrm{d}{V}_{0\mathrm{M}},$$

(22)

$${\underset{̃}{\overrightarrow{f}}}_{\mathrm{M},\mathrm{ext}}={\int }_{{\mathrm{\Gamma }}_{0\mathrm{M}}^h}\underset{̃}{N}\left({\overrightarrow{X}}_{0\mathrm{M}}\right){\overrightarrow{t}}_{0\mathrm{M}}\mathrm{d}{\mathrm{\Gamma }}_{0\mathrm{M}},$$

(23)

where $\underset{̃}{N}$ are the shape functions and $V_{0\mathrm{M}}^h$ and ${\mathrm{\Gamma }}_{0\mathrm{M}}^h$ are the discretized undeformed macroscopic volume and boundary, respectively. The integrals involved in the expressions for these column vectors are computed numerically by applying a Gauss quadrature scheme. For this, the macroscopic stress ${\mathbf{P}}_{\mathrm{M}}$ and macroscopic momentum rate ${\stackrel{.}{\overrightarrow{p}}}_{0\mathrm{M}}$ need to be evaluated at the integration points. Since these quantities are obtained from the microscopic scale, an RVE simulation has to be performed at each macroscopic integration point.

3.2. Microscopic system of equations

The microscopic system of equations is similar to the macroscopic one, with the important difference that the constitutive equations are known at the micro-scale, i.e. equations (6) to (7). By making use of equation (7), the microscopic inertia forces can be expressed in terms of the mass matrix ${\underset{\_}{\mathbf{M}}}_{\mathrm{m}}$ and the acceleration column vector ${\stackrel{··}{\overrightarrow{\underset{̃}{u}}}}_{\mathrm{m}}$ , resulting in the microscopic residual given by

$${\underset{̃}{\overrightarrow{r}}}_{\mathrm{m}}\left({\overrightarrow{u}}_{\mathrm{m}},{\stackrel{.}{\overrightarrow{u}}}_{\mathrm{m}},{\stackrel{··}{\overrightarrow{u}}}_{\mathrm{m}},t\right)={\underset{\_}{\mathbf{M}}}_{\mathrm{m}}\cdot {\stackrel{··}{\overrightarrow{\underset{̃}{u}}}}_{\mathrm{m}}+{\underset{̃}{\overrightarrow{f}}}_{\mathrm{m},\mathrm{int}}\left({\overrightarrow{u}}_{\mathrm{m}},{\stackrel{.}{\overrightarrow{u}}}_{\mathrm{m}}\right)-{\underset{̃}{\overrightarrow{f}}}_{\mathrm{m},\mathrm{ext}}\left(t\right)=\underset{̃}{\overrightarrow{0}}.$$

(24)

3.3. Time integration

At both microscopic and macroscopic scales, time integration is performed according to the implicit Newmark scheme considering a fixed time step $\mathrm{\Delta t}$ . By choosing the integration parameters $\gamma =\frac{1}{2}$ and $\beta =\frac{1}{4}$ , the average constant acceleration algorithm is obtained and the integration is unconditionally stable [35]. Newmark’s time integration relations at time step $n+1$ can be written as

$${\underset{̃}{\stackrel{··}{\overrightarrow{u}}}}_{n+1}=\frac{1}{\beta {\left(\mathrm{\Delta t}\right)}^2}\left({\underset{̃}{\overrightarrow{u}}}_{n+1}-{\underset{̃}{\overrightarrow{u}}}_n\right)-\frac{1}{\mathrm{\beta }\mathrm{\Delta }t}{\underset{̃}{\stackrel{.}{\overrightarrow{u}}}}_n+\left(1-\frac{1}{2\beta }\right){\underset{̃}{\stackrel{··}{\overrightarrow{u}}}}_n,$$

(25)

$${\underset{̃}{\stackrel{.}{\overrightarrow{u}}}}_{n+1}=\frac{\gamma }{\mathrm{\beta }\mathrm{\Delta }t}\left({\underset{̃}{\overrightarrow{u}}}_{n+1}-{\underset{̃}{\overrightarrow{u}}}_n\right)+\left(1-\frac{\gamma }{\beta }\right){\underset{̃}{\stackrel{.}{\overrightarrow{u}}}}_n+\left(1-\frac{\gamma }{2\beta }\right)\mathrm{\Delta t}{\underset{̃}{\stackrel{··}{\overrightarrow{u}}}}_n.$$

(26)

After insertion of these relations, the macroscopic and microscopic residuals (equations (20) and (24), respectively) can be expressed in terms of the displacement ${\underset{̃}{\overrightarrow{u}}}_{n+1}$ at time step $n+1$ only, as follows:

$$\underset{̃}{\overrightarrow{r}}\left({\underset{̃}{\overrightarrow{u}}}_{n+1}\right)=\underset{̃}{\overrightarrow{0}}.$$

(27)

3.4. Linearization

The resulting set of nonlinear equations is solved iteratively using the Newton–Raphson procedure. Hence, the macroscopic and microscopic residuals should be linearized. Following standard linearization techniques, the residual at iteration $i+1$ can be expressed as

$$\underset{̃}{\overrightarrow{r}}\left({\underset{̃}{\overrightarrow{u}}}_{n+1}^{i+1}\right)=\underset{̃}{\overrightarrow{r}}\left({\underset{̃}{\overrightarrow{u}}}_{n+1}^i+\mathrm{d}\underset{̃}{\overrightarrow{u}}\right)\simeq \underset{̃}{\overrightarrow{r}}\left({\underset{̃}{\overrightarrow{u}}}_{n+1}^i\right)+{\underset{\_}{\mathbf{S}}}^i\cdot \mathrm{d}\underset{̃}{\overrightarrow{u}},\mathrm{with}{\underset{\_}{\mathbf{S}}}^i={\left[\frac{\partial \underset{̃}{\overrightarrow{r}}}{\partial \underset{̃}{\overrightarrow{u}}}\right]}_{{\underset{̃}{\overrightarrow{u}}}_{n+1}^i}.$$

(28)

The expressions for the iterative macroscopic and microscopic tangent matrices ${\underset{\_}{\mathbf{S}}}_{\mathrm{M}}^i$ and ${\underset{\_}{\mathbf{S}}}_{\mathrm{m}}^i$ are obtained by linearizing the corresponding residuals, equations (20) and (24), respectively. The iterative microscopic tangent matrix ${\underset{\_}{\mathbf{S}}}_{\mathrm{m}}^i$ is given by

$${\underset{\_}{\mathbf{S}}}_{\mathrm{m}}^i=\frac{1}{\beta {\left(\mathrm{\Delta t}\right)}^2}{\underset{\_}{\mathbf{M}}}_{\mathrm{m}}+\frac{\gamma }{\mathrm{\beta }\mathrm{\Delta }t}{\underset{\_}{\mathbf{C}}}_{\mathrm{m}}^i+{\underset{\_}{\mathbf{K}}}_{\mathrm{m}}^i,$$

(29)

with tangent stiffness matrix ${\underset{\_}{\mathbf{K}}}_{\mathrm{m}}^i$ and damping matrix ${\underset{\_}{\mathbf{C}}}_{\mathrm{m}}^i$ computed in a standard manner for given (nonlinear) micro-scale constitutive relations (equations (6) and (7)).

At the macroscopic scale, the iterative tangent matrix ${\underset{\_}{\mathbf{S}}}_{\mathrm{M}}^i$ can be found by considering the iterative update of the macroscopic residual (equation (20)):

$${\underset{\_}{\mathbf{S}}}_{\mathrm{M}}^i=\frac{1}{\mathrm{d}{\underset{̃}{\overrightarrow{u}}}_{\mathrm{M}}}\left[{\int }_{V_{0\mathrm{M}}^h}\underset{̃}{N}\left({\overrightarrow{X}}_{0\mathrm{M}}\right)\mathrm{d}{\stackrel{.}{\overrightarrow{p}}}_{0\mathrm{M}}\mathrm{d}{V}_{0\mathrm{M}}+{\int }_{V_{0\mathrm{M}}^h}{\overrightarrow{\nabla }}_{0\mathrm{M}}\underset{̃}{N}\left({\overrightarrow{X}}_{0\mathrm{M}}\right)\cdot \mathrm{d}{\mathbf{P}}_{\mathrm{M}}^{\mathrm{T}}\mathrm{d}{V}_{0\mathrm{M}}\right].$$

(30)

In the most general case, the iterative change of stress $\mathrm{d}{\mathbf{P}}_{\mathrm{M}}^{\mathrm{T}}$ and momentum rate $\mathrm{d}{\stackrel{.}{\overrightarrow{p}}}_{0\mathrm{M}}$ can be related to the changes in deformation $\mathrm{d}{\mathbf{F}}_{\mathrm{M}}^{\mathrm{T}}$ and rigid body displacements $\mathrm{d}{\overrightarrow{u}}_{\mathrm{M}}$ by

$$\mathrm{d}{\mathbf{P}}_{\mathrm{M}}^{\mathrm{T}}{=}^4{\mathbf{C}}_{\mathrm{M}}^{⟨1⟩}:\mathrm{d}{\mathbf{F}}_{\mathrm{M}}^{\mathrm{T}}{+}^3{\mathbf{C}}_{\mathrm{M}}^{⟨2⟩}\cdot \mathrm{d}{\overrightarrow{u}}_{\mathrm{M}},$$

(31)

$$\mathrm{d}{\stackrel{.}{\overrightarrow{p}}}_{0\mathrm{M}}{=}^3{\mathbf{C}}_{\mathrm{M}}^{⟨3⟩}:\mathrm{d}{\mathbf{F}}_{\mathrm{M}}^{\mathrm{T}}{+}^2{\mathbf{C}}_{\mathrm{M}}^{⟨4⟩}\cdot \mathrm{d}{\overrightarrow{u}}_{\mathrm{M}},$$

(32)

where ${{}^4\mathbf{C}}_{\mathrm{M}}^{⟨1⟩},{{}^3\mathbf{C}}_{\mathrm{M}}^{⟨2⟩},{{}^3\mathbf{C}}_{\mathrm{M}}^{⟨3⟩}$ , and ${{}^2\mathbf{C}}_{\mathrm{M}}^{⟨4⟩}$ are four macroscopic constitutive tangent operators. Note that, in a homogeneous material, the stress does not depend on the rigid body displacements and the change in momentum does not depend on the deformation, i.e. ${{}^3\mathbf{C}}_{\mathrm{M}}^{⟨2⟩}={{}^3\mathbf{C}}_{\mathrm{M}}^{⟨3⟩}={}^3\mathbf{0}$ . However, in the case of a locally resonant elasto-acoustic metamaterial, this is not necessarily the case, owing to the microstructural local resonance effects. In particular, for input frequencies inside the band gap, the macroscopic stress depends on rigid body displacements, owing to inertia effects. The macroscopic constitutive tangent operators are determined from the RVE calculation. The relations for ${{}^4\mathbf{C}}_{\mathrm{M}}^{⟨1⟩},{{}^3\mathbf{C}}_{\mathrm{M}}^{⟨2⟩},{{}^3\mathbf{C}}_{\mathrm{M}}^{⟨3⟩}$ , and ${{}^2\mathbf{C}}_{\mathrm{M}}^{⟨4⟩}$ are derived in an analogous manner as for static homogenization [33], as given at the end of this section.

At both the microscopic and macroscopic scales, the iterative updates on the displacements $d\underset{̃}{\overrightarrow{u}}$ can be computed from

$${\underset{\_}{\mathbf{S}}}^i\cdot \mathrm{d}\underset{̃}{\overrightarrow{u}}=-\underset{̃}{\overrightarrow{r}}\left({\underset{̃}{\overrightarrow{u}}}_{n+1}^i\right).$$

(33)

The Newton–Raphson iterative procedure is assumed to have reached convergence when the residual $\underset{̃}{\overrightarrow{r}}\left({\underset{̃}{\overrightarrow{u}}}_{n+1}^i\right)$ normalized with respect to the external forces value ${\underset{̃}{\overrightarrow{f}}}_{\mathrm{ext}}$ is below a certain tolerance value.

3.5. Boundary conditions

At the macroscopic scale, the boundary conditions can be applied using conventional finite-element techniques. At the microscopic scale, periodic boundary conditions, as introduced in equation (11), are applied. The application of the periodic boundary conditions is standard and has been discussed in detail in the context of the quasistatic computational homogenization, e.g. by Sridhar et al. [26]. Here, it follows the same lines, in accordance with the relaxed separation of scales.

Next, the coupling between the macroscopic kinematic quantities, ${\overrightarrow{u}}_{\mathrm{M}}$ and ${\mathbf{F}}_{\mathrm{M}}$ , and the microscopic displacement field is established. The macroscopic displacement describes the amplitude of the propagating wave, given by the overall rigid body motion of the matrix domain of the RVE. Indeed, Willis [36] formulated the macroscopic displacement as the volume average of the microscopic displacements taken over the matrix domain only. An equivalent formulation was proposed by Pham et al. [23], in which the average is taken over the entire RVE and subtracted by the dynamic microfluctuations. However, a much simpler formulation [26] that works equally well is to set the microfluctuations field ${\overrightarrow{w}}_{\mathrm{m}}$ to zero in some nodes p, hereafter called “prescribed” nodes, which will be used here. Thus, from equation (4), the displacements of the prescribed nodes are given by

$${\overrightarrow{u}}_{\mathrm{m}}^{\left(p\right)}={\overrightarrow{u}}_{\mathrm{M}}+\left({\overrightarrow{X}}_{0\mathrm{m}}^{\left(p\right)}-{\overrightarrow{X}}_{0\mathrm{m}}^{\mathrm{r}}\right)\cdot \left({\mathbf{F}}_{\mathrm{M}}^{\mathrm{T}}-\mathbf{I}\right),\mathrm{with}p=1,\dots ,n_p,$$

(34)

where $n_p$ is the number of prescribed nodes, which are typically the end points of a 1D domain ( $n_p=2$ ) or the minimum required number of vertices of a rectangular or parallelepiped-shaped unit cell in 2D ( $n_p=3$ ) or 3D ( $n_p=4$ ), respectively.

3.6. Initial conditions

To apply the Newmark integration algorithm, the initial displacement and velocity of all nodes should be specified. The initial acceleration of all nodes can then be determined by making use of [37]

$${}^a\underset{\_}{\mathbf{S}}\cdot \underset{̃}{\stackrel{··}{\overrightarrow{u}}}=-\underset{̃}{\overrightarrow{r}}\mathrm{with}{}^a\underset{\_}{\mathbf{S}}=\underset{\_}{\mathbf{M}}+\mathrm{\gamma \Delta t}\underset{\_}{\mathbf{C}}+\beta {\left(\mathrm{\Delta t}\right)}^2\underset{\_}{\mathbf{K}}.$$

(35)

By setting $\underset{̃}{\overrightarrow{r}}=\underset{̃}{\overrightarrow{0}}$ and applying this relation to the initial conditions only, the initial acceleration of the free nodes can be found as

$${\underset{̃}{\stackrel{··}{\overrightarrow{u}}}}_0^{\left(f\right)}=-({{}^a\underset{\_}{\mathbf{S}}}^{\left(\mathrm{ff}\right)}{)}^{-1}\cdot {{}^a\underset{\_}{\mathbf{S}}}^{\left(\mathrm{fp}\right)}\cdot {\underset{̃}{\stackrel{··}{\overrightarrow{u}}}}_0^{\left(p\right)},$$

(36)

in which p and f denote the prescribed and all the remaining “free” nodes, respectively. This relation is applied to both the macroscopic and microscopic scales. Note that at the macro-scale the mass matrix is not available in a closed form. Therefore, it is assumed that at $t=0$ the initial mass matrix ${\left({\underset{\_}{\mathbf{M}}}_{\mathrm{M}}\right)}_0$ can be obtained by adopting ${\left({\stackrel{.}{\overrightarrow{p}}}_{0\mathrm{M}}\right)}_0={\rho }_{0\mathrm{M}}{\left({\stackrel{··}{\overrightarrow{u}}}_{\mathrm{M}}\right)}_0$ with ${\rho }_{0\mathrm{M}}$ the averaged macroscopic density, obtained using the rule of mixtures of static densities of the microscopic constituents. The initial macroscopic stiffness matrix ${\left({\underset{\_}{\mathbf{K}}}_{\mathrm{M}}\right)}_0$ and damping matrix ${\left({\underset{\_}{\mathbf{C}}}_{\mathrm{M}}\right)}_0$ are determined by averaging the microscopic quantities in the undeformed configuration.

3.7. Macroscopic stress and momentum rate

The macroscopic stress ${\mathbf{P}}_{\mathrm{M}}$ and momentum rate ${\stackrel{.}{\overrightarrow{p}}}_{0\mathrm{M}}$ must be determined after each microscopic simulation according to equations (18) and (19), respectively. The surface integral in these relations can be simplified in the case of periodic boundary conditions by making use of the anti-periodicity of the tying forces in opposite nodes [33]. In the integral, the contributions of these forces cancel each other out, leaving only the external forces applied at the prescribed nodes, yielding

$${\mathbf{P}}_{\mathrm{M}}^{\mathrm{T}}=\frac{1}{V_{0\mathrm{m}}}\sum _{i=1}^{n_p}\left({\overrightarrow{X}}_{0\mathrm{m}}^{\left(i\right)}-{\overrightarrow{X}}_{0\mathrm{m}}^{\mathrm{r}}\right)\otimes {\overrightarrow{f}}_{\mathrm{ext}}^{\left(i\right)},$$

(37)

$${\stackrel{.}{\overrightarrow{p}}}_{0\mathrm{M}}=\frac{1}{V_{0\mathrm{m}}}\sum _{i=1}^{n_p}{\overrightarrow{f}}_{\mathrm{ext}}^{\left(i\right)}.$$

(38)

Here, ${\overrightarrow{X}}_{0\mathrm{m}}^{\left(i\right)}$ are the initial position vectors of the prescribed nodes $i=1,\dots ,n_p$ and ${\overrightarrow{f}}_{\mathrm{ext}}^{\left(i\right)}$ are the reaction forces on these nodes.

3.8. Macroscopic tangents

The relations describing the macroscopic tangent tensors ${{}^4\mathbf{C}}_{\mathrm{M}}^{⟨1⟩},{{}^3\mathbf{C}}_{\mathrm{M}}^{⟨2⟩},{{}^3\mathbf{C}}_{\mathrm{M}}^{⟨3⟩}$ , and ${{}^2\mathbf{C}}_{\mathrm{M}}^{⟨4⟩}$ , as defined in equations (31) and (32), can be derived in a similar way as for the static homogenization by condensation of the microscopic stiffness [33]. Consider the converged microscopic discrete system of equations from which the dependent degrees of freedom have been eliminated by applying periodic boundary conditions. The system can be further partitioned into parts corresponding to the iterative updates of the displacements of the prescribed nodes $\mathrm{d}{\underset{̃}{u}}_{\mathrm{m}}^{\left(p\right)}$ and into parts corresponding to nodes in which the displacement is left free $\mathrm{d}{\underset{̃}{u}}_{\mathrm{m}}^{\left(f\right)}$ , as follows:

$$\left[\begin{array}{cc}{\underset{\_}{S}}_{\mathrm{m}}^{\left(\mathrm{pp}\right)}\hfill & {\underset{\_}{S}}_{\mathrm{m}}^{\left(\mathrm{pf}\right)}\hfill \\ {\underset{\_}{S}}_{\mathrm{m}}^{\left(\mathrm{fp}\right)}\hfill & {\underset{\_}{S}}_{\mathrm{m}}^{\left(\mathrm{ff}\right)}\hfill \end{array}\right]\left[\begin{array}{c}\mathrm{d}{\underset{̃}{u}}_{\mathrm{m}}^{\left(p\right)}\hfill \\ \mathrm{d}{\underset{̃}{u}}_{\mathrm{m}}^{\left(f\right)}\hfill \end{array}\right]=\left[\begin{array}{c}\mathrm{d}{\underset{̃}{f}}_{\mathrm{ext}}\hfill \\ \underset{̃}{0}\hfill \end{array}\right].$$

(39)

Note that the iterative tangent matrix ${\underset{\_}{S}}_{\mathrm{m}}$ in this expression is calculated using equation (29) by taking the stiffness matrix ${\underset{\_}{K}}^i$ at the last iteration, where convergence has been reached. From equation (39), the update of the external forces in the prescribed nodes $\mathrm{d}{\underset{̃}{f}}_{\mathrm{ext}}$ can be expressed in terms of the update of the displacements $\mathrm{d}{\underset{̃}{u}}_{\mathrm{m}}^{\left(p\right)}$ in these nodes, as

$${\underset{\_}{Ŝ}}_{\mathrm{m}}d{\underset{̃}{u}}_{\mathrm{m}}^{\left(p\right)}=d{\underset{̃}{f}}_{\mathrm{ext}},$$

(40)

with

$${\underset{\_}{Ŝ}}_{\mathrm{m}}={\underset{\_}{S}}_{\mathrm{m}}^{\left(\mathrm{pp}\right)}-{\underset{\_}{S}}_{\mathrm{m}}^{\left(\mathrm{pf}\right)}{\left({\underset{\_}{S}}_{\mathrm{m}}^{\left(\mathrm{ff}\right)}\right)}^{-1}{\underset{\_}{S}}_{\mathrm{m}}^{\left(\mathrm{fp}\right)}.$$

(41)

Equation (40) can also be rewritten in a specific vector/tensor format:

$$\sum _{j=1}^{n_p}{\widehat{\mathbf{S}}}_{\mathrm{m}}^{\left(\mathrm{ij}\right)}\cdot \mathrm{d}{\overrightarrow{u}}_{\mathrm{m}}^{\left(j\right)}=\mathrm{d}{\overrightarrow{f}}_{\mathrm{ext}}^{\left(i\right)},\mathrm{with}i=1,\dots ,n_p.$$

(42)

The macroscopic constitutive tangents ${{}^4\mathbf{C}}_{\mathrm{M}}^{⟨1⟩}$ and ${{}^3\mathbf{C}}_{\mathrm{M}}^{⟨2⟩}$ relate the iterative correction of the macroscopic stress to the macrostructural deformation and displacement variations, respectively. The expression for these tangents is derived by taking the variation of equation (37). Substitution of equation (42) for $\mathrm{d}{\overrightarrow{f}}_{\mathrm{ext}}^{\left(i\right)}$ gives

$$\mathrm{d}{\mathbf{P}}_{\mathrm{M}}^{\mathrm{T}}=\frac{1}{V_{0\mathrm{m}}}\sum _{i=1}^{n_p}\sum _{j=1}^{n_p}\left({\overrightarrow{X}}_{0\mathrm{m}}^{\left(i\right)}-{\overrightarrow{X}}_{0\mathrm{m}}^{\mathrm{r}}\right)\otimes {\widehat{\mathbf{S}}}_{\mathrm{m}}^{\left(\mathrm{ij}\right)}\cdot \mathrm{d}{\overrightarrow{u}}_{\mathrm{m}}^{\left(j\right)}.$$

(43)

After substitution of the kinematic relation for the prescribed nodal displacements $\mathrm{d}{\overrightarrow{u}}_{\mathrm{m}}^{\left(j\right)}$ , equation (34), the macroscopic constitutive tangents ${{}^4\mathbf{C}}_{\mathrm{M}}^{⟨1⟩}$ and ${{}^3\mathbf{C}}_{\mathrm{M}}^{⟨2⟩}$ can be identified by comparing the resulting relation with equation (31), i.e.

$${{}^4\mathbf{C}}_{\mathrm{M}}^{⟨1⟩}=\frac{1}{V_{0\mathrm{m}}}\sum _{i=1}^{n_p}\sum _{j=1}^{n_p}\left({\overrightarrow{X}}_{0\mathrm{m}}^{\left(i\right)}-{\overrightarrow{X}}_{0\mathrm{m}}^{\mathrm{r}}\right)\otimes {\widehat{\mathbf{S}}}_{\mathrm{m}}^{\left(\mathrm{ij}\right)}\otimes \left({\overrightarrow{X}}_{0\mathrm{m}}^{\left(j\right)}-{\overrightarrow{X}}_{0\mathrm{m}}^{\mathrm{r}}\right),$$

(44)

$${{}^3\mathbf{C}}_{\mathrm{M}}^{⟨2⟩}=\frac{1}{V_{0\mathrm{m}}}\sum _{i=1}^{n_p}\sum _{j=1}^{n_p}\left({\overrightarrow{X}}_{0\mathrm{m}}^{\left(i\right)}-{\overrightarrow{X}}_{0\mathrm{m}}^{\mathrm{r}}\right)\otimes {\widehat{\mathbf{S}}}_{\mathrm{m}}^{\left(\mathrm{ij}\right)}.$$

(45)

The remaining macroscopic constitutive tangents, ${{}^3\mathbf{C}}_{\mathrm{M}}^{⟨3⟩}$ and ${{}^2\mathbf{C}}_{\mathrm{M}}^{⟨4⟩}$ , relate the iterative correction of the macroscopic rate of momentum to the macrostructural deformation and displacement variations, respectively. The relations for these tangents are derived by taking the iterative correction form of equation (38):

$$\mathrm{d}{\stackrel{.}{\overrightarrow{p}}}_{0\mathrm{M}}=\frac{1}{V_{0\mathrm{m}}}\sum _{i=1}^{n_p}\mathrm{d}{\overrightarrow{f}}_{\mathrm{ext}}^{\left(i\right)}.$$

(46)

Analogously to these derivations, the macroscopic constitutive tangents, ${{}^3\mathbf{C}}_{\mathrm{M}}^{⟨3⟩}$ and ${{}^2\mathbf{C}}_{\mathrm{M}}^{⟨4⟩}$ , can be derived after substitution of equation (42) into equation (46):

$${{}^3\mathbf{C}}_{\mathrm{M}}^{⟨3⟩}=\frac{1}{V_{0\mathrm{m}}}\sum _{i=1}^{n_p}\sum _{j=1}^{n_p}{\widehat{\mathbf{S}}}_{\mathrm{m}}^{\left(\mathrm{ij}\right)}\otimes \left({\overrightarrow{X}}_{0\mathrm{m}}^{\left(j\right)}-{\overrightarrow{X}}_{0\mathrm{m}}^{\mathrm{r}}\right),$$

(47)

$${{}^2\mathbf{C}}_{\mathrm{M}}^{⟨4⟩}=\frac{1}{V_{0\mathrm{m}}}\sum _{i=1}^{n_p}\sum _{j=1}^{n_p}{\widehat{\mathbf{S}}}_{\mathrm{m}}^{\left(\mathrm{ij}\right)}.$$

(48)

3.9. Numerical solution procedure

The numerical implementation of the presented dynamic computational homogenization scheme leads to the solution procedure sketched in the flowchart in Figure 2. The derivation of this framework differs from previous work [23, 26] in the introduction of the iterative macroscopic tangent matrix ${\underset{\_}{\mathbf{S}}}_{\mathrm{M}}^i$ , for which the expressions of the macroscopic tangent tensors ${{}^4\mathbf{C}}_{\mathrm{M}}^{⟨1⟩},{{}^3\mathbf{C}}_{\mathrm{M}}^{⟨2⟩},{{}^3\mathbf{C}}_{\mathrm{M}}^{⟨3⟩}$ , and ${{}^2\mathbf{C}}_{\mathrm{M}}^{⟨4⟩}$ have been derived, see Sections 3.4 and 3.8.

Figure 2.

Numerical implementation of the dynamic computational homogenization model.

4. Numerical results

In this section, the dynamic computational homogenization framework is verified by comparing the results with those from direct numerical simulations.

4.1. Mechanical model system

A one-dimensional version, i.e. each node supports only one degree of freedom, being the horizontal displacement, of Liu’s locally resonant metamaterial design [4] is used to highlight the relevance of the computational homogenization framework. Although only longitudinal motions are allowed in this simplified setting, it presents both dispersion and nonlinear effects, making it sufficiently complex to assess the proposed homogenization framework.

The considered 1D unit cell (or RVE) is depicted in Figure 3(b), along with a cross-section of Liu’s unit cell (Figure 3(a)), highlighting the similarity between these models. The 1D RVE consists of two parallel parts rigidly connected at interface nodes. The lower part is composed of three material phases: an epoxy matrix (green), a silicone rubber coating (gray), and a lead core (red); the upper part consists of epoxy only. The contrast between the material properties in the lower part of the unit cell, i.e. the low stiffness of the rubber coating and the heavy mass of the lead core, ensures localized motions, which triggers local resonance effects [15]. The upper part, consisting only of the matrix material, captures the overall elastic stiffness of the locally resonant elasto-acoustic metamaterial and provides the path for wave propagation. The geometric and linear elastic material properties of the 1D RVE correspond to those considered by Sridhar et al. [26] and are given in Table 1. In addition, a small stiffness-proportional damping of the form ${\underset{\_}{C}}_{\mathrm{m}}=b{\underset{\_}{K}}_{\mathrm{m}}$ , with coefficient $b=2.242\times 10^{-6}\mathrm{s}$ , is added to the system to ensure numerical stability. The unit cell is discretized by 55 linear finite elements, as can be seen in Figure 3(b). Since low stiffness implies a lower wave speed and thus shorter wavelengths at a given frequency, a fine discretization is needed in the rubber phase, owing to its high compliance relative to the matrix phase — the element size in the rubber phase is $0.25\mathrm{mm}$ .

Figure 3.

Unit cells of locally resonant metamaterial models: (a) cross-section of Liu’s model; (b) 1D model used for the numerical simulations, including the finite element discretization. Squares indicate nodes. Both models consist of three material phases: epoxy matrix (green), silicone rubber coating (gray), and lead core (red).

Table 1.

Geometric and linear elastic material properties of the studied unit cell.

Geometrical properties			Material properties		Epoxy	Rubber	Lead
Parameter		Value, mm	Parameter
Length of RVE	$l_{\mathrm{m}}$	25.2	Density	$\rho ,{\mathrm{kg}/\mathrm{m}}^3$	$1.18\times {10}^3$	$1.30\times {10}^3$	$1.16\times {10}^4$
Thickness of rubber coating	$t_{\mathrm{c}}$	2.5	Young modulus	$E,{\mathrm{N}/\mathrm{m}}^2$	$4.35\times {10}^9$	$1.18\times {10}^5$	$4.08\times {10}^{10}$
Diameter of the lead core	d	10.0

4.1.1. Linear elastic material

For a linear elastic material behavior of all constituents, with the material properties given in Table 1, the dispersion diagram for the considered 1D unit cell is depicted in Figure 4, which has been obtained by assuming time-harmonic solutions and applying Bloch–Floquet conditions to the left and right boundary nodes [38]. It predicts a locally resonant band gap for frequencies $f_{\omega }=\omega ∕2\pi$ between $140\mathrm{Hz}$ and $310\mathrm{Hz}$ . This band gap is associated with the local resonance of the rubber-coated lead core whose natural frequency is $f_{\mathrm{het}}\approx 140\mathrm{Hz}$ . In the vicinity of the band gap, the material exhibits dispersive behavior, i.e. a nonlinear relation between the wavenumber k and the frequency $f_{\omega }$ .

Figure 4.

Dispersion diagram for 1D unit cell with the locally resonant band gap depicted by the gray shaded area: (a) real part of the dimensionless wavenumber $\Re \left(k{l}_{\mathrm{m}}\right)$ indicating the propagating part of waves; (b) imaginary part of the wavenumber $\Im \left(k\right)$ indicating the evanescent part of waves.

The homogenizability of this unit cell is in accordance with the relaxed separation of scales conditions, equation (2), as can be verified by comparing the local resonance frequency $f_{\mathrm{het}}$ with the Bragg frequency $f_{\mathrm{Bragg}}\approx c_{\mathrm{mat}}∕2l_{\mathrm{m}}$ , where $c_{\mathrm{mat}}=\sqrt{E_{\mathrm{mat}}∕{\rho }_{\mathrm{mat}}}$ is the wave speed in the matrix material and the macroscopic wavelength ${\mathrm{\Lambda }}_{\mathrm{mat}}$ is twice the unit cell size ( ${\mathrm{\Lambda }}_{\mathrm{mat}}=2l_{\mathrm{m}}$ ), indicating the limit of scale separation. For the given unit cell, $f_{\mathrm{Bragg}}=38\mathrm{kHz}\gg f_{\mathrm{het}}=140\mathrm{Hz}$ , confirming the validity of the relaxed scale separation.

4.1.2. Nonlinear material

The focus of this work is on the homogenization analysis of nonlinear locally resonant elasto-acoustic metamaterials, motivated by the fact that, for sufficiently large input excitations, especially close to the local resonance frequency, the rubber phase experiences finite strains, owing to large central mass displacements. Accordingly, the linear elastic model of the rubber must be replaced by a hyperelastic, 1D incompressible neo-Hookean constitutive model, which can be expressed as [39]

$$P=2C_{10}\left(\lambda -\frac{1}{{\lambda }^2}\right),$$

(49)

where P is the nominal (or engineering) stress, $\lambda =e+1$ is the stretch ratio, with e being the linear strain, and $C_{10}$ is a material parameter. As the neo-Hookean constitutive model recovers the linear elastic behavior for infinitesimal strains, i.e. $e\ll 1$ , a relation between $C_{10}$ and the Young modulus can be retrieved, given by $E=6C_{10}$ . The stress–strain curve for this material model is plotted together with the equivalent linear elastic model in Figure 5. Note that this constitutive model shows an asymptotic behavior at $e\to -1$ (or $\lambda \to 0$ ), which yields an infinite compressive stress and tangent stiffness. Thus, caution is required in numerically solving the nonlinear algebraic equations with the Newton–Raphson solver in the regime $e<-0.5$ .

Figure 5.

Stress–strain curve for the neo-Hookean constitutive model and the equivalent linear elastic constitutive model used to model the rubber phase.

The numerical problem to be solved consists of a structure made up of $n_{\mathrm{uc}}$ unit cells arranged in the longitudinal direction, subjected to a harmonic prescribed displacement at the left boundary node and fixed at the right end, as shown in Figure 6, and expressed as follows:

$$\mathrm{Boundary}\mathrm{conditions}:u^{\left(0\right)}\left(t\right)=Asin\left(\mathrm{\omega t}\right),u^{\left(l_{\mathrm{M}}\right)}\left(t\right)=0,\forall t\in \left[0,t_{\mathrm{end}}\right],$$

(50)

where A is the amplitude of the prescribed harmonic displacement and $\omega =2\pi f_{\omega }$ is the angular frequency of the excitation. Simulations are performed for several excitation frequencies between $50\mathrm{Hz}$ and 500 $\mathrm{Hz}$ . The initial conditions for the prescribed end nodes follow directly from equation (50), while for the remaining “free” nodes it is assumed that the system is initially at rest, i.e. ${\underset{̃}{u}}^{\left(f\right)}\left(0\right)=\underset{̃}{0}$ and ${\stackrel{.}{\underset{̃}{u}}}^{\left(f\right)}\left(0\right)=\underset{̃}{0}$ .

Figure 6.

The structural problem to be solved numerically: (a) direct numerical simulation model system; (b) computational homogenization macroscopic model systems for two different homogenization levels (h), with squares indicating nodes of quadratic finite elements and crosses indicating integration points. At each integration point, a microscopic unit cell problem is solved.

To assess the potential of the nonlinear dynamic computational homogenization, two types of simulations are carried out, i.e. wave propagation and structural dynamic analyses. In the first case, the considered structure is taken to be long enough to avoid reflections at the boundaries, enabling wave propagation analysis. In the second case, a finite structure is considered and the transient behavior after successive reflections is examined.

Both computational homogenization (CH) and direct numerical simulations (DNSs) are performed. The results from DNSs are used to assess the computational homogenization framework. In the case of DNSs, the fully discretized model is solved, as depicted in Figure 6(a). The computational homogenization simulations consist of two scales, i.e. the microscopic and macroscopic scales. At the microscopic scale, the discretized unit cell shown in Figure 3(b) is used. At the macroscopic scale, the homogenized structure with the same length $l_{\mathrm{M}}$ as the DNS model system is considered, but discretized coarsely using 1D quadratic finite elements with the three-point Gaussian quadrature for numerical integration, as shown in Figure 6(b). At each integration point of a macroscopic finite element, a microscopic unit cell problem is solved. To investigate the efficiency against the accuracy of the computational homogenization scheme, the homogenization level (h) is defined as the ratio between the total number of unit cells in the DNS model system ( $n_{\mathrm{uc}}$ ) and the total number of integration points (i.e. the number of unit cell simulations) in the macroscopic model system ( $s_{\mathrm{M}}$ ):

$$h=\frac{n_{\mathrm{uc}}}{s_{\mathrm{M}}}.$$

(51)

Thus, as h increases, the number of integration points decreases and the macroscopic mesh discretization coarsens. As a result, a reduction in computational time is expected.

The implicit Newmark time integration scheme is employed with the same settings for the multiscale and DNSs. Relatively small time steps of $\mathrm{\Delta t}=2\mathrm{\pi \omega }∕900$ , for the wave propagation analysis, and $\mathrm{\Delta t}=2\mathrm{\pi \omega }∕600$ , for the transient structural dynamic analysis, depending on the excitation angular frequency $\omega$ , are used to ensure sufficient resolution of the higher-order harmonics generated as a result of the nonlinearity.

4.2. Wave propagation analysis

For the wave propagation analysis, a locally resonant metamaterial structure consisting of $n_{\mathrm{uc}}=3360$ unit cells is considered. Thus, the total length of the structure is $l_{\mathrm{M}}=n_{\mathrm{uc}}{l}_{\mathrm{m}}\simeq 84.7\mathrm{m}$ . To avoid reflections at the boundaries, the simulations are stopped as soon as the wave reaches the right boundary. Thus, the total simulation time ( $t_{\mathrm{end}}$ ) varies with the excitation frequency.

The wave attenuation properties of locally resonant structures are usually analyzed using a transmission plot, i.e. a curve representing the relation between the transmissibility measured at a specific position and the excitation frequency. The transmissibility $\theta$ is here defined as the ratio between the output and input root-mean-square displacements, at positions $x_{\mathrm{out}}$ and $x_{\mathrm{in}}$ , respectively, as follows:

$$\theta =\frac{{ū}^{\left(x_{\mathrm{out}}\right)}}{{ū}^{\left(x_{\mathrm{in}}\right)}},$$

(52)

where the root-mean-square displacement ${ū}^{\left(x_p\right)}$ , computed over the second half of the total time interval (to exclude transient effects), is given by

$${ū}^{\left(x_p\right)}=\sqrt{\frac{1}{n_{\mathrm{\Delta t}}}\sum _{t\in \left[\frac{t_{\mathrm{end}}}{2},t_{\mathrm{end}}\right]}{\left(u^{\left(x_p\right)}\left(t\right)\right)}^2},\mathrm{with}p=″\mathrm{in}″\mathrm{or}″\mathrm{out}″,$$

and $n_{\mathrm{\Delta t}}$ the number of computed time instances within the considered time window. In this work, the input displacement is taken at $x_{\mathrm{in}}=0$ , where the harmonic displacement is prescribed, and the output response is measured at $x_{\mathrm{out}}=\frac{l_{\mathrm{M}}}{4}$ which is chosen to be far enough from the source to measure the far-field response, but close enough to account for nonlinear effects and to consider enough periods before the wave reaches the right boundary.

In Figure 7, the transmissibility curves obtained by DNSs for the linear and nonlinear metamaterials are compared. In the linear case, the transmissibility is frequency-dependent, but invariant with amplitude. The local resonance feature of the metamaterial promotes the attenuation zone observed between 100 and 300 Hz. When nonlinear local interactions of the neo-Hookean type are incorporated, the transmissibility becomes amplitude-dependent. The responses of the nonlinear metamaterial for two excitation amplitudes, $A=0.5\mathrm{mm}$ and $A=1.0\mathrm{mm}$ , are shown in Figure 7. The tension-compression asymmetry of the neo-Hookean model, with more pronounced stiffening in compression than softening in tension (see Figure 5), explains the shift of the attenuation zone toward higher frequencies with increasing amplitude.

Figure 7.

Transmissibility of the metamaterial structure based on the root-mean-square displacement at $x_{\mathrm{out}}=\frac{l_{\mathrm{M}}}{4}$ for both linear and nonlinear cases at two different excitation amplitudes computed via direct numerical simulations.

Transmissibility curves obtained by DNSs and computational homogenization for various homogenization levels are depicted in Figure 8. For the highest amplitude of excitation considered, $A=1.0\mathrm{mm}$ , which entails a significantly nonlinear response of the rubber material, a very good agreement between computational homogenization simulations and direct numerical simulations can be observed for all levels of homogenization outside the band gap region. Moreover, the band gap emerging from the locally resonant structure is captured well by all computational homogenization simulations. Computational homogenization shows very good accuracy compared with DNSs for homogenization levels h as large as 13.3. In general, the greater the amount of homogenization, the coarser the discretization and its ability to resolve the dynamics at higher frequencies decreases. The results in Figure 8 show that, as expected, the accuracy reduces with the increase in the amount of homogenization.

Figure 8.

Transmissibility of the metamaterial structure based on the macro-scale root-mean-square displacements at $x_{\mathrm{out}}=\frac{l_{\mathrm{M}}}{4}$ for the nonlinear case computed via direct numerical simulations (DNSs) and computational homogenization (CH) for several levels of homogenization h.

4.3. Transient structural dynamic analysis

In contrast with the wave propagation analysis, in practice often finite-size structures need to be analyzed where the effect of wave reflection at external boundaries and interference cannot be neglected. In this section, the applicability and efficiency of the computational homogenization framework is investigated by performing a transient analysis of a finite locally resonant elasto-acoustic metamaterial structure as defined previously, with a total length $l_{\mathrm{M}}=1120l_{\mathrm{m}}\simeq 28.2\mathrm{m}$ . The left boundary node of the metamaterial is subjected to transient excitations in the form of a sine-Gaussian pulse with a given central frequency and amplitude $A=1.0\mathrm{mm}$ .

For a pulse with central frequency $f_{\omega }=100\mathrm{Hz}$ , the time signature and frequency spectra of the applied excitation are shown in Figure 9. The transient displacement responses measured at $x_{\mathrm{meas}}=\frac{l_{\mathrm{M}}}{4}$ obtained via the computational homogenization at different levels of homogenization are compared with the response of the DNSs in Figure 10. For all levels of homogenization, the overall transient behavior is well captured. The accuracy of the computational homogenization framework compared with DNSs in space and time is quantified by computing error indicators ${\epsilon }_x$ and ${\epsilon }_t$ , respectively, defined as

$${\epsilon }_x=\frac{\sqrt{\sum _x{\left\{u_{\mathrm{CH}}^{\left(x\right)}\left(t_{\mathrm{end}}\right)-u_{\mathrm{DNS}}^{\left(x\right)}\left(t_{\mathrm{end}}\right)\right\}}^2}}{\sqrt{\sum _x{\left\{u_{\mathrm{DNS}}^{\left(x\right)}\left(t_{\mathrm{end}}\right)\right\}}^2}}\cdot 100\left[\%\right],{\epsilon }_t=\frac{\sqrt{\sum _t{\left\{u_{\mathrm{CH}}^{\left(x_{\mathrm{meas}}\right)}\left(t\right)-u_{\mathrm{DNS}}^{\left(x_{\mathrm{meas}}\right)}\left(t\right)\right\}}^2}}{\sqrt{\sum _t{\left\{u_{\mathrm{DNS}}^{\left(x_{\mathrm{meas}}\right)}\left(t\right)\right\}}^2}}\cdot 100\left[\%\right],$$

(53)

with $x_{\mathrm{meas}}=l_{\mathrm{M}}∕4$ , as before. For the considered simulation time, a good accuracy is obtained for $h=6.7$ ( ${\epsilon }_t=0.25$ %) and $h=13.3$ ( ${\epsilon }_t=11.3$ %), while for the higher homogenization level, $h=23.3$ , the error is large ( ${\epsilon }_t=91.5$ %).

Figure 9.

Transient input displacement signal with central frequency $f_{\omega }=100\mathrm{Hz}$ applied to the nonlinear metamaterial: (a) in time domain, (b) in frequency domain.

Figure 10.

Temporal displacement response of the considered nonlinear metamaterial structure with length $l_{\mathrm{M}}\simeq 28.2\mathrm{m}$ computed via direct numerical simulations (DNSs) and computational homogenization (CH) for different homogenization levels h, excitation $f_{\omega }=100\mathrm{Hz}$ , measured at $x_{\mathrm{meas}}=\frac{l_{\mathrm{M}}}{4}$ .

The dynamic response of the whole structure obtained for different homogenization levels can be compared for a particular time instance. In Figure 11(a) and (b), the dynamic responses at $t_0=0.05\mathrm{s}$ and $t_{\mathrm{end}}=0.5\mathrm{s}$ are shown, respectively. The response at $t_0=0.05\mathrm{s}$ reveals the effect of dispersion in the pulse propagation, which is expected to be strong close to the local resonance frequency, as shown in the dispersion diagram for the linear metamaterial in Figure 4. At this time instance, the pulse has not yet reached the end of the structure, and all homogenization levels capture well the structural dynamic response. As expected, for longer simulation times, in particular, at the end of the simulation, at $t=0.5\mathrm{s}$ (Figure 11(b)), the results for lower homogenization levels reveal an error that remains below 10 % ( ${\epsilon }_x=0.17$ % for $h=6.7$ and ${\epsilon }_x=7.78$ % for $h=13.3$ ), while for the higher homogenization level $h=23.3$ , the error is significantly larger ( ${\epsilon }_x=64.4$ %). The observed deviations for coarse macroscopic discretization after long simulation times can be explained by the accumulation in time of the small errors in resolving the wave, increasing with reflection at the boundary and the interference phenomenon, leading to noticeable phase deviations, as shown in the enlarged part of Figure 10 and in Figure 11(b).

Figure 11.

Displacement response along the considered nonlinear metamaterial structure with length $l_{\mathrm{M}}\simeq 28.2\mathrm{m}$ computed via direct numerical simulations (DNSs) and computational homogenization (CH) for different homogenization levels h, excitation $f_{\omega }=100\mathrm{Hz}$ , measured at: (a) $t_0=0.05\mathrm{s}$ , (b) $t_{\mathrm{end}}=0.5\mathrm{s}$ .

Similarly, a pulse displacement with central frequency $f_{\omega }=250\mathrm{Hz}$ is applied to the left end boundary of the nonlinear metamaterial structure. The time signature and frequency spectra of this excitation are shown in Figure 12. Note that the frequency content of this pulse is mostly within the band gap zone depicted, see Figures 4 and 8. Thus, a significant attenuation during the pulse propagation through the metamaterial is expected. Indeed, the dynamic responses of the structure at two time instances, $t_0=0.01\mathrm{s}$ , at the beginning of the simulation, and $t_{\mathrm{end}}=0.5\mathrm{s}$ , at the end of the simulation, are shown in Figure 13. Note that the amplitude in both cases is much smaller than the amplitude of the excitation pulse $A=1\mathrm{mm}$ . Regarding the accuracy of the homogenized solutions, as in the previous case for $f_{\omega }=100\mathrm{Hz}$ , adequate agreement is obtained for all levels of homogenization before reflection occurs, i.e. at $t_0=0.01\mathrm{s}$ (Figure 13(a)), while deviations accumulate for longer simulation time, at $t_{\mathrm{end}}=0.5\mathrm{s}$ (Figure 13(b)), especially for higher homogenization level $h=23.3$ .

Figure 12.

Transient input displacement signal with central frequency $f_{\omega }=250\mathrm{Hz}$ applied to the nonlinear metamaterial: (a) in time domain, (b) in frequency domain.

Figure 13.

Displacement response along the considered nonlinear metamaterial structure with length $l_{\mathrm{M}}\simeq 28.2\mathrm{m}$ computed via direct numerical simulations (DNSs) and computational homogenization (CH) for different homogenization levels h, excitation $f_{\omega }=250\mathrm{Hz}$ , measured at: (a) $t_0=0.01\mathrm{s}$ , (b) $t_{\mathrm{end}}=0.5\mathrm{s}$ .

In Figure 14, the transient response of the nonlinear metamaterial measured at $x_{\mathrm{meas}}=\frac{l_{\mathrm{M}}}{4}$ is shown. As expected, the response is considerably attenuated with time, as a result of the strong attenuation along the metamaterial, which reduces the amplitude of the propagating and, thus, reflected waves. Although the response of the structure for all levels of homogenization shows the same characteristic attenuation behavior, the results reveal that, as expected, the accuracy deteriorates with the increase of the amount of homogenization. Moreover, from the enlarged detail in Figure 14, one can notice that a very good agreement is achieved at the beginning of the simulation. Afterwards, deviations from the direct numerical simulations appear in the homogenized solutions, being more pronounced for the higher homogenization levels. Indeed, these deviations exhibit spurious high-frequency content, owing to numerical dispersion. To improve the accuracy of the computational homogenization, advanced numerical schemes for spatial and temporal discretization should be considered. For the macro-scale discretization, isogeometric methods, e.g. based on B-spline shape functions, promise remarkably low numerical dispersion even at low wavelengths [40]. For the time integration, high-frequency dissipative time-stepping algorithms for nonlinear elastodynamics [41–43] would include numerical dissipation to damp out spurious high-frequency contributions while being unconditionally stable.

Figure 14.

Temporal displacement response of the considered nonlinear metamaterial structure with length $l_{\mathrm{M}}\simeq 28.2\mathrm{m}$ computed via direct numerical simulations (DNSs) and computational homogenization (CH) for different homogenization levels h, excitation $f_{\omega }=250\mathrm{Hz}$ , measured at $x_{\mathrm{meas}}=\frac{l_{\mathrm{M}}}{4}$ .

The main drawback of the transient structural dynamic analyses via direct numerical simulations is the excessive computational time, especially when fine microstructures with nonlinear features are being considered. For the present transient test cases, the CPU times required for the direct numerical simulation are compared with the computational homogenization simulations in Table 2. In both test cases, for homogenization levels $h=6.7$ and $h=13.3$ , CPU time reductions of about $20\%$ and $60\%$ , respectively, are obtained while retaining sufficient accuracy, as discussed previously. Therefore, in less than half the time, computational homogenization is able to perform accurate transient dynamic analysis of a nonlinear locally resonant elasto-acoustic metamaterial while even more gain can be achieved if a somewhat larger error is tolerated (e.g. $h=23.3$ ). These results demonstrate the capability of the computational homogenization framework presented and implemented in this work to perform transient dynamic simulations for structures involving material nonlinearities at a significantly reduced cost compared with direct numerical simulations.

Table 2.

Computational time required to perform transient structural dynamic analysis in the considered nonlinear locally resonant metamaterial structure of finite size.

Computational homogenization, h	CPU time reduction w.r.t. DNSs, %
	$f_{\omega }=100\mathrm{Hz}$	$f_{\omega }=250\mathrm{Hz}$
6.7	18	24
13.3	57	59
23.3	76	76

DNS: direct numerical simulation.

5. Conclusions

In this paper, the computational homogenization framework for solving fully transient multi-scale problems proposed by Pham et al. [23] has been extended and applied to the transient dynamic analysis of a nonlinear locally resonant acoustic metamaterial. The main focus of this contribution was on the incorporation of microscopic nonlinearities in a transient computational homogenization framework. In the present case, the macroscopic tangent constitutive tensors are no longer constant, and nested sets of nonlinear equations at both scales need to be solved iteratively in combination with an implicit time integration scheme. The expressions for the macroscopic constitutive tangents have been derived as functions of the microscopic iterative tangent matrix and used to determine the iterative tangent matrix at the macroscopic scale.

Transient numerical analyses of a 1D version of Liu’s locally resonant metamaterial have been performed to evaluate the proposed computational homogenization framework and its numerical implementation. The periodic unit cell included rubber-coated inclusions, in which the rubber phase was described by the neo-Hookean material model. The performance of the computational homogenization scheme in capturing the dynamics of this nonlinear metamaterial was assessed in two situations, i.e.: (i) under free wave propagation in an infinitely long structure and (ii) transient structural dynamics of a finite-size structure. In both cases, the results were compared with those obtained from direct numerical simulations.

In the case of wave propagation analysis, for all considered macroscopic discretizations, the computational homogenization framework was able to adequately capture the local resonant band gap characteristics of the considered metamaterial. The computational homogenization showed quite good accuracy for all homogenization levels. The transient analysis of a finite-size nonlinear metamaterial demonstrated the potential of the computational homogenization framework as an efficient alternative to the conventional direct numerical simulations in nonlinear structural dynamic analysis. Accurate results were obtained with significant computational time savings. In both numerical studies, the accuracy of the computational homogenization framework deteriorated for higher levels of homogenization, as expected, but very good agreement was obtained for homogenization as high as $h=13.3$ , reducing the computational time by more than half, compared with direct numerical simulation. The ability of the computational homogenization scheme to capture the dynamics of the metamaterial structure at frequencies outside, as well as within, the local resonance band gap was demonstrated. All levels of homogenization were able to predict the highly attenuated response of the structure. The strong attenuation reveals spurious high-frequency content, to which homogenized solutions are prone. The incorporation of improved space discretization and time integration schemes to the computational homogenization framework might prevent these numerical artifacts and should be investigated further. In future work, the developed computational homogenization framework will be applied to 2D and 3D metamaterials, where the coupling between different wave polarizations in finite-size structures will become important.