# Stability criterion for waves in anisotropic solids

The equations of motion for an elastic solid are given by

\begin{align} &\nabla \cdot \boldsymbol{\sigma} + \mathbf{f} = \rho \ddot{\mathbf{u}}\\ &\boldsymbol{\sigma} = \mathbb{C}\boldsymbol{\varepsilon}\\ &\boldsymbol{\varepsilon} = \frac{1}{2}\left(\nabla \mathbf{u} + [\nabla\mathbf{u}]^T\right) \end{align}

or in index notation

\begin{align} &\sigma_{ij,j} + f_i = \rho \ddot{u_i}\\ &\sigma_{ij} = C_{ijkl}\varepsilon_{kl}\\ &\varepsilon = \frac{1}{2}(u_{i,j} + u_{j,i}) \end{align}

$\mathbf{u}$ is the displacement vector, $\mathbf{f}$ is the body force (source term), $\boldsymbol{\sigma}$ is the stress tensor, $\mathbf{\varepsilon}$ is the strain tensor, and $\mathbb{C}$ is the stiffness tensor. In the case of isotropic solids the stiffness tensor is written in terms of two different constants, for unbounded domains, the equation admits two type of waves that are uncoupled and the criterion for stability is given by the worst case of the two different cases (i.e., the one with higher speed).

For transverse isotropic materials there are 5 independent parameters that define the tensor, and 3 types of waves (2 of them are coupled). In the more general case the number of parameters is 21 and the wave are coupled.

Question: How do you find the criterion for stability in a time marching algorithm for the general case?

• How do you define "stability" and what kind of criteria are you looking for? Apr 29 '15 at 3:04
• I'm looking for stability in an explicit numerical scheme. The equivalent of the CFL condition. Apr 29 '15 at 3:05

Wave equations like this can be rewritten as a hyperbolic system of first-order conservation laws:

$$q_t + \nabla \cdot F(q) = 0.$$

The stable time step for any explicit numerical discretization depends on the CFL number, which is proportional to the maximum wave speed appearing in the problem. That speed can be found by computing the eigenvalues of the Jacobian of the flux function, $F$ and taking the largest (in absolute value).

In $d$ dimensions, $F$ has $d$ components and the proper analysis requires finding the eigenvalues of arbitrary linear combinations of those components. But for most systems, including elasticity, it is sufficient to look at the eigenvalues of each of the components of $F$.

A nice reference for this theory is LeVeque's text on finite volume methods. Elasticity is covered in detail in Chapter 22.

All the usual caveats regarding the CFL condition apply: it is a necessary but not usually sufficient condition for stability. But the sufficient condition for stability of a given discretization is generally given by the CFL condition multiplied by some constant. To find your stability criterion, you need to know both the maximum wave speed (based on the equations you're solving) and that constant (based on the discretization you're using).

• thanks for your answer. I studied the chapter in Leveque's book, and deduced the form of the matrices $A$, $B$ and $C$. And I reduced the characteristic polynomials to cubic equations, that, in general, are cumbersome to solve. To check your suggestion I simplify the equations to the transverse isotropic case, and the eigenvalues of each one are not the maxima of the speeds. I then computed the eigenvalues for linear combinations and they change depending on the direction (as described in the book). Is it necessary to find the maximum for all possible directions then? Apr 29 '15 at 22:01

In the case of an anisotropic material, the phase velocities of waves traveling through that material are determined by Christoffel's equation: $$[\rho c^2\delta_{ij} - C_{ijkl}n_jn_l][u_k]=0$$ Here, $\rho$ is the density of the material, $c$ is the propagation velocity of the wave, $\delta$ is Kronecker's delta, $n_j$ is the unit vector in the direction $j$.

It can be seen that the solution corresponds to taking the determinant of the brackets on the left, which gives an eigenvalue equation in the eigenvalues $\rho c^2$.

Thus, in the three-dimensional anisotropic case, we still have three different phase velocities for a particular direction of propagation, the largest of which has to be used for CFL analysis, in a manner similar to the way the longitudinal speed is used in an isotropic problem.

• thanks for your answer. I am aware that the Christoffel equation contains the information of the wave speed, although I will add that in order to find the maximum you need to solve for all directions (using spherical coordinates) $$\max_{\theta, \phi} \max_i{\lambda_i(\theta, \phi)}\\$$ with $\lambda_i$ eigenvalues of the Christoffel equation. My question was more related with the CFL condition for this case (I was not sure if it was just computing the maximum speed for the material). Apr 30 '15 at 19:14

