**1. Introduction**

200 Wave Processes in Classical and New Solids

pp. 12-25.

Trudeau, P. J., Whitman, R. V., and Christian, J. T. (1974), "Shear Wave Velocity and Modulus of a Marine Clay," *Journal of the Boston Society of Civil Engineers*, Vol. 61, No. 1,

Viggiani, G. (1992). Small strain stiffness of fine grained soils. Ph.D. Thesis, City University. Viggiani, G. and Atkinson, J. H. (1995a), "Stiffness of Fine-Grained Soil at Very Small

Viggiani, G. and Atkinson, J. H. (1995b), "Interpretation of Bender Element Tests,"

Viggiani, G. and Atkinson, J. H. (1997), Discussion: "Interpretation of Bender Element

Vucetic, M. and Dobry, R. (1991), "Effect of Soil Plasticity on Cyclic Response," *Journal of* 

Weaver, E., Ormsby, W. C., and Zhao, X. (2001), *Pavement Variability Study: Phase I Interim Report*, National Pooled Fund Study 2(212), Non-Nuclear Testing of Soils and Granular Bases Using the GeoGauge (Soil Stiffness Gauge) and Other Similar Devices, U.S.

Wilson, S. D. and Dietrich, R. J. (1960), "Effect of Consolidation Pressure on Elastic and Strength Properties of Clay," *Proceedings of the Soil Mechanics and Foundations Division, Research Conference on Shear Strength of Cohesive Soils*, ASCE, Boulder, CO, pp. 419-435. Woods, R. D. (1978), "Measurement of Dynamic Soil Properties," *Proceedings of the Earthquake Engineering and Soil Dynamics Specialty Conference*, ASCE, Pasadena, CA, Vol. 1, pp. 91-178. Wright, S. G., Stokoe, K. E., II, and Roesset, J. M. (1994), "SASW Measurements at Geotechnical Sites Overlaid by Water," *Dynamic Geotechnical Testing II, ASTM STP 1213*,

Wu, W., Arellano, M., Chen, D.-H., Bilyeu, J., and He., R. (1998), "Using a Stiffness Gauge as an Alternative Quality Control Device in Pavement Construction," *Texas Department of* 

Yesiller, N., Inci, G., and Miller, C. J. (2000), "Ultrasonic Testing for Compacted Clayey Soils," *Proceedings of Sessions of Geo-Denver 2000, Advances in Unsaturated Geotechnics,* 

Yu, P. and Richart, F. E., Jr. (1984), "Stress Ratio Effects on Shear Modulus of Dry Sands," *Journal of the Geotechnical Engineering Division*, ASCE, Vol. 110, No. GT3, pp. 331-345. Yun, T. S. and Santamarina, J. C. (2005), "Decementation, Softening, and Collapse: Changes in Small-Strain Shear Stiffness in Ko Loading," *Journal of Geotechnical and* 

Zen, K., Umehara, Y., and Hamada, K. (1978), "Laboratory Tests and In Situ Seismic Survey on Vibratory Shear Modulus of Clayey Soils with Various Plasticities," *Proceedings of the* 

Zeng, X. and Ni, B. (1998), "Application of Bender Elements in Measuring Gmax of Sand under Ko Condition," *Geotechnical Testing Journal*, ASTM, Vol. 21, No. 3, pp. 251-263. Zeng, X., Figueroa, J. L., and Fu, L. (2002). Measurement of base and subgrade layer stiffness using bender element technique. Proc. of the 15th ASCE Engineering Mechanics

ASCE, Geotechnical Special Publication, No. 99, Denver, CO, pp. 54-68.

*Geoenvironmental Engineering Division*, ASCE, Vol. 131, No. 3, pp. 350-358.

*5th Japanese Earthquake Engineering Symposium*, pp. 721-728.

Conference, Columbia University, New York, NY, 1-10.

Strains," *Geotechnique*, Vol. 45, No. 2, pp. 249-265.

Tests," *Geotechnique*, Vol. 47, No. 4, pp. 873-877.

*Geotechnical Engineering*, ASCE, Vol. 117, No. 1, pp. 89-107.

Department of Transportation, Federal Highway Administration.

*Geotechnique*, Vol. 45, No. 1, pp. 149-154.

Philadelphia, PA, pp. 39-57.

*Transportation Report*, Austin, TX.

Conventionally, the second-order elastodynamics equation has been employed to model waves in elastic solids:

$$
\rho \frac{\partial^2 \mathbf{w}}{\partial t^2} = \nabla \cdot \left( c^{[4]} \nabla \mathbf{w} \right) \,, \tag{1.1}
$$

where *ρ* is the density of the medium, **w** the displacement, and *c*[4] the fourth-order stiffness tensor [2]. Equation (1.1) has been derived based on the equation of motion in conjunction with the elastic constitutive equation. Equation (1.1) has been solved by the finite-difference methods, e.g., [21], and the time-domain finite-element methods, e.g., [33], for propagating waves, and the frequency-domain finite-element methods, e.g., [6], for normal mode analysis of standing waves.

Alternatively, one can use a modern numerical method [11, 12] to solve the first-order velocity-stress equations to obtain the transient solution of wave propagation in elastic solids. For example, Virieux [25] modeled waves in earth crust by using a finite-difference method to solve the velocity-stress equations. LeVeque [12, 13] simulated stress waves in isotropic, elastic solids by using a upwind scheme. Shorr [22] developed a specialized formulation of the finite-element method to simulate propagating waves. Käser and Dumbser [10] used the discontinuous Galerkin method [20] to model waves in anisotropic earth crust. Yu et al. [32] simulated non-linear stress-wave propagation in hypo-elastic media by using the space-time Conservation Element and Solution Element (CESE) method [4]. In the setting of the second-order wave equation, Eq. (1.1), ample theoretical analyses about anisotropic elasticity can be found in the literature, e.g., [2, 24]. In general, mathematical properties of the first-order velocity-stress equations must be similar to that of the second-order wave equation Eq. (1.1). However, only limited discussions, e.g., Puente et al. [18] and Yang et al. [30], about the eigen structure of the velocity-stress equations can be found in the literature.

The objective of the present chapter is to analyze the mathematical properties of the first-order velocity-stress equations. In particular, their connection to the conventional second-order

©2012 Yu et al., licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0),which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2012 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

wave equations, Eq. (1.1) is of interest. To this end, this paper will be focused on analyzing the eigen structure of the velocity-stress equations. The approach employed is identical to that demonstrated in Warming et al. [29]. To be general, the derivation will be based on the properties of a triclinic solid, which have 21 independent stiffness constants and no symmetry. As such, the result will be applicable to all anisotropic, elastic solids. Moreover, in an one-dimensional space, the Jacobian matrix of the velocity-stress equations is shown to be equivalent to the classical Christoffel matrix of the second-order wave equations. In contrast to the conventional approach, however, a harmonic solution of plane waves is not assumed. To demonstrate the capabilities of the new formulation, the CESE method is applied to solve the velocity-stress equations for modeling wave propagation in a block of beryl, a solid with hexagonal symmetry. Wave expansion from a point source is demonstrated. The calculated group velocity profile compared well with the analytical solution. Additional result of waves interacting with interfaces separating regions with different lattice orientations are also reported. All results show salient features of wave propagation.

In the remainder of this section of Introduction, the equations of elastodynamics are reviewed in order to stage the further derivation in the following sections. First, the equation of motion without the body force is considered:

$$\frac{\partial \mathbf{v}}{\partial t} = \frac{1}{\rho} \nabla \cdot \mathbf{T}\_{\prime} \tag{1.2}$$

The rest of this chapter is organized as in the following. Section 2 shows the first-order velocity-stress equations in a vector-matrix form. Section 3 proves that the velocity-stress equations are hyperbolic. Section 4 shows that the Jacobian matrix of the one-dimensional velocity-stress equations is equivalent to the classic Christoffel matrix. As a concrete example, Section 5 illustrates a special form of the velocity-stress equations for elastodynamics in beryl, a solid of hexagonal symmetry. Section 6 reports the numerical solution of wave propagation in a block of beryl. The solution is obtained by applying the CESE method to solve the velocity-stress equations. In Section 7 we provide the concluding remarks. At the end of

To proceed, a Cartesian coordinate system is employed to formulate the velocity-stress

*<sup>K</sup>*(*μ*) *<sup>∂</sup>***<sup>T</sup>** *∂xμ*

*<sup>K</sup>*(*μ*)*<sup>t</sup> <sup>∂</sup>***<sup>v</sup>** *∂xμ*

*δμ*<sup>1</sup> 000 *δμ*<sup>3</sup> *δμ*<sup>2</sup> 0 *δμ*<sup>2</sup> 0 *δμ*<sup>3</sup> 0 *δμ*<sup>1</sup> 0 0 *δμ*<sup>3</sup> *δμ*<sup>2</sup> *δμ*<sup>1</sup> 0

) must be measured based on the same Cartesian coordinates employed in formulating Eq. (2.1). This coordinate system could be arbitrary. However, the entries of *C* reported in the literature, e.g., Auld [2], are usually tabulated based on a chosen coordinate system, which is

> *<sup>A</sup>*(*μ*) *<sup>∂</sup>***<sup>u</sup>** *∂xμ*

> > *ρK*(*μ*)

�

, *A*(2), and *A*(3) can be easily reduced to their Echelon forms [23] and the result

def = (*v*1, *<sup>v</sup>*2, *<sup>v</sup>*3, *<sup>T</sup>*1, *<sup>T</sup>*2, *<sup>T</sup>*3, *<sup>T</sup>*4, *<sup>T</sup>*5, *<sup>T</sup>*6)*<sup>t</sup>*

In Eq. (2.2), *δij* is the Kronecker delta. The rank, i.e., the dimension of the column and row spaces [23], of *K*(*μ*) with *μ* = 1, 2, 3 is three. In the above equations, the value of each entry in *C*

To proceed, the velocity-stress equations, Eq. (2.1), is recast into a vector-matrix form:

3 ∑ *μ*=1

� 03×<sup>3</sup> <sup>−</sup> <sup>1</sup>

<sup>−</sup>*CK*(*μ*)*<sup>t</sup>* <sup>06</sup>×<sup>6</sup>

*∂***u** *<sup>∂</sup><sup>t</sup>* <sup>+</sup>

where 0*m*×*<sup>n</sup>* denotes a *m* × *n* null matrix. Since the rank of matrices *K*(1)

= 0,

= 0,

⎞

Velocity-Stress Equations for Wave Propagation in Anisotropic Elastic Media 203

⎠ . (2.2)

= 0. (2.3)

, *μ* = 1, 2, 3, (2.4)

, *K*(2)

, *A*(2)

, and *A*(3)

, and *K*(3) is 3,

. Matrices *A*(1)

, *A*(2), and *A*(3) is 6.

(2.1)

3 ∑ *μ*=1

3 ∑ *μ*=1

the chapter, a list of cited references and several appendices are attached.

*∂***v** *<sup>∂</sup><sup>t</sup>* <sup>−</sup> <sup>1</sup> *ρ*

*∂***T** *<sup>∂</sup><sup>t</sup>* <sup>−</sup> *<sup>C</sup>*

*K*(*μ*) def = ⎛ ⎝

aligned with the lattice orientation of the solid of interest.

*A*(*μ*) def =

shows that the rank of each of the three 9 × 9 matrices *A*(1)

**2. The model equation**

equations Eqs. (1.2) and (1.4):

where the unknown vector **u**

matrices *A*(1)

are the three Jacobian matrices:

