**2.3 Testing the scheme and the computer program 2.3.1 Stability of the scheme**

The stability of the scheme is the core problem of finite difference method which depends upon the choice of parameter *c h* <sup>1</sup> / , which is the ratio between time step and space step substantively. The choice is related to the ratio *c c* 1 2 / , i.e., the ratio between speeds of elastic longitudinal and transverse waves of the phonon field. To determine the upper bound for the ration to guarantee the stability, according to our computational practice and considering the experiences of computations for conventional materials, we choose 0.8 in all cases and results are stable.

#### **2.3.2 Accuracy test**

The stability is only a necessary condition for successful computation. We must check the accuracy of the numerical solution. This can be realized through some comparison with some well-known classical solutions (analytic as well as numerical solutions) of conventional fracture mechanics. For this purpose the material constants in the computation are chosen as 1 2 *c c* 7.34 mm/μs, 3.92 mm/ μs and 3 3 5 10 kg/m , <sup>0</sup> *p* 1 MPa which

Fig. 2. (c) Displacement component of phason field *wx* versus time

also shows the mathematical modeling of the present work is valid.

**2.3 Testing the scheme and the computer program** 

 *c h* <sup>1</sup> 

**2.3.1 Stability of the scheme** 

upon the choice of parameter

in all cases and results are stable.

**2.3.2 Accuracy test** 

and the phonon-phason coupling effect. From Fig. 2(c), in the phason field we find that the phason mode presents diffusive nature in the overall tendency, but because of influence of the phonon and phonon-phason coupling, it can also have some characters of fluctuation. So the model describes the dynamic behaviour of phonon field and phason field in deed. This

The stability of the scheme is the core problem of finite difference method which depends

substantively. The choice is related to the ratio *c c* 1 2 / , i.e., the ratio between speeds of elastic longitudinal and transverse waves of the phonon field. To determine the upper bound for the ration to guarantee the stability, according to our computational practice and

The stability is only a necessary condition for successful computation. We must check the accuracy of the numerical solution. This can be realized through some comparison with some well-known classical solutions (analytic as well as numerical solutions) of conventional fracture mechanics. For this purpose the material constants in the computation

considering the experiences of computations for conventional materials, we choose

are chosen as 1 2 *c c* 7.34 mm/μs, 3.92 mm/ μs and 3 3

/ , which is the ratio between time step and space step

0.8

5 10 kg/m , <sup>0</sup> *p* 1 MPa which

are the same with those given in classical references for conventional fracture dynamics, discussed in Fan's monograph [1] in detail. At first the comparison to the classical exact analytic solution is carried out, in this case we put 0 *w w x y* (i.e., 1 2 *KKR* 0 ) for the numerical solution. The comparison has been done with the key physical quantity dynamic stress intensity factor, which is defined by

$$K\_1(t) = \lim\_{\mathbf{x} \to \mathbf{a}\_0^+} \sqrt{\pi(\mathbf{x} - a\_0)} \sigma\_{yy}(\mathbf{x}, 0, t) \tag{10}$$

The normalized dynamic stress intensity factor can be denoted as ( )/ *static Kt K I I* , in which *static KI* is the corresponding static stress intensity factor, whose value here is taken as 0 0 *a p* **.** For the dynamic initiation of crack growth in classical fracture dynamics there is the only exact analytic solution— the Maue's solution (refer to Fan's monograph [1]), but the configuration of whose specimen is quite different from that of our specimen. Maue studied a semi-infinite crack in an infinite body, and subjected to a Heaviside impact loading at the crack surface. While our specimen is a finite size rectangular plate with a central crack, and the applied stress is at the external boundary of the specimen. Generally the Maue's model cannot describe the interaction between wave and external boundary. However, consider a very short time interval, i.e., during the period between the stress wave from the external boundary arriving at the crack tip (this time is denoted by <sup>1</sup>*t* ) and before the reflecting by external boundary stress wave emanating from the crack tip in the finite size specimen (the time is marked as <sup>2</sup>*t* ). During this special very short time interval our specimen can be seen as an "infinite specimen". The comparison given by Fig. 3 shows the numerical results are in excellent agreement with those of Maue's solution within the short interval in which the solution is valid.

Our solution corresponding to case of 0 *w w x y* is also compared with numerical

solutions of conventional crystals, e.g. Murti's solution and Chen's solutions (refer to Fan [1] and Zhu and Fan [9] for the detail), which are also shown in Fig. 3, it is evident, our solution presents very high precise.
