Waves in composites and metamaterials/Bloch waves in elastodynamics and bubbly fluids

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.

Recap

In the previous lecture we discussed Bloch wave solutions for Maxwell's equations in periodic media in the quasistatic limit. Recall the periodic medium shown in Figure 1. The lattice spacing is $\eta$.

Figure 1. Periodic medium.

The permittivity and permeability of the medium are periodic functions of the form

${\boldsymbol {\epsilon }}(\mathbf {y} +\mathbf {a} _{i})={\boldsymbol {\epsilon }}(\mathbf {y} )~;~~{\boldsymbol {\mu }}(\mathbf {y} +\mathbf {a} _{i})={\boldsymbol {\mu }}(\mathbf {y} )$

where $\mathbf {a} _{i}$ are the primitive lattice vectors. We defined periodic functions

${\begin{aligned}{\boldsymbol {\epsilon }}_{\eta }(\mathbf {x} )&:={\boldsymbol {\epsilon }}\left({\cfrac {\mathbf {x} }{\eta }}\right)~;~~{\boldsymbol {\mu }}_{\eta }(\mathbf {x} ):={\boldsymbol {\mu }}\left({\cfrac {\mathbf {x} }{\eta }}\right)~;~~\\\mathbf {D} _{\eta }(\mathbf {x} )&:=\mathbf {D} \left({\cfrac {\mathbf {x} }{\eta }}\right)~;~~\mathbf {B} _{\eta }(\mathbf {x} ):=\mathbf {B} \left({\cfrac {\mathbf {x} }{\eta }}\right)~;~~\mathbf {E} _{\eta }(\mathbf {x} ):=\mathbf {E} \left({\cfrac {\mathbf {x} }{\eta }}\right)~;~~\mathbf {H} _{\eta }(\mathbf {x} ):=\mathbf {H} \left({\cfrac {\mathbf {x} }{\eta }}\right)~.\end{aligned}}$

We then wrote Maxwell's equations (at fixed frequency) as

${\text{(1)}}\qquad {{\boldsymbol {\nabla }}\cdot \mathbf {D} _{\eta }=0~;~~{\boldsymbol {\nabla }}\cdot \mathbf {B} _{\eta }=0~;~~{\boldsymbol {\nabla }}\times \mathbf {E} _{\eta }-i\omega ~\mathbf {B} _{\eta }={\boldsymbol {0}}~;~~{\boldsymbol {\nabla }}\times \mathbf {H} _{\eta }+i\omega ~\mathbf {D} _{\eta }={\boldsymbol {0}}~.}$

We also found that the constitutive relations could be expressed as

${\text{(2)}}\qquad {\mathbf {d} _{\eta }(\mathbf {x} )={\boldsymbol {\epsilon }}_{\eta }(\mathbf {x} )\cdot \mathbf {e} _{\eta }(\mathbf {x} )~;~~\mathbf {b} _{\eta }(\mathbf {x} )={\boldsymbol {\mu }}_{\eta }(\mathbf {x} )\cdot \mathbf {h} _{\eta }(\mathbf {x} )~.}$

Next, we examined solutions to (1) that correspond to Bloch waves with wavevector $\mathbf {k}$ (possibly complex), i.e.,

$\mathbf {E} _{\eta }(\mathbf {x} )=e^{i\mathbf {k} \cdot \mathbf {x} }~\mathbf {e} _{\eta }(\mathbf {x} )~;~~\mathbf {D} _{\eta }(\mathbf {x} )=e^{i\mathbf {k} \cdot \mathbf {x} }~\mathbf {d} _{\eta }(\mathbf {x} )~;~~\mathbf {H} _{\eta }(\mathbf {x} )=e^{i\mathbf {k} \cdot \mathbf {x} }~\mathbf {h} _{\eta }(\mathbf {x} )~;~~\mathbf {B} _{\eta }(\mathbf {x} )=e^{i\mathbf {k} \cdot \mathbf {x} }~\mathbf {b} _{\eta }(\mathbf {x} )~.$

The wavelength of the Bloch waves is

$\lambda ={\cfrac {2\pi }{|{\text{Re}}(\mathbf {k} )|}}~.$

For the quasistatic limit, we look for solutions that satisfy

${\cfrac {2\pi }{|\mathbf {k} |}}\gg |\eta ~\mathbf {a} _{i}|~.$

Such a condition is satisfied if $\eta \rightarrow 0$. Following standard multiple scale analysis, we assumed that the periodic complex fields have the expansions

${\text{(3)}}\qquad {\begin{aligned}\mathbf {e} _{\eta }(\mathbf {x} )&=\mathbf {e} _{0}(\mathbf {y} )+\eta ~\mathbf {e} _{1}(\mathbf {y} )+\eta ^{2}~\mathbf {e} _{2}(\mathbf {y} )+\dots \\\mathbf {d} _{\eta }(\mathbf {x} )&=\mathbf {d} _{0}(\mathbf {y} )+\eta ~\mathbf {d} _{1}(\mathbf {y} )+\eta ^{2}~\mathbf {d} _{2}(\mathbf {y} )+\dots \\\mathbf {h} _{\eta }(\mathbf {x} )&=\mathbf {h} _{0}(\mathbf {y} )+\eta ~\mathbf {h} _{1}(\mathbf {y} )+\eta ^{2}~\mathbf {h} _{2}(\mathbf {y} )+\dots \\\mathbf {b} _{\eta }(\mathbf {x} )&=\mathbf {b} _{0}(\mathbf {y} )+\eta ~\mathbf {b} _{1}(\mathbf {y} )+\eta ^{2}~\mathbf {b} _{2}(\mathbf {y} )+\dots \end{aligned}}$

and that dependence of $\omega$ on $\eta$ and $\mathbf {k}$ has an expansion of the form