(or *c*[4]

where *K*(*μ*) with *μ* = 1, 2, 3 are 3 × 6 matrices:

where **v** is the velocity vector and *T* is the Cauchy stress tensor. *T* can be related to the strain tensor *S* by the elastic relation:

$$T = c^{[4]} S\_\prime \tag{1.3}$$

where *c*[4] is the fourth-order stiffness tensor. To proceed, Eq. (1.3) is differentiated by time *t* to yield:

$$\frac{\partial T}{\partial t} = \frac{1}{2} c^{[4]} \left( \nabla \mathbf{v} + (\nabla \mathbf{v})^{\mathbf{f}} \right) , \tag{1.4}$$

where the superscript *t* denotes transpose. The model equations include Eqs. (1.2) and (1.4). In an three-dimensional space, there are nine independent equations with the velocity and stress components as the unknowns. Next, the Voigt notation [2] is applied and the second-order tensors in Eq. (1.4) is changed to be 6-component vectors:

$$
\begin{pmatrix} T\_{11} & T\_{12} & T\_{13} \\ T\_{12} & T\_{22} & T\_{23} \\ T\_{13} & T\_{23} & T\_{33} \end{pmatrix} \rightarrow \begin{pmatrix} T\_{1} & T\_{6} & T\_{5} \\ T\_{6} & T\_{2} & T\_{4} \\ T\_{5} & T\_{4} & T\_{3} \end{pmatrix},
$$

$$
\begin{pmatrix} S\_{11} & S\_{12} & S\_{13} \\ S\_{12} & S\_{22} & S\_{23} \\ S\_{13} & S\_{23} & S\_{33} \end{pmatrix} \rightarrow \begin{pmatrix} S\_{1} & S\_{6}/2 & S\_{5}/2 \\ S\_{6}/2 & S\_{2} & S\_{4}/2 \\ S\_{5}/2 & S\_{4}/2 & S\_{3} \end{pmatrix}.
\tag{1.5}
$$

In the following, boldfaced **T** = (*T*1, *T*2, *T*3, *T*4, *T*5, *T*6)*<sup>t</sup>* and **S** = (*S*1, *S*2, *S*3, *S*4, *S*5, *S*6)*<sup>t</sup>* denote the 6-component stress and strain vectors. The non-boldfaced *T* and *S* in Eqs. (1.2) (1.3), and (1.4) denote the second-order stress and strain tensors. Aided by the Voigt notation, the elastic constitutive relation is rewritten as **T** = *C***S**, where *C* is the second-order 6 × 6 stiffness matrix. The rest of this chapter is organized as in the following. Section 2 shows the first-order velocity-stress equations in a vector-matrix form. Section 3 proves that the velocity-stress equations are hyperbolic. Section 4 shows that the Jacobian matrix of the one-dimensional velocity-stress equations is equivalent to the classic Christoffel matrix. As a concrete example, Section 5 illustrates a special form of the velocity-stress equations for elastodynamics in beryl, a solid of hexagonal symmetry. Section 6 reports the numerical solution of wave propagation in a block of beryl. The solution is obtained by applying the CESE method to solve the velocity-stress equations. In Section 7 we provide the concluding remarks. At the end of the chapter, a list of cited references and several appendices are attached.

### **2. The model equation**

2 Will-be-set-by-IN-TECH

wave equations, Eq. (1.1) is of interest. To this end, this paper will be focused on analyzing the eigen structure of the velocity-stress equations. The approach employed is identical to that demonstrated in Warming et al. [29]. To be general, the derivation will be based on the properties of a triclinic solid, which have 21 independent stiffness constants and no symmetry. As such, the result will be applicable to all anisotropic, elastic solids. Moreover, in an one-dimensional space, the Jacobian matrix of the velocity-stress equations is shown to be equivalent to the classical Christoffel matrix of the second-order wave equations. In contrast to the conventional approach, however, a harmonic solution of plane waves is not assumed. To demonstrate the capabilities of the new formulation, the CESE method is applied to solve the velocity-stress equations for modeling wave propagation in a block of beryl, a solid with hexagonal symmetry. Wave expansion from a point source is demonstrated. The calculated group velocity profile compared well with the analytical solution. Additional result of waves interacting with interfaces separating regions with different lattice orientations are

In the remainder of this section of Introduction, the equations of elastodynamics are reviewed in order to stage the further derivation in the following sections. First, the equation of motion

where **v** is the velocity vector and *T* is the Cauchy stress tensor. *T* can be related to the strain

*T* = *c*[4]

where *c*[4] is the fourth-order stiffness tensor. To proceed, Eq. (1.3) is differentiated by time *t*

where the superscript *t* denotes transpose. The model equations include Eqs. (1.2) and (1.4). In an three-dimensional space, there are nine independent equations with the velocity and stress components as the unknowns. Next, the Voigt notation [2] is applied and the second-order

> ⎛ ⎝

> ⎛ ⎝

In the following, boldfaced **T** = (*T*1, *T*2, *T*3, *T*4, *T*5, *T*6)*<sup>t</sup>* and **S** = (*S*1, *S*2, *S*3, *S*4, *S*5, *S*6)*<sup>t</sup>* denote the 6-component stress and strain vectors. The non-boldfaced *T* and *S* in Eqs. (1.2) (1.3), and (1.4) denote the second-order stress and strain tensors. Aided by the Voigt notation, the elastic constitutive relation is rewritten as **T** = *C***S**, where *C* is the second-order 6 × 6 stiffness matrix.

∇**v** + (∇**v**)*<sup>t</sup>*

*T*<sup>1</sup> *T*<sup>6</sup> *T*<sup>5</sup> *T*<sup>6</sup> *T*<sup>2</sup> *T*<sup>4</sup> *T*<sup>5</sup> *T*<sup>4</sup> *T*<sup>3</sup> �

⎞ ⎠ ,

*S*<sup>1</sup> *S*6/2 *S*5/2 *S*6/2 *S*<sup>2</sup> *S*4/2 *S*5/2 *S*4/2 *S*<sup>3</sup>

⎞ ⎠ .

∇ · *T*, (1.2)

*S*, (1.3)

, (1.4)

(1.5)

*∂***v** *<sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>1</sup> *ρ*

also reported. All results show salient features of wave propagation.

*∂T <sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>1</sup> 2 *c*[4] �

*T*<sup>11</sup> *T*<sup>12</sup> *T*<sup>13</sup> *T*<sup>12</sup> *T*<sup>22</sup> *T*<sup>23</sup> *T*<sup>13</sup> *T*<sup>23</sup> *T*<sup>33</sup>

*S*<sup>11</sup> *S*<sup>12</sup> *S*<sup>13</sup> *S*<sup>12</sup> *S*<sup>22</sup> *S*<sup>23</sup> *S*<sup>13</sup> *S*<sup>23</sup> *S*<sup>33</sup> ⎞ ⎠ →

⎞ ⎠ →

tensors in Eq. (1.4) is changed to be 6-component vectors: ⎛ ⎝

> ⎛ ⎝

without the body force is considered:

tensor *S* by the elastic relation:

to yield:

To proceed, a Cartesian coordinate system is employed to formulate the velocity-stress equations Eqs. (1.2) and (1.4):

$$\begin{aligned} \frac{\partial \mathbf{v}}{\partial t} - \frac{1}{\rho} \sum\_{\mu=1}^{3} K^{(\mu)} \frac{\partial \mathbf{T}}{\partial \mathbf{x}\_{\mu}} &= \mathbf{0}, \\\frac{\partial \mathbf{T}}{\partial t} - \mathbf{C} \sum\_{\mu=1}^{3} K^{(\mu)} \frac{\partial \mathbf{v}}{\partial \mathbf{x}\_{\mu}} &= \mathbf{0}, \end{aligned} \tag{2.1}$$

where *K*(*μ*) with *μ* = 1, 2, 3 are 3 × 6 matrices:

$$K^{(\mu)} \stackrel{\text{def}}{=} \begin{pmatrix} \delta\_{\mu 1} & 0 & 0 & 0 & \delta\_{\mu 3} \ \delta\_{\mu 2} \\ 0 & \delta\_{\mu 2} & 0 & \delta\_{\mu 3} & 0 & \delta\_{\mu 1} \\ 0 & 0 & \delta\_{\mu 3} \ \delta\_{\mu 2} & \delta\_{\mu 1} & 0 \end{pmatrix} . \tag{2.2}$$

In Eq. (2.2), *δij* is the Kronecker delta. The rank, i.e., the dimension of the column and row spaces [23], of *K*(*μ*) with *μ* = 1, 2, 3 is three. In the above equations, the value of each entry in *C* (or *c*[4] ) must be measured based on the same Cartesian coordinates employed in formulating Eq. (2.1). This coordinate system could be arbitrary. However, the entries of *C* reported in the literature, e.g., Auld [2], are usually tabulated based on a chosen coordinate system, which is aligned with the lattice orientation of the solid of interest.

To proceed, the velocity-stress equations, Eq. (2.1), is recast into a vector-matrix form:

$$\frac{\partial \mathbf{u}}{\partial t} + \sum\_{\mu=1}^{3} A^{(\mu)} \frac{\partial \mathbf{u}}{\partial x\_{\mu}} = 0. \tag{2.3}$$

where the unknown vector **u** def = (*v*1, *<sup>v</sup>*2, *<sup>v</sup>*3, *<sup>T</sup>*1, *<sup>T</sup>*2, *<sup>T</sup>*3, *<sup>T</sup>*4, *<sup>T</sup>*5, *<sup>T</sup>*6)*<sup>t</sup>* . Matrices *A*(1) , *A*(2), and *A*(3) are the three Jacobian matrices:

$$A^{(\mu)} \stackrel{\text{def}}{=} \begin{pmatrix} \mathbf{0}\_{3 \times 3} & -\frac{1}{\rho} \mathbf{K}^{(\mu)} \\ -\mathbf{C} \mathbf{K}^{(\mu)^{\mathrm{I}}} & \mathbf{0}\_{6 \times 6} \end{pmatrix}, \quad \mu = 1, 2, 3,\tag{2.4}$$

where 0*m*×*<sup>n</sup>* denotes a *m* × *n* null matrix. Since the rank of matrices *K*(1) , *K*(2) , and *K*(3) is 3, matrices *A*(1) , *A*(2), and *A*(3) can be easily reduced to their Echelon forms [23] and the result shows that the rank of each of the three 9 × 9 matrices *A*(1), *A*(2) , and *A*(3) is 6.

### **3. Hyperbolicity of the velocity-stress equations**

In this section, the hyperbolicity of the velocity-stress equations Eq. (2.3) is proved based on analyzing a linear combination of the three Jacobian matrices:

$$B \stackrel{\text{def}}{=} \sum\_{\mu=1}^{3} h\_{\mu} A^{(\mu)} = h\_1 A^{(1)} + h\_2 A^{(2)} + h\_3 A^{(3)},\tag{3.1}$$

where *In* is a *n* × *n* identity matrix and *λ* the eigenvalue. We assume that eigenvalues of *B* are either zero or non-zero. If all non-zero eigenvalues of *B* are real, all eigenvalues of *B* are real. The existence of zero eigenvalues of *B* is self-evident because rank(*B*) ≤ 6 [23]. For non-zero

det(Γ − *ρλ*<sup>2</sup> *I*3) = 0, (3.7)

Velocity-Stress Equations for Wave Propagation in Anisotropic Elastic Media 205

= *GCG<sup>t</sup>* (3.8)

*B***q**+ = *λ*+**q**+, (3.9)

<sup>+</sup>**q***v*, (3.11)

**q***v*. (3.12)

(3.10)

eigenvalues, aided by the Schur complements, Eq. (3.6) can be transformed to be

eigenvalue *λ*+, there must be a negative counterpart *λ*− = −*λ*+.

eigenvalue. To proceed, we consider the eigen problem of matrix *B*:

 − <sup>1</sup> *<sup>ρ</sup> G***q***<sup>T</sup>* −*CG<sup>t</sup>* **q***v*

**q**+ def = **q***<sup>v</sup>* **q***T* 

the associated eigenvector containing a 3-component vector **q***<sup>v</sup>* and a 6-component vector **q***T*. Aided by the definition of matrix *B*, Eq. (3.5), the eigen problem, Eq. (3.9), can be recast to

= *λ*+

 **q***<sup>v</sup>* **q***T* ,

Γ**q***<sup>v</sup>* = *ρλ*<sup>2</sup>

**<sup>q</sup>***<sup>T</sup>* <sup>=</sup> <sup>−</sup> <sup>1</sup>

*λ*+ *CG<sup>t</sup>*

where *λ*+ is a positive eigenvalue of *B* and

Γ def

is a 3 × 3 matrix. For completeness, the Schur complements are provided in Appendix 8. To proceed, we recall that the stiffness matrix *C* is symmetric and positive-definite [17]. Therefore, matrix Γ must be symmetric and positive-definite. As a result, the eigenvalues of Γ, i.e., *ρλ*2, are real and positive. Since the density *ρ* must be positive, *λ*<sup>2</sup> is positive. Therefore, the non-zero eigenvalues of *B* must be real and in positive-negative pairs, i.e., for a positive

To recap, the spectrum of *B* includes the null eigenvalue and three non-zero, positive-negative pairs of eigenvalues. Eigenvalues in pairs could be repeated such that there could be at least one and at most three pairs of distinct non-zero eigenvalues. For isotropic solids, *B* has two pairs of distinct non-zero eigenvalues. The above proof of the real spectrum of *B* is underpinned by the property of matrix *C*, which is symmetry and positive-definiteness. As will be shown in the following, this property of *C* is equally important in constructing the

To proceed, we show that matrix *B* can be diagonalized by its eigenvector matrix. Essentially, we need to find nine linearly independent eigenvectors for *B*. First, we show that there exist six linearly independent eigenvectors associated with the non-zero eigenvalues. We then show the existence of three linearly independent eigenvectors associated with the null

where

eigenvector matrix of *B*.

which gives:

where *h*1, *h*2, and *h*<sup>3</sup> are arbitrary real numbers. Without losing any generality, we let *h*1, *h*2, and *h*<sup>3</sup> be the three components of a unit vector **h**, i.e., |**h**| = 1. Equation (2.3) is hyperbolic if all eigenvalues of *B* are real, and *B* can be diagonalized by its eigenvector matrix [11].

Physical meaning of the hyperbolicity of the velocity-stress equations can be illustrated by considering a plane wave propagating in the direction of **h**, along which Eq. (2.3) is reduced to be one-dimensional. In this uni-axial direction, a new coordinate is defined as *y* def = ∑<sup>3</sup> *<sup>μ</sup>*=<sup>1</sup> *hμxμ*. Aided by the new coordinate *y*, the three spatial derivative terms in Eq. (2.3) can be combined:

$$\sum\_{\mu=1}^{3} A^{(\mu)} \frac{\partial \mathbf{u}}{\partial \mathbf{x}\_{\mu}} = \sum\_{\mu=1}^{3} A^{(\mu)} \frac{\partial \mathbf{u}}{\partial y} \frac{\partial y}{\partial \mathbf{x}\_{\mu}} = \sum\_{\mu=1}^{3} h\_{\mu} A^{(\mu)} \frac{\partial \mathbf{u}}{\partial y} = B \frac{\partial \mathbf{u}}{\partial y}. \tag{3.2}$$

Aided by Eq. (3.2), a one-dimensional version of Eq. (2.3) for modeling wave propagating in the direction of **h** becomes

$$
\frac{
\partial \mathbf{u}
}{
\partial t
} + B \frac{
\partial \mathbf{u}
}{
\partial y
} = 0.
\tag{3.3}
$$

Therefore, hyperbolicity of Eq. (2.3) implies that the eigenvalues of the Jacobian matrix *B* of Eq. (3.3) are all real and they represent different wave speeds. Moreover, Eq. (3.3) can be diagonalized so that Eq. (3.3) can be decoupled into 9 independent scalar convection equations, each of which would propagate a constant profile, i.e., a Riemann invariant, in the direction **h** by its own wave speed, i.e., the corresponding eigenvalue.

In the following, we will show that all eigenvalues of *B* are real. Then we will show that *B* is diagonalizable by its eigenvector matrix. To proceed, we introduce the following definition:

$$\mathbf{G} \stackrel{\text{def}}{=} \sum\_{\mu=1}^{3} h\_{\mu} \mathbf{K}^{(\mu)} = \begin{pmatrix} h\_1 & 0 & 0 & 0 & h\_3 \ h\_2 \\ 0 & h\_2 & 0 & h\_3 & 0 & h\_1 \\ 0 & 0 & h\_3 & h\_2 & h\_1 & 0 \end{pmatrix} . \tag{3.4}$$

Recall that |**h**| = 1. Since *h*1, *h*2, and *h*<sup>3</sup> are not all-zero, the rank of matrix *G* is 3. Aided by this definition of *G*, matrix *B* can be written as

$$B = \begin{pmatrix} \mathbf{0}\_{3 \times 3} & -\frac{1}{\rho} G \\ -CG^t & \mathbf{0}\_{6 \times 6} \end{pmatrix}. \tag{3.5}$$

Because the rank of *G* is 3, the rank of the 6 × 3 matrix *CG<sup>t</sup>* must be less than or equal to 3, such that rank(*B*) ≤ 6. The eigenvalues of *B* can be obtained by solving its characteristic equation:

$$\det(B - \lambda I\_{\mathfrak{H}}) = 0,\tag{3.6}$$

where *In* is a *n* × *n* identity matrix and *λ* the eigenvalue. We assume that eigenvalues of *B* are either zero or non-zero. If all non-zero eigenvalues of *B* are real, all eigenvalues of *B* are real. The existence of zero eigenvalues of *B* is self-evident because rank(*B*) ≤ 6 [23]. For non-zero eigenvalues, aided by the Schur complements, Eq. (3.6) can be transformed to be

$$\det(\Gamma - \rho \lambda^2 I\_3) = 0,\tag{3.7}$$

where

4 Will-be-set-by-IN-TECH

In this section, the hyperbolicity of the velocity-stress equations Eq. (2.3) is proved based on

where *h*1, *h*2, and *h*<sup>3</sup> are arbitrary real numbers. Without losing any generality, we let *h*1, *h*2, and *h*<sup>3</sup> be the three components of a unit vector **h**, i.e., |**h**| = 1. Equation (2.3) is hyperbolic if

Physical meaning of the hyperbolicity of the velocity-stress equations can be illustrated by considering a plane wave propagating in the direction of **h**, along which Eq. (2.3) is reduced to

Aided by the new coordinate *y*, the three spatial derivative terms in Eq. (2.3) can be combined:

*∂y ∂xμ* = 3 ∑ *μ*=1

Aided by Eq. (3.2), a one-dimensional version of Eq. (2.3) for modeling wave propagating in

Therefore, hyperbolicity of Eq. (2.3) implies that the eigenvalues of the Jacobian matrix *B* of Eq. (3.3) are all real and they represent different wave speeds. Moreover, Eq. (3.3) can be diagonalized so that Eq. (3.3) can be decoupled into 9 independent scalar convection equations, each of which would propagate a constant profile, i.e., a Riemann invariant, in

In the following, we will show that all eigenvalues of *B* are real. Then we will show that *B* is diagonalizable by its eigenvector matrix. To proceed, we introduce the following definition:

> ⎛ ⎝

Recall that |**h**| = 1. Since *h*1, *h*2, and *h*<sup>3</sup> are not all-zero, the rank of matrix *G* is 3. Aided by

03×<sup>3</sup> − <sup>1</sup>

−*CG<sup>t</sup>* 06×<sup>6</sup>

Because the rank of *G* is 3, the rank of the 6 × 3 matrix *CG<sup>t</sup>* must be less than or equal to 3, such that rank(*B*) ≤ 6. The eigenvalues of *B* can be obtained by solving its characteristic equation:

*ρ G*

�

*h*<sup>1</sup> 000 *h*<sup>3</sup> *h*<sup>2</sup> 0 *h*<sup>2</sup> 0 *h*<sup>3</sup> 0 *h*<sup>1</sup> 0 0 *h*<sup>3</sup> *h*<sup>2</sup> *h*<sup>1</sup> 0

*<sup>h</sup>μA*(*μ*) *<sup>∂</sup>***<sup>u</sup>**

*<sup>∂</sup><sup>y</sup>* <sup>=</sup> *<sup>B</sup> <sup>∂</sup>***<sup>u</sup>**

*<sup>∂</sup><sup>y</sup>* <sup>=</sup> 0. (3.3)

⎞

det(*B* − *λI*9) = 0, (3.6)

all eigenvalues of *B* are real, and *B* can be diagonalized by its eigenvector matrix [11].

be one-dimensional. In this uni-axial direction, a new coordinate is defined as *y*

*<sup>A</sup>*(*μ*) *<sup>∂</sup>***<sup>u</sup>** *∂y*

> *∂***u** *<sup>∂</sup><sup>t</sup>* <sup>+</sup> *<sup>B</sup> <sup>∂</sup>***<sup>u</sup>**

the direction **h** by its own wave speed, i.e., the corresponding eigenvalue.

*hμK*(*μ*) =

*B* = �

*G* def = 3 ∑ *μ*=1

this definition of *G*, matrix *B* can be written as

*<sup>h</sup>μA*(*μ*) = *<sup>h</sup>*1*A*(1) + *<sup>h</sup>*2*A*(2) + *<sup>h</sup>*3*A*(3)

, (3.1)

def = ∑<sup>3</sup>

*<sup>∂</sup><sup>y</sup>* . (3.2)

⎠ . (3.4)

. (3.5)

*<sup>μ</sup>*=<sup>1</sup> *hμxμ*.

**3. Hyperbolicity of the velocity-stress equations**

analyzing a linear combination of the three Jacobian matrices:

*<sup>B</sup>* def = 3 ∑ *μ*=1

3 ∑ *μ*=1

the direction of **h** becomes

*<sup>A</sup>*(*μ*) *<sup>∂</sup>***<sup>u</sup>** *∂xμ* = 3 ∑ *μ*=1

$$
\Gamma \stackrel{\text{def}}{=} GCG^t \tag{3.8}
$$

is a 3 × 3 matrix. For completeness, the Schur complements are provided in Appendix 8. To proceed, we recall that the stiffness matrix *C* is symmetric and positive-definite [17]. Therefore, matrix Γ must be symmetric and positive-definite. As a result, the eigenvalues of Γ, i.e., *ρλ*2, are real and positive. Since the density *ρ* must be positive, *λ*<sup>2</sup> is positive. Therefore, the non-zero eigenvalues of *B* must be real and in positive-negative pairs, i.e., for a positive eigenvalue *λ*+, there must be a negative counterpart *λ*− = −*λ*+.

To recap, the spectrum of *B* includes the null eigenvalue and three non-zero, positive-negative pairs of eigenvalues. Eigenvalues in pairs could be repeated such that there could be at least one and at most three pairs of distinct non-zero eigenvalues. For isotropic solids, *B* has two pairs of distinct non-zero eigenvalues. The above proof of the real spectrum of *B* is underpinned by the property of matrix *C*, which is symmetry and positive-definiteness. As will be shown in the following, this property of *C* is equally important in constructing the eigenvector matrix of *B*.

To proceed, we show that matrix *B* can be diagonalized by its eigenvector matrix. Essentially, we need to find nine linearly independent eigenvectors for *B*. First, we show that there exist six linearly independent eigenvectors associated with the non-zero eigenvalues. We then show the existence of three linearly independent eigenvectors associated with the null eigenvalue. To proceed, we consider the eigen problem of matrix *B*:

$$B\mathfrak{q}\_{+} = \lambda\_{+}\mathfrak{q}\_{+\prime} \tag{3.9}$$

where *λ*+ is a positive eigenvalue of *B* and

$$\mathbf{q}\_{+} \stackrel{\text{def}}{=} \begin{pmatrix} \mathbf{q}\_{\upsilon} \\ \mathbf{q}\_{\mathsf{T}} \end{pmatrix} \tag{3.10}$$

the associated eigenvector containing a 3-component vector **q***<sup>v</sup>* and a 6-component vector **q***T*. Aided by the definition of matrix *B*, Eq. (3.5), the eigen problem, Eq. (3.9), can be recast to

$$
\begin{pmatrix} -\frac{1}{\rho} G \mathbf{q}\_T \\ -C G^t \mathbf{q}\_{\overline{v}} \end{pmatrix} = \lambda\_+ \begin{pmatrix} \mathbf{q}\_{\overline{v}} \\ \mathbf{q}\_T \end{pmatrix} \lambda\_-
$$

which gives:

$$
\Gamma \mathbf{q}\_{\upsilon} = \rho \lambda\_{+}^{2} \mathbf{q}\_{\upsilon} \tag{3.11}
$$

$$\mathbf{q}\_T = -\frac{1}{\lambda\_+} \mathbf{C} \mathbf{G}^t \mathbf{q}\_v. \tag{3.12}$$

By arbitrarily specifying a value to an entry of **q***v*, **q***v* can be obtained by solving Eq. (3.11). **q***<sup>T</sup>* can then be obtained by substituting **q***<sup>v</sup>* into Eq. (3.12). As such, the eigenvector **q**+ can be fully constructed.

Essentially, the eigen problem of the 9 × 9 matrix *B* in Eq. (3.9), is reduced to a eigen problem for the 3 × 3 matrix Γ in Eq. (3.11). Since Γ is a 3 × 3, symmetric, positive-definite matrix, the spectral theorem [23] of linear algebra asserts that Γ always has three linearly independent eigenvectors, regardless whether Γ has degenerate eigenvalues or not. These three eigenvectors of Γ are denoted by **q***v*,1, **q***v*,2, and **q***v*,3. The three associated eigenvectors of *B* with positive eigenvalues are denoted by **q**+,1, **q**+,2, and **q**+,3. Aided by Eq. (3.12), **q**+,1, **q**+,2, and **q**+,3 are readily determined by **q***v*,1, **q***v*,2, and **q***v*,3. Because **q***v*,1, **q***v*,2, and **q***v*,3 are linearly independent, **q**+,1, **q**+,2, and **q**+,3 must also be linearly independent. By following exactly the same procedure, one can derive the three independent eigenvectors **q**−,1, **q**−,2, and **q**−,3 associated with the three negative eigenvalues *λ*<sup>−</sup> = −*λ*+. To proceed, since *λ*− = −*λ*+ �= *λ*+,

$$\mathbf{q}\_{\pm,i} \stackrel{\text{def}}{=} \begin{pmatrix} \mathbf{q}\_{\upsilon,i} \\ \pm \mathbf{q}\_{T,i} \end{pmatrix}, \quad i = 1,2,3, \tag{3.13}$$

*μ* = 1, 2, 3. If we let *μ* = 1, it can be straightforwardly shown that the fifth, sixth, and seventh columns of *A*(1) are all-zero. For *μ* = 2, the fourth, sixth, and eighth columns of *A*(2) are all-zero. For *μ* = 3, the fourth, fifth, and ninth columns of *A*(3) are all-zero. For completeness,

0,1 = (0, 0, 0, 1, 0, 0, 0, 0, 0)*<sup>t</sup>*

0,2 = (0, 0, 0, 0, 0, 1, 0, 0, 0)*<sup>t</sup>*

0,3 = (0, 0, 0, 0, 0, 0, 0, 1, 0)*<sup>t</sup>*

Since Eq. (2.3) is hyperbolic, we proceed to show the characteristic form of Eq. (3.3). We

*<sup>∂</sup><sup>t</sup>* <sup>+</sup> *<sup>Q</sup>*−1*BQ <sup>∂</sup>Q*−1**<sup>u</sup>**

Note that a 9 × 9 identity matrix *I*<sup>9</sup> = *QQ*−<sup>1</sup> is inserted in the second term of Eq. (3.16). Since *Q*−<sup>1</sup> is a constant matrix, it can be moved into the partial differentiation operator. Aided by

where Λ is a diagonal matrix consisting of the eigenvalues of *B*. Equation (3.17) contains nine

where *λ<sup>i</sup>* is the *i*th eigenvalue of *B* and *u*ˆ*<sup>i</sup>* the *i*th characteristic variable. If Eq. (2.3) is hyperbolic, each scalar equation in Eq. (3.18) describes convection of a characteristic variable

In this section, we show how matrix *B* is connected to the classic Christoffel matrix. Essentially, we will show that matrix Γ in Eq. (3.8) is in fact the Christoffel matrix for plane waves propagating in the direction of **h**. Recall that Γ is obtained by solving the eigen problem for non-zero eigenvalues of *B*. To proceed, two Cartesian coordinate systems are considered: (i) The global coordinate system (*x*1-*x*2-*x*3), in which the governing equations Eqs. (2.3) and

**<sup>u</sup>**<sup>ˆ</sup> def = *Q*−1**u**,

*∂***u**ˆ *<sup>∂</sup><sup>t</sup>* <sup>+</sup> <sup>Λ</sup> *<sup>∂</sup>***u**<sup>ˆ</sup>

*∂u*ˆ*<sup>i</sup>*

, **<sup>q</sup>**(3)

Velocity-Stress Equations for Wave Propagation in Anisotropic Elastic Media 207

, **<sup>q</sup>**(3)

, **<sup>q</sup>**(3)

, *μ* = 1, 2, 3, are tabulated in

,

,

. (3.15)

0,1 = (0, 0, 0, 1, 0, 0, 0, 0, 0)*<sup>t</sup>*

0,2 = (0, 0, 0, 0, 1, 0, 0, 0, 0)*<sup>t</sup>*

0,3 = (0, 0, 0, 0, 0, 0, 0, 0, 1)*<sup>t</sup>*

*<sup>∂</sup><sup>y</sup>* <sup>=</sup> 0. (3.16)

*<sup>∂</sup><sup>y</sup>* <sup>=</sup> 0, (3.17)

*<sup>∂</sup><sup>y</sup>* <sup>=</sup> 0, *<sup>i</sup>* <sup>=</sup> 1, . . . , 9, (3.18)

the eigenvectors associated with the trivial eigenvalue for *A*(*μ*)

, **<sup>q</sup>**(2)

, **<sup>q</sup>**(2)

, **<sup>q</sup>**(2)

*∂Q*−1**u**

the eigenvalue matrix Λ and the characteristic variables

Eq. (3.16) can be rewritten in the characteristic form as

*∂u*ˆ*<sup>i</sup> <sup>∂</sup><sup>t</sup>* <sup>+</sup> *<sup>λ</sup><sup>i</sup>*

**4. Coordinate transformation and the Christoffel matrix**

the following:

0,1 = (0, 0, 0, 0, 1, 0, 0, 0, 0)*<sup>t</sup>*

0,2 = (0, 0, 0, 0, 0, 1, 0, 0, 0)*<sup>t</sup>*

0,3 = (0, 0, 0, 0, 0, 0, 1, 0, 0)*<sup>t</sup>*

pre-multiply Eq. (3.3) with *Q*−<sup>1</sup> to have

**<sup>q</sup>**(1)

**<sup>q</sup>**(1)

**q**(1)

decoupled equations:

*u*ˆ*<sup>i</sup>* with a real wave speed *λi*.

are the six linearly independent eigenvectors of *B*. The corresponding non-zero eigenvalues of *B* can be written as *λ*±,*i*, *i* = 1, 2, 3, with *λ*−,*<sup>i</sup>* = −*λ*+,*i*, *i* = 1, 2, 3, and *λ*+,1 ≤ *λ*+,2 ≤ *λ*+,3. We note that it is possible for *B* to have degenerate non-zero eigenvalues. Nevertheless, all six eigenvectors associated with the non-zero eigenvalues are linearly independent. This concludes the discussions about the eigenvectors of *B* associated with the non-zero eigenvalues.

The remaining three eigenvectors of *B* are associated with the null eigenvalue. Recall that *B* is a 9 × 9 matrix with rank(*B*) ≤ 6. Thus, the nullity of *B* must be greater than or equal to 3, i.e., there are at least three linearly independent vectors spanning the null space of *B* [23]. Aided by the definition of null space and eigenvector, the basis of the null space of *B* consists of the eigenvectors associated with the null eigenvalue of *B*. As such, there must be exactly three eigenvectors associated with the null eigenvalue of *B*, because (i) there are at most nine linearly independent eigenvectors of *B*, (ii) there are at least three linearly independent eigenvectors associated with the null eigenvalues of *B*, and (iii) there are exactly six linearly independent eigenvectors associated with the non-zero eigenvalues of *B*. These three eigenvectors are denoted by **q**0,1, **q**0,2, and **q**0,3. This concludes the discussion of three linearly independent eigenvectors associated with the null eigenvalue of *B*.

Together, there are nine linearly independent eigenvectors of *B*: **q**±,*i*, **q**0,*i*, *i* = 1, 2, 3. The diagonalizing eigenvector matrix of *B* can be constructed as:

$$\mathbf{Q} \stackrel{\text{def}}{=} \left( \mathbf{q}\_{-,3} \; \mathbf{q}\_{-,2} \; \mathbf{q}\_{-,1} \; \mathbf{q}\_{0,1} \; \mathbf{q}\_{0,2} \; \mathbf{q}\_{0,3} \; \mathbf{q}\_{+,1} \; \mathbf{q}\_{+,2} \; \mathbf{q}\_{+,3} \right), \tag{3.14}$$

where *Q* is a 9 × 9 matrix composed of nine column vectors. We then perform similarity transformation Λ = *Q*−1*BQ* to diagonalize *B*, where the eigenvalue matrix Λ is composed of *λ*−,3, *λ*−,2, *λ*−,1, *λ*0, *λ*0, *λ*0, *λ*+,1, *λ*+,2, and *λ*+,3 with *λ*<sup>0</sup> = 0. This concludes the discussions of diagonalizability of *B* by the above nine linearly independent eigenvectors.

As an aside, each of Jacobian matrices *A*(*μ*) , *μ* = 1, 2, 3, is a special case of *B* with *h<sup>ν</sup>* = *δμν*, *<sup>ν</sup>* = 1, 2, 3. For these special cases, we let **<sup>q</sup>**(*μ*) ±,*i* , **<sup>q</sup>**(*μ*) 0,*<sup>i</sup>* , *<sup>i</sup>* <sup>=</sup> 1, 2, 3, be the eigenvectors of *<sup>A</sup>*(*μ*) ,

*μ* = 1, 2, 3. If we let *μ* = 1, it can be straightforwardly shown that the fifth, sixth, and seventh columns of *A*(1) are all-zero. For *μ* = 2, the fourth, sixth, and eighth columns of *A*(2) are all-zero. For *μ* = 3, the fourth, fifth, and ninth columns of *A*(3) are all-zero. For completeness, the eigenvectors associated with the trivial eigenvalue for *A*(*μ*) , *μ* = 1, 2, 3, are tabulated in the following:

$$\begin{aligned} \mathbf{q}\_{0,1}^{(1)} &= (0,0,0,0,1,0,0,0,0)^{t}, \; \mathbf{q}\_{0,1}^{(2)} = (0,0,0,1,0,0,0,0,0)^{t}, \; \mathbf{q}\_{0,1}^{(3)} = (0,0,0,1,0,0,0,0,0)^{t}, \\ \mathbf{q}\_{0,2}^{(1)} &= (0,0,0,0,0,1,0,0,0)^{t}, \; \mathbf{q}\_{0,2}^{(2)} = (0,0,0,0,0,1,0,0,0)^{t}, \; \mathbf{q}\_{0,2}^{(3)} = (0,0,0,0,1,0,0,0,0)^{t}, \\ \mathbf{q}\_{0,3}^{(1)} &= (0,0,0,0,0,0,1,0,0)^{t}, \; \mathbf{q}\_{0,3}^{(2)} = (0,0,0,0,0,0,1,0)^{t}, \; \mathbf{q}\_{0,3}^{(3)} = (0,0,0,0,0,0,0,0,1)^{t}. \end{aligned} \tag{3.15}$$

Since Eq. (2.3) is hyperbolic, we proceed to show the characteristic form of Eq. (3.3). We pre-multiply Eq. (3.3) with *Q*−<sup>1</sup> to have

$$\frac{\partial Q^{-1}\mathbf{u}}{\partial t} + Q^{-1}BQ \frac{\partial Q^{-1}\mathbf{u}}{\partial y} = 0. \tag{3.16}$$

Note that a 9 × 9 identity matrix *I*<sup>9</sup> = *QQ*−<sup>1</sup> is inserted in the second term of Eq. (3.16). Since *Q*−<sup>1</sup> is a constant matrix, it can be moved into the partial differentiation operator. Aided by the eigenvalue matrix Λ and the characteristic variables

$$\hat{\mathbf{u}} \stackrel{\text{def}}{=} \mathbf{Q}^{-1} \mathbf{u}\_{\prime}$$

Eq. (3.16) can be rewritten in the characteristic form as

6 Will-be-set-by-IN-TECH

By arbitrarily specifying a value to an entry of **q***v*, **q***v* can be obtained by solving Eq. (3.11). **q***<sup>T</sup>* can then be obtained by substituting **q***<sup>v</sup>* into Eq. (3.12). As such, the eigenvector **q**+ can be

Essentially, the eigen problem of the 9 × 9 matrix *B* in Eq. (3.9), is reduced to a eigen problem for the 3 × 3 matrix Γ in Eq. (3.11). Since Γ is a 3 × 3, symmetric, positive-definite matrix, the spectral theorem [23] of linear algebra asserts that Γ always has three linearly independent eigenvectors, regardless whether Γ has degenerate eigenvalues or not. These three eigenvectors of Γ are denoted by **q***v*,1, **q***v*,2, and **q***v*,3. The three associated eigenvectors of *B* with positive eigenvalues are denoted by **q**+,1, **q**+,2, and **q**+,3. Aided by Eq. (3.12), **q**+,1, **q**+,2, and **q**+,3 are readily determined by **q***v*,1, **q***v*,2, and **q***v*,3. Because **q***v*,1, **q***v*,2, and **q***v*,3 are linearly independent, **q**+,1, **q**+,2, and **q**+,3 must also be linearly independent. By following exactly the same procedure, one can derive the three independent eigenvectors **q**−,1, **q**−,2, and **q**−,3 associated with the three negative eigenvalues *λ*<sup>−</sup> = −*λ*+. To proceed, since

> **q**±,*<sup>i</sup>* def =

 **<sup>q</sup>***v*,*<sup>i</sup>* ±**q***T*,*<sup>i</sup>* are the six linearly independent eigenvectors of *B*. The corresponding non-zero eigenvalues of *B* can be written as *λ*±,*i*, *i* = 1, 2, 3, with *λ*−,*<sup>i</sup>* = −*λ*+,*i*, *i* = 1, 2, 3, and *λ*+,1 ≤ *λ*+,2 ≤ *λ*+,3. We note that it is possible for *B* to have degenerate non-zero eigenvalues. Nevertheless, all six eigenvectors associated with the non-zero eigenvalues are linearly independent. This concludes the discussions about the eigenvectors of *B* associated with the non-zero

The remaining three eigenvectors of *B* are associated with the null eigenvalue. Recall that *B* is a 9 × 9 matrix with rank(*B*) ≤ 6. Thus, the nullity of *B* must be greater than or equal to 3, i.e., there are at least three linearly independent vectors spanning the null space of *B* [23]. Aided by the definition of null space and eigenvector, the basis of the null space of *B* consists of the eigenvectors associated with the null eigenvalue of *B*. As such, there must be exactly three eigenvectors associated with the null eigenvalue of *B*, because (i) there are at most nine linearly independent eigenvectors of *B*, (ii) there are at least three linearly independent eigenvectors associated with the null eigenvalues of *B*, and (iii) there are exactly six linearly independent eigenvectors associated with the non-zero eigenvalues of *B*. These three eigenvectors are denoted by **q**0,1, **q**0,2, and **q**0,3. This concludes the discussion of three

Together, there are nine linearly independent eigenvectors of *B*: **q**±,*i*, **q**0,*i*, *i* = 1, 2, 3. The

where *Q* is a 9 × 9 matrix composed of nine column vectors. We then perform similarity transformation Λ = *Q*−1*BQ* to diagonalize *B*, where the eigenvalue matrix Λ is composed of *λ*−,3, *λ*−,2, *λ*−,1, *λ*0, *λ*0, *λ*0, *λ*+,1, *λ*+,2, and *λ*+,3 with *λ*<sup>0</sup> = 0. This concludes the discussions

> ±,*i* , **<sup>q</sup>**(*μ*)

**q**−,3 **q**−,2 **q**−,1 **q**0,1 **q**0,2 **q**0,3 **q**+,1 **q**+,2 **q**+,3

linearly independent eigenvectors associated with the null eigenvalue of *B*.

of diagonalizability of *B* by the above nine linearly independent eigenvectors.

diagonalizing eigenvector matrix of *B* can be constructed as:

*<sup>Q</sup>* def =

As an aside, each of Jacobian matrices *A*(*μ*)

*<sup>ν</sup>* = 1, 2, 3. For these special cases, we let **<sup>q</sup>**(*μ*)

, *i* = 1, 2, 3, (3.13)

, *μ* = 1, 2, 3, is a special case of *B* with *h<sup>ν</sup>* = *δμν*,

0,*<sup>i</sup>* , *<sup>i</sup>* <sup>=</sup> 1, 2, 3, be the eigenvectors of *<sup>A</sup>*(*μ*)

, (3.14)

,

fully constructed.

*λ*− = −*λ*+ �= *λ*+,

eigenvalues.

$$
\frac{\partial \hat{\mathbf{u}}}{\partial t} + \Lambda \frac{\partial \hat{\mathbf{u}}}{\partial y} = 0,\tag{3.17}
$$

where Λ is a diagonal matrix consisting of the eigenvalues of *B*. Equation (3.17) contains nine decoupled equations:

$$\frac{\partial \hat{u}\_i}{\partial t} + \lambda\_i \frac{\partial \hat{u}\_i}{\partial y} = 0, \quad i = 1, \dots, 9,\tag{3.18}$$

where *λ<sup>i</sup>* is the *i*th eigenvalue of *B* and *u*ˆ*<sup>i</sup>* the *i*th characteristic variable. If Eq. (2.3) is hyperbolic, each scalar equation in Eq. (3.18) describes convection of a characteristic variable *u*ˆ*<sup>i</sup>* with a real wave speed *λi*.

## **4. Coordinate transformation and the Christoffel matrix**

In this section, we show how matrix *B* is connected to the classic Christoffel matrix. Essentially, we will show that matrix Γ in Eq. (3.8) is in fact the Christoffel matrix for plane waves propagating in the direction of **h**. Recall that Γ is obtained by solving the eigen problem for non-zero eigenvalues of *B*. To proceed, two Cartesian coordinate systems are considered: (i) The global coordinate system (*x*1-*x*2-*x*3), in which the governing equations Eqs. (2.3) and (3.3) are formulated, and (ii) the local coordinate system (*x*¯1-*x*¯2-*x*¯3), in which the values of the stiffness coefficients are tabulated according to the lattice orientation of the solids, e.g., [2]. Because the values of the stiffness constants depend on the coordinate system employed, the stiffness coefficients reported in the literature need to be transformed into the global coordinate system where the model equations are formulated. One coordinate system can be obtained from the other by coordinate rotation. For completeness, the transformation rules are provided in Appendix 8. Vectors, tensors, and 6-component vectors formulated in the local coordinate system are denoted by an over-bar. The entities formulated in the global coordinate system have no bar. For completeness, Appendix 8 shows the coordinate transformation of the model equations. It turns out that the forms of the governing equations in two different coordinate systems are identical.

To connect the present formulation to the classical Christoffel matrices tabulated in the literature, we proceed to transform matrix Γ formulated in global coordinates to be Γ¯ in local coordinates. We first define the following matrices associated with the Jacobian matrices *A*(*μ*) , *μ* = 1, 2, 3:

$$\Gamma^{(\mu)} \stackrel{\text{def}}{=} K^{(\mu)} \mathcal{K}^{(\mu)^{\dagger}}, \quad \mu = 1, 2, 3,\tag{4.1}$$

which are special cases of Γ with *h<sup>ν</sup>* = *δμν*, *ν* = 1, 2, 3. Γ(1) , Γ(2) , and Γ(3) are tensors in **E**<sup>3</sup> and they can be transformed to local coordinates by using the rotation matrix *R* shown in Eq. (8.1) in Appendix 8:

$$
\bar{\Gamma}^{(\mu)} = \mathcal{R}^{\dagger} \Gamma^{(\mu)} \mathcal{R}.\tag{4.2}
$$

The algebra involved is lengthy but straightforward. Aided by Eq. (4.6) and the orthogonality

*rμ*<sup>1</sup> 000 *rμ*<sup>3</sup> *rμ*<sup>2</sup> 0 *rμ*<sup>2</sup> 0 *rμ*<sup>3</sup> 0 *rμ*<sup>1</sup> 0 0 *rμ*<sup>3</sup> *rμ*<sup>2</sup> *rμ*<sup>1</sup> 0

⎞

Velocity-Stress Equations for Wave Propagation in Anisotropic Elastic Media 209

. (4.8)

⎟⎟⎠ . (4.7)

*M*. (4.9)

*R* = *LCL*¯ *<sup>t</sup>*

⎟⎟⎠ . (4.12)

. (4.10)

⎛

⎜⎜⎝

To proceed, the relation between Γ and Γ¯ is derived in the following. Let *L* be the linear

*M* = *R<sup>t</sup>*

�

*hμK*(*μ*)

**h**¯ = *R<sup>t</sup>*

*h*<sup>1</sup> 000 ¯

By comparing Eqs. (4.10) and (4.12) to the Christoffel matrix tabulated in the literature, e.g., Auld [2], it can be seen that Γ¯ in Eq. (4.10) is indeed the Christoffel matrix. Special cases of Γ¯ in the directions along the axes of the local coordinate system, i.e., *x*¯1-, *x*¯2-, and *x*¯3-axis, are

In this section, we illustrate how the above formulation is applied to model waves in elastic solids of hexagonal symmetry. A solid of hexagonal symmetry is shown in Figure 1. The Cartesian coordinates are chosen such that the *x*3-axis is aligned with the hexagonal cylinder of the medium, and the *x*2-axis is aligned with the crystal lattice. In this coordinate system,

*hμL*(*μ*)

� 3 ∑ *μ*=1

*MCM*¯ *<sup>t</sup>*

*hμK*(*μ*)

� 3 ∑ *μ*=1

*h*3) is the unit vector of the wave propagation direction defined in the local

