Waves in composites and metamaterials/Acoustic metamaterials and negative moduli

The content of these notes is based on the lectures by Prof. Graeme W. Milton (University of Utah) given in a course on metamaterials in Spring 2007.

Acoustic Metamaterials
Recall that for isotropic materials, the stress is given by

\boldsymbol{\sigma} = \boldsymbol{\mathsf{C}}:\boldsymbol{\nabla}\mathbf{u} = \mu~[\boldsymbol{\nabla}\mathbf{u} + (\boldsymbol{\nabla}\mathbf{u})^T] + \lambda~(\boldsymbol{\nabla} \cdot \mathbf{u})~\boldsymbol{\mathit{1}} ~. $$

For a fluid, the shear modulus $$\mu = 0$$ and the Lame constant $$\lambda = \kappa - 2/3~\mu = \kappa$$, where $$\kappa$$ is the bulk modulus. Therefore, for fluids,
 * $$ \text{(1)} \qquad

\boldsymbol{\sigma} = \boldsymbol{\mathsf{C}}:\boldsymbol{\nabla}\mathbf{u} = \kappa~(\boldsymbol{\nabla} \cdot \mathbf{u})~\boldsymbol{\mathit{1}} = -p~\boldsymbol{\mathit{1}} ~. $$ The quantity
 * $$ \text{(2)} \qquad

p = -\kappa~(\boldsymbol{\nabla} \cdot \mathbf{u}) $$ is the pressure in the fluid.

Taking the Fourier transform of equations (1) and (2) we get
 * $$ \text{(3)} \qquad

\widehat{\boldsymbol{\sigma}} = -\widehat{p}~\boldsymbol{\mathit{1}} \qquad \text{and} \qquad \widehat{p} = -\widehat{\kappa}(\omega)~(\boldsymbol{\nabla} \cdot \widehat{\mathbf{u}}) ~. $$

Also recall that, for low frequency processes,
 * $$ \text{(4)} \qquad

\boldsymbol{\nabla} \cdot \widehat{\boldsymbol{\sigma}} + \omega^2~\widehat{\rho}(\omega)~\widehat{\mathbf{u}}(\omega) = \boldsymbol{0} ~. $$ For fluids, equations (4) becomes
 * $$ \text{(5)} \qquad

\boldsymbol{\nabla} \widehat{p} = \omega^2~\widehat{\rho}(\omega)~\widehat{\mathbf{u}}(\omega) ~. $$

From equations (3) and (5), we have

\boldsymbol{\nabla} \cdot \widehat{\mathbf{u}} = \boldsymbol{\nabla} \cdot \left(\cfrac{1}{\omega^2~\widehat{\rho}}~\boldsymbol{\nabla} \widehat{p}\right) = -\cfrac{\widehat{p}}{\widehat{\kappa}} ~. $$ Therefore, we get the  acoustic wave equation

{ \boldsymbol{\nabla} \cdot \left(\cfrac{1}{\widehat{\rho}}~\boldsymbol{\nabla} \widehat{p}\right) + \cfrac{\omega^2}{\widehat{\kappa}}~\widehat{p} = 0 ~. } $$ If $$\widehat{\rho}$$ and $$\widehat{\kappa}$$ only depend on $$x_2$$ and $$x_3$$ and $$\widehat{p}$$ is independent of $$x_1$$, then the three-dimensional gradient operator $$\boldsymbol{\nabla}$$ can be replaced with the two-dimensional gradient operator $$\overline{\boldsymbol{\nabla}}$$, and we get

{ \overline{\boldsymbol{\nabla}} \cdot\left(\cfrac{1}{\widehat{\rho}(\omega)}~\overline{\boldsymbol{\nabla}} \widehat{p}\right) + \cfrac{\omega^2}{\widehat{\kappa}(\omega)}~\widehat{p} = 0 } $$ and

\widehat{\mathbf{u}}(\omega) = \cfrac{1}{\omega^2~\widehat{\rho}(\omega)}~\overline{\boldsymbol{\nabla}} \widehat{p} ~. $$ Therfore, by analogy with the results from antiplane shear elasticity and TE electromagnetism, if $$\widehat{\rho}$$ and $$\widehat{\kappa}$$ are both negative, we get a negative refractive index material for acoustics.

The speed of sound in an acoustic medium is given by $$c = \sqrt{\kappa/\rho}$$. The sound speed is imaginary is either $$\kappa$$ or $$\rho$$ is negative and a material with these properties appears opaque to sound waves. However, if $$\kappa$$ and $$\rho$$ are both positive or both negative, the sound speed is real and the medium becomes transparent to acoustic waves, i.e., acoustic waves can propagate through the medium.

Let us now consider the propagation of acoustic waves across the interface shown in Figure 1. If the pressure is the same at points at equal distances from the interface as shown, the vectors $$\boldsymbol{\nabla} p$$ are reflections of each other. If $$\mathbf{n}$$ is the normal to the interface, the quantity $$\mathbf{n}\cdot\boldsymbol{\nabla} p$$ changes sign across the interface. Since $$\rho$$ also changes sign, $$\mathbf{n}\cdot\mathbf{u}$$ is continuous across the interface. Hence, the boundary condition across the interface is physical.