I will expand the answer provided by @DavidKetcheson. First the equations are rewritten as a hyperbolic system of first-order conservation laws:

$$q_t + \nabla \cdot F(q) = 0$$

or

$$q_t + A q_x + B q_y + C q_z = 0$$

Where $q$ is a state vector formed with the components of the stress tensor $(\sigma_{11}, \sigma_{22}, \sigma_{33}, \sigma_{12}, \sigma_{23}, \sigma_{13})$ and components of the velocity vector $(u, v, w)$.

$$q = \begin{pmatrix} \sigma_{11}\\ \sigma_{22}\\ \sigma_{33}\\ \sigma_{12}\\ \sigma_{23}\\ \sigma_{13}\\ u\\ v\\ w \end{pmatrix} \enspace ,$$

$$A = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 & {c}_{11} & {c}_{16} & {c}_{15}\\ 0 & 0 & 0 & 0 & 0 & 0 & {c}_{12} & {c}_{26} & {c}_{25}\\ 0 & 0 & 0 & 0 & 0 & 0 & {c}_{13} & {c}_{36} & {c}_{35}\\ 0 & 0 & 0 & 0 & 0 & 0 & {c}_{14} & {c}_{46} & {c}_{45}\\ 0 & 0 & 0 & 0 & 0 & 0 & {c}_{15} & {c}_{56} & {c}_{55}\\ 0 & 0 & 0 & 0 & 0 & 0 & {c}_{16} & {c}_{66} & {c}_{56}\\ \frac{1}{\rho} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & \frac{1}{\rho} & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & \frac{1}{\rho} & 0 & 0 & 0 \end{pmatrix} \enspace ,$$

$$B = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 & {c}_{16} & {c}_{12} & {c}_{14}\\ 0 & 0 & 0 & 0 & 0 & 0 & {c}_{26} & {c}_{22} & {c}_{24}\\ 0 & 0 & 0 & 0 & 0 & 0 & {c}_{36} & {c}_{23} & {c}_{34}\\ 0 & 0 & 0 & 0 & 0 & 0 & {c}_{46} & {c}_{24} & {c}_{44}\\ 0 & 0 & 0 & 0 & 0 & 0 & {c}_{56} & {c}_{25} & {c}_{45}\\ 0 & 0 & 0 & 0 & 0 & 0 & {c}_{66} & {c}_{26} & {c}_{46}\\ 0 & 0 & 0 & \frac{1}{\rho} & 0 & 0 & 0 & 0 & 0\\ 0 & \frac{1}{\rho} & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & \frac{1}{\rho} & 0 & 0 & 0 & 0 \end{pmatrix} \enspace ,$$

$$C = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 & {c}_{15} & {c}_{14} & {c}_{13}\\ 0 & 0 & 0 & 0 & 0 & 0 & {c}_{25} & {c}_{24} & {c}_{23}\\ 0 & 0 & 0 & 0 & 0 & 0 & {c}_{35} & {c}_{34} & {c}_{33}\\ 0 & 0 & 0 & 0 & 0 & 0 & {c}_{45} & {c}_{44} & {c}_{34}\\ 0 & 0 & 0 & 0 & 0 & 0 & {c}_{55} & {c}_{45} & {c}_{35}\\ 0 & 0 & 0 & 0 & 0 & 0 & {c}_{56} & {c}_{46} & {c}_{36}\\ 0 & 0 & 0 & 0 & 0 & \frac{1}{\rho} & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & \frac{1}{\rho} & 0 & 0 & 0 & 0\\ 0 & 0 & \frac{1}{\rho} & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix} \enspace .$$

In order to compute the speeds of the problem (as described above) we need to form the matrix $\hat{A}(n_1, n_2, n_3) = n_1 A + n_2 B + n_3 C$, where $\mathbb{n}=(n_1, n_2, n_3)$ is a unit vector and determines the direction of propagation. To find the CFL condition it is necessary to solve

$$\max_{(\theta, \phi)} \max_i \gamma_i(\theta, \phi)$$

where $(\theta, \phi)$ are spherical angles and $\gamma_i$ are the eigenvalues of the matrix $\hat{A}(\theta, \phi)$.

Based on this, and the answer provided by @DavidKetcheson it is simpler to compute the eigenvalues of the Christoffel equation, and solve the optimization problem

$$\max_{(\theta, \phi)} \max_i \lambda_i(\theta, \phi)$$

with $\lambda_i$ eigenvalues of the Christoffel equation. And the speed is just $c = \sqrt{\lambda_i/\rho}$.