Dispersion relations in a micropolar medium

We are interested in computing the dispersion relations in a homogeneous micropolar solid.

Wave propagation in micropolar solids

The equations of motion for a micropolar solid are given by [NOWACKI], [GUARIN]

\[\begin{split}\begin{align} &c_1^2 \nabla\nabla\cdot\mathbf{u}- c_2^2\nabla\times\nabla\times\mathbf{u} + K^2\nabla\times\boldsymbol{\theta} = -\omega^2 \mathbf{u} \, ,\\ &c_3^2 \nabla\nabla\cdot\boldsymbol{\theta} - c_4^2\nabla\times\nabla\times\boldsymbol{\theta} + Q^2\nabla\times\mathbf{u} - 2Q^2\boldsymbol{\theta} = -\omega^2 \boldsymbol{\theta} \, \end{align}\end{split}\]

where \(\mathbf{u}\) is the displacement vector and \(\boldsymbol{\theta}\) is the microrrotations vector, and where: \(c_1\) represents the phase/group speed for the longitudinal wave (\(P\)) that is non-dispersive as in the classical case, \(c_2\) represents the high-frequency limit phase/group speed for a transverse wave (\(S\)) that is dispersive unlike the classical counterpart, \(c_3\) represents the high-frequency limit phase/group speed for a longitudinal-rotational wave (\(LR\)) with a corkscrew-like motion that is dispersive and does not have a classical counterpart, \(c_4\) represents the high-frequency limit phase/group speed for a transverse-rotational wave (\(TR\)) that is dispersive and does not have a classical counterpart, \(Q\) represents the cut-off frequency for rotational waves appearance, and \(K\) quantifies the difference between the low-frequency and high-frequency phase/group speed for the S-wave. These parameters are defined by:

\[\begin{split}\begin{align} c_1^2 = \frac{\lambda +2\mu}{\rho},\quad &c_3^2 =\frac{\beta + 2\eta}{J},\\ c_2^2 = \frac{\mu +\alpha}{\rho},\quad &c_4^2 =\frac{\eta + \varepsilon}{J},\\ Q^2= \frac{2\alpha}{J},\quad &K^2 =\frac{2\alpha}{\rho} \, , \end{align}\end{split}\]

Dispersion relations

To identify types of propagating waves that can arise in the micropolar medium it is convenient to expand the displacement and rotation vectors in terms of scalar and vector potentials

\[\begin{split}\begin{align} \mathbf{u} &= \nabla \phi + \nabla\times\boldsymbol{\Gamma}\, ,\\ \boldsymbol{\theta} &= \nabla \tau + \nabla\times\mathbf{E}\, , \end{align}\end{split}\]

subject to the conditions:

\[\begin{split}\begin{align} &\nabla\cdot\boldsymbol{\Gamma} = 0\\ &\nabla\cdot\mathbf{E} = 0\, . \end{align}\end{split}\]

Using the above in the displacements equations of motion yields the following equations, after some manipulations

\[\begin{split}\begin{align} c_1^2 \nabla^2 \phi &= \frac{\partial^2 \phi}{\partial t^2}\, ,\\ c_3^2 \nabla^2 \tau - 2Q^2\tau &= \frac{\partial^2 \tau}{\partial t^2}\, ,\\ \begin{bmatrix} c_2^2 \nabla^2 &K^2\nabla\times\, ,\\ Q^2\nabla\times &c_4^2\nabla^2 - 2Q^2 \end{bmatrix} \begin{Bmatrix} \boldsymbol{\Gamma}\\ \mathbf{E}\end{Bmatrix} &= \frac{\partial^2}{\partial t^2} \begin{Bmatrix} \boldsymbol{\Gamma}\\ \mathbf{E}\end{Bmatrix} \, , \end{align}\end{split}\]

where we can see that the equations for the scalar potentials are uncoupled, while the ones for the vector potentials are coupled.

Writing the vector potentials as plane waves of amplitude \(\mathbf{A}\) and \(\mathbf{B}\), wave number \(\kappa\) and circular frequency \(\omega\) that propagate along the (x) axis,

\[\begin{split}\begin{align} \boldsymbol{\Gamma} &= \mathbf{A}\exp(i\kappa x - i\omega t)\\ \mathbf{E} &= \mathbf{B}\exp(i\kappa x - i\omega t)\, . \end{align}\end{split}\]

We can do these calculations using some the functions available functions in the package.

from sympy import Matrix, diff, symbols, exp, I, sqrt
from sympy import simplify, expand, solve, limit
from sympy import init_printing, pprint, factor
from continuum_mechanics.vector import lap_vec, curl, div