*h*<sup>3</sup> ¯ *h*2

*h*<sup>3</sup> 0 ¯ *h*1

*h*<sup>3</sup> are the direction cosines in the local coordination system. Aided

⎞

�

*<sup>h</sup>μK*(*μ*)*<sup>t</sup>*

�

**h**, (4.11)

of matrix *R*, Eq. (4.5) becomes

combination of *L*(1)

Γ¯ = *R<sup>t</sup>*

where **h**¯ = (¯

Γ¯(*μ*)

coordinates and ¯

Γ*R* = *R<sup>t</sup>*

*h*1, ¯ *h*2, ¯

*h*1, ¯

, *μ* = 1, 2, 3, respectively.

*L*(*μ*) =

, *L*(2)

3 ∑ *ν*=1

, and *L*(3)

Aided by Eq. (4.5), Eq. (4.8) can be expanded and becomes

3 ∑ *μ*=1

*R* = *R<sup>t</sup>*

*hμR<sup>t</sup> K*(*μ*)

Aided by Eqs. (3.8), (3.4), and (4.9), Γ can be readily transformed to Γ¯:

Next, the unit vector **h** is also transformed into the local coordinates:

by Eqs. (4.7) and (4.11), matrix *L* defined in Eq. (4.8) can be expanded as

⎛

¯

0 ¯ *h*<sup>2</sup> 0 ¯

0 0 ¯ *h*<sup>3</sup> ¯ *h*<sup>2</sup> ¯ *h*<sup>1</sup> 0

⎜⎜⎝

*L* =

� 3 ∑ *μ*=1

*L* =

*GCG<sup>t</sup>*

*h*2, and ¯

**5. Elastic solids with hexagonal symmetry**

:

*rμνK*(*ν*) =

*L* def = 3 ∑ *μ*=1

<sup>Γ</sup>¯(*μ*) has the three eigenvectors ¯**q**(*μ*) *<sup>v</sup>*,*<sup>i</sup>* , *<sup>i</sup>* <sup>=</sup> 1, 2, 3. **<sup>q</sup>**(*μ*) *<sup>v</sup>*,*<sup>i</sup>* and ¯**q**(*μ*) *<sup>v</sup>*,*<sup>i</sup>* are the same vector defined in the global and local coordinate systems, respectively. The algebraic relation between the Γ(*μ*) and Γ¯(*μ*) in Eq. (4.2) needs to be further elaborated. To proceed, Bond's matrix [2] is used to transform the stiffness matrix from the local coordinates to the global coordinates:

$$\mathbf{C} = M\bar{\mathbf{C}}M^{\dagger}.\tag{4.3}$$

Details of the transformation are provided in Appendix 8. For a wide range of media, entries of *C*¯ in the local coordinate system are widely available in the literature, e.g., [2]. Substituting Eqs. (4.1) and (4.3) into Eq. (4.2) gives:

$$\bar{\Gamma}^{(\mu)} = L^{(\mu)} \bar{\mathbb{C}} L^{(\mu)^{\dagger}}, \quad \mu = 1, 2, 3,\tag{4.4}$$

where

$$L^{(\mu)} \stackrel{\text{def}}{=} \mathbb{R}^{\mathfrak{t}} K^{(\mu)} M \tag{4.5}$$

is a 3 × 6 matrix with its entries as the components of the direction cosine corresponding to *xμ*-axis. To proceed, Appendix 8 shows that:

$$K^{(\mu)}M = R\sum\_{\nu=1}^{3} r\_{\mu\nu} K^{(\nu)}.\tag{4.6}$$

The algebra involved is lengthy but straightforward. Aided by Eq. (4.6) and the orthogonality of matrix *R*, Eq. (4.5) becomes

$$L^{(\mu)} = \sum\_{\nu=1}^{3} r\_{\mu\nu} K^{(\nu)} = \begin{pmatrix} r\_{\mu 1} & 0 & 0 & 0 & r\_{\mu 3} \ r\_{\mu 2} \\ 0 & r\_{\mu 2} & 0 & r\_{\mu 3} & 0 & r\_{\mu 1} \\ 0 & 0 & r\_{\mu 3} \ r\_{\mu 2} & r\_{\mu 1} & 0 \end{pmatrix}. \tag{4.7}$$

To proceed, the relation between Γ and Γ¯ is derived in the following. Let *L* be the linear combination of *L*(1) , *L*(2) , and *L*(3) :

$$L \stackrel{\text{def}}{=} \sum\_{\mu=1}^{3} h\_{\mu} L^{(\mu)}.\tag{4.8}$$

Aided by Eq. (4.5), Eq. (4.8) can be expanded and becomes

8 Will-be-set-by-IN-TECH

(3.3) are formulated, and (ii) the local coordinate system (*x*¯1-*x*¯2-*x*¯3), in which the values of the stiffness coefficients are tabulated according to the lattice orientation of the solids, e.g., [2]. Because the values of the stiffness constants depend on the coordinate system employed, the stiffness coefficients reported in the literature need to be transformed into the global coordinate system where the model equations are formulated. One coordinate system can be obtained from the other by coordinate rotation. For completeness, the transformation rules are provided in Appendix 8. Vectors, tensors, and 6-component vectors formulated in the local coordinate system are denoted by an over-bar. The entities formulated in the global coordinate system have no bar. For completeness, Appendix 8 shows the coordinate transformation of the model equations. It turns out that the forms of the governing equations in two different

To connect the present formulation to the classical Christoffel matrices tabulated in the literature, we proceed to transform matrix Γ formulated in global coordinates to be Γ¯ in local coordinates. We first define the following matrices associated with the Jacobian matrices

*CK*(*μ*)*<sup>t</sup>*

they can be transformed to local coordinates by using the rotation matrix *R* shown in Eq. (8.1)

the global and local coordinate systems, respectively. The algebraic relation between the Γ(*μ*) and Γ¯(*μ*) in Eq. (4.2) needs to be further elaborated. To proceed, Bond's matrix [2] is used to

*C* = *MCM*¯ *<sup>t</sup>*

Details of the transformation are provided in Appendix 8. For a wide range of media, entries of *C*¯ in the local coordinate system are widely available in the literature, e.g., [2]. Substituting

*CL*¯ (*μ*)*<sup>t</sup>*

is a 3 × 6 matrix with its entries as the components of the direction cosine corresponding to

3 ∑ *ν*=1

*rμνK*(*ν*)

*K*(*μ*)

*L*(*μ*) def = *R<sup>t</sup>*

*K*(*μ*)*M* = *R*

Γ(*μ*)

*<sup>v</sup>*,*<sup>i</sup>* and ¯**q**(*μ*)

Γ¯(*μ*) = *R<sup>t</sup>*

*<sup>v</sup>*,*<sup>i</sup>* , *<sup>i</sup>* <sup>=</sup> 1, 2, 3. **<sup>q</sup>**(*μ*)

transform the stiffness matrix from the local coordinates to the global coordinates:

Γ¯(*μ*) = *L*(*μ*)

, *μ* = 1, 2, 3, (4.1)

*R*. (4.2)

. (4.3)

, *μ* = 1, 2, 3, (4.4)

*M* (4.5)

. (4.6)

*<sup>v</sup>*,*<sup>i</sup>* are the same vector defined in

, and Γ(3) are tensors in **E**<sup>3</sup> and

, Γ(2)

Γ(*μ*) def

which are special cases of Γ with *h<sup>ν</sup>* = *δμν*, *ν* = 1, 2, 3. Γ(1)

= *K*(*μ*)

coordinate systems are identical.

<sup>Γ</sup>¯(*μ*) has the three eigenvectors ¯**q**(*μ*)

Eqs. (4.1) and (4.3) into Eq. (4.2) gives:

*xμ*-axis. To proceed, Appendix 8 shows that:

*A*(*μ*)

, *μ* = 1, 2, 3:

in Appendix 8:

where

$$L = \sum\_{\mu=1}^{3} h\_{\mu} \mathbb{R}^{t} K^{(\mu)} M = \mathbb{R}^{t} \left( \sum\_{\mu=1}^{3} h\_{\mu} K^{(\mu)} \right) M. \tag{4.9}$$

Aided by Eqs. (3.8), (3.4), and (4.9), Γ can be readily transformed to Γ¯:

$$\bar{\Gamma} = R^\dagger \Gamma R = R^\dagger G C G^\dagger R = R^\dagger \left( \sum\_{\mu=1}^3 h\_{\mu} K^{\{\mu\}} \right) M \bar{C} M^\dagger \left( \sum\_{\mu=1}^3 h\_{\mu} K^{\{\mu\}^\dagger} \right) R = L \bar{C} L^\dagger. \tag{4.10}$$

Next, the unit vector **h** is also transformed into the local coordinates:

$$
\tilde{\mathbf{h}} = \mathcal{R}^t \mathbf{h},
\tag{4.11}
$$

where **h**¯ = (¯ *h*1, ¯ *h*2, ¯ *h*3) is the unit vector of the wave propagation direction defined in the local coordinates and ¯ *h*1, ¯ *h*2, and ¯ *h*<sup>3</sup> are the direction cosines in the local coordination system. Aided by Eqs. (4.7) and (4.11), matrix *L* defined in Eq. (4.8) can be expanded as

$$L = \begin{pmatrix} \bar{h}\_1 & 0 & 0 & 0 & \bar{h}\_3 \ \bar{h}\_2 \\ 0 & \bar{h}\_2 & 0 & \bar{h}\_3 & 0 & \bar{h}\_1 \\ 0 & 0 & \bar{h}\_3 \ \bar{h}\_2 \ \bar{h}\_1 & 0 \end{pmatrix}. \tag{4.12}$$

By comparing Eqs. (4.10) and (4.12) to the Christoffel matrix tabulated in the literature, e.g., Auld [2], it can be seen that Γ¯ in Eq. (4.10) is indeed the Christoffel matrix. Special cases of Γ¯ in the directions along the axes of the local coordinate system, i.e., *x*¯1-, *x*¯2-, and *x*¯3-axis, are Γ¯(*μ*) , *μ* = 1, 2, 3, respectively.

#### **5. Elastic solids with hexagonal symmetry**

In this section, we illustrate how the above formulation is applied to model waves in elastic solids of hexagonal symmetry. A solid of hexagonal symmetry is shown in Figure 1. The Cartesian coordinates are chosen such that the *x*3-axis is aligned with the hexagonal cylinder of the medium, and the *x*2-axis is aligned with the crystal lattice. In this coordinate system,

**Figure 1.** The Cartesian coordinates and the lattice structure of solids of hexagonal symmetry.

9 of the 21 stiffness constants are non-zero and 5 of them, i.e., *c*11, *c*12, *c*13, *c*<sup>33</sup> and *c*44, are independent. The stiffness matrix is

$$\mathbf{C} = \begin{pmatrix} c\_{11} \ c\_{12} \ c\_{13} & 0 & 0 & 0 \\ c\_{12} \ c\_{11} \ c\_{13} & 0 & 0 & 0 \\ c\_{13} \ c\_{13} \ c\_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & c\_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & c\_{44} & 0 \\ 0 & 0 & 0 & 0 & 0 & c\_{66} \end{pmatrix} . \tag{5.1}$$

can be directly taken from the literature. To proceed, we transform the model equations from the local coordinates (*x*¯1, *x*¯2, *x*¯3) to the global coordinates (*x*1, *x*2, *x*3). In general, the rotation can be defined by the three Euler angles [1]: *α*, *β*, and *γ*. The coordinate transformation matrix

> cos *γ* cos *β* cos *α* − sin *γ* sin *α* cos *γ* cos *β* sin *α* + sin *γ* cos *α* − cos *γ* sin *β* − sin *γ* cos *β* cos *α* − cos *γ* sin *α* − sin *γ* cos *β* sin *α* + cos *γ* cos *α* sin *γ* sin *β*

With *rμν* denoting the (*μ*, *ν*) entry of *R*, we have *∂x*¯*μ*/*∂x<sup>ν</sup>* = *rμν* with *μ*, *ν* = 1, 2, 3. Aided by

solids of hexagonal symmetry, Eq. (5.4) is reduced to be one-dimensional for modeling wave

As shown in Fig. 2, two angles are needed to specify a particular direction in the three-dimensional space. The Euler angles are chosen to be (*θ*, −*φ*, 0), where 0 ≤ *θ* ≤ 2*π* and −*π*/2 ≤ *φ* ≤ *π*/2. The transformation matrix *R* of a given pair of angles (*θ*, −*φ*) is

Aided by the above coordinate transformation, the Jacobian matrix *A*¯(1) in Eq. (5.5) becomes

� 03×<sup>3</sup> *A*¯ *<sup>v</sup> A*¯ *<sup>T</sup>* 06×<sup>6</sup>

cos *φ* cos *θ* cos *φ* sin *θ* sin *φ* − sin *θ* cos *θ* 0 − sin *φ* cos *θ* − sin *φ* sin *θ* cos *φ*

�

⎞

⎟⎟⎟⎟⎟⎟⎠ ,

− cos *φ* cos *θ* 000 − sin *φ* − cos *φ* sin *θ* 0 − cos *φ* sin *θ* 0 − sin *φ* 0 − cos *φ* cos *θ* 0 0 − sin *φ* − cos *φ* sin *θ* − cos *φ* cos *θ* 0

9 ∑ *j*=1 *a*¯ (1) *ij ∂uj ∂x*¯1

3 ∑ *μ*=1

9 ∑ *j*=1 *a*¯ (*μ*) *ij*

*∂uj ∂x*¯*μ*

Velocity-Stress Equations for Wave Propagation in Anisotropic Elastic Media 211

*ij* . To connect the present formulation to the Christoffel matrices for

*<sup>∂</sup><sup>t</sup>* <sup>+</sup>

*<sup>∂</sup><sup>t</sup>* <sup>+</sup>

sin *β* cos *α* sin *β* sin *α* cos *β*

⎞ ⎠ .

= 0, *i* = 1, 2, . . . , 9, (5.4)

= 0, *i* = 1, 2, . . . , 9, (5.5)

⎠ . (5.6)

⎞ ⎠ ,

⎞

, (5.7)

*R*, Eq. (8.1), can be readily specified:

the chain rule, the model equation Eq. (2.3) becomes

<sup>=</sup> 0, or *<sup>∂</sup>ui*

<sup>=</sup> 0, or *<sup>∂</sup>ui*

⎛ ⎝

−*c*<sup>11</sup> cos *φ* cos *θ* −*c*<sup>12</sup> cos *φ* sin *θ* −*c*<sup>13</sup> sin *φ* −*c*<sup>12</sup> cos *φ* cos *θ* −*c*<sup>11</sup> cos *φ* sin *θ* −*c*<sup>13</sup> sin *φ* −*c*<sup>13</sup> cos *φ* cos *θ* −*c*<sup>13</sup> cos *φ* sin *θ* −*c*<sup>33</sup> sin *φ*

−*c*<sup>66</sup> cos *φ* sin *θ* −*c*<sup>66</sup> cos *φ* cos *θ* 0

0 −*c*<sup>44</sup> sin *φ* −*c*<sup>44</sup> cos *φ* sin *θ* −*c*<sup>44</sup> sin *φ* 0 −*c*<sup>44</sup> cos *φ* cos *θ*

*<sup>A</sup>*¯(*μ*) *<sup>∂</sup>***<sup>u</sup>** *∂x*¯*μ*

(*ν*)

*∂x*¯1

*R*(*θ*, −*φ*, 0) =

*R*(*α*, *β*, *γ*) = ⎛ ⎝

> *∂***u** *<sup>∂</sup><sup>t</sup>* <sup>+</sup>

(*μ*) *ij* <sup>=</sup> <sup>∑</sup><sup>3</sup>

*<sup>A</sup>*¯ *<sup>v</sup>* <sup>=</sup> <sup>1</sup> *ρ*

*A*¯ *<sup>T</sup>* =

⎛ ⎝

⎛

⎜⎜⎜⎜⎜⎜⎝

where *a*¯

3 ∑ *μ*=1

*<sup>ν</sup>*=<sup>1</sup> *rμνa*

*<sup>∂</sup><sup>t</sup>* <sup>+</sup> *<sup>A</sup>*¯(1) *<sup>∂</sup>***<sup>u</sup>**

*<sup>A</sup>*¯(1) = *<sup>r</sup>*11*A*(1) + *<sup>r</sup>*12*A*(2) + *<sup>r</sup>*13*A*(3) =

propagation along the *x*¯1 axis:

*∂***u**

where *c*<sup>66</sup> = (*c*<sup>11</sup> − *c*12)/2. Aided by Eq. (5.1), the Jacobian matrix *A*(1) , i.e., Eq. (2.4) in the model equation Eq. (2.3) become

$$A^{(1)} = \begin{pmatrix} \mathbf{0}\_{3 \times 3} & -\frac{1}{\rho} K^{(1)} \\ -\mathcal{K}^{(1)^\dagger} & \mathbf{0}\_{6 \times 6} \end{pmatrix}, \quad \text{where} \quad -\mathcal{K}^{(1)^\dagger} = \begin{pmatrix} -c\_{11} & 0 & 0 \\ -c\_{12} & 0 & 0 \\ -c\_{13} & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & -c\_{44} \\ 0 & -c\_{66} & 0 \end{pmatrix}, \tag{5.2}$$

Similarly, for *A*(2) and *A*(3) , the related matrices are

$$-\mathcal{L}K^{(2)}{}^{t} = \begin{pmatrix} 0 & -c\_{12} & 0 \\ 0 & -c\_{11} & 0 \\ 0 & -c\_{13} & 0 \\ 0 & 0 & -c\_{44} \\ 0 & 0 & 0 \\ -c\_{66} & 0 & 0 \end{pmatrix} \quad \text{and} \quad -\mathcal{L}K^{(3)}{}^{t} = \begin{pmatrix} 0 & 0 & -c\_{13} \\ 0 & 0 & -c\_{13} \\ 0 & 0 & -c\_{33} \\ 0 & -c\_{44} & 0 \\ -c\_{44} & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}. \tag{5.3}$$

