Dynamic Analysis of Cracked Octagonal Quasicrystals

Abstract

We focus on the dynamic fracture problem of octagonal quasicrystals by applying a rectangular sample with a Griffith crack which is often used in classical elastic media based on the method of finite difference. This paper mainly is to investigate the variation of phonon, phason fields, and stress singularity around the crack tip including the stress intensity factor. In addition, the moving boundary due to the crack propagation has also been treated by introducing an additional condition for determining solution. The influence of wave propagation and diffusion in the dynamic process is also discussed in detail. Through comparing the results of octagonal quasicrystals with the results of crystal, this paper proclaims the influence of phonon-phason coupling in dynamic fracture problem of octagonal quasicrystals which should not be neglected.

1. Introduction

Quasicrystals are a form of solids different from both crystals and glassy solids, which possess quasiperiodic long-range translational symmetry and noncrystallographic rotational symmetry [1]. Due to the quasiperiodicity, quasicrystals have a special type of elastic degrees of freedom, termed as phason degrees of freedom apart from phonon degrees which exist too for crystals. Under given load conditions, they are giving rise to two displacement fields 𝐮(𝐫,𝑡) for phonon field and 𝐰(𝐫,𝑡) for phason field, which also yield the elastic strain tensors 𝜀_𝑖𝑗 and 𝑤_𝑖𝑗, respectively [2–4],𝜀_𝑖𝑗≡12𝜕𝜕𝑥𝑢^𝑖_𝑗+𝜕𝑢^𝑗𝜕𝑥_𝑖,𝑤^𝑖𝑗≡𝜕𝑤^𝑖𝜕𝑥_𝑗𝑤_𝑖𝑗≠𝑤_𝑗𝑖.(1.1) In linear elasticity the stress tensors are related to the strain tensors obtained by [2]𝜎_𝑖𝑗=𝐶_{𝑖𝑗𝑘𝑙}𝜀_𝑘𝑙+𝑅_{𝑖𝑗𝑘𝑙}𝑤_𝑘𝑙,𝐻_𝑖𝑗=𝑅_{𝑘𝑙𝑖𝑗}𝜀_𝑘𝑙+𝐾_{𝑖𝑗𝑘𝑙}𝑤_𝑘𝑙,(1.2) where 𝑢_𝑖, 𝑤_𝑖 are phonon and phason displacements, 𝜎_𝑖𝑗 and 𝜀_𝑖𝑗 phonon stresses and strains, 𝐻_𝑖𝑗 and 𝑤_𝑖𝑗 phason stresses and strains, 𝐶_{𝑖𝑗𝑘𝑙}, 𝐾_{𝑖𝑗𝑘𝑙}, and 𝑅_{𝑖𝑗𝑘𝑙} the phonon, phason, phonon-phason coupling elastic constants, respectively. But for the elasto-/hydro-dynamics, there are different theoretical points of view, for example, those being put forward by Bak [5, 6] and by Lubensky et al [7]. Here we adopt the viewpoint raised by Fan [2]:𝜌𝜕²𝑢_𝑖𝜕𝑡²=𝜕𝜎^𝑖𝑗𝜕𝑥_𝑗1,Γ_𝑤𝜕𝑤_𝑖𝜕𝑡=𝜕𝜕𝐻𝑥^𝑖𝑗_𝑗,(1.3) in which Γ_𝑤=1/𝜅 denotes the kinetic coefficient of phason field.

For this subject, Zhu and Fan and Wang et al. have used the icosahedral Al-Pd-Mn quasicrystal and the decagonal Al-Ni-Co quasicrystal to simulate dynamic behaviour of quasicrystals [8, 9]. However, little attention has been devoted to the octagonal quasicrystals which occupy a very important position in this solids. Of course, there are many theoretical and experimental results reported up to now [10–33]. Their results could be better convinced if they have considered the dynamic behaviour of octagonal quasicrystals.

In the following sections, we will reveal the dynamic fracture behaviour of octagonal quasicrystals by introducing a rectangular sample with a Griffith crack based on the method of finite difference. In the end, we explain physical significance of the results obtained in the paper. Of course, further studies of the quasicrystals should be consummated by the future results. It is hoped that the dynamics question of quasicrystals will be observed by experiments.

2. Basic Formulas and Simplified Model