${\text{(4)}}\qquad \omega =\omega _{\eta }^{j}(\mathbf {k} )=\omega _{0}^{j}(\mathbf {y} )+\eta ~\omega _{1}^{j}(\mathbf {y} )+\eta ^{2}~\omega _{2}^{j}(\mathbf {y} )+\dots$

Defining

${\boldsymbol {\nabla }}_{y}:=\left({\frac {\partial }{\partial y_{1}}},{\frac {\partial }{\partial y_{2}}},{\frac {\partial }{\partial y_{3}}}\right)$

we then showed that we got the equations in the quasistatic limit

${\text{(5)}}\qquad {{\boldsymbol {\nabla }}_{y}\cdot \mathbf {d} _{0}(\mathbf {y} )=0~;~~{\boldsymbol {\nabla }}_{y}\cdot \mathbf {b} _{0}(\mathbf {y} )=0~;~~{\boldsymbol {\nabla }}_{y}\times \mathbf {e} _{0}(\mathbf {y} )={\boldsymbol {0}}~;~~{\boldsymbol {\nabla }}_{y}\times \mathbf {h} _{0}(\mathbf {y} )={\boldsymbol {0}}}$

and the constitutive equations

${\text{(6)}}\qquad {\mathbf {d} _{0}(\mathbf {y} )={\boldsymbol {\epsilon }}(\mathbf {y} )\cdot \mathbf {e} _{0}(\mathbf {y} )~;~~\mathbf {b} _{0}(\mathbf {y} )={\boldsymbol {\mu }}(\mathbf {y} )\cdot \mathbf {h} _{0}(\mathbf {y} )~.}$

Note that the quasistatic equations (5), (6) have the same form as the Maxwell equations (1) and (2).

We also found that the volume averaged fields ($\langle \bullet \rangle$) satisfy the equations

${\text{(7)}}\qquad {\langle {\boldsymbol {\nabla }}_{y}\cdot \mathbf {d} _{1}(\mathbf {y} )\rangle =0~;~~\langle {\boldsymbol {\nabla }}_{y}\cdot \mathbf {b} _{1}(\mathbf {y} )\rangle =0~;~~\langle {\boldsymbol {\nabla }}_{y}\times \mathbf {e} _{1}(\mathbf {y} )\rangle ={\boldsymbol {0}}~;~~\langle {\boldsymbol {\nabla }}_{y}\times \mathbf {h} _{1}(\mathbf {y} )\rangle ={\boldsymbol {0}}~.}$

along with the necessary conditions

${\text{(8)}}\qquad {\mathbf {k} \times \langle \mathbf {e} _{0}(\mathbf {y} )\rangle -\omega _{0}^{j}~\langle \mathbf {b} _{0}(\mathbf {y} )\rangle ={\boldsymbol {0}}~;~~\mathbf {k} \times \langle \mathbf {h} _{0}(\mathbf {y} )\rangle +\omega _{0}^{j}~\langle \mathbf {d} _{0}(\mathbf {y} )\rangle ={\boldsymbol {0}}~.}$

Effective Permittivity and Permeability

Let ${\boldsymbol {\epsilon }}_{\text{eff}}$ be an effective permittivity tensor associated with the field ${\boldsymbol {\epsilon }}(\mathbf {y} )$ and let ${\boldsymbol {\mu }}_{\text{eff}}$ be an effective permeability tensor associated with the field ${\boldsymbol {\mu }}(\mathbf {y} )$. These effective tensors are defined through

${\text{(9)}}\qquad \langle \mathbf {d} _{0}\rangle ={\boldsymbol {\epsilon }}_{\text{eff}}\cdot \langle \mathbf {e} _{0}\rangle ~;~~\langle \mathbf {b} _{0}\rangle ={\boldsymbol {\mu }}_{\text{eff}}\cdot \langle \mathbf {h} _{0}\rangle ~.$

Plugging (9) into (8) gives

${\text{(10)}}\qquad \mathbf {k} \times \langle \mathbf {e} _{0}\rangle -\omega _{0}^{j}~{\boldsymbol {\mu }}_{\text{eff}}\cdot \langle \mathbf {h} _{0}\rangle ={\boldsymbol {0}}~;~~\mathbf {k} \times \langle \mathbf {h} _{0}\rangle +\omega _{0}^{j}~{\boldsymbol {\epsilon }}_{\text{eff}}\cdot \langle \mathbf {e} _{0}\rangle ={\boldsymbol {0}}~.$

Eliminate $\langle \mathbf {h} _{0}\rangle$ from (10) to get

${\cfrac {1}{\omega _{0}^{j}}}~\mathbf {k} \times \left\{{\boldsymbol {\mu }}_{\text{eff}}^{-1}\cdot [\mathbf {k} \times \langle \mathbf {e} _{0}\rangle ]\right\}+\omega _{0}^{j}~{\boldsymbol {\epsilon }}_{\text{eff}}\cdot \langle \mathbf {e} _{0}\rangle ={\boldsymbol {0}}$

or,

${\text{(11)}}\qquad {\boldsymbol {\epsilon }}_{\text{eff}}^{-1}\cdot \left(\mathbf {k} \times \left\{{\boldsymbol {\mu }}_{\text{eff}}^{-1}\cdot [\mathbf {k} \times \langle \mathbf {e} _{0}\rangle ]\right\}\right)+(\omega _{0}^{j})^{2}~\langle \mathbf {e} _{0}\rangle ={\boldsymbol {0}}~.$

Define the matrix ${\boldsymbol {A}}_{\text{eff}}(\mathbf {k} )$ such that for all vectors $\mathbf {v}$

${\text{(12)}}\qquad {\boldsymbol {A}}_{\text{eff}}(\mathbf {k} )\cdot \mathbf {v} ={\boldsymbol {\epsilon }}_{\text{eff}}^{-1}\cdot \left(\mathbf {k} \times \left\{{\boldsymbol {\mu }}_{\text{eff}}^{-1}\cdot [\mathbf {k} \times \mathbf {v} ]\right\}\right)~.$

