**Author details**

16 Will-be-set-by-IN-TECH

**Figure 6.** The computational domain is divided into three regions. The orientations of the solids in the central, left and right regions are *θ* = 0◦ and *φ* = 0◦, *θ* = 0◦ and *φ* = 60◦, and *θ* = 0◦ and *φ* = 30◦,

elastic solids. The eigen structure of the equations has been thoroughly studied by analyzing

hyperbolicity of the equations has been rigorously proved by showing the real spectrum of matrix *B* and diagonalizability of *B*. *B* has a degenerate null eigenvalue. For non-zero eigenvalues, *B* has at least one and at most three positive-negative pairs of real eigenvalues. Moreover, *B* has nine linearly independent eigenvectors, which can be used to diagonalize *B*. The eigenvalues represent the wave speeds, and the eigenvectors are related to the polarization of the waves. Since the first-order velocity-stress equations are hyperbolic, the Cauchy problem is well-posed [9]. For numerical simulation, the proven hyperbolicity justifies the use of modern finite-volume methods, originally developed for solving conservation laws, to solve the velocity-stress equations. For upwind methods [12] and discontinuous Galerkin methods [18], information about eigenvalues and eigenvectors of matrix *B* are critically important. The derived characteristic relations of the governing equations can be directly used to derive the Riemann solver and the flux function as the building blocks of the methods. Moreover, by clearly defining the local and global coordinate systems, we also show that matrix *B* is directly connected to the classic Christoffel matrix. Therefore, nine coupled first-order velocity-stress equations are shown to be directly related to the classic second-order elastodynamic equations. To demonstrate the application of the velocity-stress equations, we apply the space-time CESE method to solve the equations for modeling wave propagation in a block of beryl. The calculated energy profiles compare well with the analytical solution. We also simulate wave propagation in a heterogeneous solid composed of three blocks of beryl with different lattice orientations. Numerical results show

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

. The

the composed Jacobian matrix of the first-order equations, i.e., *B* = ∑<sup>3</sup>

respectively.

complex wave/interface interactions.

Sheng-Tao John Yu, Yung-Yu Chen and Lixiang Yang *The Department of Mechanical Engineering, The Ohio State University, Columbus, OH 43210, USA*

**Figure 7.** Waves propagation in a heterogeneous domain. The total energy density (*e*/ max(*e*)) are plotted at different times: (a) *t* = 72 *μ*s, (b) *t* = 96 *μ*s, (c) *t* = 108 *μ*s, (d) *t* = 180 *μ*s. Δ*t* = 60 ns, and CFL number = 0.95.

## **8. Appendix**

#### **A: The Schur complements**

Consider matrix *M*, which can be divided into 4 sub matrices:

$$M = \begin{pmatrix} A & B \\ C & D \end{pmatrix}'$$

where *A*, *B*, *C*, and *D* are *p* × *p*, *p* × *q*, *q* × *p*, and *q* × *q* matrices, respectively. If *A*−<sup>1</sup> and *D*−<sup>1</sup> exist, then matrices *D* − *CA*−1*B* and *A* − *BD*−1*C* are the Schur complements of *A* and *D* [15]. In this case, matrix *M* can be factored out as

$$
\begin{pmatrix} A \ B \\ C \ D \end{pmatrix} = \begin{pmatrix} I\_p & 0\_{p \times q} \\ CA^{-1} & I\_q \end{pmatrix} \begin{pmatrix} A & B \\ 0\_{q \times p} \ D - CA^{-1}B \end{pmatrix} / 2
$$

or

$$
\begin{pmatrix} A \ B \\ C \ D \end{pmatrix} = \begin{pmatrix} A - BD^{-1} \mathcal{C} \ B \\ 0\_{q \times p} & D \end{pmatrix} \begin{pmatrix} I\_p & 0\_{p \times q} \\ D^{-1} \mathcal{C} & I\_q \end{pmatrix} \dots$$

This factorization is useful to calculate det(*M*).

#### **B: Rotation of Cartesian coordinate system**

This appendix provides the transformation relations for the global and local coordinate systems. Let **e**1, **e**2, and **e**<sup>3</sup> be the natural basis of the global coordinate system, and let **e**¯1, ¯**e**2, and ¯**e**<sup>3</sup> be the natural basis of the local coordinate system. An arbitrary vector **b** can be represented by using either the global coordinates or the local coordinates as: **b** = ∑3 *<sup>μ</sup>*=<sup>1</sup> *<sup>b</sup>μ***e***<sup>μ</sup>* <sup>=</sup> <sup>∑</sup><sup>3</sup> *<sup>μ</sup>*=<sup>1</sup> ¯ *bμ***e**¯*μ*. Let **b** denote the vector defined in the global coordinates and **b**¯ defined in the local coordinates. The transformation between **b** and **b**¯ is written as: **b** = *R***b**¯, where