The above Jacobian matrices in the model equations are valid only if the local coordinates (*x*¯1, *x*¯2, *x*¯3) are employed, in which the entries of the stiffness matrix such as that shown in Eq. (5.1) can be directly taken from the literature. To proceed, we transform the model equations from the local coordinates (*x*¯1, *x*¯2, *x*¯3) to the global coordinates (*x*1, *x*2, *x*3). In general, the rotation can be defined by the three Euler angles [1]: *α*, *β*, and *γ*. The coordinate transformation matrix *R*, Eq. (8.1), can be readily specified:

10 Will-be-set-by-IN-TECH

x3

**Figure 1.** The Cartesian coordinates and the lattice structure of solids of hexagonal symmetry.

*C* =

*ρK*(1)

⎛

⎜⎜⎜⎜⎜⎜⎝

where *c*<sup>66</sup> = (*c*<sup>11</sup> − *c*12)/2. Aided by Eq. (5.1), the Jacobian matrix *A*(1)

�

, the related matrices are

⎞

⎟⎟⎟⎟⎟⎟⎠

9 of the 21 stiffness constants are non-zero and 5 of them, i.e., *c*11, *c*12, *c*13, *c*<sup>33</sup> and *c*44, are

*c*<sup>11</sup> *c*<sup>12</sup> *c*<sup>13</sup> 000 *c*<sup>12</sup> *c*<sup>11</sup> *c*<sup>13</sup> 000 *c*<sup>13</sup> *c*<sup>13</sup> *c*<sup>33</sup> 000 000 *c*<sup>44</sup> 0 0 0000 *c*<sup>44</sup> 0 00000 *c*<sup>66</sup>