Then (11) can be written as

${\text{(13)}}\qquad {{\boldsymbol {A}}_{\text{eff}}(\mathbf {k} )\cdot \langle \mathbf {e} _{0}\rangle +(\omega _{0}^{j})^{2}~\langle \mathbf {e} _{0}\rangle ={\boldsymbol {0}}~.}$

The eigenvalue problem (13) may be used to determine the dispersion relations $\omega _{0}^{j}$ and also possible values of $\langle \mathbf {e} _{0}\rangle$ associated with each Block wave mode. Once we know $\langle \mathbf {e} _{0}\rangle$ we can then find $\langle \mathbf {h} _{0}\rangle$.

Note that the eigenvalue of ${\boldsymbol {A}}_{\text{eff}}$ is zero when $\mathbf {e} _{0}=\mathbf {k}$. So, in practice, it is necessary to examine only the other two eigenvalues.

Summary:

If the average field $\langle \mathbf {e} _{0}\rangle$ is to be a Bloch wave solution of wavevector $\mathbf {k}$ and frequency $\omega _{0}^{j}$, then $\langle \mathbf {e} _{0}\rangle$ must be an eigenvector of the matrix ${\boldsymbol {A}}_{\text{eff}}(\mathbf {k} )$ and $(\omega _{0}^{j})^{2}$ must be the associated eigenvalue. More general solutions can be obtained by superposing different Block wave solutions.

Isotropic constituents

When the permittivity and permeability of the medium are isotropic, equation (12) becomes

${\boldsymbol {A}}_{\text{eff}}(\mathbf {k} )\cdot \mathbf {v} ={\cfrac {1}{\epsilon _{\text{eff}}~\mu _{\text{eff}}}}\cdot \left(\mathbf {k} \times \mathbf {k} \times \mathbf {v} \right)$

or,

${\text{(14)}}\qquad {\boldsymbol {A}}_{\text{eff}}(\mathbf {k} )\cdot \mathbf {v} ={\cfrac {1}{\epsilon _{\text{eff}}~\mu _{\text{eff}}}}\cdot \left[(\mathbf {k} \cdot \mathbf {k} )~\mathbf {v} -(\mathbf {k} \cdot \mathbf {v} )~\mathbf {k} \right]~.$

Since $\mathbf {k} \cdot \langle \mathbf {e} _{0}\rangle =0$, equation (13) becomes

${\text{(15)}}\qquad {\cfrac {1}{\epsilon _{\text{eff}}~\mu _{\text{eff}}}}~(\mathbf {k} \cdot \mathbf {k} )+(\omega _{0}^{j})^{2}=0~.$

If we assume that

$\mathbf {k} =\left({\cfrac {2\pi }{\lambda }}+{\cfrac {i}{\delta }}\right)~\mathbf {n}$

where $\lambda$ is the wavelength, $\delta$ is the attenuation length, and $\mathbf {n}$ is a real unit vector, equation (15) gives

${\cfrac {1}{\epsilon _{\text{eff}}~\mu _{\text{eff}}}}~\left({\cfrac {2\pi }{\lambda }}+{\cfrac {i}{\delta }}\right)^{2}+(\omega _{0}^{j})^{2}=0~.$

Therefore for a given frequency, the wavelength and attenuation length are determined by the complex permittivity $\epsilon _{\text{eff}}$ and the complex permeability $\mu _{\text{eff}}$.

Elastic wave propagation in the quasistatic limit

Let us now perform the exercise for linear elastic materials. Recall that the momentum equation at fixed frequency is given by

${\boldsymbol {\nabla }}\cdot {\boldsymbol {\sigma }}(\mathbf {x} )=-\omega ^{2}~\rho (\mathbf {x} )~\mathbf {u} (\mathbf {x} )$

where ${\boldsymbol {\sigma }}$ is the stress, $\rho$ is the mass density, and $\mathbf {u}$ is the displacement. For infinitesimal deformations, the strain is related to the displacement by

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

Also, the stress is related to the strain for linear elastic materials by

${\boldsymbol {\sigma }}(\mathbf {x} )={\boldsymbol {\mathsf {C}}}(\mathbf {x} ):{\boldsymbol {\epsilon }}(\mathbf {x} )~.$

For a periodic medium, define

${\boldsymbol {\sigma }}_{\eta }(\mathbf {x} )={\boldsymbol {\sigma }}(\mathbf {x} /\eta )~;~~{\boldsymbol {\epsilon }}_{\eta }(\mathbf {x} )={\boldsymbol {\epsilon }}(\mathbf {x} /\eta )~;~~\mathbf {u} _{\eta }(\mathbf {x} )=\mathbf {u} (\mathbf {x} /\eta )~;~~\rho _{\eta }(\mathbf {x} )=\rho (\mathbf {x} /\eta )~;~~{\boldsymbol {\mathsf {C}}}_{\eta }(\mathbf {x} )={\boldsymbol {\mathsf {C}}}(\mathbf {x} /\eta )~.$

Then the governing equations can be written as

${\text{(16)}}\qquad {\boldsymbol {\nabla }}\cdot {\boldsymbol {\sigma }}_{\eta }(\mathbf {x} )=-\omega ^{2}~\rho _{\eta }(\mathbf {x} )~\mathbf {u} _{\eta }(\mathbf {x} )~;~~{\boldsymbol {\epsilon }}_{\eta }(\mathbf {x} )={\frac {1}{2}}\left\{{\boldsymbol {\nabla }}\mathbf {u} _{\eta }(\mathbf {x} )+[{\boldsymbol {\nabla }}\mathbf {u} _{\eta }(\mathbf {x} )]^{T}\right\}~;~~{\boldsymbol {\sigma }}_{\eta }(\mathbf {x} )={\boldsymbol {\mathsf {C}}}_{\eta }(\mathbf {x} ):{\boldsymbol {\epsilon }}_{\eta }(\mathbf {x} )~.$

