A Bloch-based procedure for dispersion analysis of lattices with periodic time-varying properties

TL;DR

This paper introduces a Bloch-based method for analyzing the dispersion properties of lattices with periodically varying parameters, enabling efficient prediction of wave behavior in dynamic, time-modulated media.

Contribution

The paper presents a novel eigenvalue-based approach combining Bloch theorem and harmonic balance for dispersion analysis of time-varying lattices, applicable across multiple physical systems.

Findings

Effectively predicts dispersion in spring-mass superlattices

Demonstrates broad applicability to elastic media

Shows potential for non-reciprocal wave control

Abstract

We present a procedure for the systematic estimation of the dispersion properties of linear discrete systems with periodic time-varying coefficients. The approach relies on the analysis of a single unit cell, making use of Bloch theorem along with the application of a harmonic balance methodology over an imposed solution ansatz. The solution of the resulting eigenvalue problem is followed by a procedure that selects the eigen-solutions corresponding to the ansatz, which is a plane wave defined by a frequency-wavenumber pair. Examples on spring-mass superlattices demonstrate the effectiveness of the method at predicting the dispersion behavior of linear elastic media. The matrix formulation of the problem suggests the broad applicability of the proposed technique. Furthermore, it is shown how dispersion can inform about the dynamic behavior of time-modulated finite lattices. The technique can be extended to multiple areas of physics, such as acoustic, elastic and electromagnetic systems, where periodic time-varying material properties may be used to obtain non-reciprocal wave propagation.

Equations

$$\mathbf{M}{\ddot{\mathbf{u}}}_n+{\mathbf{K}}^{(l)}(t){\mathbf{u}}_{n-1}+\mathbf{K}(t){\mathbf{u}}_n+{\mathbf{K}}^{(r)}(t){\mathbf{u}}_{n+1}=0.$$

$$\mathbf{K}(t)=\mathbf{K}(t+T_m).$$

$${\mathbf{K}}^{(l)}(t)$$

$$\mathbf{K}(t)$$

$${\mathbf{K}}^{(r)}(t)$$

$${\mathbf{u}}_n(t)=\mathbf{a}(t)e^{i(n\kappa {\lambda }_m+\omega t)},$$

$$\mathbf{a}(t)=\sum _{p=-\infty }^{\infty }{\mathbf{a}}_p{e}^{ip{\omega }_{m}t}.$$

$${\mathbf{u}}_{n-1}(t)=e^{-i\mu }{\mathbf{u}}_n(t),{\mathbf{u}}_{n+1}(t)=e^{i\mu }{\mathbf{u}}_n(t),$$

$$\mathbf{M}{\ddot{\mathbf{u}}}_n(t)+\left({\mathbf{K}}^{(l)}(t)e^{-i\mu }+\mathbf{K}(t)+{\mathbf{K}}^{(r)}(t)e^{i\mu }\right){\mathbf{u}}_n(t)=0,$$

$$\mathbf{M}{\ddot{\mathbf{u}}}_n(t)+\hat{\mathbf{K}}(\mu ,t){\mathbf{u}}_n(t)=0.$$

$$\hat{\mathbf{K}}(\mu ,t)={\mathbf{K}}^{(l)}(t)e^{-i\mu }+\mathbf{K}(t)+{\mathbf{K}}^{(r)}(t)e^{i\mu }.$$

$$-(\omega +p{\omega }_m{)}^2\mathbf{M}+\sum _{q=-\infty }^{\infty }{\hat{\mathbf{K}}}_q(\mu ){\mathbf{a}}_{p-q}=0.$$

$$-(\omega +p{\omega }_m{)}^2\mathbf{M}+\sum _{q=-P+p}^{P+p}{\hat{\mathbf{K}}}_q(\mu ){\mathbf{a}}_{p-q}=0.$$

$${\lambda }_{r,p}={\omega }_r+p{\omega }_m$$

$${\mathbf{\alpha }}_{r,p}=[{\mathbf{a}}_{-P}^{(p)},...,{\mathbf{a}}_0^{(p)},...,{\mathbf{a}}_{+P}^{(p)}{]}_r^T.$$

$${\mathbf{u}}_n^{(r)}(t)=e^{in\mu }\sum _{p=-P}^{+P}\left(\sum _{q=-P}^{+P}{\mathbf{a}}_q^{(p)}{|}_r{e}^{i({\omega }_r+q{\omega }_m)t}\right)e^{ip{\omega }_{m}t}$$

$${\mathbf{u}}_n^{(r)}(t)=e^{i(n\mu +{\omega }_{r}t)}\sum _{p=-P}^{+P}\left(\sum _{q=-P}^{+P}{\mathbf{a}}_q^{(p)}{|}_r{e}^{i(p+q){\omega }_{m}t}\right)$$

$${\mathbf{u}}_n^{(r)}(t)\approx e^{i(n\mu +{\omega }_{r}t)}\sum _{p=-P}^{+P}\left(\sum _{q=-P}^{+P}{\mathbf{a}}_q^{(p)}{|}_r{\delta }_{p+q,0}\right)$$

