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

HAL
open science

Nonlinear dynamics of two-dimensional granular periodic lattices

Diala Bitar, Nadim Kacem, Noureddine Bouhaddi

► To cite this version:

Diala Bitar, Nadim Kacem, Noureddine Bouhaddi. Nonlinear dynamics of two-dimensional granular periodic lattices. ISMA2016, International Conference on Noise and Vibration Engineering, Sep 2016, Louvain, Belgium. hal-01512532

HAL Id: hal-01512532

https://hal.science/hal-01512532

Submitted on 23 Apr 2017

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Nonlinear dynamics of two-dimensional granular periodic lattices

D. Bitar¹, N. Kacem¹, N. Bouhaddi¹

¹ FEMTO-ST Institute, UMR 6174, Department of Applied Mechanics, University of Franche-Comté, UBFC, 24 rue de l'Épitaphe F-25000 Besançon, France.

e-mail: diala.bitar@femto-st.fr

## Abstract

We study wave propagation in a damped driven discrete nonlinear two-dimensional (2D) periodic structure. The considered system is composed of identical spherical particles and subjected to harmonic horizontal base excitation, where its collective nonlinear dynamics is caused by the Hertzian contact law and the applied compressive load. The coupled equations governing the nonlinear vibrations of the considered system have been solved using an analytical-numerical solving procedure, which allows a detailed study of complex model interactions of the periodic structure. The multiple scales method coupled with standing wave decomposition were used to transform this nonlinear system into a set of complex differential equations, which has been numerically solved. Numerical simulations have been performed in order to analyze the stability, the modal interactions and the bifurcation topologies resulting from the collective dynamics of particles in contact, in a 2D configuration.

## 1 Introduction

The effect of structural periodicity on wave propagation has been extensively studied in a wide array of fields, such as vibrations of spring-mass systems, arrays of coupled pendulums, electrostatically coupled microbeams, electromagnetic waves in photonic crystals and sound waves in phononic crystals. Several researches were based on studying wave propagation in two-dimensional linear periodic structures. For instance, Leamy [1] et al. detailed an exact wave-based approach for characterizing wave propagation in two-dimensional periodic lattices, while Ruzzene et al. [2] evaluated the dynamic behavior of 2D cellular structures with the focus on the effect of the dynamics of the propagation of elastic waves within the structure. Moreover, Zou et al. [3] presented a detailed calculations of the dispersion relations of bulk waves propagating in two-dimensional (2D) piezoelectric composite structures and showed that the first bandgaps could be controlled according to need. In the field of structural mechanics, Langley [4] et al. considered a beam grillage constructed from strips of aluminum with bolted joints, while in mechanical vibrations Zhou et al. [5] extended the study of wave finite element method to the two-dimensional periodic beam grillage to calculate the wave dispersion characteristics.

Wave propagation in nonlinear periodic media is spare and as most of the structures are inherently nonlinear it is important to investigate the effect of such nonlinearities on the wave propagation characteristics. Nonlinear periodic structures may exhibit a variety of rich dynamic properties such as multi-stability, modal interactions, bifurcation topology transfer, existence of highly stable localized solutions, solitary waves and variation in wave speeds, propagation direction related to wave amplitude and nonlinearity and so on. For instance, in the field of acoustics, Narisetti et al. [6] developed a perturbation approach for two-dimensional mono-atomic lattice of identical masses connected by nonlinear spring stiffness of type Duffing between two adjacent masses to evaluate the influence of nonlinearities on the location of band gaps, group velocity
magnitudes and the direction of energy propagation. In the field of photonics and phononics, Maldovan et al. [7] demonstrated theoretically the simultaneous localization of photons and phonons in the same spatial region by introducing lattice defects in a periodic array of dielectric/elastic material that exhibits gaps for both electromagnetic and elastic waves. In Biomedical, Sheldon et al. [8] modeled the muscular layer of the uterus as a 2D periodic structures of cells coupled through resistors that represents gaps junctions, where they demonstrated that the spatial heterogeneity enhances and modulates excitability in a mathematical model of the muscular layer of the uterus. Moreover, in optics, Kantner et al. [9] showed that arbitrary stable spatio-temporal periodic patterns can be created in 2D lattices of coupled oscillators with inhomogeneous coupling delay.

