Waves in composites and metamaterials/Elastodynamics and electrodynamics

< Waves in composites and metamaterials

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.

Dissipation

Recall from the previous lecture that the average rate of work done in a cycle of oscillation of material with frequency dependent mass is


  \begin{align}
     \mathcal{P} & = \cfrac{\omega}{2~\pi}
       \int_{0}^{2\pi/\omega}~\mathbf{F}(t)\cdot\mathbf{V}(t)~\text{d}t  \\
       & = \cfrac{\text{Re}(\widehat{\mathbf{F}})\cdot\text{Re}(\widehat{\mathbf{V}}) +
                  \text{Im}(\widehat{\mathbf{F}})\cdot\text{Im}(\widehat{\mathbf{V}})}{2} \\
       & = \omega~[\text{Re}(\widehat{\mathbf{V}})\cdot\text{Im}[\boldsymbol{M}(\omega)]\cdot\text{Re}(\widehat{\mathbf{V}}) + 
                   \text{Im}(\widehat{\mathbf{V}})\cdot\text{Im}[\boldsymbol{M}(\omega)]\cdot\text{Im}(\widehat{\mathbf{V}})].
  \end{align}

This quadratic form will be non-negative for all choices of \widehat{\mathbf{V}} if and only if \text{Im}(\boldsymbol{M}(\omega)) is positive semidefinite for all real \omega>0. Therefore, a restriction on the behavior of such materials is that


  {
  \text{Im}[\boldsymbol{M}(\omega)] \ge 0 ~.
  }

Similarly, for electrodynamics, the average power dissipated into heat is given by

 \text{(1)} \qquad 
  \mathcal{P} = \cfrac{\omega}{2~\pi}
       \int_{0}^{2\pi/\omega}~\left[\mathbf{E}(t)\cdot\frac{\partial \mathbf{D}(t)}{\partial t} + 
         \mathbf{H}(t)\cdot\frac{\partial \mathbf{B}(t)}{\partial t}\right]~dt ~.

In this case, the quantity \mathbf{E} is equivalent to the voltage and the quantity rate of change of electrical displacement \partial\mathbf{D}/\partial t is equivalent to the current (recall that in electrostatics the power is given by \mathcal{P} = V~I). In addition, we also have a contribution due to magnetic induction.

Let us assume that the fields can be expressed in harmonic form, i.e.,


  \mathbf{E}(t) = \text{Re}[\widehat{\mathbf{E}}~e^{-i\omega t}] ~;~~
  \mathbf{H}(t) = \text{Re}[\widehat{\mathbf{H}}~e^{-i\omega t}]

or equivalently as


  \mathbf{E}(t) = \text{Re}(\widehat{\mathbf{E}})~\cos(\omega t) + 
                \text{Im}(\widehat{\mathbf{E}})~\sin(\omega t) ~;~~
  \mathbf{H}(t) = \text{Re}(\widehat{\mathbf{H}})~\cos(\omega t) + 
                \text{Im}(\widehat{\mathbf{H}})~\sin(\omega t)~.

Also, recall that,


  \frac{\partial \mathbf{D}}{\partial t} = \boldsymbol{\nabla} \times \mathbf{H} = -i\omega~\boldsymbol{\epsilon}\cdot\mathbf{E}(t) 
  \qquad \text{and} \qquad
  \frac{\partial \mathbf{B}}{\partial t} = -\boldsymbol{\nabla} \times \mathbf{E} = -i\omega~\boldsymbol{\mu}\cdot\mathbf{H}(t) ~.

Therefore, for real \omega and real \mathcal{P}, we can write equation (1) as (with the substitution z = \omega t),


  \begin{align}
  \mathcal{P} = \cfrac{\omega}{2~\pi}
       \int_{0}^{2\pi}~ & 
       \left[\text{Re}(\widehat{\mathbf{E}})~\cos z + \text{Im}(\widehat{\mathbf{E}})~\sin z\right]\cdot
       \left\{\text{Im}(\boldsymbol{\epsilon})\cdot
       \left[\text{Re}(\widehat{\mathbf{E}})~\cos z + \text{Im}(\widehat{\mathbf{E}})~\sin z\right]\right\} + \\
        &
       \left[\text{Re}(\widehat{\mathbf{H}})~\cos z + \text{Im}(\widehat{\mathbf{H}})~\sin z\right]\cdot
       \left\{\text{Im}(\boldsymbol{\mu})\cdot
       \left[\text{Re}(\widehat{\mathbf{H}})~\cos z + \text{Im}(\widehat{\mathbf{H}})~\sin z\right]\right\}~\text{d}z
  \end{align}