$$\dot{\mathbf{u}}(n,t)=\sum _{\omega }\sum _{\mu }{\mathbf{c}}^{(\omega ,\mu )}{e}^{i(\mu n-\omega t)}.$$

$${\mathbf{c}}^{(\omega ,\mu )}=⟨\dot{\mathbf{u}}(n,t),e^{i(\mu n-\omega t)}⟩=\frac{1}{TN}{\int }_0^T\sum _{n=1}^N\dot{\mathbf{u}}(n,t)e^{i(\mu n-\omega t)}dt.$$

$$\Phi (\omega ,\mu )=\sum _{r=1}^R\frac{1}{2}m{\left(c_r^{(\omega ,\mu )}\right)}^2=\frac{m}{2(TN{)}^2}\sum _{r=1}^R{\left|{\int }_0^T\sum _{n=1}^N{\dot{u}}_r(n,t)e^{i(\mu n-\omega t)}dt\right|}^2$$

$$\Phi (\omega ,\mu )=\frac{m}{2(TN{)}^2}\sum _{r=1}^R{\left|\Delta t\sum _{j=0}^{T/\Delta t}\sum _{n=1}^N{\dot{u}}_r(n,j\Delta t)e^{i(\mu n-\omega j\Delta t)}\right|}^2$$

$$m{\ddot{u}}_r+k_{r-1}(t)(u_r-u_{r-1})+k_r(t)(u_r-u_{r+1})+k_{g,r}(t)u_r=0.$$

$${\ddot{u}}_r+{\omega }_0^2(1+{\beta }_m{\Phi }_{r-1}(t))(u_r-u_{r-1})+{\omega }_0^2(1+{\beta }_m{\Phi }_r(t))(u_r-u_{r+1})+{\omega }_0^2{\gamma }_g(1+{\beta }_{g,m}{\Phi }_{g,r}(t))u_r=0.$$

$${\beta }_m=\frac{k_m}{k_0},{\beta }_{g,m}=\frac{k_{g,m}}{k_g},{\gamma }_g=\frac{k_g}{k_0}$$

$$k_r(t)=k_0(1+{\beta }_m\cos \left({\omega }_{m}t+{\varphi }_r\right)),$$

$$\begin{array}{rl}{\hat{\mathbf{K}}}_q(\mu )=& {\gamma }_g{\omega }_0^2\mathbf{I}{\delta }_{q,0}+{\omega }_0^2\left[\begin{array}{cccccc}2& -1& 0& 0& \dots & -e^{-i\mu }\\ -1& 2& -1& 0& \dots & 0\\ 0& -1& 2& -1& \dots & 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ -e^{i\mu }& 0& 0& 0& \dots & 2\end{array}\right]{\delta }_{q,0}\\ +& {\omega }_0^2\frac{{\beta }_m}{2}\left[\begin{array}{ccccc}{e}^{iq{\varphi }_R}+e^{qi{\varphi }_1}& -e^{iq{\varphi }_1}& 0& \dots & -e^{iq{\varphi }_R}{e}^{-i\mu }\\ -e^{iq{\varphi }_1}& e^{iq{\varphi }_1}+e^{iq{\varphi }_2}& -e^{iq{\varphi }_2}& \dots & 0\\ 0& -e^{iq{\varphi }_2}& e^{iq{\varphi }_2}+e^{iq{\varphi }_3}& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ -e^{iq{\varphi }_R}{e}^{i\mu }& 0& 0& \dots & e^{iq{\varphi }_{R-1}}+e^{qi{\varphi }_R}\end{array}\right]{\delta }_{q,\pm 1}\end{array}$$

$$k_{g,r}(t)=k_{g,0}(1+{\beta }_{g,m}\cos \left({\omega }_{m}t+{\varphi }_r\right)),$$

$${\hat{\mathbf{K}}}_q(\mu )$$

$$k_{r}(t)=k_{0}(1+\beta_{m}\big{\{}2H\big{[}\cos\left(\omega_{m}t+\phi_{s}\right)\big{]}-1\big{\}}),$$

$${\hat{\mathbf{K}}}_q(\mu )$$

Full text

1 Introduction

Periodic structures are finding increasing applications in diverse fields of science and engineering due to their unique dynamic properties [1]. One growing area of research is controlling the direction of wave propagation using periodic arrays of active materials or meta-materials. During the last decades, a number of researchers have achieved unidirectional wave propagation by spatio-temporal modulation of material properties. For example, Oliner and Cassedy [2] investigated modulated electric circuits and obtained expressions for dispersion curves using continued fractions. Using their work, other authors [3, 4, 5] studied time-varying electromagnetic media with the objective of investigating wave dispersion. Similarly, in mechanical systems, Lurie [6] studied the dynamics of time varying composites and derived their effective properties. Extending this work, Shui et al. [7, 8] recently demonstrated asymmetrical wave properties in spatiotemporal composites. Furthermore, Wright and Cobbold [9, 10] extended multiple scattering theory to phononic crystals with time-varying properties. Using this concept, Engheta and coworkers [11, 12] demonstrated one-way phonon transport in waveguides by a spatio-temporal modulation in graphene. Another example in this line of work is that of Li et al. [13], whereby numerical simulations show tunable wave propagation in helicoidal phononic crystals with coupling between in-plane torsional and out-of-plane longitudinal waves. By modulating the out-of-plane stiffness in space and time using external control, they alter the in-plane wave propagation. Similarly, Wang et al. [14] investigated both numerically and experimentally a tunable acoustic unidirectional device that provide asymmetrical transmission properties.