Hertzian chains of spherical particles represent a famous example of nonlinear discrete structure, where the nonlinearity is caused by the contact law and the applied compressive load. The dynamic behavior one-dimensional (1D) chains of spheres have been widely studied [10]. Heterogeneous, 1D strongly compressed spheres chain has been used to generate solitary wave [11, 12]. It was observed that the wave speed changes with increase in amplitude at certain frequency leading to small changes in the band edges [13]. In addition, Theocharis et al. [14] describe the dynamic behavior of nonlinear periodic phononic structures, along with how such structures can be utilized to affect the propagation of mechanical waves. However, the study of 2D granular lattices has not been widely studied. Recently, Leonard et al. [15] studied the stress wave properties in a basic 2D granular crystals. Aware about the richness and the complexity of the collective dynamics occurring between localized modes in periodic nonlinear arrays [16], we study the nonlinear dynamics of 2D spherical particles in contact under compression.

In this context, we investigate wave propagation in a discrete nonlinear two-dimensional damped granular periodic lattice subject to an external harmonic excitation with dynamics governed by a Hertzian contact model. The coupled equations governing the nonlinear vibrations of the considered system have been solved using an analytical-numerical solving procedure, which allows a detailed study of complex model interactions of the periodic structure. The multiple scales method coupled with standing wave decomposition were used to transform this nonlinear system into a set of complex algebraic equations, which has been numerically solved. Numerical simulations have been performed in order to analyze the stability, the modal interactions and the bifurcation topologies resulting from the collective dynamics of particles in contact, in a 2D configuration.

## 2 Equivalent Mass-Spring system

Consider a 2D periodic structure of identical unit cells subjected to harmonic horizontal base excitation. Each unit cell consist of spherical particle coupled with two uniform beams and subjected to compressive forces applied to all ends (Figures 1(a) 1(b))). Hence, an equivalent spring-mass system can be constructed as shown in Figure 1(c). The displacement of the i, jth sphere from its equilibrium position in the initially compressed chain is defined by $u_{i,j}$. The periodic system is subjected to harmonic horizontal base excitation $\ddot{y}(\tau)$. For dynamical displacements with small amplitudes relative to those due to static load, where a power series expansion was calculated in order to obtain the following i, jth particle equivalent equation of motion.

$$ 
\begin{aligned}
& m_{eq}\ddot{u}_{i,j} + k_{eq}u + c_x(-\dot{u}_{i-1,j} + 2\dot{u}_{i,j} - \dot{u}_{i+1,j}) + c_y(-\dot{u}_{i,j-1} + 2\dot{u}_{i,j} - \dot{u}_{i,j+1}) \\
& \qquad -k_x(u_{i-1,j} - 2u_{i,j} + u_{i+1,j}) - k_y(u_{i,j+1} - 2u_{i,j} + u_{i,j+1}) \\
& \qquad -k_{2x}[(u_{i-1,j} - u_{i,j})^2 - (u_{i,j} - u_{i+1,j})^2] - k_{2y}[(u_{i,j-1} - u_{i,j})^2 - (u_{i,j} - u_{i,j+1})^2] \\
& \qquad +k_{3x}[(u_{i-1,j} - u_{i,j})^3 - (u_{i,j} - u_{i+1,j})^3] + k_{3y}[(u_{i,j-1} - u_{i,j})^3 - (u_{i,j} - u_{i,j+1})^3] \\
&= -m_{eq}\ddot{y} \\
&= Y_0\cos(\Omega\tau)
\end{aligned}
\quad (1) $$
<figure><img src="image_3.png" /><figcaption>Figure 1: (a) 2D periodic nonlinear lattice of spheres in contact under compression and subjected to horizontal base excitation. (b) A view of the system from the top. (c) An equivalent Mass-Spring system.</figcaption></figure>