A1, A2, A3, B1, B2, B3 = symbols("A1 A2 A3 B1 B2 B3")
kappa, omega, t, x = symbols("kappa omega t x")
c1, c2, c3, c4, K, Q = symbols("c1 c2 c3 c4 K Q", positive=True)

We define the vector potentials \(\boldsymbol{\Gamma}\) and \(\mathbf{E}\).

Gamma = Matrix([A1, A2, A3]) * exp(I*kappa*x - I*omega*t)
E = Matrix([B1, B2, B3]) * exp(I*kappa*x - I*omega*t)

And compute the equations using the vector operators. Namely, the Laplace (vector.lap_vec()) and curl (vector.curl()) operators.

eq1 = c2**2 * lap_vec(Gamma) + K**2*curl(E) - Gamma.diff(t, 2)
eq2 = Q**2 * curl(Gamma) + c4**2*lap_vec(E) - 2*Q**2*E - E.diff(t, 2)
eq1 = simplify(eq1/exp(I*kappa*x - I*omega*t))
eq2 = simplify(eq2/exp(I*kappa*x  - I*omega*t))
eq = eq1.col_join(eq2)

We can compute the matrix for the system using .jacobian()

M = eq.jacobian([A1, A2, A3, B1, B2, B3])
M

\[\begin{split}\left[\begin{matrix}- c_{2}^{2} \kappa^{2} + \omega^{2} & 0 & 0 & 0 & 0 & 0\\0 & - c_{2}^{2} \kappa^{2} + \omega^{2} & 0 & 0 & 0 & - i K^{2} \kappa\\0 & 0 & - c_{2}^{2} \kappa^{2} + \omega^{2} & 0 & i K^{2} \kappa & 0\\0 & 0 & 0 & - 2 Q^{2} - c_{4}^{2} \kappa^{2} + \omega^{2} & 0 & 0\\0 & 0 & - i Q^{2} \kappa & 0 & - 2 Q^{2} - c_{4}^{2} \kappa^{2} + \omega^{2} & 0\\0 & i Q^{2} \kappa & 0 & 0 & 0 & - 2 Q^{2} - c_{4}^{2} \kappa^{2} + \omega^{2}\end{matrix}\right]\end{split}\]

And, we are interested in the determinant of the matrix \(M\).

factor(M.det())

\[\left(c_{2} \kappa - \omega\right) \left(c_{2} \kappa + \omega\right) \left(2 Q^{2} + c_{4}^{2} \kappa^{2} - \omega^{2}\right) \left(- K^{2} Q^{2} \kappa^{2} + 2 Q^{2} c_{2}^{2} \kappa^{2} - 2 Q^{2} \omega^{2} + c_{2}^{2} c_{4}^{2} \kappa^{4} - c_{2}^{2} \kappa^{2} \omega^{2} - c_{4}^{2} \kappa^{2} \omega^{2} + \omega^{4}\right)^{2}\]

The roots for this polynomial (in \(\omega^2\)) represent the dispersion relations.

disps = solve(M.det(), omega**2)
for disp in disps:
    display(disp)

\[c_{2}^{2} \kappa^{2}\]

\[2 Q^{2} + c_{4}^{2} \kappa^{2}\]

\[Q^{2} + \frac{c_{2}^{2} \kappa^{2}}{2} + \frac{c_{4}^{2} \kappa^{2}}{2} - \frac{1}{2} \sqrt{4 K^{2} Q^{2} \kappa^{2} + 4 Q^{4} - 4 Q^{2} c_{2}^{2} \kappa^{2} + 4 Q^{2} c_{4}^{2} \kappa^{2} + c_{2}^{4} \kappa^{4} - 2 c_{2}^{2} c_{4}^{2} \kappa^{4} + c_{4}^{4} \kappa^{4}}\]

\[Q^{2} + \frac{c_{2}^{2} \kappa^{2}}{2} + \frac{c_{4}^{2} \kappa^{2}}{2} + \frac{1}{2} \sqrt{4 K^{2} Q^{2} \kappa^{2} + 4 Q^{4} - 4 Q^{2} c_{2}^{2} \kappa^{2} + 4 Q^{2} c_{4}^{2} \kappa^{2} + c_{2}^{4} \kappa^{4} - 2 c_{2}^{2} c_{4}^{2} \kappa^{4} + c_{4}^{4} \kappa^{4}}\]

References

[NOWACKI]

Nowacki, W. (1986). Theory of asymmetric elasticity. pergamon Press, Headington Hill Hall, Oxford OX 3 0 BW, UK, 1986.

[GUARIN]

Guarín-Zapata, N., Gomez, J., Valencia, C., Dargush, G. F., & Hadjesfandiari, A. R. (2020). Finite element modeling of micropolar-based phononic crystals. Wave Motion, 92, 102406.