, where <sup>−</sup> *CK*(1)*<sup>t</sup>*

and <sup>−</sup> *CK*(3)*<sup>t</sup>*

The above Jacobian matrices in the model equations are valid only if the local coordinates (*x*¯1, *x*¯2, *x*¯3) are employed, in which the entries of the stiffness matrix such as that shown in Eq. (5.1)

⎞

⎟⎟⎟⎟⎟⎟⎠

=

=

⎛

⎜⎜⎜⎜⎜⎜⎝

⎛

⎜⎜⎜⎜⎜⎜⎝

x2

. (5.1)

, i.e., Eq. (2.4) in the

⎞

⎟⎟⎟⎟⎟⎟⎠

⎞

⎟⎟⎟⎟⎟⎟⎠

, (5.2)

. (5.3)

x1

independent. The stiffness matrix is

model equation Eq. (2.3) become

� 03×<sup>3</sup> <sup>−</sup> <sup>1</sup>

<sup>−</sup>*CK*(1)*<sup>t</sup>* <sup>06</sup>×<sup>6</sup>

*A*(1) =

Similarly, for *A*(2) and *A*(3)

<sup>−</sup>*CK*(2)*<sup>t</sup>*

=

⎛

⎜⎜⎜⎜⎜⎜⎝

$$\begin{pmatrix} \mathcal{R}(\mathfrak{a}, \mathfrak{f}, \gamma) \\\\ \begin{pmatrix} \cos \gamma \cos \beta \cos \mathfrak{a} - \sin \gamma \sin \mathfrak{a} & \cos \gamma \cos \beta \sin \mathfrak{a} + \sin \gamma \cos \mathfrak{a} & -\cos \gamma \sin \beta \\\ -\sin \gamma \cos \beta \cos \mathfrak{a} - \cos \gamma \sin \mathfrak{a} & -\sin \gamma \cos \beta \sin \mathfrak{a} + \cos \gamma \cos \mathfrak{a} & \sin \gamma \sin \beta \\\ \sin \beta \cos \mathfrak{a} & \sin \beta \sin \mathfrak{a} & \cos \beta \end{pmatrix}. \\ \end{pmatrix}$$