Expoliting this periodic modulation of properties, numerous works have demonstrated non-reciprocal wave propagation. For example, Fleury et al. [15] achieved one way transport in acoustic channels using a circulator-like device. In a similar line of work, Alù and coworkers [16] demonstrated nonreciprocal excitation in electric circuits by time varying capacitance: a small modulation causes the counter-propagating modes to have distinct frequencies, which results in time-reversal symmetry breaking. Asymptotic analysis is employed to predict frequency splitting for a $3$ port system. Deymier and coworkers [17] recently demonstrated unidirectional wave propagation numerically, modulating the spring stiffness periodically in space and time in a one dimensional spring-mass chain. The authors obtained asymmetrical dispersion diagrams by extracting the data from extensive numerical simulations using the spectral energy density (SED) method [18].

Most of the studies above rely on numerical simulations, which are time consuming, particularly for higher dimensional structures and are not convenient for design and optimization studies. A systematic procedure for obtaining dispersion information of lattices with periodic time-varying properties is therefore of interest for design and analysis of time-varying systems. In this paper, we present a technique for obtaining the dispersion diagrams of discretized linear elastic media with space-time stiffness variability. The method consists of four steps: the appropriate selection of an extended unit cell, a suitable choice of the solution ansatz, the application of Bloch wave analysis and a filtering of the obtained eigen-solutions to identify the frequencies associated with the fundamental plane wave solution.

In periodic time-varying lattices, the size of a unit cell is based not only on the geometry and time independent material properties, but also needs to takes into account the size of the stiffness variability within the system [1, 19, 20]. In order to make use of a suitable ansatz, we draw upon the Hill-determinant method [21], which is widely used to obtain criteria for the stability regions of linear systems with time varying coefficients. Through superposition principle, the motion of a particle is expressed by a plane wave whose amplitude is periodic. Then, we employ a Bloch wave based analysis over the extended unit cell to obtain the eigenvalues for each wavenumber. These three steps were followed by Trainiti and Ruzzene [22] in the application of this analytical method to the analysis of continuous mechanical systems with periodic time-varying coefficients. However, the multiplicity of frequencies in the Hill-determinant method causes eigenvalues whose associated eigenvectors represent a different frequency in displacements to arise, leading to non-physical curves in their dispersion diagram. Our final step applies a filtering procedure using the eigenvectors such that the dispersion diagram reflects only the possible wavenumber-frequency couples.

The outline of this paper is as follows: Section 2 presents the developed procedure for computing the dispersion diagram for spring-mass chains with periodic time varying stiffness. Section 3 illustrates its implementation in the analysis of three distinct systems, and the results are validated with numerical solutions. Sec. 4 shows the usefulness of dispersion diagrams in predicting the broken time-reversal behavior of infinite and finite systems with time varying coefficients. Finally, Sec. 5 summarizes our work and presents potential future research directions.

2 Theoretical background:

2.1 Time varying periodic lattice configuration

We consider a one-dimensional (1D) structure that is characterized by periodically modulated properties. The structure is governed by linear interactions, defined by stiffness coefficients that vary in time, and by constant inertia coefficients. The stiffness modulation is expressed as a traveling wave propagating with veloctiy $v_{m}=\lambda_{m}/T_{m}$, where $\lambda_{m}$ and $T_{m}$ respectively denote the spatial wavelength and temporal period of the modulation. Thus, at any given instant of time, it is possible to describe the structure as the assembly of unit cells that are identified by one spatial modulation period $\lambda_{m}$. Based on this description, and the assumption that the structure has been discretized through a Finite Element procedure, we may express the generalized equation of motion for the $n$-th cell of the assembly as

[TABLE]

where $\bm{M},\bm{K}$ and $\bm{u}_{n}$ denote the mass, stiffness matrices and a vector of degrees of freedom of the unit cell $n$, while $\bm{K}^{(l)},\bm{K}^{(r)}$ describe the interactions of this unit cell with its neighboring unit cells. For simplicity of notation, but without loss of generality, it is assumed that the mass is lumped so that there is no inertial coupling across the cells.

The stiffness matrices are all expressed as a periodic functions of time with period $T_{m}$ so that the following relation holds

[TABLE]

Accordingly, each of the matrices in Eqn. (1) can be expanded in terms of their Fourier series and expressed as

[TABLE]

where $\bm{K}_{q}$, $\bm{K}_{q}^{(l)}$ and $\bm{K}_{q}^{(r)}$ are the corresponding matrix coefficients and $\omega_{m}=2\pi/T_{m}$ is the frequency associated with the temporal modulation.

2.2 Spatio-temporal plane wave expansion for the estimation of dispersion