We focus on the solution for the plane elasticity of the octagonal quasicrystals that is, the atom arrangement along the 𝑧-direction is periodic and along the 𝑥−𝑦 plane is quasiperiodic. For the octagonal quasicrystals, we have the elastic constants matrix [𝐂𝐊𝐑]:[𝐂𝐊𝐑]=⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣𝐿+𝑅𝐿002𝑀𝐿+−𝐿00𝑅2𝑀𝑀𝑀000𝑀𝑀000−𝐾𝑅00𝑅₁−𝑅00𝑅−−00𝑅𝑅𝑅𝑅00𝑅−𝑅00𝐾₂𝐾₂0000−𝑅−𝑅0𝐾0₁𝐾₁00+𝐾₂+𝐾₃𝐾₃00𝑅𝑅00𝐾₃𝐾₁+𝐾₂+𝐾₃⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦.(2.1) In this case, the stress-strain relations are reduced to (2.2) from (2.1),𝜎𝜎_{𝑥𝑥𝑦𝑦}=𝐿𝜀_𝑥𝑥+𝜀_𝑦𝑦+2𝑀𝜀_𝑥𝑥+𝑅𝑤_𝑥𝑥+𝑤_𝑦𝑦,𝜎_𝑥𝑦=𝐿𝜀_𝑥𝑥+𝜀_𝑦𝑦+2𝑀𝜀_𝑦𝑦−𝑅𝑤_𝑥𝑥+𝑤_𝑦𝑦,𝐻_𝑥𝑥=𝜎_𝑦𝑥=2𝑀𝜀_𝑥𝑦+𝑅𝑤_𝑦𝑥−𝑤_𝑥𝑦,𝐻_𝑦𝑦=𝐾₁𝑤_𝑥𝑥+𝐾₂𝑤_𝑦𝑦+𝑅𝜀_𝑥𝑥−𝜀_𝑦𝑦,𝐻_𝑥𝑦=𝐾₁𝑤_𝑦𝑦+𝐾₂𝑤_𝑥𝑥+𝑅𝜀_𝑥𝑥−𝜀_𝑦𝑦,𝐻_𝑦𝑥=𝐾₁+𝐾₂+𝐾₃𝑤_𝑥𝑦+𝐾₃𝑤_𝑦𝑥−2𝑅𝜀_𝑥𝑦,=𝐾₁+𝐾₂+𝐾₃𝑤_𝑦𝑥+𝐾₃𝑤_𝑥𝑦+2𝑅𝜀_𝑥𝑦,(2.2) where 𝐿=𝐶₁₂, 𝑀=(𝐶₁₁−𝐶₁₂)/2 are the phonon elastic constants, 𝐾₁, 𝐾₂, and 𝐾₃ are the phason elastic constants, and 𝑅 is phonon-phason coupling elastic constant.

Substituting (2.2) into (1.3), we obtain the equations of motion of octagonal quasicrystals given by displacement components as follows:𝜕²𝑢_𝑥𝜕𝑡²=𝑐²¹𝜕²𝑢_𝑥𝜕𝑥²+𝑐²¹−𝑐²²𝜕²𝑢_𝑦𝜕𝑥𝜕𝑦+𝑐²²𝜕²𝑢_𝑥𝜕𝑦²+𝑐²³𝜕²𝑤_𝑥𝜕𝑥²+2𝜕²𝑤_𝑦𝜕𝑥𝜕𝑦−𝜕²𝑤_𝑥𝜕₂𝜕𝑦²,𝑢_𝑦𝜕𝑡²=𝑐²²𝜕²𝑢_𝑦𝜕𝑥²+𝑐²¹−𝑐²²𝜕²𝑢_𝑥𝜕𝑥𝜕𝑦+𝑐²¹𝜕²𝑢_𝑦𝜕𝑦²+𝑐²³𝜕²𝑤_𝑦𝜕𝑥²−2𝜕²𝑤_𝑥𝜕𝑥𝜕𝑦−𝜕²𝑤_𝑦𝜕𝑤_𝑥𝜕𝑦²,𝜕𝑡=𝑑²¹𝜕²𝑤_𝑥𝜕𝑥²+𝑑²²𝜕²𝑤_𝑥𝜕𝑦²+𝑑²²−𝑑²¹𝜕²𝑤_𝑦𝜕𝑥𝜕𝑦+𝑑²³𝜕²𝑢_𝑥𝜕𝑥²−2𝜕²𝑢_𝑦𝜕𝑥𝜕𝑦−𝜕²𝑢_𝑥𝜕𝑤_𝑦𝜕𝑦²,𝜕𝑡=𝑑²¹𝜕²𝑤_𝑦𝜕𝑥²+𝑑²²𝜕²𝑤_𝑦𝜕𝑦²+𝑑²²−𝑑²¹𝜕²𝑤_𝑥𝜕𝑥𝜕𝑦+𝑑²³𝜕²𝑢_𝑦𝜕𝑥²+2𝜕²𝑢_𝑥𝜕𝑥𝜕𝑦−𝜕²𝑢_𝑦𝜕𝑦²,(2.3) in which𝑐₁=𝐿+2𝑀𝜌,𝑐²=𝑀𝜌,𝑐³=𝑅𝑑₁𝜌,=𝐾₁𝜅,𝑑²=𝐾₁+𝐾₂+𝐾₃𝜅,𝑑³=𝑅𝜅,(2.4) note that constants 𝑐₁, 𝑐₂, and 𝑐₃ have the meaning of elastic wave speeds, while 𝑑²¹, 𝑑²², and 𝑑²³ do not represent wave speed; they are diffusive coefficients.

To ensure the uniform configuration along the periodic direction physically and geometrically, we consider such model that a rectangular octagonal quasicrystal containing a Griffith crack in the center is subjected a dynamic tensile stress on both upper and lower boundaries. Assume that the Griffith crack penetrates through the solid along the periodic direction. As shown in Figure 1(a), the size of the specimen is that ED equals to 2𝐿 and EF equals to 2𝐻, and a is the length of the crack, which is constant, and 𝑝(𝑡) is the dynamic load, which varies with time. This represents the condition of the initiation of crack growth. As shown in Figure 1(b), the conditions are almost the same to that in Figure 1(a) except that the length of the crack 𝑎(𝑡) is a function of time instead of constant and the dynamic load 𝑝(𝑡)becomes constant.

(a)

(b)

Figure 1 

Sample of cracked quasicrystal under tension: (a) Crack stable with a dynamic load, (b) Crack fast propagation with a constant load.