With *rμν* denoting the (*μ*, *ν*) entry of *R*, we have *∂x*¯*μ*/*∂x<sup>ν</sup>* = *rμν* with *μ*, *ν* = 1, 2, 3. Aided by the chain rule, the model equation Eq. (2.3) becomes

$$\frac{\partial \mathbf{u}}{\partial t} + \sum\_{\mu=1}^{3} \bar{A}^{(\mu)} \frac{\partial \mathbf{u}}{\partial \bar{\mathbf{x}}\_{\mu}} = 0, \quad \text{or} \quad \frac{\partial u\_i}{\partial t} + \sum\_{\mu=1}^{3} \sum\_{j=1}^{9} \bar{a}\_{ij}^{(\mu)} \frac{\partial u\_j}{\partial \bar{\mathbf{x}}\_{\mu}} = 0, \quad i = 1, 2, \dots, 9,\tag{5.4}$$

where *a*¯ (*μ*) *ij* <sup>=</sup> <sup>∑</sup><sup>3</sup> *<sup>ν</sup>*=<sup>1</sup> *rμνa* (*ν*) *ij* . To connect the present formulation to the Christoffel matrices for solids of hexagonal symmetry, Eq. (5.4) is reduced to be one-dimensional for modeling wave propagation along the *x*¯1 axis:

$$\frac{\partial \mathbf{u}}{\partial t} + \bar{A}^{(1)} \frac{\partial \mathbf{u}}{\partial \tilde{x}\_1} = 0, \quad \text{or} \quad \frac{\partial u\_i}{\partial t} + \sum\_{j=1}^9 \bar{a}\_{ij}^{(1)} \frac{\partial u\_j}{\partial \tilde{x}\_1} = 0, \quad i = 1, 2, \dots, 9,\tag{5.5}$$

As shown in Fig. 2, two angles are needed to specify a particular direction in the three-dimensional space. The Euler angles are chosen to be (*θ*, −*φ*, 0), where 0 ≤ *θ* ≤ 2*π* and −*π*/2 ≤ *φ* ≤ *π*/2. The transformation matrix *R* of a given pair of angles (*θ*, −*φ*) is

