nwhm/instability

The homogenized equations for a linear 1D medium lead to dispersion relations with unstable modes. This happens in both the LY and the FC approach. Manuel has a proof that the instability appears whenever both \(Z\) and \(c\) vary in space (at least for layered media).

What is the cause of the instability? It can be shown that the unstable modes have large wavenumber – at least \(O(\delta)\). These modes cannot possibly be accurately captured by the homogenized equations, so it is unsurprising that they are represented in an unstable way. We can run stable spectral simulations of the linear equations by just keeping the grid coarse enough so that the unstable modes aren’t resolved. For nonlinear simulations, we will probably need to filter.

Perhaps the unstable modes correspond to those whose group velocity is larger than the maximum physical velocity in the medium? We should check.

Similar behavior is observed in models for crystals - see CLL

Regarding stable modes and wave lengths

From the symmetry section, we consider the combined PDE representing the linear homogenized PDE for \(\sigma\):

\(K_h^{-1}\rho_m\sigma_{tt}-\sigma_{xx}=\delta^2\left(\alpha_1+\gamma_1\right)\sigma_{xxxx}+O\left(\delta^4\right),\)

where:

\(\alpha_1+\gamma_1=\frac{\left(Z_A^2-Z_B^2\right)^2}{192K_m^2\rho_m^2}.\)

Now, to get the dispersion relation assume \(\sigma(x,t)=\sigma_0e^{i(kx-\omega t)}\) and plug it back into the previous PDE to get the dispersion relation:

\(\omega^2-\hat{c}^2k^2+\hat{c}^2\delta^2(\alpha_1+\gamma_1)k^4=0,\)

where \(\hat{c}=\sqrt{\frac{K_h}{\rho_m}}\). The dispersion relation can also be written as:

\(\omega=\pm\hat{c}k\sqrt{1-\delta^2(\alpha_1+\gamma_1)k^2}.\)

To have a purely wave solution and, moreover, to have a stable solution we need \(1-\delta^2(\alpha_1+\gamma_1)k^2\geq 0\), which implies:

\(k^2\leq\frac{1}{\delta^2(\alpha_1+\gamma_1)},\)

or, by using the expression for \(\alpha_1+\gamma_1\) from symmetry section:

\(k\leq\frac{1}{\delta}\frac{\sqrt{192}K_m\rho_m}{(Z_A^2-Z_B^2)}.\)

So by reducing the impedance contrast we get more stable modes. This can also be understood in terms of wavelength:

\(\lambda\geq\delta\frac{2\pi(Z_A^2-Z_B^2)}{\sqrt{192} K_m^2\rho_m^2},\)

which means that anything with a wavelength smaller is unstable. We can define a critical wave number \(k_*\) and a corresponding critical wave length \(\lambda_*\) as the largest wave number and smallest wave length to have the solution stable. They are given by:

\(k_*=\frac{1}{\delta}\frac{\sqrt{192}K_m\rho_m}{(Z_A^2-Z_B^2)}.\)

and

\(\lambda_*=\delta\frac{2\pi(Z_A^2-Z_B^2)}{\sqrt{192} K_m^2\rho_m^2}.\)

Regarding the velocity group

Considering the dispersion relation \(\omega=\omega(k)\) we can get the group velocity by:

\(v_g=\frac{\partial\omega}{\partial k}=\hat{c}\frac{1-2\delta^2(\alpha_1+\gamma_1)k^2}{\sqrt{1-\delta^2(\alpha_1+\gamma_1)k^2}}\rightarrow-\infty\) as \(k\rightarrow k_*\)

This predicts that the group velocity goes to zero and then becomes negative. Is this accurate, or an artifact of the homogenization process? We could test it using appropriate wavepacket initial conditions.