As shown in Figures 1(a) and 1(b), only the upper right quarter needs to consider due to the symmetry. Referring to the upper right part and considering a fix grips case, the following boundary conditions should be satisfied:𝜎𝑢_𝑥=0,𝜎_𝑦𝑥=0,𝑤_𝑥=0,𝐻_𝑦𝑥=0on𝑥=0for0≤𝑦≤𝐻,_𝑥𝑥𝜎=0,𝜎_𝑦𝑥=0,𝐻_𝑥𝑥=0,𝐻_𝑦𝑥=0on𝑥=𝐿for0≤𝑦≤𝐻,_𝑦𝑦𝜎=𝑝(𝑡),𝜎_𝑥𝑦=0,𝐻_𝑦𝑦=0,𝐻_𝑥𝑦=0on𝑦=𝐻for0≤𝑥≤𝐿,_𝑦𝑦𝑢=0,𝜎_𝑥𝑦=0,𝐻_𝑦𝑦=0,𝐻_𝑥𝑦=0on𝑦=0for0≤𝑥≤𝑎(𝑡),_𝑦=0,𝜎_𝑥𝑦=0,𝑤_𝑦=0,𝐻_𝑥𝑦=0on𝑦=0for𝑎(𝑡)≤𝑥≤𝐿,(2.5) in which 𝑝(𝑡) is the dynamic load.

In addition, the initial conditions including the initial displacement and velocity are equal to zeroes:𝑤𝑢_𝑥(𝑥,𝑦,𝑡)∣_𝑡=0=0,𝑢_𝑦(𝑥,𝑦,𝑡)∣_𝑡=0=0,_𝑥𝜕𝑢_𝑥(𝑥,𝑦,𝑡)∣_𝑡=0=0,𝑤_𝑦(𝑥,𝑦,𝑡)∣_𝑡=0=0,𝜕𝑡(𝑥,𝑦,𝑡)∣^𝑡=0=0,𝜕𝑢^𝑦𝜕𝑡(𝑥,𝑦,𝑡)∣^𝑡=0=0.(2.6)

In the following section, we let the related parameters 2𝐻=40mm and 2𝐿=20mm, the initial length 𝑎₀=2.4mm of the crack, and elastic constants can be found in Table 1.

Table 1 

The related parameters used in numerical computation [24, 26, 28] (unit: 10¹² dyn/cm²).

3. Demonstrations of the Finite Difference Method and the Results

The final governing equation (2.3) along with boundary conditions equation (2.5) and initial conditions equation (2.6) are very complicated, analytic solution for the boundary-initial value problem is not available at present, which has to be solved by numerical method. Here we extend the method of finite difference of Shmuely and Alterman [22] scheme for analyzing crack problem for conventional engineering materials to quasicrystalline materials. A grid is imposed on the upper right of the specimen shown in Figure 2, the grid is extended beyond the half step.

Figure 2 

Finite difference scheme quarter specimen.

Taking the finite difference scheme proposed by [22], (2.3) and the conditions equations (2.5) and (2.6) can be reduced to certain difference equations to solve. The details on numerical computation are listed in the Appendix. To check the validity of the model and the accuracy of the algorithm, at first we give a numerical example for phason displacement. Figure 3 depicts displacement component of phason field 𝑤_𝑥 versus time. According to the present model phasons represent diffusion, the computational results indicate this is true, and the component 𝑤_𝑥 is identical to the fundamental solution of pure diffusion in mathematical physics, that is,𝑤∼1√𝑡−𝑡₀𝑒^{−(𝑥−𝑥0)2/Γ𝑤(𝑡−𝑡0)},(3.1) where 𝑡 is the time 𝑡₀ is a special value of 𝑡, 𝑥 is a distance, and 𝑥₀ is a special value of 𝑥. Though there is small fluctuation around the fundamental solution due to phonon-phason coupling, this checking confirms that the model is reasonable and the finite difference method and its computer implementation are effective and of high accuracy. In addition the phonons present the character of wave propagation this is natural and easily understood, and need not to give additional discussion and graphic illustration.

Figure 3 

Displacement component of phason field 𝑤𝑥 versus time and comparison with the fundamental solution of pure diffusion equation in mathematical physics.

In the “phase” of dynamic initiation of crack growth, the specimen with stationary crack assuming 𝑎(𝑡)=constant is subjected to a rapidly varying load 𝑝(𝑡)=𝑝𝑓(𝑡), where 𝑓(𝑡) is taken as the Heaviside function. The stresses at the crack tip present singularity of order 𝑟^−1/2, where 𝑟 denotes the distance from the crack tip, so that we can obtain the dynamic stress intensity factor for phonon stress field such as𝐾_𝐼(𝑡)=lim_𝑥→𝑎+√2𝜋(𝑥−𝑎)𝜎_𝑦𝑦(𝑥,0,𝑡)(3.2) and the normalized dynamics stress intensity factor; that is, 𝐾_𝐼(𝑡)/√𝜋𝑎𝑝 (√𝜋𝑎𝑝 denotes the static stress intensity factor of a Griffith crack) versus time is depicted in Figure 4.

Figure 4 

Normalized S.I.F. for phonon stress field versus time.

In Figure 4, there are two curves in the figure: one represents quasicrystal, that is, 𝑅/𝑀=0.01, and the other describes periodic crystals corresponding to 𝑅/𝑀=0. Because of the phonon-phason coupling effect, the mechanical properties of the quasicrystals are obviously different from the classical crystals. But the phonon wave propagation plays dominating role in dynamics of quasicrystals.

The follows are the cases that the specimen with a stationary crack are subjected to two most common types of pulse dynamic load. The pulses possess different period 𝑇₀, respectively.