$$R = \begin{pmatrix} r\_{11} \ r\_{12} \ r\_{13} \\ r\_{21} \ r\_{22} \ r\_{23} \\ r\_{31} \ r\_{32} \ r\_{33} \end{pmatrix} \tag{8.1}$$

To transform the 6-component stress vector **T**, we use the Bond matrix [2]. Aided by matrix

*r*21*r*<sup>31</sup> *r*22*r*<sup>32</sup> *r*23*r*<sup>33</sup> *r*22*r*<sup>33</sup> + *r*32*r*<sup>23</sup> *r*21*r*<sup>33</sup> + *r*23*r*<sup>31</sup> *r*21*r*<sup>32</sup> + *r*22*r*<sup>31</sup> *r*11*r*<sup>31</sup> *r*12*r*<sup>32</sup> *r*13*r*<sup>33</sup> *r*12*r*<sup>33</sup> + *r*13*r*<sup>32</sup> *r*11*r*<sup>33</sup> + *r*13*r*<sup>31</sup> *r*11*r*<sup>32</sup> + *r*12*r*<sup>31</sup> *r*11*r*<sup>21</sup> *r*12*r*<sup>22</sup> *r*13*r*<sup>23</sup> *r*12*r*<sup>23</sup> + *r*13*r*<sup>22</sup> *r*11*r*<sup>23</sup> + *r*13*r*<sup>21</sup> *r*11*r*<sup>22</sup> + *r*12*r*<sup>21</sup>

The application of *M* to the stress vector is **T** = *M***T**¯ , where **T** is the 6-component stress vector defined in the global coordinates and **T**¯ in local coordinates. Similarly, the Bond matrix for strain *N* [2] is used to transform the 6-component strain vector **S** as **S** = *N***S**¯. Simple

equations in the global and local coordinate systems are stated as **T** = *C***S** and **T**¯ = *C*¯**S**¯. Aided

This appendix shows that the velocity-stress equations, Eq. (1.4), remain in the same form when formulated in different Cartesian coordinate systems. To proceed, we consider an arbitrary tensor *E*. In the following, we show that ∇ · *E* = ∇ · ¯ (*ER*), where ∇· ¯ is the divergence operator defined in the local coordinate system, and *R* is the rotation matrix *R* defined in

> 3 ∑ *k*=1

*∂eij ∂x*¯*<sup>k</sup>*

*rjk* =

*<sup>i</sup>* . (8.5)

*TR*) = <sup>1</sup> *ρ* ∇ · ¯ *T*¯.

3 ∑ *k*=1

*∂ ∂x*¯*<sup>k</sup>*

∇ · *T*. (8.6)

*SR*. (8.7)

⎛ ⎝ 3 ∑ *j*=1 *eijrjk*

⎞ ⎠

by the Bond's matrices *M* and *N*, *C*¯ can be related to *C* in the following equation:

**T** = *M***T**¯ = *M*(*C*¯**S**¯) = *MC*¯(*N*−1**S**)=(*MCM*¯ *<sup>t</sup>*

<sup>13</sup> 2*r*12*r*<sup>13</sup> 2*r*11*r*<sup>13</sup> 2*r*11*r*<sup>12</sup>

⎞

⎟⎟⎟⎟⎟⎟⎠

. (8.4)

<sup>23</sup> 2*r*22*r*<sup>23</sup> 2*r*21*r*<sup>23</sup> 2*r*21*r*<sup>22</sup>

<sup>33</sup> 2*r*32*r*<sup>33</sup> 2*r*31*r*<sup>33</sup> 2*r*31*r*<sup>32</sup>

. Finally, we verify Eq. (4.3). The elastic constitutive

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

)**S** = *C***S**.

*R*, the Bond transformation matrix *M* is

*r*2 <sup>11</sup> *<sup>r</sup>*<sup>2</sup>

*r*2 <sup>21</sup> *<sup>r</sup>*<sup>2</sup>

*r*2 <sup>31</sup> *<sup>r</sup>*<sup>2</sup>

manipulation shows that *N*−<sup>1</sup> = *M<sup>t</sup>*

**C: Rotation of governing equations**

3 ∑ *j*=1

= �

*∂eij ∂xj* = 3 ∑ *j*=1

∇ · ¯ (*ER*)

Aided by Eq. (8.5), Eq. (8.6) can be rewritten as:

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

�

(∇ · *E*)*<sup>i</sup>* =

<sup>12</sup> *<sup>r</sup>*<sup>2</sup>

<sup>22</sup> *<sup>r</sup>*<sup>2</sup>

<sup>32</sup> *<sup>r</sup>*<sup>2</sup>