Let us examine Bloch wave solutions to these equations of the form

${\text{(17)}}\qquad {\boldsymbol {\sigma }}_{\eta }(\mathbf {x} )=e^{i~\mathbf {k} \cdot \mathbf {x} }~{\widehat {\boldsymbol {\sigma }}}_{\eta }(\mathbf {x} )~;~~{\boldsymbol {\epsilon }}_{\eta }(\mathbf {x} )=e^{i~\mathbf {k} \cdot \mathbf {x} }~{\widehat {\boldsymbol {\epsilon }}}_{\eta }(\mathbf {x} )~;~~\mathbf {u} _{\eta }(\mathbf {x} )=e^{i~\mathbf {k} \cdot \mathbf {x} }~{\widehat {\mathbf {u} }}_{\eta }(\mathbf {x} )~.$

Plug (17) into (16) to get

${\begin{aligned}{\boldsymbol {\nabla }}\cdot \left(e^{i~\mathbf {k} \cdot \mathbf {x} }~{\widehat {\boldsymbol {\sigma }}}_{\eta }(\mathbf {x} )\right)&=-\omega ^{2}~\rho _{\eta }(\mathbf {x} )~e^{i~\mathbf {k} \cdot \mathbf {x} }~{\widehat {\mathbf {u} }}_{\eta }(\mathbf {x} )\\e^{i~\mathbf {k} \cdot \mathbf {x} }~{\widehat {\boldsymbol {\epsilon }}}_{\eta }(\mathbf {x} )&={\frac {1}{2}}\left\{{\boldsymbol {\nabla }}\left(e^{i~\mathbf {k} \cdot \mathbf {x} }~{\widehat {\mathbf {u} }}_{\eta }(\mathbf {x} )\right)+\left[{\boldsymbol {\nabla }}\left(e^{i~\mathbf {k} \cdot \mathbf {x} }~{\widehat {\mathbf {u} }}_{\eta }(\mathbf {x} )\right)\right]\right\}\\e^{i~\mathbf {k} \cdot \mathbf {x} }~{\widehat {\boldsymbol {\sigma }}}_{\eta }(\mathbf {x} )&={\boldsymbol {\mathsf {C}}}_{\eta }(\mathbf {x} ):\left(e^{i~\mathbf {k} \cdot \mathbf {x} }~{\widehat {\boldsymbol {\epsilon }}}_{\eta }(\mathbf {x} )\right)\end{aligned}}$

or,

${\begin{aligned}&e^{i~\mathbf {k} \cdot \mathbf {x} }~\left({\boldsymbol {\nabla }}\cdot {\widehat {\boldsymbol {\sigma }}}_{\eta }+i~{\widehat {\boldsymbol {\sigma }}}_{\eta }~\cdot \mathbf {k} \right)=-e^{i~\mathbf {k} \cdot \mathbf {x} }~\omega ^{2}~\rho _{\eta }~{\widehat {\mathbf {u} }}_{\eta }\\&e^{i~\mathbf {k} \cdot \mathbf {x} }~{\widehat {\boldsymbol {\epsilon }}}_{\eta }=e^{i~\mathbf {k} \cdot \mathbf {x} }~\left[{\cfrac {i}{2}}\left({\widehat {\mathbf {u} }}_{\eta }\otimes \mathbf {k} +\mathbf {k} \otimes {\widehat {\mathbf {u} }}_{\eta }\right)+{\frac {1}{2}}~\left({\boldsymbol {\nabla }}{\widehat {\mathbf {u} }}_{\eta }+[{\boldsymbol {\nabla }}{\widehat {\mathbf {u} }}_{\eta }]^{T}\right)\right]\\&e^{i~\mathbf {k} \cdot \mathbf {x} }~{\widehat {\boldsymbol {\sigma }}}_{\eta }=e^{i~\mathbf {k} \cdot \mathbf {x} }~{\boldsymbol {\mathsf {C}}}_{\eta }:{\widehat {\boldsymbol {\epsilon }}}_{\eta }\end{aligned}}$

or,

${\text{(18)}}\qquad {\begin{aligned}&{\boldsymbol {\nabla }}\cdot {\widehat {\boldsymbol {\sigma }}}_{\eta }+i~{\widehat {\boldsymbol {\sigma }}}_{\eta }~\cdot \mathbf {k} +\omega ^{2}~\rho _{\eta }~{\widehat {\mathbf {u} }}_{\eta }=0\\&{\widehat {\boldsymbol {\epsilon }}}_{\eta }={\cfrac {i}{2}}\left({\widehat {\mathbf {u} }}_{\eta }\otimes \mathbf {k} +\mathbf {k} \otimes {\widehat {\mathbf {u} }}_{\eta }\right)+{\frac {1}{2}}~\left({\boldsymbol {\nabla }}{\widehat {\mathbf {u} }}_{\eta }+[{\boldsymbol {\nabla }}{\widehat {\mathbf {u} }}_{\eta }]^{T}\right)\\&{\widehat {\boldsymbol {\sigma }}}_{\eta }={\boldsymbol {\mathsf {C}}}_{\eta }:{\widehat {\boldsymbol {\epsilon }}}_{\eta }~.\end{aligned}}$

Equations (18) have solutions if the frequency $\omega$ takes one of a discrete set of values, $\omega _{\eta }^{j}(\mathbf {k} )$, $j=1,2,\dots$. Let us assume that the fields have perturbation expansions in powers of $\eta$, i.e,