From Figure 5, it can be easily concluded that S.I.F. is influenced greatly by the shape of the load of the pulse imposed. Also it is interesting that there is a time of lag between the variables and the load getting extremum. The reason for that is the influence of stress wave propagation as well as the effect of refraction and reflection. It is obviously seen that the shape of these sets of waves in the initial is like sawtooth. there will be of a bit fluctuations in the process when the wave reaches the second peak value. From the above analysis it can be concluded that the curve of S.I.F. is influenced greatly by the shape of the load of the pulse imposed, and the shapes of both are similar.

Figure 5 

Normalized S.I.F. versus time for different wave length of zigzag pulse dynamic loads.

From Figure 6, we conclude that the curved line of S.I.F. has relationships with the energy of the relative load. The larger the time of the pulse continues, the more mildly the S.I.F. varies also in the same conditions S.I.F. is larger on average in the long run.

Figure 6 

Normalized S.I.F. versus time for different wave length of rectangular pulse dynamic loads.

Now we focus on the discussion for the “phase” of fast crack propagation. To explore the inertia effect caused by the fast crack propagation, the specimen are designed under the action of constant load 𝑝(𝑡)=𝑝 rather than time-varying load, but the crack grows with high speed. The problem for fast crack propagation is a nonlinear problem, because one part of the boundaries’ crack is of unknown length beforehand. For this moving boundary problem, we must give additional condition for determining solution. That is, we must give a criterion checking crack propagation or crack arrest at the growing crack tip. This criterion can be imposed in different ways, for example, the critical stress criterion or critical energy criterion. The stress criterion is used in this paper: 𝜎_𝑦𝑦<𝜎_𝑐 represents crack arrest, 𝜎_𝑦𝑦=𝜎_𝑐 represents critical state, and 𝜎_𝑦𝑦>𝜎_𝑐 represents crack propagation. Here 𝜎_𝑐=450 Mpa is adopted as this criterion. The detailed introduction is cited from [2].

The crack velocity for quasicrystals and periodic crystals is constructed in Figure 7, from which, we can see that the velocity in quasicrystals is lower than that of the periodic crystals. The simulation reveals the dominating role of highly coordinate atomic environment as structure-intrinsic obstacles for crack propagation. For explaining this phenomenon we can explain a specific crack propagation mechanism: because the phonon-phason coupling makes the quasicrystals different from periodic crystals, the crack tip propagation velocity maybe hindered by phonon-phason coupling effect. Though the term seems like that introduced by [32], the meaning of them is different.

Figure 7 

The crack velocity versus normalized crack size with different coupling constants with 𝑝=1 Mpa.

Next we will explore the velocity under different loads in quasicrystal. The above described procedure was conducted, keeping the same geometry and material constants. With various loads, the relation between velocity and crack growth is depicted in Figure 8. It is understandable, while the load is increasing, the time to reach the stress criterion is short; so that the velocity becomes fast.

Figure 8 

Variation of crack velocity versus crack growth size for different load levels.

For the profile of growing crack shown in Figure 9, presents roughness when load level is growing, but there is no experimental observation for fast crack propagation though Ebert et al. had made some observation by scanning tunneling microscopy for quasistatic crack growth [33]. Because the fast propagation and quasistatic crack growth belong to two different regimes, the comparison cannot be done.

Figure 9 

Normalized crack growth size (𝑎−𝑎0)/𝑎0 of crack tip versus time for different load levels.

In Figure 10, we can see that the fluctuation of normalized S.I.F. is very large when the criterion 𝜎_𝑐=450 Mpa, that is, because the surface of the crack is not flat.

Figure 10 

Normalized S.I.F. for quasicrystal in the criterion 𝜎𝑐 = 450 Mpa.

4. Conclusion and Discussion

The dynamic fracture process of octagonal quasicrystals occupies an important position in dynamic response of quasicrystals. Of course, many analytic solutions which are related to the static mechanics of octagonal quasicrystals have been obtained. But obtaining the analytic solution of dynamic fracture process of octagonal quasicrystals is very difficult. In this paper, we simulate dynamic fracture process of octagonal quasicrystals. Because the analytic solution about it is not available, we develop the finite difference procedure based on argument of Bak as well as argument of Lubensky et al., which can also be seen as a collaborating model of wave propagation and motion of diffusion. Numerical results reveal the validity of wave propagation behavior of phonon field and behavior of motion of diffusion of phason field; the interaction between phonons and phasons is also recorded. The dominant physical quantities (i.e., stress intensity factors) have been discussed in detail. The formulism is applied to analysis of dynamic initiation of crack growth and crack fast propagation for two-dimensional octagonal quasicrystals of point group 8 mm the displacement and stress fields around the tip of stationary and propagating cracks are revealed. The results show that the phonon-phason coupling plays an important role in dynamic fracture behaviour of octagonal quasicrystals, which is different from crystals.

Appendix

The Details of Finite Difference Scheme