where $m_{eq}$ is the effective mass of the two-beams-sphere system, $k_{eq}$ is the equivalent structural beams linear stiffness, $c_x = \xi_x m_{eq} \omega_0$ and $c_y = \xi_y m_{eq} \omega_0$ are the damping coefficients in the x and y directions respectively. $\omega_0^2 = \frac{k_{eq}+2k_x+2k_y}{m_{eq}}$ is the natural frequency where $k_x$ and $k_y$ represent the linear coupling spring constants. In addition, $k_{2x}$, $k_{2y}$ and $k_{3x}$, $k_{3y}$ are respectively the coupling quadratic and cubic parameters in the x and y directions and $Y_0 = m_{eq} \Omega^2 y_0$ the harmonic based excitation amplitude ($y(\tau) = y_0 \cos(\Omega \tau)$).

For boundary conditions, we suppose that we have fixed identical systems in both directions, where we define $2(N + M + 1)$ extra variables and set them to zero as: $u_{0,j} = u_{N+1,j} = u_{i,0} = u_{0,M+1} = 0$.

## 2.1 Normalized equations

It proves convenient to define the following scaled dimensionless variables:
$$ \tau = \omega_0 t \qquad x_{i,j} = \frac{u_{i,j}}{u_D} \qquad (2) $$

where $u_D = \frac{Y_0}{2(c_x+c_y)\omega_0}$ is the dynamic displacement of the associated linear system while neglecting the linear coupling with $Q = \frac{m_{eq}\omega_0}{2(c_x+c_y)}$. By replacing these variable into the equations of motion and dividing by $\frac{Y_0\omega_0 m_{eq}}{2(c_x+c_y)}$, we obtain the following nondimensional system:

$$ \begin{aligned} & \ddot{x}_{i,j} + x_i + \xi_x(-\dot{x}_{i-1,j} + 2\dot{x}_{i,j} - \dot{x}_{i+1,j}) + \xi_y(-\dot{x}_{i,j-1} + 2\dot{x}_{i,j} - \dot{x}_{i,j+1}) \\ & -R_x(x_{i-1,j} - 2x_{i,j} + x_{i+1,j}) - R_y(x_{i,j-1} - 2x_{i,j} + x_{i,j+1}) \\ & -\beta_{2x}[(x_{i-1,j} - x_{i,j})^2 - (x_{i,j} - x_{i+1,j})^2] - \beta_{2y}[(x_{i,j-1} - x_{i,j})^2 - (x_{i,j} - x_{i,j+1})^2] \\ & +\beta_{3x}[(x_{i-1,j} - x_{i,j})^3 - (x_{i,j} - x_{i+1,j})^3] + \beta_{3y}[(x_{i,j-1} - x_{i,j})^3 - (x_{i,j} - x_{i,j+1})^3] \\ & = -2(\xi_x + \xi_y) \cos(\frac{\Omega}{\omega_0}t) \end{aligned} \qquad (3) $$

The parameters appearing in Equation (3) are:

$$ \begin{aligned} \xi_x &= \frac{c_x}{m_{eq}\omega_0}, & \xi_y &= \frac{c_y}{m_{eq}\omega_0}, & R_x &= \frac{k_x}{k_{eq}+2k_x+2k_y}, & R_y &= \frac{k_x}{k_{eq}+2k_x+2k_y}, \\ \beta_{2x} &= \frac{k_{2x}Y_0}{2(c_x+c_y)\omega_0^3 m_{eq}}, & \beta_{2y} &= \frac{k_{2y}Y_0}{2(c_x+c_y)\omega_0^3 m_{eq}}, \\ \beta_{3x} &= \frac{k_{3x}Y_0^2}{4(c_x+c_y)\omega_0^4 m_{eq}}, & \text{and} & \beta_{3y} &= \frac{k_{3y}Y_0^2}{4(c_x+c_y)\omega_0^4 m_{eq}} \end{aligned} \qquad (4) $$

## 2.2 Linear study

Equation (3) describes a system of Duffing oscillators linearly coupled by springs $k_x$ and $k_y$ in the x and y directions respectively and subjected to harmonic horizontal base excitation. To determine the natural frequencies and their eigenvectors, Equation (3) can be written in matrix form as