Expanding out, and using the fact that


  \int_{0}^{2\pi} \cos^2 z~\text{d}z = \int_{0}^{2\pi} \sin^2 z~\text{d}z = \pi \quad
  \text{and} \quad \int_{0}^{2\pi} \sin(2 z)~\text{d}z = 0

we have,

 \text{(2)} \qquad 
  \mathcal{P} = \cfrac{\omega}{2}
       \left[\text{Re}(\widehat{\mathbf{E}})\cdot\text{Im}(\boldsymbol{\epsilon})\cdot\text{Re}(\widehat{\mathbf{E}}) + 
             \text{Im}(\widehat{\mathbf{E}})\cdot\text{Im}(\boldsymbol{\epsilon})\cdot\text{Im}(\widehat{\mathbf{E}}) + 
             \text{Re}(\widehat{\mathbf{H}})\cdot\text{Im}(\boldsymbol{\mu})\cdot\text{Re}(\widehat{\mathbf{H}}) + 
             \text{Im}(\widehat{\mathbf{H}})\cdot\text{Im}(\boldsymbol{\mu})\cdot\text{Im}(\widehat{\mathbf{H}})\right] ~.

Since \omega \ge 0 and the power \mathcal{P} \ge 0, the quadratic forms in equation (2) require that


  \text{Im}(\boldsymbol{\epsilon}) > 0  \qquad \text{and} \qquad \text{Im}(\boldsymbol{\mu}) > 0 ~.

Note that if the permittivity is expressed as


   \boldsymbol{\epsilon} = \boldsymbol{\epsilon}_0 + i\cfrac{\boldsymbol{\sigma}}{\omega}

the requirement \text{Im}(\boldsymbol{\epsilon}) > 0 implies that the conductivity \boldsymbol{\sigma} > 0. Therefore, if the conductivity is greater than zero, there will be dissipation.

Brief introduction to elastodynamics

A concise introduction to the theory of elasticity can be found in Atkin80. In this section, we consider the linear theory of elasticity for infinitesimal strains and small displacements.

Consider the body (\Omega) shown in Figure~1. Let \Gamma be a subpart of the body (in the interior of \Omega or sharing a part of the surface of \Omega). Postulate the existence of a force \mathbf{t}(\mathbf{x},\mathbf{n}) per unit area on the surface of \Gamma where \mathbf{n} is the outward unit normal to the surface of \Gamma. Then \mathbf{t} is the force exerted on \Gamma by the material outside \Gamma or by surface tractions.

Figure 1. Illustration of the concept of stress.

From the balance of forces on a small tetrahedron (\Gamma), we can show that \mathbf{t}(\mathbf{x},\mathbf{n}) is linear in \mathbf{n}. Therefore,


  \mathbf{t} = \boldsymbol{\sigma}\cdot\mathbf{n}

where \boldsymbol{\sigma} is a second-order tensor called the stress tensor.

Since the tetrahedron cannot rotate at infinite velocity as its size goes to zero (conservation of angular momentum), we can show that the stress tensor is symmetric, i.e.,

 
  \boldsymbol{\sigma} = \boldsymbol{\sigma}^T ~.

In particular, for a fluid,


  \boldsymbol{\sigma} = - p~\boldsymbol{\mathit{1}}

where p is the pressure.

Let us assume that the stress depends only on the strain (and not on strain gradients or strain rates), where the strain is defined as

 \text{(3)} \qquad 
  \boldsymbol{\epsilon} = \frac{1}{2}\left[\boldsymbol{\nabla}\mathbf{u} + (\boldsymbol{\nabla}\mathbf{u})^T\right] ~.