The dispersion relations for the considered time-varying structure can be estimated by seeking for a plane wave solution with modulated amplitude, which is expressed as

[TABLE]

where $\kappa$ is the wavenumber and $\bm{a}(t)=\bm{a}(t+T_{m})$ is a periodic amplitude function in time. The frequencies in the displacement amplitude $\bm{a}(t)$ depend on the stiffness modulation frequency $\omega_{m}=2\pi/T_{m}$ and it can be expressed as a Fourier series in the form

[TABLE]

Based on the solution form expressed by Eqn. (4), the following relations hold

[TABLE]

where $\mu=\kappa\lambda_{m}$ is the non-dimensional wavenumber. Substituting into Eqn. (1) gives

[TABLE]

which may be written as

[TABLE]

with

[TABLE]

Next, substituting Eqn. (5) into Eqn. (7) and performing harmonic balance, by collecting the terms with frequency $\omega+p\omega_{m}$ we reach

[TABLE]

Choosing a truncation order $P$ for the displacement amplitude $\bm{a}(t)$, i.e., $\bm{a}_{p}=\mathbf{0}$ for $|p|>P$, Eqn. (9) reduces to a system of $(2P+1)R$ equations, with $R$ being the number of degrees of freedom in a unit cell. This system of equations defines a quadratic eigenvalue problem that can be solved in terms of frequency $\omega$ for assigned values of $\mu$. Specifically, $\mu$ spans the reciprocal lattice space, limited to $\mu\in[-\pi,\,\,\,+\pi]$ which defines the First Brillouin zone. It is here considered instructive to analyze the structure of the eigenvalue problem in Eqn. (9). By doing a change of variable $p\to p+r\omega_{m}$, the resulting expression may be written as $-((\omega+r\omega_{m})+p\omega_{m})^{2}\bm{M}+\sum_{q=-\infty}^{\infty}\bm{\hat{K}}_{q}(\mu)\bm{a}_{(p+r)-q}=\mathbf{0}$. Comparing it with Eqn. (9), we observe that if $(\omega,\bm{a}_{p})$ is an eigensolution of the infinite dimensional eigenvalue problem, then $(\omega+r\omega_{m},\bm{a}_{p+r})$ is also a solution.

The solution of the following quadratic eigenvalue problem with truncated terms

[TABLE]

thus leads to $R\times(2P+1)$ eigenvalues that are of the general form

[TABLE]

with $r=1,..,R$ and $p=-P,...,+P$. The associated eigvenvectors can be expressed as

[TABLE]

Equation (10) shows how the eigenvalues are clustered into $R$ groups of $2P+1$ values centered at frequencies $\omega_{r}$ and separated by integer multiples of the modulation frequency $\omega_{m}$. Thus, the solution approach introduces the challenge of identifying the $R$ dispersion branches associated with the sought plane wave solution from the $R(2P+1)$ frequencies obtained from the eigenvalue problem. This challenge is here addressed through a procedure that identifies the plane wave branches based on a weighting factor corresponding to the magnitude of the eigenvector components associated with the fundamental plane wave term.

To illustrate the rationale of this weighting procedure, we consider the plane wave solution imposed in Eqns. (4) and (5). The $r$-th family of solutions ($2P+1$ solutions) associated with this plane wave may be expressed as

[TABLE]

which can be re-organized to read

[TABLE]

Thus, the fundamental plane wave term has a magnitude that is identified for $p+q=0$. It is reasonable to expect that such fundamental is the leading term in the expansion such that

[TABLE]

where $\delta_{p+q,0}$ is the Kronecker delta, i.e. $\delta_{i,j}=1$ if $i=j$, or $\delta_{i,j}=0$ if $i\neq j$. According to the above equation, the dispersion branch corresponding to the fundamental plane wave can be effectively obtained weighting each branch by the magnitude of the fundamental component and applying a thresholding value which filters the branches, thereby avoiding their plotting if the associated eigenvector magnitude is low compared to that of the fundamental. The effectiveness of this procedure in tracking the branch corresponding to the plane wave of interest is illustrated through the examples presented in the next section.

2.3 Spectral energy density method

We now briefly describe the spectral energy density method for obtaining the dispersion diagrams as used by other studies [17]. In this method, displacements are obtained by numerical integration. Full numerical simulations of the system velocity field under a range of excitation frequencies are needed, and the results are used to construct the dispersion diagrams.

Consider a finite lattice having $R$ degrees of freedom per unit cell and a total length of $N$ unit cells, i.e. a total of $RxN$ degrees of freedom. For the algorithm, first the velocity field of each simulation $\dot{\bm{u}}(n,t)$ is projected into the orthogonal Fourier basis $e^{i(\mu n-\omega t)}$

[TABLE]

[TABLE]

where $T$ is the total integration time. The spectral energy density is a measure of the energy content of the each plane wave with specific frequency $\omega$ and wavenumber $\mu$, by considering the kinetic energy of the masses, $(1/2)m\lVert\bm{c}^{(\omega,\mu)}\rVert^{2}$. Given Eqn. (14), the SED $\Phi(\omega,\mu)$ is obtained as