$$R(\theta, -\phi, 0) = \begin{pmatrix} \cos\phi\cos\theta & \cos\phi\sin\theta & \sin\phi \\ -\sin\theta & \cos\theta & 0 \\ -\sin\phi\cos\theta & -\sin\phi\sin\theta \cos\phi \end{pmatrix} . \tag{5.6}$$

Aided by the above coordinate transformation, the Jacobian matrix *A*¯(1) in Eq. (5.5) becomes

$$\begin{aligned} \bar{A}^{(1)} &= r\_{11}A^{(1)} + r\_{12}A^{(2)} + r\_{13}A^{(3)} = \left(\frac{\mathbf{0}\_{3\times 2}}{\bar{A}\_{T}}\mathbf{\bar{0}}\_{6\times 6}^{2}\right), \\ \bar{A}\_{\upsilon} &= \frac{1}{\rho} \begin{pmatrix} -\cos\phi\cos\theta & 0 & 0 & 0 & -\sin\phi & -\cos\phi\sin\theta \\ 0 & -\cos\phi\sin\theta & 0 & -\sin\phi & 0 & -\cos\phi\cos\theta \\ 0 & 0 & -\sin\phi - \cos\phi\sin\theta & -\cos\phi\cos\theta & 0 \end{pmatrix}, \\ \bar{A}\_{T} &= \begin{pmatrix} -c\_{11}\cos\phi\cos\theta & -c\_{12}\cos\phi\sin\theta & -c\_{13}\sin\phi \\ -c\_{12}\cos\phi\cos\theta & -c\_{11}\cos\phi\sin\theta & -c\_{13}\sin\phi \\ 0 & -c\_{44}\cos\phi\cos\theta & -c\_{44}\cos\phi\sin\theta \\ 0 & -c\_{44}\sin\phi & -c\_{44}\cos\phi\sin\theta \\ -c\_{46}\cos\phi\sin\theta & 0 & -c\_{44}\cos\phi\cos\theta \\ -c\_{66}\cos\phi\sin\theta & -c\_{66}\cos\phi\cos\theta & 0 \end{pmatrix}, \end{aligned}$$