Here \mathbf{u}(\mathbf{x},t) is the displacement field. Note that a gradient of the displacement field is used to define the strain because rigid body motions should not affect \boldsymbol{\sigma} and a rigid body rotation gives zero strains (for small displacements).

Assume that \boldsymbol{\sigma} depends linearly on \boldsymbol{\epsilon} so that


  \boldsymbol{\sigma}(\mathbf{x},t) = \int \text{d}\mathbf{x}'~
    \left[\int \boldsymbol{\mathsf{K}}_{\varepsilon}(\mathbf{x},\mathbf{x}',t' - t):\boldsymbol{\epsilon}(\mathbf{x},t')~\text{d}t'\right] ~.

Note that this assumption ignores preexisting internal stresses such as those found in prestressed concrete. If the material can be approximated as being local, then

 \text{(4)} \qquad 
  \boldsymbol{\sigma}(\mathbf{x},t) = \int \boldsymbol{\mathsf{K}}_{\varepsilon}(\mathbf{x},t' - t):\boldsymbol{\epsilon}(\mathbf{x},t')~\text{d}t' ~.

Taking the Fourier transforms of equation (4), we get

 \text{(5)} \qquad 
  \widehat{\boldsymbol{\sigma}}(\mathbf{x},\omega) = \boldsymbol{\mathsf{C}}(\mathbf{x},\omega):\widehat{\boldsymbol{\epsilon}}(\mathbf{x},\omega)

where


  \boldsymbol{\mathsf{C}}(\mathbf{x},\omega) = \int \boldsymbol{\mathsf{K}}_{\varepsilon}(\mathbf{x},\tau)~e^{-i\omega\tau}~\text{d}\tau ~;~~
  \tau = t' - t ~.

In index notation, equation (5) can be written as

 
  \widehat{\sigma}_{ij} = C_{ijkl}~\widehat{\epsilon}_{kl} ~.

Causality implies that stresses at time t can only depend on strains of previous times, i.e., if t' \le t or \tau \le 0. Therefore,


  \boldsymbol{\mathsf{K}}_{\varepsilon}(\mathbf{x},\tau) = \boldsymbol{\mathsf{0}} \qquad \text{if} \qquad \tau > 0 ~.

This in turn implies that the integral converges only if \text{Im}(\omega) > 0, i.e., \boldsymbol{\mathsf{C}}(\mathbf{x},\omega) is analytic when \text{Im}(\omega) > 0.

In the absence of body forces, the equation of motion of the body can be written as

 \text{(6)} \qquad 
  \boldsymbol{\nabla} \cdot \boldsymbol{\sigma} = \rho(\mathbf{x})~\frac{\partial^2 \mathbf{u}}{\partial t^2}

where \rho is the mass density, \boldsymbol{\nabla} \cdot \boldsymbol{\sigma} is the internal force per unit volume, and \partial^2 \mathbf{u}/\partial t^2 is the acceleration. Hence, this is just the expression of Newton's second law for continuous systems.

For a material which has a frequency dependent mass, equation (6) may be written as

 \text{(7)} \qquad 
  \boldsymbol{\nabla} \cdot \boldsymbol{\sigma} = \int_{-\infty}^{\infty} \rho(\mathbf{x}, t - t')~\frac{\partial^2 \mathbf{u}}{\partial t^2}~\text{d}t'

where causality implies that if t' > t then \rho = 0.

Taking the Fourier transform of equation (7), we get

 \text{(8)} \qquad 
  \boldsymbol{\nabla} \cdot \widehat{\boldsymbol{\sigma}}(\mathbf{x},\omega) = 
     -\omega^2~\widehat{\rho}(\mathbf{x},\omega)~\widehat{\mathbf{u}}(\mathbf{x},\omega)~.

Substituting equation (5) into equation (8) we get

 \text{(9)} \qquad 
  \boldsymbol{\nabla} \cdot (\boldsymbol{\mathsf{C}}:\widehat{\boldsymbol{\epsilon}}) + \omega^2~\widehat{\rho}~\widehat{\mathbf{u}} = 0~.

Also, taking the Fourier transform of equation (3), we have

 \text{(10)} \qquad 
  \widehat{\boldsymbol{\epsilon}} = \frac{1}{2}\left[\boldsymbol{\nabla} \widehat{\mathbf{u}} + (\boldsymbol{\nabla} \widehat{\mathbf{u}})^T\right] ~.

Since \widehat{\boldsymbol{\sigma}} and \widehat{\boldsymbol{\epsilon}} are symmetric, we must have


  C_{ijkl} = C_{jikl} = C_{ijlk} ~.

Because of this symmetry, we can replace \widehat{\boldsymbol{\epsilon}} by \boldsymbol{\nabla} \widehat{\mathbf{u}} in equation (9) to get

 
  \boldsymbol{\nabla} \cdot (\boldsymbol{\mathsf{C}}:\boldsymbol{\nabla \widehat{\mathbf{u}})} + \omega^2~\widehat{\rho}~\widehat{\mathbf{u}} = 0~.

Dropping the hats, we then get the wave equation for elastodynamics

 \text{(11)} \qquad 
  {
  \boldsymbol{\nabla} \cdot (\boldsymbol{\mathsf{C}}:\boldsymbol{\nabla \mathbf{u})} + \omega^2~\rho~\mathbf{u} = 0~.
  }

Antiplane shear

Let us now consider the case of antiplane shear. Assume that \boldsymbol{\mathsf{C}} is isotropic, i.e.,


  \boldsymbol{\mathsf{C}}:\boldsymbol{\nabla}\mathbf{u} = \mu~[\boldsymbol{\nabla}\mathbf{u} + (\boldsymbol{\nabla}\mathbf{u})^T] + \lambda~\text{tr}(\boldsymbol{\nabla}\mathbf{u})~\boldsymbol{\mathit{1}}

where \mu is the shear modulus and \lambda is the Lame modulus.

Let us assume that \mu and \lambda are independent of x_1, i.e.,


  \mu \equiv \mu(x_2,x_3) \qquad \text{and} \qquad
  \lambda \equiv \lambda(x_2,x_3) ~.

Let us look for a solution with u_2 = u_3 = 0 and u_1 independent of x_1, i.e., u_1 \equiv u_1(x_2,x_3). This is an out of plane mode of deformation.

Then, noting that \text{tr}(\boldsymbol{\nabla}\mathbf{u}) = \boldsymbol{\nabla} \cdot \mathbf{u}, we have


  \text{tr}(\boldsymbol{\nabla}\mathbf{u}) = \frac{\partial u_i}{\partial x_i} = 0 ~.

Therefore,


  [\boldsymbol{\mathsf{C}}:\boldsymbol{\nabla}\mathbf{u}]_{ij} = 
    \mu~\left[\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right]

or


  \boldsymbol{\mathsf{C}}:\boldsymbol{\nabla}\mathbf{u} = \mu
    \begin{bmatrix}
       0 & \frac{\partial u_1}{\partial x_2} & \frac{\partial u_1}{\partial x_3} \\
       \\
       \frac{\partial u_1}{\partial x_2} & 0 & 0 \\
       \\
       \frac{\partial u_1}{\partial x_3} & 0 & 0 
    \end{bmatrix} ~.

Therefore


  \boldsymbol{\nabla} \cdot (\boldsymbol{\mathsf{C}}:\boldsymbol{\nabla\mathbf{u})} = 
    \begin{bmatrix}
       \frac{\partial }{\partial x_2}\left(\mu~ \frac{\partial u_1}{\partial x_2}\right) +
           \frac{\partial }{\partial x_3}\left(\mu~ \frac{\partial u_1}{\partial x_3}\right) \\
       0 \\
       0
    \end{bmatrix} ~.

Plugging into the wave equation (11) we get


   \frac{\partial }{\partial x_2}\left(\mu~ \frac{\partial u_1}{\partial x_2}\right) +
       \frac{\partial }{\partial x_3}\left(\mu~ \frac{\partial u_1}{\partial x_3}\right) + \omega^2~\rho~u_1
    = 0

or (using the two-dimensional gradient operator \overline{\boldsymbol{\nabla}})

 \text{(12)} \qquad 
   {
   \overline{\boldsymbol{\nabla}} \cdot(\mu~\overline{\boldsymbol{\nabla} u_1)} + \omega^2~\rho~u_1 = 0 ~. 
   }

TM and TE modes in electromagnetism

Let us now consider the TM (transverse magnetic field) and TE (transverse electric field) modes in electromagnetism and look for parallels with antiplane shear in elastodynamics.

Recall the Maxwell equations (with hats dropped)

 \text{(13)} \qquad 
  \boldsymbol{\nabla} \times \mathbf{E} = i\omega\mu\mathbf{H} ~;~~
  \boldsymbol{\nabla} \times \mathbf{H} = -i\omega\epsilon\mathbf{E} ~.

Assume that \mu and \epsilon are scalars which are independent of x_1, i.e., \mu \equiv \mu(x_2,x_3) and \epsilon \equiv \epsilon(x_2,x_3).

For the TE case, we look for solutions with E_2 = E_3 = 0 and E_1 independent of x_1, i.e., E_1 \equiv E_1(x_2, x_3).

Then,


  \boldsymbol{\nabla} \times \mathbf{E} = \left[0, \frac{\partial E_1}{\partial x_3}, - \frac{\partial E_1}{\partial x_2}\right]~.

This implies that


  \mathbf{H} = \left[ 0, \cfrac{1}{i\omega\mu}~\frac{\partial E_1}{\partial x_3}, 
    - \cfrac{1}{i\omega\mu}~\frac{\partial E_1}{\partial x_2}\right]~.

Therefore,


  \boldsymbol{\nabla} \times \mathbf{H} = \left[ 
    -\frac{\partial }{\partial x_3}\left(\cfrac{1}{i\omega\mu}~\frac{\partial E_1}{\partial x_3}\right)
    -\frac{\partial }{\partial x_2}\left(\cfrac{1}{i\omega\mu}~\frac{\partial E_1}{\partial x_2}\right),
    0, 0\right]~.

or,


  \boldsymbol{\nabla} \times \mathbf{H} = \cfrac{i}{\omega}\left[ 
    \frac{\partial }{\partial x_2}\left(\cfrac{1}{\mu}~\frac{\partial E_1}{\partial x_2}\right) +
    \frac{\partial }{\partial x_3}\left(\cfrac{1}{\mu}~\frac{\partial E_1}{\partial x_3}\right),
    0, 0\right] 
    = \cfrac{i}{\omega} \left[\overline{\boldsymbol{\nabla}} \cdot\left(\cfrac{1}{\mu(x_2,x_3)}
        \overline{\boldsymbol{\nabla}} E_1\right), 0, 0\right] ~.

Plugging into equation (13) we get the TE equation


  {
  \overline{\boldsymbol{\nabla}} \cdot\left(\cfrac{1}{\mu(x_2,x_3)}\overline{\boldsymbol{\nabla}} E_1\right) + 
    \omega^2~\epsilon(x_2,x_3)~E_1 = 0 ~.
  }

This equation has the same form as equation (12).

More generally, if


  \boldsymbol{\mu} = \boldsymbol{\mu}(x_2,x_3) = 
    \begin{bmatrix}
      \mu_{11} & 0 & 0 \\ 0 & \mu_{22} & \mu_{23} \\ 0 & \mu_{23} & \mu_{33}
    \end{bmatrix} = 
    \begin{bmatrix}
      \mu_{11} & \left[\mathsf{0}\right] \\ \left[\mathsf{0}\right] & \left[\mathsf{M}\right]
    \end{bmatrix} 
    \quad \text{where} \quad
    \left[\mathsf{M}\right] = 
    \begin{bmatrix}
      \mu_{22} & \mu_{23} \\ \mu_{23} & \mu_{33}
    \end{bmatrix} \equiv \boldsymbol{M}

and


  \boldsymbol{\epsilon} = \boldsymbol{\epsilon}(x_2,x_3) = 
    \begin{bmatrix}
      \epsilon_{11} & \left[\mathsf{0}\right] \\ \left[\mathsf{0}\right] & \left[\mathsf{N}\right]
    \end{bmatrix} 
    \quad \text{where} \quad
    \left[\mathsf{N}\right] = 
    \begin{bmatrix}
      \epsilon_{22} & \epsilon_{23} \\ \epsilon_{23} & \epsilon_{33}
    \end{bmatrix} \equiv \boldsymbol{N}

we get the TE equation


  {
  \overline{\boldsymbol{\nabla}} \cdot\left[\left(\boldsymbol{R}_{\perp}\cdot\boldsymbol{M}^{-1}\cdot\boldsymbol{R}_{\perp}^T\right)\cdot
       \overline{\boldsymbol{\nabla}} E_1\right] + \omega^2~\epsilon_{11}~E_1~\mathbf{1} = \mathbf{0} 
  }
  \qquad \text{where} \qquad
  \boldsymbol{R}_{\perp} \equiv \left[\mathsf{R}\right]_{\perp} = 
    \begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix} ~.

Similarly, there is a TM equation with H_2 = H_3 = 0 of the form


  {
  \overline{\boldsymbol{\nabla}} \cdot\left[\left(\boldsymbol{R}_{\perp}\cdot\boldsymbol{N}^{-1}\cdot\boldsymbol{R}_{\perp}^T\right)
       \cdot\overline{\boldsymbol{\nabla}} H_1\right] + \omega^2~\mu_{11}~H_1~\mathbf{1} = \mathbf{0}
  }

which for the isotropic case reduces to


  {
  \overline{\boldsymbol{\nabla}} \cdot\left(\cfrac{1}{\epsilon(x_2,x_3)}\overline{\boldsymbol{\nabla}} H_1\right) + 
    \omega^2~\mu(x_2,x_3)~H_1 = 0 ~.
  }

The general solution independent of x_1 is a superposition of the TE and TM solutions. This can be seen by observing that the Maxwell equations decouple under these conditions and a general solution can be written as


  (E_1,E_2,E_3) = (E_1, 0, 0) + (0, E_2, E_3)

where the first term represents the TE solution. We can show that the second term represents the TM solution by observing that


  \boldsymbol{\nabla} \times (0, E_2, E_3) = 
    \left[\frac{\partial E_3}{\partial x_2} - \frac{\partial E_2}{\partial x_3},0,0\right]

implying that H_2 = H_3 = 0 which is the TM solution.

A resonant structure

Consider the periodic geometry shown in Figure 2. The matrix material has a high value of shear modulus (\mu) while the split-ring shaped region has a low shear modulus or is a void. The material inside the ring has the same shear modulus as the matrix material and is connected to the matrix by a thin ligament. The system is subjected to a displacement u_1 in the x_1 direction (parallel to the axis of each cylindrical split ring).

Figure 2. A periodic geometry containing split hollow cylinders of soft material in a matrix of stiff material. The x_1 direction is parallel to the axis of each cylinder.

Clearly, each periodic component of the system behaves like a mass attached to a spring. This is a resonant structure and the effective density \rho_{11}^{*}(\omega) can be negative. A detailed treatment of the problem can be found in Movchan04. Note that the governing equation for this problem is


   \overline{\boldsymbol{\nabla}} \cdot(\mu~\overline{\boldsymbol{\nabla}} u_1) + \omega^2~\rho~u_1 = 0 ~.

Let us compare this problem with the TM case where H_1 is the out of plane magnetic induction. The governing equation now is


   \overline{\boldsymbol{\nabla}} \cdot(\cfrac{1}{\epsilon}~\overline{\boldsymbol{\nabla}} H_1) + \omega^2~\mu~H_1 = 0 ~.

If the value of 1/\epsilon in the region of the void (ring) is small and hence \epsilon is large (which implies that the conductivity \sigma is large), analogy with the equation of elastodynamics implies that the effective permeability \mu_{11}^{*} can be negative for this material.

References

This article is issued from Wikiversity - version of the Wednesday, September 07, 2011. The text is available under the Creative Commons Attribution/Share Alike but additional terms may apply for the media files.