[TABLE]

where $r$ is the position of the particle within the unit cell. We performed the time integration numerically using the trapezoidal rule

[TABLE]

where $\Delta t$ is the grid spacing of the numerical integration.

Last, Eqn. (17) is evaluated for each $(\omega,\kappa)$ pair and the result is a contour plot. The dispersion diagrams are given by the high contour values of $\Phi$, indicating the allowance of plane waves with that specified frequency/wavenumber pair.

3 Verification of the proposed method: Dispersion analysis of spring-mass lattices

The procedure explained in Section 2 is illustrated for the case of the spring-mass lattice shown in Fig. 1. The masses of the lattice are considered constant in time, and are all equal to $m$. The masses are connected to their adjacent neighbors by springs of time-modulated constant $k_{r}(t)=k_{0}+k_{m}\Phi_{r}(t)$, and to the ground by springs that may be also time-modulated, i.e. $k_{g,r}(t)=k_{g_{0}}+k_{g_{m}}\Phi_{g,r}(t)$. According to the derivations in the previous section, the procedure considers periodic time modulations, so that the two functions $\Phi(t),\Phi_{m}(t)$ are both periodic of period $T_{m}$, i.e. $\Phi_{r}(t)=\Phi_{r}(t+T_{m}),\Phi_{g,r}(t)=\Phi_{g,r}(t+T_{m})$.

The equation of motion for the $r$-th mass is

[TABLE]

A characteristic frequency $\omega_{0}^{2}=k_{0}/m$ is introduced, so that Eqn. (18) can be rewritten in non-dimensional form as follows

[TABLE]

where

[TABLE]

We note that the above expression corresponds to a system of Mathieu-Hill equations, whose stability with respect to the modulation parameters has been extensively investigated [21]. The method presented here applies when the parameters lie in the stable regions. In our examples, the modulation properties are small compared to the constant order terms and we assume that solutions lie within the stable regions.

Three examples of stiffness modulated systems are studied. In the first two examples, we analyze the dispersion of one-dimensional lattices having a harmonic modulation of their properties in space and time. In the last example, we study a spring-mass chain with a traveling square-wave modulation, thereby demonstrating the applicability of the method when there are multiple modulation frequencies.

3.1 Harmonic stiffness modulation

A first example considers a stiffness modulation that is imposed along the length of the chain with spatial wavelength $\lambda_{m}=Ra$, where $a$ is the distance between neighboring masses, while $R$ denotes the number of masses in a spatial modulation wavelength $\lambda_{m}$. Figure 1 displays a schematic of the chain along with the spring stiffness modulation.

We first consider a periodic modulation of the springs between adjacent masses, so that the stiffness of the $r$-th spring is

[TABLE]

where $\phi_{r}=2\pi r/R$ is the spatial phase shift related to the position $r$ of the considered spring. Also, the ground stiffnesses are considered constant in time, therefore in Eqn. (19) $\beta_{g,m}=0$. A unit cell is identified that contains $R$ degrees of freedom, whose behavior is governed by a system of equations of the kind of Eqn. (1), subsequently transformed into Eqn. (9), with $\bm{M}=\bm{I}$, where $\bm{I}\in\mathbb{R}^{R\times R}$ is the identity matrix, while

[TABLE]

We study a chain with modulation wavelength $\lambda_{m}=3a$ ($R=3$ masses in each unit cell) and analyze its dispersion properties estimated from the solution of Eqn. (9). Given the considered form of modulation, which is harmonic, the solution is conducted by considering $P=1$ terms, which is found to lead to results that are in agreement with those estimated through numerical simulations and for higher order expansions. Results are presented for the following set of parameters ${\Omega}_{m}={\omega_{m}}/{\omega_{0}}=0.2$, $\beta_{m}=0.15$, $\gamma_{g}=1$.

The solution is conducted by imposing $\mu\in[-\pi,\,\,\pi]$ which provides $R\times(2P+1)=9$ frequency branches. These branches are plotted in terms of the non-dimensional frequency $\Omega=\omega/\omega_{0}$ in Fig. 2(a). Figure 2(a) displays $3$ families of curves that appear parallel and separated by $\Omega_{m}$. These branches cross for certain wavenumber values and fold at the boundaries of the first Brillouin zone. The dispersion diagram corresponds to the branch associated with the fundamental plane wave component among the $3$ families of curves. To identify this branch, the filtering method described in the previous section is applied, whereby the plotting of each curve is based on a weighting parameter defined by the amplitude of the eigenvector corresponding to the fundamental and a thresholding procedure. The result of this process is shown in Fig. 2(b), which shows only the single fundamental branch. As expected the lattice is characterized by a high-pass cut-off at $\Omega=1$ for $\mu=0$ which is defined by the ground spring stiffness, and by a low-pass cut off for $\mu=\pm\pi$. Of significance is the asymmetry about the $\mu=0$ axis, which is generated by the modulation. This leads to asymmetric bandgaps seen in Fig. 2(b), which result in the ability of the lattice to support one directional wave motion. For frequencies within one of these bandgaps, frequency is only defined for either a positive or negative wavenumber value which corresponds to waves propagating in the forward or backward direction only. For example, for $\Omega\in[2.05,2.15]$ the band diagram displays only a branch with negative slope, i.e., negative group velocity and thus excitation of the lattice within this frequency range will result in a wave propagating only in the $-x$ direction. Thus, lattices with periodic time-varying properties are characterized by non-reciprocal behavior for frequencies belonging to non-symmetric bandgaps.