Here we extend the method of finite difference of Shmuely and Alterman scheme for analyzing crack problem for conventional engineering materials to quasicrystalline materials. A grid is imposed on the upper right of the specimen shown in Figure 2. For convenience, the mesh size ℎ is taken to be the same in both 𝑥 and 𝑦 directions. The grid is extended beyond the half step by adding four special grid lines 𝑥=−ℎ/2, 𝑥=𝐿+ℎ/2, 𝑦=−ℎ/2, and 𝑦=𝐻+ℎ/2 which form the grid boundaries. Denoting the time step by 𝜏 and using central difference approximations, the finite difference formulations for the above text are𝑢_𝑥(𝑥,𝑦,𝑡+𝜏)=+2𝑢_𝑥𝜏(𝑥,𝑦,𝑡)−𝑢_𝑥(𝑥,𝑦,𝑡−𝜏)ℎ𝑐¹²𝑢_𝑥+𝜏(𝑥+ℎ,𝑦,𝑡)−2𝑢_𝑥(𝑥,𝑦,𝑡)+𝑢_𝑥(𝑥−ℎ,𝑦,𝑡)ℎ²𝑐²¹−𝑐²²−𝑢𝑢_𝑦𝑦(𝑥+ℎ,𝑦+ℎ,𝑡)−𝑢_𝑦(𝑥+ℎ,𝑦−ℎ,𝑡)+𝜏(𝑥−ℎ,𝑦+ℎ,𝑡)+𝑢_𝑦(𝑥−ℎ,𝑦−ℎ,𝑡)ℎ𝑐²²𝑢_𝑥+𝜏(𝑥,𝑦+ℎ,𝑡)−2𝑢_𝑥(𝑥,𝑦,𝑡)+𝑢_𝑥(𝑥,𝑦−ℎ,𝑡)ℎ𝑐³²𝑤_𝑥+2𝜏(𝑥+ℎ,𝑦,𝑡)−2𝑤_𝑥(𝑥,𝑦,𝑡)+𝑤_𝑥(𝑥−ℎ,𝑦,𝑡)ℎ²𝑐²³+𝑤𝑤_𝑦𝑦(𝑥+ℎ,𝑦+ℎ,𝑡)−𝑤_𝑦(𝑥+ℎ,𝑦−ℎ,𝑡)−𝑤_𝑦(𝑥−ℎ,𝑦+ℎ,𝑡)−𝜏(𝑥−ℎ,𝑦−ℎ,𝑡)ℎ𝑐³²𝑤_𝑥𝑢(𝑥,𝑦+ℎ,𝑡)−2𝑤_𝑥(𝑥,𝑦,𝑡)+𝑤_𝑥(𝑥,𝑦−ℎ,𝑡),(A.1a)_𝑦(𝑥,𝑦,𝑡+𝜏)=2𝑢_𝑦(𝑥,𝑦,𝑡)−𝑢_𝑦(𝑥,𝑦,𝑡−𝜏)+𝜏ℎ𝑐²²×𝑢_𝑦(𝑥+ℎ,𝑦,𝑡)−2𝑢_𝑦(𝑥,𝑦,𝑡)+𝑢_𝑦(𝑥−ℎ,𝑦,𝑡)+𝜏2ℎ²𝑐²¹×𝑢_𝑥−𝑐²²+𝜏(𝑥+ℎ,𝑦+ℎ,𝑡)−𝑢_𝑥(𝑥+ℎ,𝑦−ℎ,𝑡)−𝑢_𝑥(𝑥−ℎ,𝑦+ℎ,𝑡)+𝑢_𝑥(𝑥−ℎ,𝑦−ℎ,𝑡)ℎ𝑐¹²𝑢_𝑦(𝑥,𝑦+ℎ,𝑡)−2𝑢_𝑦(𝑥,𝑦,𝑡)+𝑢_𝑦(𝑥,𝑦−ℎ,𝑡)+𝜏ℎ𝑐³²×𝑤_𝑦(𝑥+ℎ,𝑦,𝑡)−2𝑤_𝑦(𝑥,𝑦,𝑡)+𝑤_𝑦(𝑥−ℎ,𝑦,𝑡)−2𝜏2ℎ²𝑐²³×𝑤+_𝑥𝑤(_𝑥𝑥+ℎ,𝑦+ℎ,𝑡)−𝑤_𝑥(𝑥+ℎ,𝑦−ℎ,𝑡)−𝑤_𝑥(𝑥−ℎ,𝑦+ℎ,𝑡)(𝑥−ℎ,𝑦−ℎ,𝑡)−𝜏ℎ𝑐³²𝑤×𝑤_𝑦(𝑥,𝑦+ℎ,𝑡)−2𝑤_𝑦(𝑥,𝑦,𝑡)+𝑤_𝑦(𝑥,𝑦−ℎ,𝑡),(A.1b)_𝑥(𝑥,𝑦,𝑡+𝜏)=𝑤_𝑥(𝑥,𝑦,𝑡)−𝑑²³𝜏ℎ²+𝑑₂₁𝑢_𝑥(𝑥+ℎ,𝑦,𝑡)−2𝑢_𝑥(𝑥,𝑦,𝑡)+𝑢_𝑥(𝑥−ℎ,𝑦,𝑡)𝜏ℎ²+𝑑₂₂𝑤_𝑥(𝑥+ℎ,𝑦,𝑡)−2𝑤_𝑥(𝑥,𝑦,𝑡)+𝑤_𝑥(𝑥−ℎ,𝑦,𝑡)𝜏ℎ²−1𝑤_𝑥(𝑥,𝑦+ℎ,𝑡)−2𝑤_𝑥(𝑥,𝑦,𝑡)+𝑤_𝑥(𝑥,𝑦−ℎ,𝑡)2𝑑²³𝜏ℎ²+𝑢𝑢_𝑦𝑦(𝑥+ℎ,𝑦+ℎ,𝑡)−𝑢_𝑦(𝑥+ℎ,𝑦−ℎ,𝑡)−𝑢_𝑦(𝑥−ℎ,𝑦+ℎ,𝑡)(𝑥−ℎ,𝑦−ℎ,𝑡)+𝑑²²−𝑑²¹𝜏(2ℎ)²×𝑤+_𝑦𝑤(𝑥_𝑦+ℎ,𝑦+ℎ,𝑡)−𝑤_𝑦(𝑥+ℎ,𝑦−ℎ,𝑡)−𝑤_𝑦(𝑥−ℎ,𝑦+ℎ,𝑡)−𝑑₂₃(𝑥−ℎ,𝑦−ℎ,𝑡)𝜏ℎ²𝑤𝑢_𝑥(𝑥,𝑦+ℎ,𝑡)−2𝑢_𝑥(𝑥,𝑦,𝑡)+𝑢_𝑥(𝑥,𝑦−ℎ,𝑡),(A.1c)_𝑦(𝑥,𝑦,𝑡+𝜏)=𝑤_𝑦(𝑥,𝑦,𝑡)−𝑑²³𝜏ℎ²+𝑑₂₁𝑢_𝑦(𝑥+ℎ,𝑦,𝑡)−2𝑢_𝑦(𝑥,𝑦,𝑡)+𝑢_𝑦(𝑥−ℎ,𝑦,𝑡)𝜏ℎ²+𝑑₂₂𝑤_𝑦(𝑥+ℎ,𝑦,𝑡)−2𝑤_𝑦(𝑥,𝑦,𝑡)+𝑤_𝑦(𝑥−ℎ,𝑦,𝑡)𝜏ℎ²−1𝑤_𝑦(𝑥,𝑦+ℎ,𝑡)−2𝑤_𝑦(𝑥,𝑦,𝑡)+𝑤_𝑦(𝑥,𝑦−ℎ,𝑡)2𝑑²³𝜏ℎ²+𝑢𝑢_𝑥𝑥(𝑥+ℎ,𝑦+ℎ,𝑡)−𝑢_𝑥(𝑥+ℎ,𝑦−ℎ,𝑡)−𝑢_𝑥(𝑥−ℎ,𝑦+ℎ,𝑡)(𝑥−ℎ,𝑦−ℎ,𝑡)+𝑑²²−𝑑²¹𝜏(2ℎ)²×𝑤+_𝑥𝑤(_𝑥𝑥+ℎ,𝑦+ℎ,𝑡)−𝑤_𝑥(𝑥+ℎ,𝑦−ℎ,𝑡)−𝑤_𝑥(𝑥−ℎ,𝑦+ℎ,𝑡)−𝑑₂₃(𝑥−ℎ,𝑦−ℎ,𝑡)𝜏ℎ²𝑢_𝑦(𝑥,𝑦+ℎ,𝑡)−2𝑢_𝑦(𝑥,𝑦,𝑡)+𝑢_𝑦(𝑥,𝑦−ℎ,𝑡).(A.1d)