$$ \underbrace{\ddot{X} + C\dot{X} + K_L X}_{\text{Linear part}} + \underbrace{F_{NL}(X)}_{\text{Nonlinear part}} = F(t), \qquad (5) $$

where $X = [x_{1,1}, x_{1,2}, \dots, x_{1,M}, \dots \ x_{i,1}, x_{i,2}, \dots, x_{i,M}, \dots]$ is the N×M-dimensional variable vector, C the damping matrix, $F(t) = -2(\xi_x + \xi_y) \cos(\frac{\Omega}{\omega_0}t)[1, \dots, 1]^T$ the excitation vector, $K_{NL}(X)$ the nonlinear stiffness vector and

$$ K_L = \begin{pmatrix} K_1 & -K_2 & & & & 0 \\ -K_2 & K_1 & -K_2 & & & 0 \\ & \ddots & \ddots & & & 0 \\ 0 & & -K_2 & K_1 & -K_2 & 0 \\ & & & -K_2 & -K_2 & K_1 \end{pmatrix} \qquad (6) $$

where

$$ K_1 = \begin{pmatrix} 1 + 2R_x + 2R_y & -R_y & & & & 0 \\ -R_y & 1 + 2R_x + 2R_y & -R_y & & & 0 \\ & \ddots & \ddots & & & 0 \\ 0 & & -R_y & 1 + 2R_x + 2R_y & -R_y & 0 \\ & & & -R_y & 1 + 2R_x + 2R_y & 0 \end{pmatrix} \qquad (7) $$
and $K_2 = R_x I_M$. The eigenvalues $\omega_{i,j}$ are represented in the following form for few coupled oscillators in contact:

$$ 
\begin{aligned}
N = 1, M = 1: & \omega_{1,1} = \sqrt{1 + 2R_x + 2R_y} \\
N = 2, M = 1: & \omega_{1,1} = \sqrt{1 + R_x + 2R_y} \\
& \omega_{2,1} = \sqrt{1 + 3R_x + 2R_y}
\end{aligned}
\qquad (8) $$

We may express all normal frequencies relative to the same nondimensional reference frequency which is 1, so that

$$ \omega_{i,j} = \sqrt{(1 + \lambda_i R_x + \lambda_j R_y)} \quad (i = 1, \dots, N \text{ and } j = 1, \dots, M), \qquad (9) $$

We suppose that $k_x \ll k_{eq} + 2k_x + 2k_y$ and $k_y \ll k_{eq} + 2k_x + 2k_y$. Consequently, $R_x$ and $R_y \ll 1$ and:

$$ \omega_{i,j} \approx 1 + \frac{1}{2}\lambda_i R_x + \frac{1}{2}\lambda_j R_y \quad (i = 1, \dots, N \text{ and } j = 1, \dots, M). \qquad (10) $$

This assumption leads to the creation of linear closed modes which permits to study the effect of the mode localization on the collective dynamics. In the following section, we will proceed to solve the normalized equivalent differential System (3) using a perturbation technique. Before doing so, we should note that the quadratic nonlinearities are much smaller than the Duffing coupling parameters, thereby their nonlinear effects are considered negligible.

## 2.3 Analytical-numerical approach

In this paper, we use the multiple scales method as an approximate analytical solution to solve nonlinear differential systems. The main advantages of the present approach are the capacity of handling weakly coupled nonlinear systems, which permits to visualize all physical responses branches and their properties in terms of modal interactions and bifurcation topology transfer.

We introduce the parameters listed in Equation (11) and set the external frequency an amount $\varepsilon\omega_0\Omega_D$ away from the resonant frequency, whereby they contribute to the equations of amplitude.

$$ \xi_x = \frac{1}{2}\varepsilon\eta_x \quad \xi_x = \frac{1}{2}\varepsilon\eta_x \quad R_x = \frac{1}{2}\varepsilon\Gamma_x \quad R_y = \frac{1}{2}\varepsilon\Gamma_y \quad \text{and} \quad 2(\xi_x + \xi_y) = \varepsilon^{3/2}g \qquad (11) $$

