Stability criterion for waves in anisotropic solids

Question

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?

Answer 1

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).

Answer 2

In the case of an anisotropic material, the phase velocities of waves traveling through that material are determined by Christoffel's equation: \begin{equation} [\rho c^2\delta_{ij} - C_{ijkl}n_jn_l][u_k]=0 \end{equation} 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.

Answer 3

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}$.