The displacements at mesh points located at the special lines are determined by satisfying the boundary conditions; we obtain, respectively, for points on the grid lines 𝑥=−ℎ/2 and 𝑥=𝐿+ℎ/2𝑢_𝑥⎛⎜⎜⎝−ℎ𝐿+2ℎ2,𝑦,𝑡⎞⎟⎟⎠=𝑢^𝑥⎛⎜⎜⎝ℎ𝐿−2ℎ2,𝑦,𝑡⎞⎟⎟⎠±12𝑑²¹𝑐²¹−2𝑐²²+𝑐²³𝑑²³𝑐²¹𝑑²¹−𝑐²³𝑑²³×⎡⎢⎢⎣𝑢_𝑦⎛⎜⎜⎝ℎ𝐿−2ℎ2,𝑦+ℎ,𝑡⎞⎟⎟⎠−𝑢^𝑦⎛⎜⎜⎝ℎ𝐿−2ℎ∓12,𝑦−ℎ,𝑡⎞⎟⎟⎠⎤⎥⎥⎦2𝑐²³𝑑²²−2𝑑²¹−𝑑²⁴𝑐²¹𝑑²¹−𝑐²³𝑑²³×⎡⎢⎢⎣𝑤_𝑦⎛⎜⎜⎝ℎ𝐿−2ℎ2,𝑦+ℎ,𝑡⎞⎟⎟⎠−𝑤^𝑦⎛⎜⎜⎝ℎ𝐿−2ℎ𝑤2,𝑦−ℎ,𝑡⎞⎟⎟⎠⎤⎥⎥⎦,(A.2a)_𝑥⎛⎜⎜⎝−ℎ𝐿+2ℎ2,𝑦,𝑡⎞⎟⎟⎠=𝑤^𝑥⎛⎜⎜⎝ℎ𝐿−2ℎ2,𝑦,𝑡⎞⎟⎟⎠±𝑑²³𝑐²²−𝑐²¹𝑐²¹𝑑²¹−𝑐²³𝑑²³×⎡⎢⎢⎣𝑢_𝑦⎛⎜⎜⎝ℎ𝐿−2ℎ2,𝑦+ℎ,𝑡⎞⎟⎟⎠−𝑢^𝑦⎛⎜⎜⎝ℎ𝐿−2ℎ±12,𝑦−ℎ,𝑡⎞⎟⎟⎠⎤⎥⎥⎦2𝑐²¹𝑑²²−𝑑²¹−𝑑²⁴−𝑐²³𝑑²³𝑐²¹𝑑²¹−𝑐²³𝑑²³×⎡⎢⎢⎣𝑤_𝑦⎛⎜⎜⎝ℎ𝐿−2ℎ2,𝑦+ℎ,𝑡⎞⎟⎟⎠−𝑤^𝑦⎛⎜⎜⎝ℎ𝐿−2ℎ𝑢2,𝑦−ℎ,𝑡⎞⎟⎟⎠⎤⎥⎥⎦,(A.2b)_𝑦⎛⎜⎜⎝−ℎ𝐿+2ℎ2,𝑦,𝑡⎞⎟⎟⎠=𝑢^𝑦⎛⎜⎜⎝ℎ𝐿−2ℎ±12,𝑦,𝑡⎞⎟⎟⎠2⎡⎢⎢⎣𝑢^𝑥⎛⎜⎜⎝ℎ𝐿−2ℎ2,𝑦+ℎ,𝑡⎞⎟⎟⎠−𝑢^𝑥⎛⎜⎜⎝ℎ𝐿−2ℎ∓12,𝑦−ℎ,𝑡⎞⎟⎟⎠⎤⎥⎥⎦2𝑐²³𝑑²²+𝑑²⁴𝑐²²𝑑²²−𝑐²³𝑑²³⎡⎢⎢⎣𝑤_𝑥⎛⎜⎜⎝ℎ𝐿−2ℎ2,𝑦+ℎ,𝑡⎞⎟⎟⎠−𝑤^𝑥⎛⎜⎜⎝ℎ𝐿−2ℎ𝑤2,𝑦−ℎ,𝑡⎞⎟⎟⎠⎤⎥⎥⎦,(A.2c)_𝑦⎛⎜⎜⎝−ℎ𝐿+2ℎ2,𝑦,𝑡⎞⎟⎟⎠=𝑤^𝑦⎛⎜⎜⎝ℎ𝐿−2ℎ2,𝑦,𝑡⎞⎟⎟⎠±12𝑐²³𝑑²³+𝑐²²𝑑²⁴𝑐²²𝑑²²−𝑐²³𝑑²³×⎡⎢⎢⎣𝑤_𝑥⎛⎜⎜⎝ℎ𝐿−2ℎ2,𝑦+ℎ,𝑡⎞⎟⎟⎠−𝑤^𝑥⎛⎜⎜⎝ℎ𝐿−2ℎ2,𝑦−ℎ,𝑡⎞⎟⎟⎠⎤⎥⎥⎦,(A.2d)where (A.2a) and (A.2b) related to 𝑥=−ℎ/2 are not valid. From the first condition of (2.5), at 𝑥=0, 𝑢_𝑥=0 and 𝑤_𝑥=0. To satisfy the condition the displacements 𝑢_𝑥 and 𝑤_𝑥 at 𝑥=−ℎ/2 are approximated by𝑢_𝑥𝑥,−ℎ𝑤2,𝑡=−𝑢^𝑥𝑥,ℎ2,𝑡_𝑥,𝑥,−ℎ2,𝑡=−𝑤^𝑥𝑥,ℎ2,𝑡.(A.3) On the grid line 𝑦=−ℎ/2 and 𝑦=𝐻+ℎ/2, we obtain𝑢_𝑥⎛⎜⎜⎝𝑥,−ℎ𝐿+2ℎ2,𝑡⎞⎟⎟⎠=𝑢^𝑥⎛⎜⎜⎝𝑥,−ℎ𝐿−2ℎ±12,𝑡⎞⎟⎟⎠2⎡⎢⎢⎣𝑢^𝑦⎛⎜⎜⎝𝑥+ℎ,ℎ𝐿−2ℎ2,𝑡⎞⎟⎟⎠−𝑢^𝑦⎛⎜⎜⎝𝑥−ℎ,ℎ𝐿−2ℎ±12,𝑡⎞⎟⎟⎠⎤⎥⎥⎦2𝑐²³𝑑²¹−𝑑²³𝑐²²𝑑²¹−𝑐²³𝑑²²⎡⎢⎢⎣𝑤_𝑦⎛⎜⎜⎝𝑥+ℎ,ℎ𝐿−2ℎ2,𝑡⎞⎟⎟⎠−𝑤^𝑦⎛⎜⎜⎝𝑥−ℎ,ℎ𝐿−2ℎ𝑤2,𝑡⎞⎟⎟⎠⎤⎥⎥⎦,(A.4a)_𝑥⎛⎜⎜⎝𝑥,−ℎ𝐿+2ℎ2,𝑡⎞⎟⎟⎠=𝑤^𝑥⎛⎜⎜⎝𝑥,−ℎ𝐿−2ℎ±12,𝑡⎞⎟⎟⎠2𝑐²³𝑑²²−𝑐²²𝑑²³𝑐²²𝑑²¹−𝑐²³𝑑²²⎡⎢⎢⎣𝑤_𝑦⎛⎜⎜⎝𝑥+ℎ,ℎ𝐿−2ℎ2,𝑡⎞⎟⎟⎠−𝑤^𝑦⎛⎜⎜⎝𝑥−ℎ,ℎ𝐿−2ℎ𝑢2,𝑡⎞⎟⎟⎠⎤⎥⎥⎦,(A.4b)_𝑦⎛⎜⎜⎝𝑥,−ℎ𝐿+2ℎ2,𝑡⎞⎟⎟⎠=𝑢^𝑦⎛⎜⎜⎝𝑥,ℎ𝐿−2ℎ2,𝑡⎞⎟⎟⎠±12𝑐²³𝑑²²+𝑑²¹𝑐²¹−2𝑐²²𝑐²¹𝑑²¹−𝑐²³𝑑²²×⎡⎢⎢⎣𝑢_𝑥⎛⎜⎜⎝𝑥+ℎ,ℎ𝐿−2ℎ2,𝑡⎞⎟⎟⎠−𝑢^𝑥⎛⎜⎜⎝𝑥−ℎ,ℎ𝐿−2ℎ±12,𝑡⎞⎟⎟⎠⎤⎥⎥⎦2𝑐²³𝑑²³−𝑑²¹𝑐²¹𝑑²¹−𝑐²³𝑑²²⎡⎢⎢⎣𝑤_𝑥⎛⎜⎜⎝𝑥+ℎ,ℎ𝐿−2ℎ2,𝑡⎞⎟⎟⎠−𝑤^𝑥⎛⎜⎜⎝𝑥−ℎ,ℎ𝐿−2ℎ𝑤2,𝑡⎞⎟⎟⎠⎤⎥⎥⎦,(A.4c)_𝑦⎛⎜⎜⎝𝑥,−ℎ𝐿+2ℎ2,𝑡⎞⎟⎟⎠=𝑤^𝑦⎛⎜⎜⎝𝑥,ℎ𝐿−2ℎ±𝑑2,𝑡₂₂⎞⎟⎟⎠𝑐²¹−𝑐²²𝑐²¹𝑑²¹−𝑐²³𝑑²²⎡⎢⎢⎣𝑢_𝑥⎛⎜⎜⎝𝑥+ℎ,ℎ𝐿−2ℎ2,𝑡⎞⎟⎟⎠−𝑢^𝑥⎛⎜⎜⎝𝑥−ℎ,ℎ𝐿−2ℎ±12,𝑡⎞⎟⎟⎠⎤⎥⎥⎦2𝑐²¹𝑑²³−𝑐²³𝑑²²𝑐²¹𝑑²¹−𝑐²³𝑑²²⎡⎢⎢⎣𝑤_𝑥⎛⎜⎜⎝𝑥+ℎ,ℎ𝐿−2ℎ2,𝑡⎞⎟⎟⎠−𝑤^𝑥⎛⎜⎜⎝𝑥−ℎ,ℎ𝐿−2ℎ2,𝑡⎞⎟⎟⎠⎤⎥⎥⎦,(A.4d)in which, (A.4c) and (A.4d) related to 𝑦=−ℎ/2 is valid only along the crack surface, namely, only for 𝑥≤𝑎−ℎ/2 at 𝑦=0, in which the crack terminates. From the last condition of (2.5), at 𝑦=0 and the ahead of the crack, 𝑢_𝑦=0,𝑤_𝑦=0. To satisfy this condition the displacements 𝑢_𝑦 and 𝑤_𝑦 at 𝑦=−ℎ/2 are approximated by𝑢_𝑦𝑥,−ℎ𝑤2,𝑡=−𝑢^𝑦𝑥,ℎ2,𝑡_𝑦,𝑥,−ℎ2,𝑡=−𝑤^𝑦𝑥,ℎ2,𝑡.(A.5)