${\text{(19)}}\qquad {\begin{aligned}&{\widehat {\boldsymbol {\sigma }}}_{\eta }(\mathbf {x} )={\boldsymbol {\sigma }}_{0}(\mathbf {y} )+\eta ~{\boldsymbol {\sigma }}_{1}(\mathbf {y} )+\eta ^{2}~{\boldsymbol {\sigma }}_{2}(\mathbf {y} )+\dots \\&{\widehat {\boldsymbol {\epsilon }}}_{\eta }(\mathbf {x} )={\boldsymbol {\epsilon }}_{0}(\mathbf {y} )+\eta ~{\boldsymbol {\epsilon }}_{1}(\mathbf {y} )+\eta ^{2}~{\boldsymbol {\epsilon }}_{2}(\mathbf {y} )+\dots \\&{\widehat {\mathbf {u} }}_{\eta }(\mathbf {x} )=\mathbf {u} _{0}(\mathbf {y} )+\eta ~\mathbf {u} _{1}(\mathbf {y} )+\eta ^{2}~\mathbf {u} _{2}(\mathbf {y} )+\dots \end{aligned}}$

where the functions ${\boldsymbol {\sigma }}_{i}$, ${\boldsymbol {\epsilon }}_{i}$, $\mathbf {u} _{i}$ are periodic. Let us also assume that

$\omega =\omega _{\eta }^{j}(\mathbf {k} )=\omega _{0}^{j}(\mathbf {y} )+\eta ~\omega _{1}^{j}(\mathbf {y} )+\eta ^{2}~\omega _{2}^{j}(\mathbf {y} )+\dots$

Plugging these expansions into (18) and using the definition

${\boldsymbol {\nabla }}_{y}:=\left({\frac {\partial }{\partial y_{1}}},{\frac {\partial }{\partial y_{2}}},{\frac {\partial }{\partial y_{3}}}\right)\qquad \implies {\boldsymbol {\nabla }}={\cfrac {1}{\eta }}~{\boldsymbol {\nabla }}_{y}\quad {\text{since}}\quad \mathbf {y} ={\cfrac {\mathbf {x} }{\eta }}$

we get

${\text{(20)}}\qquad {\begin{aligned}&{\cfrac {1}{\eta }}{\boldsymbol {\nabla }}_{y}\cdot ({\boldsymbol {\sigma }}_{0}+\eta ~{\boldsymbol {\sigma }}_{1}+\dots )+i~({\boldsymbol {\sigma }}_{0}+\eta ~{\boldsymbol {\sigma }}_{1}+\dots )~\cdot \mathbf {k} +\rho _{\eta }~(\omega _{0}^{j}+\eta ~\omega _{1}^{j}+\dots )^{2}~(\mathbf {u} _{0}+\eta ~\mathbf {u} _{1}+\dots )=0\\&{\boldsymbol {\epsilon }}_{0}+\eta ~{\boldsymbol {\epsilon }}_{1}+\dots ={\cfrac {i}{2}}\left[(\mathbf {u} _{0}+\eta ~\mathbf {u} _{1}+\dots )\otimes \mathbf {k} +\mathbf {k} \otimes (\mathbf {u} _{0}+\eta ~\mathbf {u} _{1}+\dots )\right]+\\&\qquad \qquad \qquad \qquad {\cfrac {1}{2\eta }}~\left({\boldsymbol {\nabla }}_{y}(\mathbf {u} _{0}+\eta ~\mathbf {u} _{1}+\dots )+[{\boldsymbol {\nabla }}_{y}(\mathbf {u} _{0}+\eta ~\mathbf {u} _{1}+\dots )]^{T}\right)\\&{\boldsymbol {\sigma }}_{0}+\eta ~{\boldsymbol {\sigma }}_{1}+\dots ={\boldsymbol {\mathsf {C}}}_{\eta }:({\boldsymbol {\epsilon }}_{0}+\eta ~{\boldsymbol {\epsilon }}_{1}+\dots )~.\end{aligned}}$

Collecting terms containing $1/\eta$ from (20) we get

${\begin{aligned}&{{\boldsymbol {\nabla }}_{y}\cdot {\boldsymbol {\sigma }}_{0}=0}{\text{(21)}}\qquad \\&{{\boldsymbol {\nabla }}_{y}\mathbf {u} _{0}+[{\boldsymbol {\nabla }}_{y}\mathbf {u} _{0}]^{T}=0}~.{\text{(22)}}\qquad \end{aligned}}$

Collecting terms of order $\eta ^{0}$ from (20) we get

${\begin{aligned}&{{\boldsymbol {\nabla }}_{y}\cdot {\boldsymbol {\sigma }}_{1}+i~{\boldsymbol {\sigma }}_{0}\cdot \mathbf {k} +\rho _{\eta }~(\omega _{0}^{j})^{2}~\mathbf {u} _{0}=0}{\text{(23)}}\qquad \\&{{\boldsymbol {\epsilon }}_{0}={\cfrac {i}{2}}~(\mathbf {u} _{0}\otimes \mathbf {k} +\mathbf {k} \otimes \mathbf {u} _{0})+{\frac {1}{2}}~({\boldsymbol {\nabla }}\mathbf {u} _{1}+[{\boldsymbol {\nabla }}\mathbf {u} _{1}]^{T})}{\text{(24)}}\qquad \\&{{\boldsymbol {\sigma }}_{0}={\boldsymbol {\mathsf {C}}}_{\eta }:{\boldsymbol {\epsilon }}_{0}}{\text{(25)}}\qquad ~.\end{aligned}}$

From (22) we see that $\mathbf {u} _{0}(\mathbf {y} )$ must be linear in $\mathbf {y}$. We also know that (from our definition), $\mathbf {u} _{0}(\mathbf {y} )$ is periodic in $\mathbf {y}$. A function that is both linear and periodic must be constant. Therefore,

$\mathbf {u} _{0}(\mathbf {y} )={\widehat {\mathbf {u} }}_{0}=~{\text{constant}}~.$

Then we can write (24) in the form