Since this situation is analgous to the one we observed earlier for electromagnetism, the effect of a slab of negative $$\kappa,\rho$$ material will be to translate a source of acoustic waves and the slab will act like a perfect lens.

Negative Elastic Moduli
Early work on negative elastic moduli can be found in Lakes01. More recent work on ultrasonic metamaterials can be found in Fang06. Consider the array of Helmholtz resonators shown in Figure 2 (an example of such a resonator in everyday life is a soda bottle which resonates when you blow over the top of the neck.)  The resonator can be thought of as a spring-mass system where the air inside the cavity acts as a spring and the water in the narrow neck acts as a mass. There is some frequency at which the spring-mass system resonates.

A model of the Helmholtz resonator
A model of the Helmholtz resonator is shown in Figure 3. For simplicity, we assume that each cavity has a square cross section as do the piston arms. The cross sectional area is assumed to be $$d^2$$. The air in each cavity is model with a spring of complex spring constant $$K$$ and the water in the neck is modeled as a rigid body of mass $$m$$. The piston is filled with a compressible fluid with complex bulk modulus $$\kappa$$.

Let the force applied on the pistons be $$F(t)$$. Then the pressure in the fluid is
 * $$ \text{(6)} \qquad

p(t) = \cfrac{F(t)}{d^2} ~. $$ Since the fluid transmits the pressure to the mass and the area of cross section of the cavity is $$d^2$$, the force applied by the fluid on the mass is also $$F(t)$$.

Let $$f(t)$$ be the force that the spring applies on the mass. By symmetry, the same forces is applied to both cavities shown in Figure 3. Let $$x(t)$$ denote the position of the piston and let $$y(t)$$ be the position of the mass.

Assume harmonic time dependence of the quantities $$F(t)$$, $$f(t)$$, $$x(t)$$, and $$y(t)$$. Then
 * $$ \text{(7)} \qquad

F(t) = \text{Re}(\widehat{F}~e^{-i\omega t}) ~; f(t) = \text{Re}(\widehat{f}~e^{-i\omega t}) ~; x(t) = \text{Re}(\widehat{x}~e^{-i\omega t}) ~; y(t) = \text{Re}(\widehat{y}~e^{-i\omega t}) ~. $$

Assume that $$p = 0$$ at $$t = 0$$. Then, we also have
 * $$ \text{(8)} \qquad

p(t) = \text{Re}(\widehat{p}~e^{-i\omega t}) ~. $$ Newton's law then implies that
 * $$ \text{(9)} \qquad

m~\cfrac{d^2 y}{d t^2} = F(t) - f(t) ~. $$ Substituting equations (7) into equation (9), we get
 * $$ \text{(10)} \qquad

-m~\omega^2~\widehat{y} = \widehat{F} - \widehat{f} ~. $$ Also, from Hooke's law
 * $$ \text{(11)} \qquad

f(y) = K~y(t) \qquad \implies \qquad \widehat{f} = K~\widehat{y} ~. $$ Combining equations (10) and (11) we get
 * $$ \text{(12)} \qquad

\widehat{F} = (K - m~\omega^2)~\widehat{y} ~. $$ Also, from equation (6) we have
 * $$ \text{(13)} \qquad

\widehat{p} = \cfrac{\widehat{F}}{d^2} ~. $$ Combining equations (13) and (12) we have
 * $$ \text{(14)} \qquad

\widehat{p} = \cfrac{(K - m~\omega^2)}{d^2}~\widehat{y} ~. $$ The change in volume of the fluid due to the motions of the pistons and the masses is given by
 * $$ \text{(15)} \qquad

\Delta V(t) = 2~d^2~[y(t) - x(t)]~. $$ Now, from the constitutive relation for the fluid,
 * $$ \text{(16)} \qquad

p(t) = -\kappa~\cfrac{\Delta V(t)}{V_0} = -\kappa~\cfrac{2~d^2~[y(t) - x(t)]}{V_0} ~. $$ where $$\kappa$$ is the complex bulk modulus of the fluid and $$V_0$$ is the initial volume.

Substituting equations (7) and (8) into equation (16) we get
 * $$ \text{(17)} \qquad

\widehat{p} = \cfrac{2~\kappa~d^2~(\widehat{x} - \widehat{y})}{V_0} ~. $$ From equations (14) and (17), eliminating $$\widehat{y}$$, we get

\widehat{p} = \cfrac{2~\kappa~d^2}{V_0} \left(\widehat{x} - \cfrac{d^2}{K - m~\omega^2}~\widehat{p}\right) ~. $$ Solving for $$\widehat{p}$$, we have

\widehat{p} = \cfrac{2~\kappa~d^2}{V_0} \left[1 + \cfrac{2~\kappa~d^4}{V_0(K - m~\omega^2)}\right]^{-1} ~\widehat{x} ~. $$ Therefore,
 * $$ \text{(18)} \qquad