Appendix 8. Consider the *i*th component of ∇ · *E* with *i* = 1, 2, 3:

3 ∑ *k*=1 *∂eij ∂x*¯*<sup>k</sup>*

Next, apply coordinate transformation to the equation of motion, Eq. (1.2):

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

∇ · ¯ (*TR*) = <sup>1</sup>

Then perform the coordinate transformation for constitutive equations, Eq. (1.3):

*Rt*

*ρ* ∇ · ¯ (*R<sup>t</sup>*

*c*[4]

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

*∂x*¯*<sup>k</sup> ∂xj* = 3 ∑ *j*=1

The above equation shows that ∇ · *E* = ∇ · ¯ (*ER*) because *R* is a constant spatial tensor.

⎛

⎜⎜⎜⎜⎜⎜⎝

*M* def =

is the rotation matrix and it is orthonormal, i.e., *R*−<sup>1</sup> = *R<sup>t</sup>* or ∑<sup>3</sup> *<sup>k</sup>*=<sup>1</sup> *rikrjk* <sup>=</sup> <sup>∑</sup><sup>3</sup> *<sup>k</sup>*=<sup>1</sup> *rkirkj* = *δij*, *i*, *j* = 1, 2, 3. The natural bases of the global and local coordinate systems are also related through *R* defined in Eq. (8.1):

$$\mathbf{e}\_{l} = \sum\_{j=1}^{3} r\_{i\bar{j}} \mathbf{\bar{e}}\_{j\prime} \quad \mathbf{\bar{e}}\_{l} = \sum\_{j=1}^{3} r\_{j\bar{l}} \mathbf{e}\_{\bar{j}\prime} \quad i = 1, 2, 3. \tag{8.2}$$

To proceed, consider the transformation of a tensor *E*:

$$E\mathbf{b} = R(\bar{E}\bar{\mathbf{b}}) = R(\bar{E}R^{\dagger}R\bar{\mathbf{b}}) = (R\bar{E}R^{\dagger})(R\bar{\mathbf{b}}) = (R\bar{E}R^{\dagger})\mathbf{b} \implies E = R\bar{E}R^{\dagger} \tag{8.3}$$

where *E* is defined in the global coordinates, and *E*¯ in the local coordinates.

To transform the 6-component stress vector **T**, we use the Bond matrix [2]. Aided by matrix *R*, the Bond transformation matrix *M* is

$$M \stackrel{\text{def}}{=} \begin{pmatrix} r\_{11}^2 & r\_{12}^2 & r\_{13}^2 & 2r\_{12}r\_{13} & 2r\_{11}r\_{13} & 2r\_{11}r\_{12} \\ r\_{21}^2 & r\_{22}^2 & r\_{23}^2 & 2r\_{22}r\_{23} & 2r\_{21}r\_{23} & 2r\_{21}r\_{22} \\ r\_{31}^2 & r\_{32}^2 & r\_{33}^2 & 2r\_{32}r\_{33} & 2r\_{31}r\_{33} & 2r\_{31}r\_{32} \\ r\_{21}r\_{31} & r\_{22}r\_{32} & r\_{23}r\_{33} & r\_{22}r\_{33} + r\_{32}r\_{23} & r\_{21}r\_{33} + r\_{23}r\_{31} & r\_{21}r\_{32} + r\_{22}r\_{31} \\ r\_{11}r\_{31} & r\_{12}r\_{32} & r\_{13}r\_{33} & r\_{12}r\_{33} + r\_{13}r\_{32} & r\_{11}r\_{33} + r\_{13}r\_{31} & r\_{11}r\_{32} + r\_{12}r\_{31} \\ r\_{11}r\_{21} & r\_{12}r\_{21} & r\_{13}r\_{23} & r\_{12}r\_{23} + r\_{13}r\_{22} & r\_{11}r\_{23} + r\_{13}r\_{21} & r\_{11}r\_{22} + r\_{12}r\_{21} \end{pmatrix} . \tag{8.4}$$

The application of *M* to the stress vector is **T** = *M***T**¯ , where **T** is the 6-component stress vector defined in the global coordinates and **T**¯ in local coordinates. Similarly, the Bond matrix for strain *N* [2] is used to transform the 6-component strain vector **S** as **S** = *N***S**¯. Simple manipulation shows that *N*−<sup>1</sup> = *M<sup>t</sup>* . Finally, we verify Eq. (4.3). The elastic constitutive equations in the global and local coordinate systems are stated as **T** = *C***S** and **T**¯ = *C*¯**S**¯. Aided by the Bond's matrices *M* and *N*, *C*¯ can be related to *C* in the following equation:

$$\mathbf{T} = M\mathbf{\bar{T}} = M(\mathbf{\bar{C}}\bar{\mathbf{S}}) = M\mathbf{\bar{C}}(N^{-1}\mathbf{S}) = (M\mathbf{\bar{C}}M^{\dagger})\mathbf{S} = C\mathbf{S}.$$

#### **C: Rotation of governing equations**

18 Will-be-set-by-IN-TECH

� *A B C D*

where *A*, *B*, *C*, and *D* are *p* × *p*, *p* × *q*, *q* × *p*, and *q* × *q* matrices, respectively. If *A*−<sup>1</sup> and *D*−<sup>1</sup> exist, then matrices *D* − *CA*−1*B* and *A* − *BD*−1*C* are the Schur complements of *A* and *D* [15].

� ,

� � *A B*

0*q*×*<sup>p</sup> D* − *CA*−1*B*

� � *Ip* 0*p*×*<sup>q</sup> D*−1*C Iq*

*bμ***e**¯*μ*. Let **b** denote the vector defined in the global coordinates and **b**¯

⎞

)(*R***b**¯)=(*RER*¯ *<sup>t</sup>*

� ,

� .

⎠ (8.1)

*<sup>k</sup>*=<sup>1</sup> *rkirkj* = *δij*,

, (8.3)

*<sup>k</sup>*=<sup>1</sup> *rikrjk* <sup>=</sup> <sup>∑</sup><sup>3</sup>

*rji***e***j*, *i* = 1, 2, 3. (8.2)

)**b** ⇒ *E* = *RER*¯ *<sup>t</sup>*