${\text{(26)}}\qquad {\boldsymbol {\epsilon }}_{0}={\cfrac {i}{2}}~({\widehat {\mathbf {u} }}_{0}\otimes \mathbf {k} +\mathbf {k} \otimes {\widehat {\mathbf {u} }}_{0})+{\frac {1}{2}}~({\boldsymbol {\nabla }}\mathbf {u} _{1}+[{\boldsymbol {\nabla }}\mathbf {u} _{1}]^{T})\qquad \implies {{\boldsymbol {\epsilon }}_{0}={\frac {1}{2}}~\left({\boldsymbol {\nabla }}_{y}{\widehat {\mathbf {u} }}(\mathbf {y} )+[{\boldsymbol {\nabla }}_{y}{\widehat {\mathbf {u} }}(\mathbf {y} )]^{T}\right)}$

where

${\widehat {\mathbf {u} }}(\mathbf {y} )=i~(\mathbf {k} \cdot \mathbf {y} )~{\widehat {\mathbf {u} }}_{0}+\mathbf {u} _{1}(\mathbf {y} )~.$

Note that equations (21), (26), and (25) have a form similar to that of the elasticity equations in the absence of inertial and body forces.

The complex effective elasticity tensor may be defined via the relation

${\text{(27)}}\qquad {\langle {\boldsymbol {\sigma }}_{0}\rangle ={\boldsymbol {\mathsf {C}}}_{\text{eff}}:\langle {\boldsymbol {\epsilon }}_{0}\rangle ~.}$

Taking the average of (26) over the periodic cell, we get

$\langle {\boldsymbol {\epsilon }}_{0}\rangle ={\cfrac {i}{2}}~({\widehat {\mathbf {u} }}_{0}\otimes \mathbf {k} +\mathbf {k} \otimes {\widehat {\mathbf {u} }}_{0})+{\frac {1}{2}}~\left(\langle {\boldsymbol {\nabla }}\mathbf {u} _{1}\rangle +\langle [{\boldsymbol {\nabla }}\mathbf {u} _{1}]^{T}\rangle \right)$

or

${\text{(28)}}\qquad {\langle {\boldsymbol {\epsilon }}_{0}\rangle ={\cfrac {i}{2}}~({\widehat {\mathbf {u} }}_{0}\otimes \mathbf {k} +\mathbf {k} \otimes {\widehat {\mathbf {u} }}_{0})~.}$

Also, taking the average of (23) we get

$\langle {\boldsymbol {\nabla }}_{y}\cdot {\boldsymbol {\sigma }}_{1}\rangle +i~\langle {\boldsymbol {\sigma }}_{0}\rangle \cdot \mathbf {k} +\langle \rho _{\eta }\rangle ~(\omega _{0}^{j})^{2}~{\widehat {\mathbf {u} }}_{0}=0$

or,

${\text{(29)}}\qquad {i~\langle {\boldsymbol {\sigma }}_{0}\rangle \cdot \mathbf {k} +\langle \rho _{\eta }\rangle ~(\omega _{0}^{j})^{2}~{\widehat {\mathbf {u} }}_{0}=0~.}$

Plugging (28) into (27) gives

${\text{(30)}}\qquad \langle {\boldsymbol {\sigma }}_{0}\rangle ={\cfrac {i}{2}}~{\boldsymbol {\mathsf {C}}}_{\text{eff}}:({\widehat {\mathbf {u} }}_{0}\otimes \mathbf {k} +\mathbf {k} \otimes {\widehat {\mathbf {u} }}_{0})~.$

Plugging (30) into (29) gives

$-{\frac {1}{2}}~\left[{\boldsymbol {\mathsf {C}}}_{\text{eff}}:({\widehat {\mathbf {u} }}_{0}\otimes \mathbf {k} +\mathbf {k} \otimes {\widehat {\mathbf {u} }}_{0})\right]\cdot \mathbf {k} +\langle \rho _{\eta }\rangle ~(\omega _{0}^{j})^{2}~{\widehat {\mathbf {u} }}_{0}=0~.$

From the minor symmetry of the tensor ${\boldsymbol {\mathsf {C}}}_{\text{eff}}$, i.e.,

$C_{ijkl}^{\text{eff}}=C_{ijlk}^{\text{eff}}$

we have

${\boldsymbol {\mathsf {C}}}_{\text{eff}}:({\widehat {\mathbf {u} }}_{0}\otimes \mathbf {k} )={\boldsymbol {\mathsf {C}}}_{\text{eff}}:(\mathbf {k} \otimes {\widehat {\mathbf {u} }}_{0})~.$

Hence

$\left[{\boldsymbol {\mathsf {C}}}_{\text{eff}}:(\mathbf {k} \otimes {\widehat {\mathbf {u} }}_{0})\right]\cdot \mathbf {k} =\langle \rho _{\eta }\rangle ~(\omega _{0}^{j})^{2}~{\widehat {\mathbf {u} }}_{0}$

or,

${{\boldsymbol {A}}_{\text{eff}}(\mathbf {k} )~{\widehat {\mathbf {u} }}_{0}=(\omega _{0}^{j})^{2}~{\widehat {\mathbf {u} }}_{0}}$

where, for any vector $\mathbf {v}$,

${{\boldsymbol {A}}_{\text{eff}}(\mathbf {k} )~\mathbf {v} ={\cfrac {1}{\langle \rho _{\eta }\rangle }}~\left[{\boldsymbol {\mathsf {C}}}_{\text{eff}}:(\mathbf {k} \otimes \mathbf {v} )\right]\cdot \mathbf {k} ~.}$

The quantity ${\boldsymbol {A}}_{\text{eff}}$ is the effective acoustic tensor. The dispersion relations may then be calculated in a manner similar to that for electromagnetism.

Effective properties of bubbly fluid

So far we have found that, in the quasistatic limit, the Bloch solutions satisfy equations which are directly analogous to the electromagnetic or elasticity equations but with complex fields and complex effective tensors. This implies that if we have a formula for the effective tensor which is valid for real tensors, then (by analytic continuation) we can use the same formula when the tensors take complex values.