**Figure 2.** Rotations of the Cartesian coordinate system. (a) Rotation with respect to the *x*<sup>3</sup> axis. (b) Rotation with respect to the *x*¯2 axis. (c) Overall rotation.

Matrix *A*¯(1) is a 9 × 9 matrix, there are nine eigenvalues, which are calculated by solving the polynomial:

$$\det(A^{\tilde{1}(1)} - \lambda I\_{\theta}) = \det\begin{pmatrix} -\lambda I\_3 & \bar{A}\_v \\ \bar{A}\_T & -\lambda I\_6 \end{pmatrix} = 0\_\lambda$$

where *In* is a rank *n* identity matrix. Based on the Schur complement [15, pp. 475], the non-trivial eigenvalues of *A*¯(1) can be obtained by solving a simpler equation

$$\det(\bar{A}\_{\upsilon}\bar{A}\_{T} - \kappa^{2}I\_{3}) = 0,\tag{5.8}$$

Equation (5.10) are identical to the Christoffel equations for hexagonal solids as shown in the

Velocity-Stress Equations for Wave Propagation in Anisotropic Elastic Media 213

<sup>3</sup> are positive and real. As such, the non-trivial eigenvalues of *<sup>A</sup>*¯(1) are <sup>±</sup>*κ<sup>i</sup>* with *<sup>i</sup>* <sup>=</sup> 1, 2, and 3. Moreover, *κ*1, *κ*2, and *κ*<sup>3</sup> are the three wave speeds in the hexagonal solids. Simple manipulation also shows that the other three trivial eigenvalues of *A*¯(1) are null. Therefore,

In this section, we apply the space-time Conservation Element and Solution Element (CESE) method [4] to solve the velocity-stress equations for modeling wave propagation in a block of beryl, a elastic solid of hexagonal symmetry. The CESE method is a modern numerical method for time accurate solutions of linear and non-linear hyperbolic partial differential equations. The CESE method treats space and time in a unified way in calculating the space-time flux. The method is simple, robust, and its operational count is comparable to that of a second-order central difference scheme. The backbone of the CESE method is the *a* scheme [4], which is neutrally stable without artificial dissipation for solving a scalar convection equation with a constant coefficient. For more complex wave equations, the *a* scheme is extended with added artificial damping, including the *a*-*α*-*�* scheme [4] for shock capturing and the *c*-*τ* scheme [5] for Courant-Friedrichs-Lewy (CFL) number insensitive calculations. To model multi-dimensional problems, the CESE method uses unstructured meshes [27] with triangles and tetrahedra as the basic elements in two- and three-dimensional spaces. Extensions to use quadrilaterals and hexahedra are also available [36]. In the past, the CESE method has been successfully applied to solve a wide range of various fluid dynamics problems, including aero-acoustics [14], cavitations [19], complex shock waves [8], detonations [26], magnetohydrodynamics (MHD)[34, 35], etc. Moreover, the method has been applied to model electromagnetic waves [28] and nonlinear stress waves in isotropic solids [3, 31, 32]. In the

present chapter, the following results are obtained by using the *a*-*α*-*�* scheme.

The two-dimensional velocity-stress equations are employed to model wave propagation in a block of beryl, which is a solid of hexagonal symmetry. The material properties are taken from [16]: *ρ* = 2.7 g/cm3, *c*<sup>11</sup> = 26.94 × 1011 dynes/cm2, *c*<sup>12</sup> = 9.61 × 1011 dynes/cm2, *c*<sup>13</sup> = 6.61 × 1011 dynes/cm2, *c*<sup>33</sup> = 23.63 × 1011 dynes/cm2, and *c*<sup>44</sup> = 6.53 × 1011 dynes/cm2. The computational domain is −1 m ≤ *x*(1) ≤ 1 m and −1 m ≤ *x*(2) ≤ 1 m. A unstructured mesh composed of triangular elements is used. The square domain is divided into 2.2 millions of triangular elements. The maximum length of the element edge is 2.68 × 10−<sup>3</sup> m. The initial condition is a two-dimensional step function **v** = **a**H(*σ* − |**x**|), where **a** represents a constant initial velocity to initiate the expanding wave from the origin. The constant *σ* in the initial

The two-dimensional code is parallelized by domain decomposition. First, a graph of element connectivity is built based on the unstructured mesh. The connectivity graph is processed by METIS [7] to partition the domain. As shown in Fig. 3, the overall spatial domain is decomposed into 16 sub-domains for parallel computing. Numerical calculations for elements in each sub-domain is distributed to a workstation as a part of a networked cluster. In all cases, Δ*t* = 65 × 10−<sup>9</sup> s and the maximum CFL number *ν* is about 0.95. Three cases of wave propagation are calculated: (1) calculation of the group velocity profile to assess

<sup>1</sup>, *<sup>κ</sup>*<sup>2</sup> <sup>2</sup>, and

Shown in Eq. (5.8), matrix *A*¯ *vA*¯ *<sup>T</sup>* is real and symmetric. Thus, its three eigenvalues *κ*<sup>2</sup>

**6. Modeling wave propagation in a block of beryl**

literature, e.g., Auld [2].

all nine eigenvalues are real.

condition is chosen to be 5 × 10−<sup>3</sup> m.

*κ*2

where

$$
\bar{A}\_{\upsilon}\bar{A}\_{T} = \frac{1}{\rho} \begin{pmatrix}
\Gamma\_{11} \; \Gamma\_{12} \; \Gamma\_{13} \\
\Gamma\_{21} \; \Gamma\_{22} \; \Gamma\_{23} \\
\Gamma\_{31} \; \Gamma\_{32} \; \Gamma\_{33}
\end{pmatrix} /
$$

and

$$\begin{aligned} \Gamma\_{11} &= c\_{11}\cos^2\phi\cos^2\theta + c\_{44}\sin^2\phi + c\_{66}\cos^2\phi\sin^2\theta, \\ \Gamma\_{22} &= c\_{11}\cos^2\phi\sin^2\theta + c\_{44}\sin^2\phi + c\_{66}\cos^2\phi\cos^2\theta, \\ \Gamma\_{33} &= c\_{33}\sin^2\phi + c\_{44}\cos^2\phi\sin^2\theta + c\_{44}\cos^2\phi\cos^2\theta. \\ \Gamma\_{12} &=\Gamma\_{21} = c\_{12}\cos^2\phi\cos\theta\sin\theta + c\_{66}\cos^2\phi\cos\theta\sin\theta, \\ \Gamma\_{13} &=\Gamma\_{31} = c\_{13}\cos\phi\sin\phi\cos\theta + c\_{44}\cos\phi\sin\phi\cos\theta, \\ \Gamma\_{23} &=\Gamma\_{32} = c\_{13}\cos\phi\sin\phi\sin\theta + c\_{44}\cos\phi\sin\phi\sin\theta. \end{aligned} \tag{5.9}$$

To proceed, according to Eq. (5.6), we let **h**¯ = (¯ *h*1, ¯ *h*2, ¯ *h*3) be the unit vector parallel to the *x*¯1 axis in the local coordinate system, i.e., ¯ *h*<sup>1</sup> = *r*<sup>11</sup> = cos *φ* cos *θ*, ¯ *<sup>h</sup>*<sup>2</sup> <sup>=</sup> *<sup>r</sup>*<sup>12</sup> <sup>=</sup> cos *<sup>φ</sup>* sin *<sup>θ</sup>*, and ¯ *h*<sup>3</sup> = *r*<sup>13</sup> = sin *φ*. Aided by the definition of **h**¯, Eq. (5.9) can be rewritten as

$$\begin{aligned} \Gamma\_{11} &= c\_{11}\bar{h}\_1^2 + c\_{44}\bar{h}\_3^2 + c\_{66}\bar{h}\_2^2, \\ \Gamma\_{22} &= c\_{11}\bar{h}\_2^2 + c\_{44}\bar{h}\_3^2 + c\_{66}\bar{h}\_1^2, \\ \Gamma\_{33} &= c\_{44}\bar{h}\_1^2 + c\_{44}\bar{h}\_2^2 + c\_{33}\bar{h}\_3^2, \\ \Gamma\_{12} &= \Gamma\_{21} = (c\_{12} + c\_{66})\bar{h}\_1\bar{h}\_2, \\ \Gamma\_{13} &= \Gamma\_{31} = (c\_{13} + c\_{44})\bar{h}\_1\bar{h}\_3, \\ \Gamma\_{23} &= \Gamma\_{32} = (c\_{13} + c\_{44})\bar{h}\_2\bar{h}\_3. \end{aligned} \tag{5.10}$$

Equation (5.10) are identical to the Christoffel equations for hexagonal solids as shown in the literature, e.g., Auld [2].

Shown in Eq. (5.8), matrix *A*¯ *vA*¯ *<sup>T</sup>* is real and symmetric. Thus, its three eigenvalues *κ*<sup>2</sup> <sup>1</sup>, *<sup>κ</sup>*<sup>2</sup> <sup>2</sup>, and *κ*2 <sup>3</sup> are positive and real. As such, the non-trivial eigenvalues of *<sup>A</sup>*¯(1) are <sup>±</sup>*κ<sup>i</sup>* with *<sup>i</sup>* <sup>=</sup> 1, 2, and 3. Moreover, *κ*1, *κ*2, and *κ*<sup>3</sup> are the three wave speeds in the hexagonal solids. Simple manipulation also shows that the other three trivial eigenvalues of *A*¯(1) are null. Therefore, all nine eigenvalues are real.