In constructing the approximation ((A.2a)-(A.2b)–(A.5)) we follow a method proposed by Alterman and Rotenberg. According to this method, derivatives perpendicular to the boundary are proposed by uncentered differences and derivatives parallel to the boundary by centered difference. The real boundary can be considered as located at a distance of half the mesh size from the grid boundaries. The four grid corners require a special treatment. Difference methods of handling the discontinuities at such points have been proposed in the past. Here we found that satisfactory results are obtained when the displacements sought are extrapolated from those given along both sides of the corner in question. Accordingly, the components 𝑢_𝑥, 𝑢_𝑦, 𝑤_𝑥, 𝑤_𝑦 at (−ℎ/2,−ℎ/2) are given by𝑢_𝑥𝑢_𝑦−2ℎ,−ℎ2,𝑡=𝑢^𝑥𝑢_𝑦ℎ2,−2ℎ,𝑡+𝑢^𝑥𝑢_𝑦−0.5⎡⎢⎢⎣𝑢_𝑥−2ℎ,ℎ2,𝑡𝑢_𝑦32ℎ,−2ℎ,𝑡𝑢^𝑥𝑢_𝑦𝑤−2ℎ,32ℎ,𝑡_𝑥⎤⎥⎥⎦,𝑤_𝑦−2ℎ,−2ℎ,𝑡=𝑤^𝑥𝑤_𝑦ℎ2,−2ℎ,𝑡+𝑤^𝑥𝑤_𝑦−0.5⎡⎢⎢⎣𝑤_𝑥−2ℎ,ℎ2,𝑡𝑤_𝑦32ℎ,−2ℎ,𝑡+𝑤^𝑥𝑤_𝑦−2ℎ,32ℎ,𝑡⎤⎥⎥⎦.(A.6)

Similar expressions are used for deriving the displacement components at (−ℎ/2,𝐻+ℎ/2), (𝐿+ℎ/2,𝐿+ℎ/2), and (𝐿+ℎ/2,−ℎ/2). By following relevant stability criterion of the scheme, the computation is always stable and achieves high exactness. The detailed introduction can be seen in [2].

Acknowledgments

The authors are very grateful to the reviewers for their helpful suggestions and patient assistance. The work is supported by the National Natural Science Foundation of China (Grant no. 10672022 and 10972035).