Let us examine one such situation which arises in bubbly fluids. We can assume bubbly fluids to be an assembly of coated spheres. For such a geometry, we have Hashin's relations Hashin62,Hashin62a when the material properties are real. Let us now discuss how the Hashin relations are obtained.

Hashin's relations for assemblages of coated spheres

Consider the coated sphere shown in Figure~2. The sphere has a permittivity $\epsilon _{1}$ and the coating has a permittivity $\epsilon _{2}$. Let the coated sphere be embedded in a matrix with permittivity $\epsilon _{0}$.

Figure 2. Coated sphere in a matrix.

For certain choices of the permittivities, current will tends to flow around the coated sphere while for other choices current will be attracted toward the coated sphere. Therefore we expect that there might be a situation in which the permittivities are such that the coated sphere will have no effect on a current flowing through the matrix. In such a situation, we could continue to add spheres and completely fill space (except for a set of measure zero) without affecting the field as can be seen in Figure 3. Clearly, in that case, the effective permittivity of the assemblage of coated spheres is equal to that of the matrix.

Figure 3. Assemblage of coated spheres.

In determining the effective properties of an assemblage of coated spheres we make the following assumptions:

1. the coated spheres do not overlap the boundary of the unit cell.
2. fluxes and potentials at the boundary remain unaltered due to addition of coated spheres.

The goal is to find a matrix whose properties are such that when a coated sphere is added to it, the fields in the matrix remain unaltered.

Effective Permittivity

Consider the single coated sphere shown in Figure 2. Let us assume that the sphere is centered at the origin. We look for a solution to the time-independent Maxwell's equations

${\boldsymbol {\nabla }}\times \mathbf {E} =0~;~~{\boldsymbol {\nabla }}\cdot \mathbf {D} =0~;~~\mathbf {D} =\epsilon ~\mathbf {E}$

with potentials

${\begin{aligned}\varphi _{1}(x)&=a_{1}~r~\cos \theta &\qquad {\text{in the core}}\\\varphi _{2}(x)&=\left(a_{2}~r+{\cfrac {b_{2}}{r^{2}}}\right)~\cos \theta &\qquad {\text{in the coating}}\\\varphi _{\text{eff}}(x)&=a_{\text{eff}}~r~\cos \theta &\qquad {\text{in the effective medium.}}\end{aligned}}$

Then the electric field is given by

${\begin{aligned}\mathbf {E} _{1}={\boldsymbol {\nabla }}\varphi _{1}(x)&=a_{1}~[\cos \theta ~\mathbf {e} _{r}-\sin \theta ~\mathbf {e} _{\theta }]\\\mathbf {E} _{2}={\boldsymbol {\nabla }}\varphi _{2}(x)&=\left(a_{2}~r-{\cfrac {2b_{2}}{r^{3}}}\right)~\cos \theta ~\mathbf {e} _{r}-\left(a_{2}~r+{\cfrac {b_{2}}{r^{3}}}\right)~\sin \theta ~\mathbf {e} _{\theta }\\\mathbf {E} _{\text{eff}}={\boldsymbol {\nabla }}\varphi _{\text{eff(x)}}&=a_{\text{eff}}~[\cos \theta ~\mathbf {e} _{r}-\sin \theta ~\mathbf {e} _{\theta }]~.\end{aligned}}$

The potentials satisfy Laplace's equations

$\nabla ^{2}\varphi _{1}=\nabla ^{2}\varphi _{2}=\nabla ^{2}\varphi _{\text{eff}}=0$

and we only need to match the boundary conditions at the interfaces to get a solution, i.e.,

${\text{(31)}}\qquad a_{1}~r_{c}=a_{2}~r_{c}+{\cfrac {b_{2}}{r_{c}^{2}}}~;~~a_{\text{eff}}~r_{e}=a_{2}~r_{e}+{\cfrac {b_{2}}{r_{e}^{2}}}~.$

Continuity of the tangential component of the electric field at the interface implies that

${\text{(32)}}\qquad \epsilon _{1}~a_{1}=\epsilon _{2}~\left[a_{2}-{\cfrac {2b_{2}}{r_{c}^{3}}}\right]~;~~\epsilon _{\text{eff}}~a_{\text{eff}}=\epsilon _{2}~\left[a_{2}-{\cfrac {2b_{2}}{r_{e}^{3}}}\right]~.$

Combining (31) and (32) gives

${\cfrac {b_{2}}{r_{c}^{3}}}\left[1+{\cfrac {3~\epsilon _{2}}{\epsilon _{1}-\epsilon _{2}}}\right]={\cfrac {b_{2}}{r_{e}^{3}}}\left[1+{\cfrac {3~\epsilon _{2}}{\epsilon _{\text{eff}}-\epsilon _{2}}}\right]~.$

Defining the volume fraction $f$ as

$f_{1}=1-f_{2}={\cfrac {r_{c}^{3}}{r_{e}^{3}}}$

leads to the following expression for $\epsilon _{\text{eff}}$:

${\epsilon _{\text{eff}}=\epsilon _{2}+{\cfrac {3~f_{1}~\epsilon _{2}(\epsilon _{1}-\epsilon _{2})}{3~\epsilon _{2}+f_{2}~(\epsilon _{1}-\epsilon _{2})}}~.}$

Effective Bulk Modulus

In this case, consider the coated sphere shown in Figure~2 with the difference that each material now has two properties- the bulk and shear moduli. Let the bulk and shear moduli of the sphere be $\kappa _{1}$ and $\mu _{1}$, let those of the coating be $\kappa _{2}$ and $\mu _{2}$, and let the moduli of the effective medium (the matrix in the figure) be $\kappa _{\text{eff}}$ and $\mu _{\text{eff}}$, respectively. As before, the radius of the sphere is $r_{c}$ and the outer boundary of the coating has a radius $r_{e}$.

The governing equations are