We now compare the dispersion curves obtained using our Bloch based procedure with the dispersion data extracted from explicit numerical simulations using the SED method. The SED results are obtained by computing the transient response of the system to harmonic excitation at frequency $\Omega$, which is imposed as a prescribed displacement on one of the masses. The dispersion diagram is obtained through the excitation over a range of frequencies with each frequency excited individually. The velocity field obtained in each simulation is then used to compute the spectral energy density $\Phi(\Omega,\mu)$, which is shown as a contour plot to trace the dispersion diagram. In our calculation, we considered a closed system consisting of $70$ unit cells, whereby Born-Karman boundary conditions are enforced to avoid edge reflections [23]. The results, displayed in Fig. 3, show the excellent agreement between the SED method and the procedure presented in this paper.

The next example considers the modulation of the ground stiffness, which is expressed as

[TABLE]

while the inter-mass stiffnesses are kept constant in time, i.e. $\beta_{m}=0$. The corresponding stiffness coefficients are given by

[TABLE]

The considered lattice consists of unit cells with $R=3$ masses. Results are obtained for ${\Omega}_{m}=0.2$, $\beta_{g,m}=0.15$, $\gamma_{g}=2$. For this set of parameters, the obtained dispersion branches are shown in Fig. 3(a), while the results obtained with the SED method are presented in Fig. 3(b). The results, displayed in Fig. 3, show the excellent agreement between the SED method and the procedure presented in this paper.

3.2 Square stiffness modulation

Having demonstrated the effectiveness of our procedure in predicting dispersion diagrams for chains with harmonic modulation of stiffness, we now present a more complex example where the modulation is a traveling square wave and is thus composed of multiple harmonic components. The case of stiffness modulation imposed on the inter-mass springs is investigated for illustration purposes. The spring constant of the $s$-th spring is expressed as

[TABLE]

where $H$ is the Heaviside function. A schematic of the modulated spring-mass chain is shown in Fig. 4.

Similar to the harmonic modulation, the stiffness modulation square wave travels with phase velocity $V_{m}=\omega_{m}\lambda_{m}/2\pi$. The stiffness coefficients in Eqn. (9) are now given by

[TABLE]

where

[TABLE]

The dispersion diagram of a lattice with modulation wavelength $\lambda_{m}=3a$ ($R=3$) and stiffness parameters ${\Omega}_{m}=0.25$, $\beta_{m}=0.05$, $\gamma_{g}=1$ is shown in Fig. 5(a), while the SED calculation is presented in Fig. 5(b). The number of terms used in the series expressing the stiffness modulation of the structure in Eqn. (26) is $2P+1=31$ ($P=15$). Since the square wave is a superposition of multiple frequencies, a larger number of terms than in the harmonic stiffness modulation case are required to obtain an accurate dispersion diagram for this system. The convergence rate of the Fourier series with respect to the number of terms is slow due to the Gibbs phenomenon when representing the square wave function in the Fourier basis.

The square-wave modulated system also exhibits asymmetric dispersion behavior illustrating that the system violates time reversal symmetry. Note that the dispersion diagram is different from the harmonically modulated chain and displays several gaps in each branch instead of one due to the presence of multiple frequencies in the modulation. A comparison of Figs. 5(a) and 5(b) shows that the proposed method predicts the gaps associated with stiffness modulation accurately. Finally, we note that the weighting process shows some traces of branches at low frequencies and low wavenumbers. There are two parallel branches to the actual dispersion curve and these arise due to the corresponding components of the eigenvector having a non-negligible amplitude. These additional branches may filtered out through a more aggressive threshold, which would need to be adapted in the case of multiple component modulations. Despite these minor traces, good agreement in Fig. 5 indicates that the Bloch-based analysis accurately predicts the dispersion diagrams where the modulation involves multiple frequencies.

4 Non-reciprocal behavior of time-modulated lattices

Upon observation of the time reversal symmetry breaking properties of modulated lattices, we now illustrate their implications on wave motion and dynamic response of both infinite and finite lattices. The discussion first illustrates the effects of the modulation parameters on the observed asymmetric gaps and relates them to the bandgaps observed in space-periodic, time-invariant systems. Next, numerical simulation of the transient response of modulated lattices illustrate the one-directional wave properties of this class of lattices, which are predicted by the analysis of dispersion. Finally, the study of a finite lattice shows how frequency modulation leads to the non-reciprocal steady-state response, which can be also predicted from observations made on the dispersion diagrams.

4.1 Analysis of symmetry and asymmetry of the dispersion diagram