*M* =

� *Ip* 0*p*×*<sup>q</sup> CA*−<sup>1</sup> *Iq*

> � *A* − *BD*−1*C B* 0*q*×*p D*

This appendix provides the transformation relations for the global and local coordinate systems. Let **e**1, **e**2, and **e**<sup>3</sup> be the natural basis of the global coordinate system, and let **e**¯1, ¯**e**2, and ¯**e**<sup>3</sup> be the natural basis of the local coordinate system. An arbitrary vector **b** can be represented by using either the global coordinates or the local coordinates as: **b** =

defined in the local coordinates. The transformation between **b** and **b**¯ is written as: **b** = *R***b**¯,

*i*, *j* = 1, 2, 3. The natural bases of the global and local coordinate systems are also related

3 ∑ *j*=1

*r*<sup>11</sup> *r*<sup>12</sup> *r*<sup>13</sup> *r*<sup>21</sup> *r*<sup>22</sup> *r*<sup>23</sup> *r*<sup>31</sup> *r*<sup>32</sup> *r*<sup>33</sup>

*R* =

*rij***e**¯*j*, ¯**e***<sup>i</sup>* =

*R***b**¯)=(*RER*¯ *<sup>t</sup>*

where *E* is defined in the global coordinates, and *E*¯ in the local coordinates.

is the rotation matrix and it is orthonormal, i.e., *R*−<sup>1</sup> = *R<sup>t</sup>* or ∑<sup>3</sup>

3 ∑ *j*=1

**e***<sup>i</sup>* =

To proceed, consider the transformation of a tensor *E*:

*E***b** = *R*(*E*¯**b**¯) = *R*(*ER*¯ *<sup>t</sup>*

⎛ ⎝

**8. Appendix**

or

∑3

where

*<sup>μ</sup>*=<sup>1</sup> *<sup>b</sup>μ***e***<sup>μ</sup>* <sup>=</sup> <sup>∑</sup><sup>3</sup>

**A: The Schur complements**

In this case, matrix *M* can be factored out as

� *A B C D*

> � *A B C D*

This factorization is useful to calculate det(*M*).

**B: Rotation of Cartesian coordinate system**

*<sup>μ</sup>*=<sup>1</sup> ¯

through *R* defined in Eq. (8.1):

� =

> � =

Consider matrix *M*, which can be divided into 4 sub matrices:

This appendix shows that the velocity-stress equations, Eq. (1.4), remain in the same form when formulated in different Cartesian coordinate systems. To proceed, we consider an arbitrary tensor *E*. In the following, we show that ∇ · *E* = ∇ · ¯ (*ER*), where ∇· ¯ is the divergence operator defined in the local coordinate system, and *R* is the rotation matrix *R* defined in Appendix 8. Consider the *i*th component of ∇ · *E* with *i* = 1, 2, 3:

$$\begin{split} (\nabla \cdot E)\_i &= \sum\_{j=1}^3 \frac{\partial e\_{ij}}{\partial x\_j} = \sum\_{j=1}^3 \sum\_{k=1}^3 \frac{\partial e\_{ij}}{\partial \bar{x}\_k} \frac{\partial \bar{x}\_k}{\partial x\_j} = \sum\_{j=1}^3 \sum\_{k=1}^3 \frac{\partial e\_{ij}}{\partial \bar{x}\_k} r\_{jk} = \sum\_{k=1}^3 \frac{\partial}{\partial \bar{x}\_k} \left( \sum\_{j=1}^3 e\_{ij} r\_{jk} \right) \\ &= \left( \nabla \cdot (ER) \right)\_i. \end{split} \tag{8.5}$$

The above equation shows that ∇ · *E* = ∇ · ¯ (*ER*) because *R* is a constant spatial tensor.

Next, apply coordinate transformation to the equation of motion, Eq. (1.2):

$$R^t \frac{\partial \mathbf{v}}{\partial t} = \frac{1}{\rho} R^t \nabla \cdot \mathbf{T}. \tag{8.6}$$

Aided by Eq. (8.5), Eq. (8.6) can be rewritten as:

$$\frac{\partial \bar{\mathbf{v}}}{\partial t} = \frac{1}{\rho} \mathbf{R}^t \bar{\nabla} \cdot (TR) = \frac{1}{\rho} \bar{\nabla} \cdot (\mathbf{R}^t TR) = \frac{1}{\rho} \bar{\nabla} \cdot \bar{T}.$$

Then perform the coordinate transformation for constitutive equations, Eq. (1.3):

$$R^tTR = R^t c^{|4|} SR.\tag{8.7}$$

Recast *c*[4] in Eq. (8.7) into the index form:

$$\begin{split} \bar{T}\_{\rm mn} &= r\_{\rm im} T\_{\rm ij} r\_{\rm jn} = r\_{\rm im} c\_{\rm ijkl} \mathbf{S}\_{\rm kl} r\_{\rm jn} = r\_{\rm im} c\_{\rm ijkl} r\_{\rm ko} \bar{\mathbf{S}}\_{\rm op} r\_{\rm lp} r\_{\rm jn} \\ &= r\_{\rm im} r\_{\rm jn} r\_{\rm ko} r\_{\rm lp} c\_{\rm ijkl} \bar{\mathbf{S}}\_{\rm op} = \mathfrak{c}\_{\rm mnop} \bar{\mathbf{S}}\_{\rm op} \end{split} \tag{8.8}$$

[6] Kacimi, A. E. & Laghrouche, O. [2009]. Numerical modelling of elastic wave scattering in frequency domain by the partition of unity finite element method, *International Journal*

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

[7] Karypis, G. & Kumar, V. [1998]. A fast and high quality multilevel scheme for partitioning

[8] Kim, C., Yu, S. J. & Zhang, Z. [2004]. Cavity flow in scramjet engine by Space-Time

[9] Kreiss, H. O. [1973]. *Methods for the approximate solution of time dependent problems*, Global Atmospheric Research Programme - WMO-ICSU Joint Organizing Committee. [10] Käser, M. & Dumbser, M. [2006]. An arbitrary high-order discontinuous galerkin method for elastic waves on unstructured meshes; i. the two-dimensional isotropic case with

[11] Kulikovskii, A., Pogorelov, N. & Semenov, A. Y. [2000]. *Mathematical Aspects of Numerical*

[12] LeVeque, R. J. [2002a]. *Finite-Volume Methods for Hyperbolic Problems*, Cambridge texts in

[13] LeVeque, R. J. [2002b]. Finite-volume methods for non-linear elasticity in heterogeneous

[14] Loh, C., Hultgren, L. & Chang, S. [2001]. Wave computation in compressible flow using space-time conservation element and solution element method, *AIAA Journal*

[15] Meyer, C. D. [2000]. *Matrix Analysis and Applied Linear Algebra*, Society for Industrial and

[16] Musgrave, M. J. P. [1954]. On the propagation of elastic waves in aeolotropic media. II. media of hexagonal symmetry, *Proceedings of the Royal Society of London. Series A,*

[17] Musgrave, M. J. P. [1970]. *Crystal Acoustics; Introduction to the Study of Elastic Waves and*

[18] Puente, J. d. l., Käser, M., Dumbser, M. & Igel, H. [2007]. An arbitrary high-order discontinuous galerkin method for elastic waves on unstructured meshes; IV. anisotropy,

[19] Qin, J. R., Yu, S. J. & Lai, M. [2001]. Direct calculations of cavitating flows in fuel delivery pipe by the Space-Time CESE method, *Journal of Fuels and Lubricants, SAE Transaction*

[20] Reed, W. H. & Hill, T. R. [1973]. Triangular mesh methods for the neutron transport

[21] Saenger, E. H., Gold, N. & Shapiro, S. A. [2000]. Modeling the propagation of elastic

[23] Strang, G. [2006]. *Linear Algebra and Its Applications*, 4th ed edn, Thomson, Brooks/Cole,

[24] Ting, T. C. [1996]. *Anisotropic Elasticity: Theory and Applications*, Oxford University Press,

[25] Virieux, J. [1986]. P-SV wave propagation in heterogeneous media: Velocity-stress

[26] Wang, B., He, H. & Yu, S. J. [2005]. Direct calculation of wave implosion for detonation

equation, *Technical Report LA-UR-73-479*, Los Alamos Scientific Laboratory.

waves using a modified finite-difference grid, *Wave Motion* 31(1): 77–92.

[22] Shorr, B. F. [2004]. *The Wave Finite Element Method*, 1 edn, Springer.

*for Numerical Methods in Engineering* 77(12): 1646–1669.

irregular graphs, *SIAM Journal on Scientific Computing* 20(1): 359–392.

conservation and solution element method, *AIAA Journal* 42(5): 912–919.

external source terms, *Geophysical Journal International* 166(2): 855–877.

media, *International Journal for Numerical Methods in Fluids* 40(1-2): 93–104.

*Solution of Hyperbolic Systems*, 1 edn, Chapman and Hall/CRC.

applied mathematics, Cambridge University Press, Cambridge.

39(5): 794–801.

108: 1720–1725.

Belmont, Calif.

New York.

Applied Mathematics, Philadelphia.

*Mathematical and Physical Sciences* 226(1166): 356–366.

*Vibrations in Crystals*, Holden-day, San Francisco.

*Geophysical Journal International* 169(3): 1210–1228.

finite-difference method, *Geophysics* 51(4): 889–901.

initiation, *AIAA Journal* 43(10): 2157–2169.

where *cijkl* and *c*¯*mnop* are the components of the fourth-order stiffness tensor *c*[4] and *c*¯ [4] , respectively. For conciseness, Eq. (8.8) uses Einstein summation convention. By using the direct notation, Eq. (8.8) can be recast as *T*¯ = *c*¯ [4] *S*¯.

## **D: Supplemental algebraic relations for the direction cosine matrices**

This appendix provides the algebraic relations to facilitate the derivation of Eq. (4.6). Aided by Eqs. (2.2) and (8.4), we expand *K*(*μ*)*M*, *μ* = 1, 2, 3, as:

*K*(1) *M* = ⎛ ⎝ *r*2 <sup>11</sup> *<sup>r</sup>*<sup>2</sup> <sup>12</sup> *<sup>r</sup>*<sup>2</sup> <sup>13</sup> 2*r*12*r*<sup>13</sup> 2*r*11*r*<sup>13</sup> 2*r*11*r*<sup>12</sup> *r*11*r*<sup>21</sup> *r*12*r*<sup>22</sup> *r*13*r*<sup>23</sup> *r*12*r*<sup>23</sup> + *r*13*r*<sup>22</sup> *r*11*r*<sup>23</sup> + *r*13*r*<sup>21</sup> *r*11*r*<sup>22</sup> + *r*12*r*<sup>21</sup> *r*11*r*<sup>31</sup> *r*12*r*<sup>32</sup> *r*13*r*<sup>33</sup> *r*12*r*<sup>33</sup> + *r*13*r*<sup>32</sup> *r*11*r*<sup>33</sup> + *r*13*r*<sup>31</sup> *r*11*r*<sup>32</sup> + *r*12*r*<sup>31</sup> ⎞ ⎠ , (8.9) *K*(2)*M* = ⎛ ⎝ *r*11*r*<sup>21</sup> *r*12*r*<sup>22</sup> *r*13*r*<sup>23</sup> *r*12*r*<sup>23</sup> + *r*13*r*<sup>22</sup> *r*11*r*<sup>23</sup> + *r*13*r*<sup>21</sup> *r*11*r*<sup>22</sup> + *r*12*r*<sup>21</sup> *r*2 <sup>21</sup> *<sup>r</sup>*<sup>2</sup> <sup>22</sup> *<sup>r</sup>*<sup>2</sup> <sup>23</sup> 2*r*22*r*<sup>23</sup> 2*r*21*r*<sup>23</sup> 2*r*21*r*<sup>22</sup> *r*21*r*<sup>31</sup> *r*22*r*<sup>32</sup> *r*23*r*<sup>33</sup> *r*22*r*<sup>33</sup> + *r*32*r*<sup>23</sup> *r*21*r*<sup>33</sup> + *r*23*r*<sup>31</sup> *r*21*r*<sup>32</sup> + *r*22*r*<sup>31</sup> ⎞ ⎠ , (8.10) *K*(3) *M* = ⎛ ⎝ *r*11*r*<sup>31</sup> *r*12*r*<sup>32</sup> *r*13*r*<sup>33</sup> *r*12*r*<sup>33</sup> + *r*13*r*<sup>32</sup> *r*11*r*<sup>33</sup> + *r*13*r*<sup>31</sup> *r*11*r*<sup>32</sup> + *r*12*r*<sup>31</sup> *r*21*r*<sup>31</sup> *r*22*r*<sup>32</sup> *r*23*r*<sup>33</sup> *r*22*r*<sup>33</sup> + *r*32*r*<sup>23</sup> *r*21*r*<sup>33</sup> + *r*23*r*<sup>31</sup> *r*21*r*<sup>32</sup> + *r*22*r*<sup>31</sup> *r*2 <sup>31</sup> *<sup>r</sup>*<sup>2</sup> <sup>32</sup> *<sup>r</sup>*<sup>2</sup> <sup>33</sup> 2*r*32*r*<sup>33</sup> 2*r*31*r*<sup>33</sup> 2*r*31*r*<sup>32</sup> ⎞ ⎠ . (8.11)

By comparing Eqs. (8.9), (8.10), and (8.11) with Eq. (4.7), Eq. (4.6) is reached.

## **9. References**


[6] Kacimi, A. E. & Laghrouche, O. [2009]. Numerical modelling of elastic wave scattering in frequency domain by the partition of unity finite element method, *International Journal for Numerical Methods in Engineering* 77(12): 1646–1669.

20 Will-be-set-by-IN-TECH

*o p* = *c*¯*mnopS*¯

*o prl prjn*

*o p*, (8.8)

⎞

⎞

⎞

⎠ , (8.9)

⎠ , (8.10)

⎠ . (8.11)

[4] ,

*mn* = *rimTijrjn* = *rimcijklSklrjn* = *rimcijklrkoS*¯

where *cijkl* and *c*¯*mnop* are the components of the fourth-order stiffness tensor *c*[4] and *c*¯

**D: Supplemental algebraic relations for the direction cosine matrices**

respectively. For conciseness, Eq. (8.8) uses Einstein summation convention. By using the

This appendix provides the algebraic relations to facilitate the derivation of Eq. (4.6). Aided

*r*11*r*<sup>21</sup> *r*12*r*<sup>22</sup> *r*13*r*<sup>23</sup> *r*12*r*<sup>23</sup> + *r*13*r*<sup>22</sup> *r*11*r*<sup>23</sup> + *r*13*r*<sup>21</sup> *r*11*r*<sup>22</sup> + *r*12*r*<sup>21</sup> *r*11*r*<sup>31</sup> *r*12*r*<sup>32</sup> *r*13*r*<sup>33</sup> *r*12*r*<sup>33</sup> + *r*13*r*<sup>32</sup> *r*11*r*<sup>33</sup> + *r*13*r*<sup>31</sup> *r*11*r*<sup>32</sup> + *r*12*r*<sup>31</sup>

*r*11*r*<sup>21</sup> *r*12*r*<sup>22</sup> *r*13*r*<sup>23</sup> *r*12*r*<sup>23</sup> + *r*13*r*<sup>22</sup> *r*11*r*<sup>23</sup> + *r*13*r*<sup>21</sup> *r*11*r*<sup>22</sup> + *r*12*r*<sup>21</sup>

*r*21*r*<sup>31</sup> *r*22*r*<sup>32</sup> *r*23*r*<sup>33</sup> *r*22*r*<sup>33</sup> + *r*32*r*<sup>23</sup> *r*21*r*<sup>33</sup> + *r*23*r*<sup>31</sup> *r*21*r*<sup>32</sup> + *r*22*r*<sup>31</sup>

*r*11*r*<sup>31</sup> *r*12*r*<sup>32</sup> *r*13*r*<sup>33</sup> *r*12*r*<sup>33</sup> + *r*13*r*<sup>32</sup> *r*11*r*<sup>33</sup> + *r*13*r*<sup>31</sup> *r*11*r*<sup>32</sup> + *r*12*r*<sup>31</sup> *r*21*r*<sup>31</sup> *r*22*r*<sup>32</sup> *r*23*r*<sup>33</sup> *r*22*r*<sup>33</sup> + *r*32*r*<sup>23</sup> *r*21*r*<sup>33</sup> + *r*23*r*<sup>31</sup> *r*21*r*<sup>32</sup> + *r*22*r*<sup>31</sup>

[1] Arfken, G. [2005]. *Mathematical methods for physicists*, 6th ed. / edn, Elsevier, Boston.

[3] Cai, M., Yu, S. J. & Zhang, M. [2006]. Theoretical and numerical solutions of linear and nonlinear elastic waves in a thin rod, *AIAA 2006-4778, 42nd AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, Sacramento, Califonia, July 9-12, 2006*, Sacramento,

[4] Chang, S. [1995]. The method of space-time conservation element and solution element - a new approach for solving the Navier-Stokes and euler equations, *J. Comput. Phys.*

[5] Chang, S. & Wang, X. [2003]. Multi-Dimensional courant number insensitive CE/SE euler solvers for applications involving highly nonuniform meshes, *39th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit*, Huntsville, Alabama.

[2] Auld, B. A. [1990]. *Acoustic Fields and Waves in Solids*, 2nd edn, R.E. Krieger.

<sup>13</sup> 2*r*12*r*<sup>13</sup> 2*r*11*r*<sup>13</sup> 2*r*11*r*<sup>12</sup>

<sup>23</sup> 2*r*22*r*<sup>23</sup> 2*r*21*r*<sup>23</sup> 2*r*21*r*<sup>22</sup>

<sup>33</sup> 2*r*32*r*<sup>33</sup> 2*r*31*r*<sup>33</sup> 2*r*31*r*<sup>32</sup>

[4] *S*¯.

= *rimrjnrkorl pcijklS*¯

Recast *c*[4] in Eq. (8.7) into the index form:

*T*¯

direct notation, Eq. (8.8) can be recast as *T*¯ = *c*¯

*K*(1) *M* = ⎛ ⎝

*r*2 <sup>11</sup> *<sup>r</sup>*<sup>2</sup>

*r*2 <sup>21</sup> *<sup>r</sup>*<sup>2</sup>

*r*2 <sup>31</sup> *<sup>r</sup>*<sup>2</sup>

*K*(2)*M* = ⎛ ⎝

*K*(3) *M* = ⎛ ⎝

**9. References**

Califonia.

119(2): 295–324.

by Eqs. (2.2) and (8.4), we expand *K*(*μ*)*M*, *μ* = 1, 2, 3, as:

<sup>12</sup> *<sup>r</sup>*<sup>2</sup>

<sup>22</sup> *<sup>r</sup>*<sup>2</sup>

<sup>32</sup> *<sup>r</sup>*<sup>2</sup>

By comparing Eqs. (8.9), (8.10), and (8.11) with Eq. (4.7), Eq. (4.6) is reached.


22 Will-be-set-by-IN-TECH

[27] Wang, X. & Chang, S. [1999]. A 2D Non-Splitting unstructured triangular mesh euler solver based on the Space-Time conservation element and solution element method,

[28] Wang, X., Chen, C. & Liu, Y. [2002]. The space-time CE/SE method for solving maxwell's equations in time-domain, *Antennas and Propagation Society International Symposium, 2002.*

[29] Warming, R. F., Beam, R. M. & Hyett, B. J. [1975]. Diagonalization and simultaneous symmetrization of the Gas-Dynamic matrices, *Mathematics of Computation*

[30] Yang, L., Chen, Y. & Yu, S. J. [2011]. Velocity-Stress equations for waves in solids of hexagonal symmetry solved by the Space-Time CESE method, *ASME Journal of Vibration*

[31] Yang, L., Lowe, R. L., Yu, S. J. & Bechtel, S. E. [2010]. Numerical solution by the CESE method of a First-Order hyperbolic form of the equations of dynamic nonlinear elasticity,

[32] Yu, S. J., Yang, L., Lowe, R. L. & Bechtel, S. E. [2010]. Numerical simulation of linear and nonlinear waves in hypoelastic solids by the CESE method, *Wave Motion* 47(3): 168–182. [33] Zhang, J. & Gao, H. [2009]. Elastic wave modelling in 3-D fractured media: an explicit

[34] Zhang, M., Yu, S. J., Lin, S., Chang, S. & Blankson, I. [2004]. Solving magnetohydrodynamic equations without special treatment for Divergence-Free

[35] Zhang, M., Yu, S. J., Lin, S. H., Chang, S. & Blankson, I. [2006]. Solving the MHD equations by the space-time conservation element and solution element method, *Journal*

[36] Zhang, Z., Yu, S. T. J. & Chang, S. [2002]. A Space-Time conservation element and solution element method for solving the two- and Three-Dimensional unsteady euler equations using quadrilateral and hexahedral meshes, *Journal of Computational Physics*

*Computational Fluid Dynamics Journal* 8(2): 309–325.

*ASME Journal of Vibrations and Acoustics* 132(5): 051003.

magnetic field, *AIAA Journal* 42(12): 2605–2608.

*of Computational Physics* 214(2): 599–617.

approach, *Geophysical Journal International* 117(3): 1233–1241.

*IEEE*, Vol. 1, pp. 164–167 vol.1.

*and Acoustics* 133(2): 021001.

29(132): 1037–1045.

222 Wave Processes in Classical and New Solids

175(1): 168–199.