${\boldsymbol {\nabla }}\cdot {\boldsymbol {\sigma }}=0~;~~{\boldsymbol {\epsilon }}={\frac {1}{2}}\left[{\boldsymbol {\nabla }}\mathbf {u} +({\boldsymbol {\nabla }}\mathbf {u} )^{T}\right]~;~~{\boldsymbol {\sigma }}=\lambda ~{\text{tr}}{\boldsymbol {\epsilon }}~{\boldsymbol {\mathit {1}}}+2~\mu ~{\boldsymbol {\epsilon }}~;~~\kappa =\lambda +{\cfrac {2}{3}}~\mu ~.$

Let the matrix be subject to a hydrostatic state of stress, i.e.,

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

Let us look for a solution that does not perturb this field when the coated sphere is added to the matrix.

Therefore, we look for a solution with a radial displacement field

${\begin{aligned}\mathbf {u} _{1}(\mathbf {x} )&=a_{1}~r~\mathbf {e} _{r}&\qquad {\text{in the core}}\\\mathbf {u} _{2}(\mathbf {x} )&=\left(a_{2}~r+{\cfrac {b_{2}}{r^{2}}}\right)~\mathbf {e} _{r}&\qquad {\text{in the coating}}\\\mathbf {u} _{\text{eff}}(\mathbf {x} )&=a_{\text{eff}}~r~\mathbf {e} _{r}&\qquad {\text{in the effective medium}}\end{aligned}}$

From the continuity of displacements at the interfaces, we have

$a_{1}~r_{c}=a_{2}~r_{c}+{\cfrac {b_{2}}{r_{c}^{2}}}~;~~a_{\text{eff}}~r_{e}=a_{2}~r_{e}+{\cfrac {b_{2}}{r_{e}^{2}}}~.$

And the continuity of radial tractions ${\boldsymbol {\sigma }}\cdot \mathbf {e} _{r}$ across the interface implies that

$\kappa _{1}~a_{1}=\kappa _{2}~a_{2}-{\cfrac {4~b_{2}~\mu _{2}}{3~r_{c}^{3}}}~;~~\kappa _{\text{eff}}~a_{\text{eff}}=\kappa _{2}~a_{2}-{\cfrac {4~b_{2}~\mu _{2}}{3~r_{e}^{3}}}~;~~$

From the continuity relations and using the definition

$f_{1}=1-f_{2}={\cfrac {r_{c}^{3}}{r_{e}^{3}}}$

we can show that the effective bulk modulus is given by

${\text{(33)}}\qquad {\kappa _{\text{eff}}=\kappa _{2}+{\cfrac {f_{1}}{{\cfrac {1}{\kappa _{1}-\kappa _{2}}}+{\cfrac {f_{2}}{\kappa _{2}+{\cfrac {4}{3}}~\mu _{2}}}}}~.}$

Bubbly Fluid

For the assemblage of coated spheres, when the bulk and shear moduli of the two phases are real, the effective bulk modulus is given by equation (33). We can then use the Bloch wave solutions to show that the same result holds when the moduli are complex.

Returning to the bubbly fluid problem, suppose that phase 1 is gas and phase 2 is water. The mixture is then a bubbly fluid. If we assume that water is incompressible, then $\kappa _{2}\rightarrow \infty$. Hence from (33) we have

${\text{(34)}}\qquad \kappa _{\text{eff}}\approx {\cfrac {\kappa _{1}+{\cfrac {4}{3}}~f_{2}~\mu _{2}}{f_{1}}}~.$

Also, assume that $\kappa _{1}$ (air) is independent of frequency.

Now consider a plane shear wave propagating into a bubbly fluid. Let the frequency of the wave be $\omega$ and let it be real. Let the wave be spatially attenuated with a complex wavevector $\mathbf {k}$. Then the associated strain field is given by

${\boldsymbol {\epsilon }}(\mathbf {x} ,\omega )=e^{i~\mathbf {k} \cdot \mathbf {x} }~{\widehat {\boldsymbol {\epsilon }}}(\omega )~.$

The stress field is given by (assuming a Newtonian fluid)

${\boldsymbol {\sigma }}(\mathbf {x} ,\omega )=-i~\omega ~\eta _{\mu }~e^{i~\mathbf {k} \cdot \mathbf {x} }~{\widehat {\boldsymbol {\epsilon }}}(\omega )=-i~\omega ~\eta _{\mu }~{\boldsymbol {\epsilon }}(\mathbf {x} ,\omega )$

where $\eta _{\mu }$ is the shear viscosity of water and is assume to be independent of frequency. Therefore, we can think of water as having a complex shear modulus, i.e.,

$\mu _{2}=-i~\omega ~\eta _{\mu }~.$

Plugging this in (34) gives

$\kappa _{\text{eff}}\approx {\cfrac {\kappa _{1}}{f_{1}}}-i~\eta _{\kappa }^{\text{eff}}~\omega$

where the effective bulk viscosity is defined as

$\eta _{\kappa }^{\text{eff}}\approx {\cfrac {4~f_{2}~\eta _{\mu }}{3~f_{1}}}~.$

When $f_{1}\rightarrow 0$, we get

$\eta _{\kappa }^{\text{eff}}\rightarrow {\cfrac {4~\eta _{\mu }}{3~f_{1}}}~.$

Therefore, the shear viscosity of water has been converted into the bulk viscosity of the bubbly fluid. This is the reason that sound is damped strongly in bubbly fluids.

References

• Z. Hashin. The elastic moduli of heterogeneous materials. J. Appl. Mech., 29:143--150, 1962.
• Z. Hashin and S. Shtrikman. A variational approach to the theory of the effective magnetic permeability of multiphase materials. J. Appl. Phys., 33(10):3125--3131, 1962.
• G. W. Milton. Theory of Composites. Cambridge University Press, New York, 2002.
• M. I. Hussein. Reduced Bloch mode expansion for periodic media band structure calculations. Proc. R. Soc. A, 465:2825-2848, 2009.