The analysis of symmetries and their lack of thereof in the dispersion diagram due to spatio-temporal modulations is conducted with reference to the spring-mass lattice with harmonic ground stiffness modulation, whose parameters and stiffness properties are described in Section 3 and Eqn. (22). For reference purposes we first consider the case where only spatial modulation exists, i.e. when $\Omega_{m}=0$. This leads to a time-invariant periodic system, which is known to feature a series of bandgaps whose frequency location and width are defined by the amplitude and period of the modulation. In this case, the period is $\lambda_{m}=3a$, while the modulation amplitude is $\beta_{m}=0.2$. The dispersion diagram for this lattice, shown in Fig. 6(a), shows the presence of two bandgaps, shaded in red and blue for visualization purposes, which are commonly expected to occur for a system with periodic spatial variations in mechanical properties. The diagrams are plotted over an extended wavenumber range, $\mu\in[-3\pi,+3\pi]$, to illustrate its periodicity as well as to highlight its symmetry with respect to the $\mu=0$ axis. In the presence of a modulation that occurs both in space and in time, defined in this case by a modulation frequency $\Omega_{m}=0.2$, the symmetry of the dispersion diagrams is broken as shown in Figure 6(b). Two asymmetric gaps occur at different frequencies for waves traveling to the left or to the right. The difference in the central frequency of the gaps, $\Delta\Omega=\Omega_{m}=0.2$, is thus exactly equal to the value of the modulation frequency. This statement, which was proven analytically in [22] for continouus media, can be explained intuitively by observing that the impedance mismatch encountered by the propagating wave travels at velocity $v_{m}=\omega_{m}/\kappa_{m}$, so the rate at which the same wave perceives it is different when the wave travels to the left or to the right, by a difference that is the velocity of the modulation wave. Another perspective to this observation is to apply a change in variable $\bar{x}=x-(\omega_{m}/\kappa_{m})t$ to the space-time modulated lattice, which recovers a space-only modulated lattice, but in the new coordinate system where the frequency of any traveling wave is shifted by a quantity $\omega_{m}$. One interesting observation is that the impedance mismatch remains the same, therefore the gaps have the same width, but occur at shifted frequencies. Similar observations are reported by Cassedy and Oliner [2] in continuous systems subjected to space-time periodic modulation, where a proportional relation between the difference in frequency in the gaps ($\Delta\Omega$) and the phase velocity of the modulation wave $v_{m}$ was also derived.

4.2 One-directional wave motion in modulated lattices

The implications of asymmetric dispersion properties in terms of wave motion can be illustrated through numerical simulations conducted on lattices with the considered spatio-temporal modulations. These lattices violate mechanical reciprocity by allowing one-directional wave propagation. Consider a lattice with harmonic modulation of the first-neighbor interaction springs as discussed in Sec. 3.1, having the properties $\Omega_{m}=0.2$, $\beta_{m}=0.15$, $\gamma_{g}=1$. The dispersion diagram of this lattice exhibits an asymmetric bandgap in the branch with positive slope for $\Omega\in[2.05,2.15]$, see Fig. 2(b). At frequencies in this bandgap, only modes with negative group velocities are allowed, which means that the lattice allows propagation in the $-x$ direction. The response of the lattice is evaluated for transient excitation over one mass at a frequency inside this bandgap, $\Omega=2.1$. For the numerical simulations, a lattice with $70$ unit cells and a modulation period $\lambda_{m}=3a$ is considered, which corresponds to a total of $210$ masses. Reflections at the edges are avoided by applying Born-Karman boundary conditions. Figure 7 displays the evolution of displacement along the chain against the normalized time $\tau=\omega_{0}t$. In the simulations, the center mass is excited for a duration of $\omega_{0}t=200$. The total simulation time is $\omega_{0}t=600$. For the integration, we used Verlet algorithm [24] with a time step $\Delta t=2\times 10^{-8}$. The resulting displacement field is strongly asymmetric, as most of the energy of wave motion is supported by waves that travel along the $-x$ direction. Some waves with frequency different from $\omega$ are also excited and they travel in both directions with group velocities different from the primary wave. Their amplitude is small compared to the primary wave, as can be observed in Fig. 7. Note that the dispersion diagrams only predict the wavenumbers and group velocities of waves traveling with the excitation frequency $\omega$ and the transient response of the chain is indeed consistent with the dispersion analysis.

4.2.1 Steady state response of a finite modulated lattice: acoustic circulator

Finally, we consider the free out-of-plane vibration response of a three mass system arranged in a circle as illustrated in Fig. 8. The masses are all identical and equal to $m$. Without loss of generality, the ground spring stiffness is set to $k_{g}=0$.

The spring stiffness has two components: a constant part and a time varying modulation component. A traveling wave modulation of frequency $\omega_{m}$ and amplitude $k_{m}$ is applied to the springs connecting the masses and their stiffness is $k_{r}=k_{0}+k_{m}\cos(\omega_{m}t+2\pi r/3)$. The governing equations of this system read $\bm{M}\ddot{\bm{u}}+\bm{\hat{K}}\bm{u}=\mathbf{0}$. They can be written in dimensionless form as $\ddot{\bm{u}}+\bm{K}\bm{u}=\mathbf{0}$, where the stiffness matrix $\bm{K}$ is

