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

12 Will-be-set-by-IN-TECH

x¯1

x 3

x¯3

det( ¯ *<sup>A</sup>*(1) <sup>−</sup> *<sup>λ</sup>I*9) = det

non-trivial eigenvalues of *A*¯(1) can be obtained by solving a simpler equation

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

*h*<sup>3</sup> = *r*<sup>13</sup> = sin *φ*. Aided by the definition of **h**¯, Eq. (5.9) can be rewritten as Γ<sup>11</sup> = *c*<sup>11</sup> ¯

Γ<sup>22</sup> = *c*<sup>11</sup> ¯

Γ<sup>33</sup> = *c*<sup>44</sup> ¯

*ρ*

⎛ ⎝

Γ<sup>11</sup> = *c*<sup>11</sup> cos<sup>2</sup> *φ* cos2 *θ* + *c*<sup>44</sup> sin2 *φ* + *c*<sup>66</sup> cos<sup>2</sup> *φ* sin2 *θ*, Γ<sup>22</sup> = *c*<sup>11</sup> cos<sup>2</sup> *φ* sin2 *θ* + *c*<sup>44</sup> sin2 *φ* + *c*<sup>66</sup> cos2 *φ* cos<sup>2</sup> *θ*, Γ<sup>33</sup> = *c*<sup>33</sup> sin2 *φ* + *c*<sup>44</sup> cos2 *φ* sin2 *θ* + *c*<sup>44</sup> cos2 *φ* cos<sup>2</sup> *θ*. Γ<sup>12</sup> = Γ<sup>21</sup> = *c*<sup>12</sup> cos<sup>2</sup> *φ* cos *θ* sin *θ* + *c*<sup>66</sup> cos<sup>2</sup> *φ* cos *θ* sin *θ*, Γ<sup>13</sup> = Γ<sup>31</sup> = *c*<sup>13</sup> cos *φ* sin *φ* cos *θ* + *c*<sup>44</sup> cos *φ* sin *φ* cos *θ*, Γ<sup>23</sup> = Γ<sup>32</sup> = *c*<sup>13</sup> cos *φ* sin *φ* sin *θ* + *c*<sup>44</sup> cos *φ* sin *φ* sin *θ*.

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

*<sup>h</sup>*<sup>2</sup> <sup>=</sup> *<sup>r</sup>*<sup>12</sup> <sup>=</sup> cos *<sup>φ</sup>* sin *<sup>θ</sup>*, and ¯

*h*2 <sup>1</sup> <sup>+</sup> *<sup>c</sup>*<sup>44</sup> ¯ *h*2 <sup>3</sup> <sup>+</sup> *<sup>c</sup>*<sup>66</sup> ¯ *h*2 2,

*h*2 <sup>2</sup> <sup>+</sup> *<sup>c</sup>*<sup>44</sup> ¯ *h*2 <sup>3</sup> <sup>+</sup> *<sup>c</sup>*<sup>66</sup> ¯ *h*2 1,

*h*2 <sup>1</sup> <sup>+</sup> *<sup>c</sup>*<sup>44</sup> ¯ *h*2 <sup>2</sup> + *<sup>c</sup>*<sup>33</sup> ¯ *h*2 3,

Γ<sup>12</sup> = Γ<sup>21</sup> = (*c*<sup>12</sup> + *c*66)¯

Γ<sup>13</sup> = Γ<sup>31</sup> = (*c*<sup>13</sup> + *c*44)¯

Γ<sup>23</sup> = Γ<sup>32</sup> = (*c*<sup>13</sup> + *c*44)¯

*h*<sup>1</sup> = *r*<sup>11</sup> = cos *φ* cos *θ*, ¯

*h*1 ¯ *h*2,

*h*1 ¯ *h*3,

*h*2 ¯ *h*3.

x 1

(b)

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

where *In* is a rank *n* identity matrix. Based on the Schur complement [15, pp. 475], the

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

φ

x <sup>2</sup>x¯<sup>2</sup>

� −*λI*<sup>3</sup> *<sup>A</sup>*¯ *<sup>v</sup> A*¯ *<sup>T</sup>* −*λI*<sup>6</sup>

Γ<sup>11</sup> Γ<sup>12</sup> Γ<sup>13</sup> Γ<sup>21</sup> Γ<sup>22</sup> Γ<sup>23</sup> Γ<sup>31</sup> Γ<sup>32</sup> Γ<sup>33</sup> x1

� = 0,

det(*A*¯ *vA*¯ *<sup>T</sup>* − *κ*<sup>2</sup> *I*3) = 0, (5.8)

⎞ ⎠ , x3

θ

*h*3) be the unit vector parallel to the *x*¯1

(c)

x¯3

x2

(5.9)

(5.10)

<sup>x</sup>¯<sup>1</sup> <sup>x</sup>¯<sup>2</sup>

φ

x1

polynomial:

where

and

x3x3

θ

x1

(a)

x2

Rotation with respect to the *x*¯2 axis. (c) Overall rotation.