We calculate the responses of the spheres lattice subjected to harmonic horizontal excitation using secular perturbation theory, by expanding $x_{i,j}(t)$ as a sum of standing-wave modes with slowly varying amplitude [17].

$$ 
\begin{aligned}
x_{i,j}(t) = \varepsilon^{1/2} \sum_{r=1}^{N} \sum_{p=1}^{M} (A_{r,p}(T) \sin(nq_r) \sin(dq_p) e^{i\omega_{r,p}t} + \text{c.c.}) \\
+ \varepsilon^{3/2} x_{i,j}^{(1)}(t) + \dots, \quad i = 1, \dots, N \text{ and } j = 1, \dots, M
\end{aligned}
\qquad (12) $$

where $T = \varepsilon t$ is a slow time variable, that authorizes the complex amplitude $A_{r,p}(T)$ to vary slowly in time. Since we proposed fixed boundary conditions, the possible wave components $q_r, d_p$ can be given as
$$ \left\{ \begin{array}{ll} q_r = \frac{r\pi}{N+1}, & r = 1, \dots, N \\ d_p = \frac{p\pi}{M+1}, & p = 1, \dots, M \end{array} \right. \qquad (13) $$

After replacing the proposed solution (12) into the equation of motion, we can get at the order of $\epsilon^{3/2}$, $N \times M$ equations of the form:

$$ \ddot{x}_{i,j}^{(1)} + x_{i,j}^{(1)} = \sum_{r} \sum_{p} (r, p^{\text{th}} \text{secular term}) e^{i\omega_{r,p}t} + \text{other terms} \quad (14) $$

We must eliminate the $N \times M$ secular terms so that $x_{i,j}(t)$ remains finite which allows us to determine the equations for the slowly varying amplitudes $A_{r,p}(T)$. To extract the equations for these amplitudes, we make use of the orthogonality of the modes, by multiplying the $r, p^{\text{th}}$ secular term by $\sin(nq_r) \sin(dq_p)$ and summing over n and m. We also express all normal frequencies relative to 1, so that:

$$ \omega_{r,p} = 1 + \epsilon\Omega_{r,p} \qquad (15) $$

We find that the equation of the $r, p^{\text{th}}$ amplitude $A_{r,p}(T)$ is given by:

$$ \begin{aligned} & 2i\omega_{r,p}\frac{dA_{r,p}}{dT} + 2i\omega_{r,p}(\eta_x\sin^2(\frac{q_r}{2}) + \eta_y\sin^2(\frac{d_p}{2}))A_{r,p} + 2(\Gamma_x\sin^2(\frac{q_r}{2}) + \Gamma_y\sin^2(\frac{d_p}{2}))A_{r,p} \\ & +3\beta_{3x}\sin(\frac{q_r}{2})\sum_{i,k,l}\sum_{j,s,o}\sin(\frac{q_i}{2})\sin(\frac{q_k}{2})\sin(\frac{q_l}{2})A_{i,j}A_{k,s}A_{l,o}^*\Delta_{ikl;r}^{(1)}\Delta_{js;p}^{(2)}e^{i(\Omega_{i,j}+\Omega_{k,s}-\Omega_{l,o}-\Omega_{r,p})T} \\ & +3\beta_{3y}\sin(\frac{d_p}{2})\sum_{i,k,l}\sum_{j,s,o}\sin(\frac{d_p}{2})\sin(\frac{d_r}{2})\sin(\frac{d_s}{2})A_{i,j}A_{k,s}A_{l,o}^*\Delta_{ikl;r}^{(2)}\Delta_{js;p}^{(1)}e^{i(\Omega_{i,j}+\Omega_{k,s}-\Omega_{l,o}-\Omega_{r,p})T} \\ & = \frac{2}{(N+1)(M+1)}g\sum_{n}\sum_{m}\sin(nq_r)\sin(md_p)e^{i(\Omega_D-\Omega_{r,p})T} \end{aligned} \qquad (16) $$

In solving Equations (16), we write the periodic steady state solution in the following form:

$$ A_{r,p} = a_{r,p}e^{i[\Omega_D - \Omega_{r,p}]T} \qquad (17) $$

Substituting Equation (17) into (16), we obtain the required Equation (5) for the complex amplitudes $a_{r,p}$.

$$ \begin{aligned} & -2\omega_{r,p}(\Omega_D - \Omega_{r,p})a_{r,p} + 2i\omega_{r,p}(\eta_x\sin^2(\frac{q_r}{2}) + \eta_y\sin^2(\frac{d_p}{2}))a_{r,p} + 2(\Gamma_x\sin^2(\frac{q_r}{2}) + \Gamma_y\sin^2(\frac{d_p}{2}))a_{r,p} \\ & +3\beta_{3x}\sin(\frac{q_r}{2})\sum_{i,k,l}\sum_{j,s,o}\sin(\frac{q_i}{2})\sin(\frac{q_k}{2})\sin(\frac{q_l}{2})a_{i,j}a_{k,s}a_{l,o}^*\Delta_{ikl;r}^{(1)}\Delta_{js;p}^{(2)} \\ & +3\beta_{3y}\sin(\frac{d_p}{2})\sum_{i,k,l}\sum_{j,s,o}\sin(\frac{d_p}{2})\sin(\frac{d_r}{2})\sin(\frac{d_s}{2})a_{i,j}a_{k,s}a_{l,o}^*\Delta_{ikl;r}^{(2)}\Delta_{js;p}^{(1)} \\ & = \frac{2}{(N+1)(M+1)}g\sum_{n}\sum_{m}\sin(nq_r)\sin(md_p) \end{aligned} \qquad (18) $$

These complex algebraic equations which define the time independent mode amplitudes are the main result of the perturbation analysis applied to the normalized equivalent differential system. In order to study the collective dynamics of the weakly coupled 2D periodic structure, we start by solving the algebraic system for a set of design parameters listed in Table 1 which satisfies the modes localization assumption given by Equation (10).
<table><caption>Table 1: Design parameters for the corresponding periodic structure depicted in Figure 1(c)</caption><thead><tr><th>Configuration</th><th>η<sub>x</sub></th><th>η<sub>y</sub></th><th>Γ<sub>x</sub></th><th>Γ<sub>y</sub></th><th>β<sub>3x</sub></th><th>β<sub>3y</sub></th><th>g</th></tr></thead><tbody><tr><td>1</td><td>8.10<sup>-3</sup></td><td>8.10<sup>-3</sup></td><td>16.10<sup>-3</sup></td><td>16.10<sup>-3</sup></td><td>29.10<sup>-3</sup></td><td>29.10<sup>-3</sup></td><td>16.10<sup>-3</sup></td></tr><tr><td>2</td><td>8.10<sup>-3</sup></td><td>8.10<sup>-3</sup></td><td>16.10<sup>-3</sup></td><td>16.10<sup>-3</sup></td><td>15.10<sup>-3</sup></td><td>15.10<sup>-3</sup></td><td>16.10<sup>-3</sup></td></tr><tr><td>3</td><td>8.10<sup>-3</sup></td><td>4.10<sup>-3</sup></td><td>16.10<sup>-3</sup></td><td>10.10<sup>-3</sup></td><td>5.10<sup>-3</sup></td><td>19.10<sup>-3</sup></td><td>12.10<sup>-3</sup></td></tr><tr><td>4</td><td>8.10<sup>-3</sup></td><td>18.10<sup>-4</sup></td><td>16.10<sup>-3</sup></td><td>7.10<sup>-3</sup></td><td>2.10<sup>-3</sup></td><td>20.10<sup>-3</sup></td><td>98.10<sup>-4</sup></td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Figure 2: Response intensities of two weakly coupled oscillators under harmonic base external excitation as a function of the detuning parameter $\Omega_D$, for the first set of design parameters listed in Table 1.</figcaption></figure>

## 3 Results and discussions

Before starting our investigations we should note the following statements. First, the second member of the Equation (18) is proportional to the sum of standing waves, which is null for all i,j modes where at least one position is even. For a single direction ($N = 1$, or $M = 1$), the general periodic structure is reduced to an array of coupled oscillators, where a detailed study has been carried out [16].

Particularly, for two coupled complex algebraic equations, we used Mathematica to solve them numerically, where the stability has been performed based on the eigenvalues of the Jacobian matrix of the differential System (16) for each point. Figure 2 displays the intensity responses of two weakly coupled oscillators ($N = 1, M = 2$), as a function of the detuning parameter $\Omega_D$ for the first design parameters listed in Table 1. Black, orange and gray curves indicate the Single Mode (SM), Double Mode (DM) and the unstable solution branches. SM solutions are generated by the null trivial solution $a_{1,2} = 0$, where $a_{1,1}$ corresponds to a single forced Duffing oscillator. Remarkably, there are additional branches labeled DM, involving the excitation of both modes collectively. With two weakly coupled oscillators we can obtain up to four stable solutions for a given frequency. Hence, for 2D periodic structure we expect additional complex phenomena not obtained for one dimensional array.

Therefore, we study the case of 2 × 2 weakly coupled oscillators, for several design parameters. Starting by the case where we consider the same linear and nonlinear coupling parameters in both directions. Figure 4 shows the squares of the amplitudes of the different modes as a function of the detuning parameters for the second design parameters in Table 1. $DM_{i,j}$ represent the Double Mode solution branches generated by two different modes and $QM$ is the Quadruple Mode involved by the excitation of all modes collectively.

Note that the modes localization have interesting effects on the response intensities with additional multi-modal solutions that share the same frequency range with bifurcation topology transfer. In addition, as we
choose identical coupling parameters, $a_{1,2}$ and $a_{2,1}$ become symmetric and interact in the same manner with $a_{1,1}$. Moreover, the Triple Mode (TM) solution branches are missing in this case. In fact, as by replacing one of the non-excited modes $a_{1,2}$, $a_{2,1}$ or $a_{2,2}$ by 0 into any of their equations, we obtain a null trivial solution, which corresponds to the case of a DM.

<figure><img src="image_5.png" /><figcaption>Figure 3: Response intensities of 2 × 2 weakly coupled oscillators under harmonic base external excitation as a function of the detuning parameter $\Omega_D$, for the second set of design parameters listed in table 1.</figcaption></figure>

By choosing different linear and nonlinear coupling parameters, we calculate the response intensities for the third design parameters and plot them in figure 4. Remarkably, the DM solution branches are being separated, with additional narrow multimodal solution branches. In addition by decreasing the coupling parameters in the y direction, the $DM_{1,1;1,2}$ and QM move away from the resonant branch in $|A_{1,1}|^2$ and their frequency range decreases.

For the fourth set of design parameters, the coupling in the y direction is highly decreased. The responses intensities plotted in Figure 5 show that the small solution branches join together in one $DM_{1,1;1,2}$ branch heading in the direction of the non-resonant branch of $|A_{1,1}|^2$, While the response intensity of $|A_{2,1}|^2$ takes an elliptical form.
<figure><img src="image_6.png" /><figcaption>Figure 4: Response intensities of 2 × 2 weakly coupled oscillators under harmonic base external excitation as a function of the detuning parameter $\Omega_D$, for the third set of design parameters listed in table 1.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 5: Response intensities of 2 × 2 weakly coupled oscillators under harmonic base external excitation as a function of the detuning parameter $\Omega_D$, for the fourth set of design parameters listed in table 1.</figcaption></figure>
# 4 Conclusion

To conclude, a 2D periodic nonlinear structure composed of coupled two-beams-Mass structures under harmonic horizontal base excitation has been considered. The equations of motions of the equivalent Mass-Spring system has been normalized and solved using an analytical-numerical approach based on secular perturbation theory coupled with standing waves decomposition. Intensity responses for few coupled oscillators has been calculated for several design parameters that satisfies the modes localization condition. Complex nonlinear phenomena with high number of stable multimodal solution branches has been observed. Considering different linear and nonlinear coupling parameters in the different directions, the intensity responses topology can be modified with additional multimodal solution branches.

# Acknowledgements

This project has been performed in cooperation with the Labex ACTION program (contract ANR-11-LABX-01-01).

# References

[1] M. J. Leamy, *Exact wave-based bloch analysis procedure for investigating wave propagation in two-dimensional periodic lattices*, Journal of Sound and Vibration, 331(7):1580–1596, 2012.

[2] F. Soranna, M. Ruzzene, F. Scarpa, *Wave beaming effects in two-dimensional cellular structures*, Smart Materials and Structures, 12(3):363372, 2003.

[3] B. Liang, J. Cheng, X. Zou, Q. Chen, *Control of the elastic wave bandgaps in two-dimensional piezoelectric periodic structures*, Smart Materials and Structures, 17(1):363372, 2008.

[4] R. S. Langley, *The response of two-dimensional periodic structures to point harmonic forcing*, Journal of Sound and Vibration, 197(4):447–469, 1996.

[5] M. N. Ichchou, A. M. Zine, C. W. Zhou, J. P. Lainé, *Multi-scale modelling for two-dimensional periodic structures using a combined mode/wave based approach*, Computers and Structures, 154:145–162, 2015.

[6] M. J. Leamy, R. K. Narisetti, M. Ruzzene, *A perturbation approach for analyzing dispersion and group velocities in two-dimensional nonlinear periodic lattices*, Journal of Vibration and Acoustics, 133(6):0739–3717, 2011.

[7] M. Maldovan, E. L. Thomas, *Simultaneous localization of photons and phonons in two-dimensional periodic structures*, Applied Physics Letters, 88(25), 2006.

[8] R. E. Sheldon, M. Baghdadi, C. McCloskey, A. M. Blanks, A. Shmygol, H. A. van den Berg, *Spatial heterogeneity enhances and modulates excitability in a mathematical model of the myometrium*, Journal of The Royal Society Interface, 10(86), 2013.

[9] S. Yanchuk, M. Kantner, E. Schöll, *Delay-induced patterns in a two-dimensional lattice of coupled oscillators*, Scientific Reports, 5:8522, 2015.

[10] V. F. Nesterenko, *Propagation of nonlinear compression pulses in granular media*, Journal of Applied Mechanics and Technical Physics, 24(5):733–743, 1983.

[11] V.F. Nesterenko, *Dynamics of Heterogeneous Materials*, Dynamics of Heterogeneous Materials, Springer, 2001.
[12] J. Yang, C. Silvestro, D. Khatri, L. De Nardo, C. Daraio, *Interaction of highly nonlinear solitary waves with linear elastic media*, Phys. Rev. E, 83:046606, Apr 2011.

[13] N. Boechler, C. Daraio, R. K. Narisetti, M. Ruzzene, M. J. Leamy, *Analytical and experimental analysis of bandgaps in nonlinear on dimensional periodic structures*, In Tsung-Tsong Wu and Chien-Ching Ma, editors, *IUTAM Symposium on Recent Advances of Acoustic Waves in Solids*, volume 26 of *IUTAM Bookseries*, pages 209–219, Springer Netherlands, 2010.

[14] G.Theocharis, N.Boechler, C.Daraio, *Nonlinear Periodic Phononic Structures and Granular Crystals*, pages 217–251, Springer Berlin Heidelberg, Berlin, Heidelberg, 2013.

[15] A. Leonard, F. Fraternali, C. Daraio, *Directional wave propagation in a highly nonlinear square packing of spheres*, Experimental Mechanics, 53(3):327–337, 2013.

[16] D. Bitar, N. Kacem, N. Bouhaddi, M. Collet, *Collective dynamics of periodic nonlinear oscillators under simultaneous parametric and external excitations*, Nonlinear Dynamics, 82(1-2):749–766, 2015.

[17] R. Lifshitz, M. C. Cross, *Response of parametrically driven nonlinear coupled oscillators with application to micromechanical and nanomechanical resonator arrays*, Physical review B, 67, 2003.