[TABLE]

and $c.c$ stands for complex conjugate. To solve the system of equations, we follow a procedure similar to that of Alu and coworkers [16]. The following change of variables is introduced to diagonalize the above system

[TABLE]

Note that we chose a basis that is not orthogonal because, as demonstrated below, the solution is expressed conveniently in the above basis in the presence of modulation. Using the transformation in Eqn. (28) leads to the following equation for $\bm{a}=[a_{1}\;a_{2}\;a_{3}]^{T}$

[TABLE]

where

[TABLE]

Note that the $a_{1}$ mode becomes uncoupled from the $a_{2}$ and $a_{3}$ modes. To solve for these modes, we assume a solution of the form

[TABLE]

where $\omega$ is an unknown frequency and $A_{2},A_{3}$ are scalar constants. Substituting the above equation into the system of equations (29) results in a homogeneous system of equations and the condition for the existence of non-trivial solutions leads to

[TABLE]

Solving the above system leads to the frequencies

[TABLE]

Substituting values $\beta_{m}=0.09$ and $\Omega_{m}=0.5$ for the modulation parameters, we find that $A_{2}\gg A_{3}$ for the mode $\Omega_{1}$ and $A_{3}\gg A_{2}$ for the other mode $\Omega_{2}$. Thus transforming back to the original displacement basis using the inverse of the transform in Eqn. (28), and considering only the term dominant between $A_{2}$ and $A_{3}$ in each of the modes, the solution is expressed as

[TABLE]

Note that the frequencies split as a consequence of modulation from $\Omega=\sqrt{3}$ and the $a_{2}$ and $a_{3}$ modes are called clockwise and counter-clockwise propagating modes, respectively, based on the sign of the group velocity.

We now interpret the above results using the dispersion analysis of an infinite chain of masses subjected to the same modulation. The normalized modulation stiffness $\beta_{m}$ and frequency $\Omega_{m}$ are set to $0.09$ and $0.5$, respectively. Figure 9(a) displays the dispersion diagram for a chain having a modulation wavelength $\lambda_{m}=3a$. Similar to the previous case, there are $3$ branches in the dispersion diagram as a unit cell has three degrees of freedom. Also, the presence of modulation breaks the symmetry about the $\mu=0$ axis, where $\mu=3a\kappa$. Figure 9(b) displays a zoomed-in view of the dispersion diagram where the two branches intersect close to $\mu=0$. We observe that the branches do not intersect at $\mu=0$ as a result of modulation.

The steady state free vibration response corresponds to a constraint $u_{i}=u_{i+3}$ imposed on an infinite chain. This constraint is satisfied for wavenumbers $\mu$ which have $e^{i\mu}=1$. This condition is equivalent to $\mu/2\pi=n$, where $n$ is an integer and is equivalent to $\mu=0$ in the First Brillouin zone. Figure 9 displays that there are 3 values of frequency $\Omega$ at wavenumber $\mu=0$. One frequency is $\Omega=0$, while there are two non-zero frequencies. In the absence of modulation, these frequencies would be $\Omega=\sqrt{3}$. Modulation results in a split of frequencies, and the two frequencies are derived in Eqn. (32). From the dispersion diagram, we infer that the wave with frequency $\Omega_{2}-\Omega_{m}/2=1.7285$ is a right traveling wave while the wave with frequency $\Omega_{1}+\Omega_{m}/2=1.7347$ is a left traveling wave. The dispersion analysis thus predicts accurately the frequency splitting due to stiffness modulation in a 3-mass system. This frequency splitting is the basis for designing non-reciprocal acoustic circulators. When the system is excited at the frequency $\Omega=1.7315$, the combination of the $a_{2}$ and $a_{3}$ modes results in small ratio $A_{3}/A_{2}$, as demonstrated by Alu and coworkers [16], see frequency response function in Fig. 10 (damping coefficient was set to c=0.01). This phenomena can be exploited to build a nonreciprocal mechanical $3$-port device.

5 Conclusions

The paper presents a Bloch-based method to study the dispersion properties of systems with periodic time-varying coefficients. The method allows us to obtain the dispersion diagram of the lattice by solving a quadratic eigenvalue problem over the unit cell, whose length depends on the spatial modulation wavelength. The applicability of our method is demonstrated for both harmonic modulation and modulation having multiple frequencies, and it is validated with dispersion diagrams obtained from numerical simulations using SED method. Our method shows excellent agreement with the dispersion diagrams obtained from full numerical simulations for all the considered cases and predicts unidirectional wave propagation at certain frequencies. Potential future research directions include extending our formulation to higher dimensions and wave propagation in other physical systems, as well as designing lattices with specific nonreciprocal mechanical wave propagation characteristics.

6 Acknowledgements

The authors are indebted to the US Army Research Office (Grant number W911NF1210460), the US Air Force Office of Scientific Research (Grant number FA9550-13-1-0122), the University Carlos III de Madrid and the Ministerio de Ciencia e Innovación de España (Project DPI/2014-57989-P) for the financial support.