\widehat{F} = d^2~\widehat{p} = \cfrac{2~\kappa~d^4}{V_0} \left[1 + \cfrac{2~\kappa~d^4}{V_0(K - m~\omega^2)}\right]^{-1} ~\widehat{x} ~. $$ By definition, the Young's modulus relates the stress to the strain. For our model this implies that

\cfrac{\widehat{F}}{d^2} = Y~\cfrac{\widehat{\Delta L}}{L_0} $$ where $$Y$$ is a complex Young's modulus, $$\Delta L(t)$$ is the change in length and $$L_0$$ is the initial length of the region between the two pistons. Therefore,
 * $$ \text{(19)} \qquad

\widehat{F} = \cfrac{d^2}{L_0}~Y~(2~\widehat{x}) ~. $$ From equations (18) and (19) we can deduce an expression for the Young's modulus of the form
 * $$ \text{(20)} \qquad

{ Y(\omega) = \cfrac{\kappa~d^2~L_0}{V_0} \left[1 + \cfrac{2~\kappa~d^4}{V_0(K - m~\omega^2)}\right]^{-1}~. } $$ If in particular $$\kappa$$ and $$K$$ are real (purely elastic springs with not damping) and positive, then a plot of $$Y$$ as a function of $$\omega$$ has the form shown in Figure 4. So it will appear that the system will have a negative Young's modulus for frequencies higher than $$\sqrt{K/m}$$ !

In fact, several other negative modulus materials can be envisaged.

Why are negative moduli a surprise?
It is common for people to use a generalized Maxwell model in linear viscoelasticity (see Figure 5). The question is: what is the relation between $$F$$ and $$u$$?

Recall that for a single Maxwell element ($$j$$), the displacements in the spring and the dashpot are given by

u_1(t) = \cfrac{f_j(t)}{K_j} \qquad \text{and} \qquad \cfrac{d^2 u_2(t)}{d t^2} = \cfrac{f_j(t)}{\eta_j} ~. $$ Fourier transforming these equations gives

\widehat{u}_1(\omega) = \cfrac{\widehat{f}_j(\omega)}{K_j} \qquad \text{and} \qquad i\omega~\widehat{u}_2(\omega) = \cfrac{\widehat{f}_j(\omega)}{\eta_j} ~. $$ The total displacement of the Maxwell element is

\widehat{u}_j(\omega) = \widehat{u}_1 + \widehat{u}_2 = \left[\cfrac{1}{K_j} -\cfrac{i}{\omega~\eta_j}\right]~\widehat{f}_j(\omega)~. $$ Dropping the hats, we have

u_j = \left(\cfrac{1}{K_j} - \cfrac{i}{\omega~\eta_j}\right)~f_j ~. $$ Therefore, from the balance of forces (for the generalized Maxwell model)

F = f_0 + \sum f_j = \left[K_0 + \sum_j \left(\cfrac{1}{K_j} - \cfrac{i}{\omega~\eta_j}\right)^{-1}\right]~u ~. $$ The quantity

K_{\text{eff}} = K_0 + \sum_j \left(\cfrac{1}{K_j} - \cfrac{i}{\omega~\eta_j}\right)^{-1} = K_0 + \sum_j \cfrac{K_j}{1 - i/\omega\tau_j} $$ is the effective spring constant for the model ($$\tau_j = \eta_j/K_j$$ is called the relaxation time).

Similarly, for the Young's modulus and converting the sum into an integral, we have

Y_{\text{eff}}(\omega) = Y_0 + \int_0^{\infty} \cfrac{H(\tau)}{1 - i/\omega\tau}~\text{d}\tau $$ where $$\tau$$ is the relaxation spectrum, $$H(\tau) \ge 0$$, and $$Y_0 \ge 0$$. Therefore,

Y_{\text{eff}}(\omega) = Y_0 + \int_0^{\infty} \cfrac{H(\tau)~(1 + i/\omega\tau)}{1 + 1/\omega^2\tau^2}~\text{d}\tau ~. $$ If we discretize the spectrum we get

Y_{\text{eff}}(\omega) = Y_0 + \sum_j \cfrac{H_j~(1 + i/\omega\tau_j)}{1 + 1/\omega^2\tau_j^2} $$ where $$H_j \ge 0$$, $$\tau_j > 0$$, and $$Y_0 \ge 0$$. This implies that $$\text{Im}(Y_{\text{eff}}(\omega)) \ge 0$$ and $$\text{Re}(Y_{\text{eff}}(\omega)) \ge 0$$ for all $$\omega$$. Clearly, viscoelastic models cannot represent negative elastic moduli and can therefore fail badly when used to model materials such as Helmholtz resonators.

In fact, for viscoelastic material models, we find that $$Y(\omega)$$ is analytic in the entire complex plane except for isolated poles (at intervals of $$\omega\tau$$) in the negative imaginary axis. This is contrary to what we observe for the frequency dependent models we have dealt with so far for which causality forces the moduli to be analytic only the upper half complex $$\omega$$ plane.