To proceed, according to Eq. (5.6), we let **h**¯ = (¯

axis in the local coordinate system, i.e., ¯

x2

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 condition is chosen to be 5 × 10−<sup>3</sup> m.

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

**Figure 3.** Two-dimensional mesh (2.2 million triangular elements). (a) close look at the mesh around origin. (b) decomposed 16 sub-domains.

(a) (b)

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

solid. The orientation of the solid is set to be *<sup>θ</sup>* = <sup>0</sup>◦ and *<sup>φ</sup>* = <sup>90</sup>◦, to have *<sup>x</sup>*(1) aligned to *<sup>x</sup>*3, and *x*(2) aligned to *x*2. The group velocities have the analytical solution. The group velocities of pure-shear (SH), quasi-longitudinal (qL), and quasi-shear (qS) waves in a block of beryl are plotted in Fig. 4(a) [16]. The simulation result, shown in Fig. 4(b), compares well with the

In Case 2, the orientation of the computational plane is in the direction of *θ* = 20◦ and *φ* = 70◦. Similar to that in Case 1, the initial condition is a velocity source at the origin with *σ* = 0.00 m. Figure 5 shows two snapshots of the calculated energy profiles at *t* = 26 and 130 *μ*s. With the non-reflecting boundary condition implemented on the four boundaries of the computational

For Case 3, we calculate wave propagation in a domain composed of three blocks of beryl with three different lattice orientations. The partition of the domain is shown in Fig. 6. Shown in Fig. 7(a), an initial, cylindrical wave is generated at the origin. This is because of the lattice orientation of the central block such that an isotropic wave expansion occurs. This wave expands and interacts with the interfaces separating the neighbouring blocks of beryl with different lattice orientations. Figure 7 shows four snapshots of the calculated energy profiles at different times. Due to complex wave/interface interactions, all wave modes, including SH, qL and qS are excited at the interface. As a results, refracted waves at different speeds can be seen in Figs. 7(b) and 7(c). The original circular wave front is fractured after passing through the interfaces, as shown in Fig. 7(d). Complex wave features have been successfully

In this chapter, we have presented the first-order velocity-stress equations as an alternative to the conventional second-order wave equations for modeling wave propagation in anisotropic

domain, waves propagate out of the domain without spurious reflection.

**Figure 5.** Plots for the normalized total energy density (*e*/ max(*e*)) at two times: (a) *t* = 26 *μ*s, (b)

*t* = 130 *μ*s.

analytical solution.

captured by the numerical results.

**7. Concluding remarks**

**Figure 4.** Comparison between the analytical solutions of the group velocities (in SH, qS, and qL polarization) and the calculated energy profiles for beryl at *t* = 91 *μ*s, including: (a) The analytical solution, and (b) The calculated energy density profiles.

numerical accuracy, (2) waves propagation in a plane of an arbitrary orientation to assess the non-reflective boundary condition, and (3) waves propagation in a plane composed of three blocks of beryl in different lattice orientations.

In Case 1, two-dimensional numerical simulations are performed to calculate the density of the total energy, which is the summation of the kinetic energy and the strain energy normalized by area. For elastic solids, the group velocity is identical to the energy velocity [2]. Here, we consider wave propagation in 100 (*x*2-*x*3) or 010 (*x*1-*x*3) plane of the hexagonal

**Figure 5.** Plots for the normalized total energy density (*e*/ max(*e*)) at two times: (a) *t* = 26 *μ*s, (b) *t* = 130 *μ*s.

solid. The orientation of the solid is set to be *<sup>θ</sup>* = <sup>0</sup>◦ and *<sup>φ</sup>* = <sup>90</sup>◦, to have *<sup>x</sup>*(1) aligned to *<sup>x</sup>*3, and *x*(2) aligned to *x*2. The group velocities have the analytical solution. The group velocities of pure-shear (SH), quasi-longitudinal (qL), and quasi-shear (qS) waves in a block of beryl are plotted in Fig. 4(a) [16]. The simulation result, shown in Fig. 4(b), compares well with the analytical solution.

In Case 2, the orientation of the computational plane is in the direction of *θ* = 20◦ and *φ* = 70◦. Similar to that in Case 1, the initial condition is a velocity source at the origin with *σ* = 0.00 m. Figure 5 shows two snapshots of the calculated energy profiles at *t* = 26 and 130 *μ*s. With the non-reflecting boundary condition implemented on the four boundaries of the computational domain, waves propagate out of the domain without spurious reflection.

For Case 3, we calculate wave propagation in a domain composed of three blocks of beryl with three different lattice orientations. The partition of the domain is shown in Fig. 6. Shown in Fig. 7(a), an initial, cylindrical wave is generated at the origin. This is because of the lattice orientation of the central block such that an isotropic wave expansion occurs. This wave expands and interacts with the interfaces separating the neighbouring blocks of beryl with different lattice orientations. Figure 7 shows four snapshots of the calculated energy profiles at different times. Due to complex wave/interface interactions, all wave modes, including SH, qL and qS are excited at the interface. As a results, refracted waves at different speeds can be seen in Figs. 7(b) and 7(c). The original circular wave front is fractured after passing through the interfaces, as shown in Fig. 7(d). Complex wave features have been successfully captured by the numerical results.
