Shock and Vibration
Volume 2018 (2018), Article ID 5894356, 8 pages
https://doi.org/10.1155/2018/5894356

Research Article

Hole Defects Affect the Dynamic Fracture Behavior of Nearby Running Cracks

R. S. Yang,1,2 C. X. Ding,¹ L. Y. Yang,¹ P. Xu,¹ and C. Chen¹

¹School of Mechanics and Civil Engineering, China University of Mining and Technology (Beijing), Beijing 100083, China
²State Key Laboratory for Geomechanics and Deep Underground Engineering (Beijing), Beijing 100083, China

Correspondence should be addressed to C. X. Ding; dingcx91@sina.com

Received 24 November 2017; Revised 8 March 2018; Accepted 12 March 2018; Published 18 April 2018

Academic Editor: Marco Alfano

Copyright © 2018 R. S. Yang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Effects of defects on the dynamic fracture behavior of engineering materials cannot be neglected. Using the experimental system of digital laser dynamic caustics, the effects of defects on the dynamic fracture behavior of nearby running cracks are studied. When running cracks propagate near to defects, the crack path deflects toward the defect; the degree of deflection is greater for larger defect diameters. When the running crack propagates away from the defect, the degree of deflection gradually reduces and the original crack path is restored. The intersection between the caustic spot and the defect is the direct cause of the running crack deflection; the intersection area determines the degree of deflection. In addition, the defect locally inhibits the dynamic stress intensity factor of running cracks when they propagate toward the defect and locally promotes the dynamic stress intensity factor of running cracks when they propagate away from the defect.

1. Introduction

Defects such as cracks and holes inevitably exist in most engineering materials. Developments in fracture mechanics and material mechanics have shown that small defects have a significant impact on the mechanical properties of materials [1–3]. Material fracture follows a process of crack initiation, propagation, and arrest under external loading. Defects such as cracks and holes in materials have clear effects on this process. The fracture behavior of defective materials under dynamic loading is often more complex than that under static loading [4] and is important in many engineering applications, such as automobile collisions, gravity piling, and rock blasting. Therefore, fracture behavior of defective materials under dynamic loading is the focus of much current research.

Many experimental studies have been carried out on the dynamic mechanical behavior of materials and structures [5–13]. Zhang and Zhao [5] carried out a semidisc dynamic three-point bending test with impact loading. Their results showed that dynamic loading rate has a great influence on fracture toughness and arrest toughness of materials. Cortet et al. [6] found attractive and repulsive phenomena of two opposite-running cracks. Koivisto et al. [7] and Fender et al. [8] conducted further studies on such phenomenon and revealed the cause of these crack interactions using a simple geometric model. Yang et al. [9] noted that when a running crack traverses hole-like defects with different diameters, the stress intensity factor and propagation velocity of the running crack change proportionally.

Theoretical calculations and numerical simulations have also been studied in depth [14–20]. The hybrid displacement technique and the finite element method (FEM) were used by Soh and Yang to describe the effects of material defects on dynamic behavior of cracks [14]. Wang et al. [15] developed a distinct lattice spring model to study the stress intensity factor in crack tips when traversing a manufactured empty hole, achieving results that were in good agreement with laboratory experiments. Wu and Wong [16] adopted the numerical manifold method to investigate the effect of microcracks on the mechanical properties of granite under different loading rates. They found that the dynamic strength and deformation properties of granite improved with increasing loading rate.

The dynamic fracture of defective material is mainly affected by three factors: dynamic loading strain rate, interactions between running cracks and the primary crack, and interactions between running cracks and defects. Studies on interactions between running cracks and defects typically focus on cases where the running crack traverses the defect directly. However, the size of defects in actual engineering materials is often small. Running cracks produced by dynamic fracture of the material generally propagate areas near the defect, rather than directly traversing it. The effects of such indirect interactions on the dynamic properties of materials are worthy of further study. Therefore, in the present study, a laboratory experiment combined with dynamic three-point bending and the optical testing method of caustics is performed. The effects of defects with different diameters on the dynamic fracture behavior of nearby running cracks are studied. A numerical simulation using extended FEM is conducted, which provides results that are in good agreement with the experimental results.

2. Experimental Setup

2.1. Experimental Principle and System

Materials deform when subjected to stress, becoming thinner under tensile stress and thicker under compressive stress. For transparent materials without external stress, vertical incident light is transmitted directly, leading to a uniform bright area on the reference plane. Running cracks occurs in structures under certain loading conditions. Under these conditions, the optical properties of the stressed area at the crack tip change, causing vertical incident light to refract. The corresponding area on the reference plane therefore becomes a shaded area. This is the basic principle of caustics formation. The shaded area on the reference plane is called the caustic spot, and the boundary of the caustic spot represents the caustics.

Figure 1 shows a diagram of the caustics at a mode Ι running crack tip under impact loading. The maximum diameter of caustic spot is , which can be used to characterize the stress field of the crack tip, to a certain extent. The caustics equation is

where and are the polar radius and polar angle of an arbitrary point near the crack tip in a polar coordinate system, ; is the distance between the specimen plane and the reference plane; is the optical stress constant of specimen material.

Figure 1: Diagram of caustics at a mode Ι crack tip.

Based on the principle of caustics, an experimental system of digital laser dynamic caustics is designed and set up. The beam path of the system is shown in Figure 2. In the experimental system, the green laser beam from a pump laser device (DWGD; 300–1500-mW; Beijing Laserwave Optoelectronics Technology Co., Ltd.) passes through a beam expander. The divergent beam is converged into a parallel beam through field lens 1 and then focused through field lens 2 onto a high-speed camera (Fastcam-SA5). The impact loading experiments are carried out in the parallel beam between the two field lenses. Specimen fracture and crack propagation are recorded by the high-speed camera and then exported to a computer for later processing. The laser device emits a green laser with excellent optical coherence. The maximum shooting speed of the high-speed camera is 10⁶ fps. In this experiment, the impact loading is caused by free fall of a 1.2 kg drop hammer released from a height of 0.5 m. Shooting speed of the high-speed camera is set to 10⁵ fps; thus, the interval between two successive images is 10 μs.

Figure 2: Beam path of the experimental system of digital laser dynamic caustics.

2.2. Specimen Specifications and Parameters

The specimen material is polymethyl methacrylate (PMMA), 200 mm × 50 mm × 5 mm. The dynamic elasticity modulus and dynamic Poisson’s ratio of PMMA are 6.1 GPa and 0.31, respectively. A diagram of the specimen after laser cutting is shown in Figure 3. A vertical 5 mm precrack is processed at the midpoint on the bottom edge of the sample; the upper endpoint of the precrack is denoted as O. A hole with radius is cut in center of the specimen to simulate a defect. The eccentric distance between the left edge of the defect and the axis of the specimen is 3 mm. The vertical distance between O and the center of the defect is 20 mm. The impact position of the drop hammer is midpoint P of the upper edge of the sample. The distance between the two supports at bottom of the specimen is 140 mm.

Figure 3: Diagram of the specimen.

Five sets of test specimens are designed for different values of : S1 ( mm, without defect, used as a reference specimen), S2 ( mm), S3 ( mm), S4 ( mm), and S5 ( mm). Loading conditions and other external factors are the same during all experiments.

3. Impact Fracture

3.1. Fracture Form

Under impact loading, the running crack initiates at the endpoint of the precrack O and extends toward the impact position P. The maximum vertical displacement of the running crack is  mm (the distance between O and P). The crack path of S1 is straight and the cross section is smooth. In contrast, the fracture forms of the other four specimens show large differences, as shown in Figure 4. In areas near to the defect, the running crack deflects toward the defect. The degree of deflection is greater for larger defect diameters. Figure 5 shows the relationship between deflecting displacement and vertical displacement (15 mm ≤  mm). The deflecting displacement increases gradually as the running crack approaches the defect, reaching a maximum when vertical displacement is 20 mm (at the same horizontal position as center of the defect). When the running crack moves away from the defect, the deflecting displacement gradually reduces and the original propagation path is restored. Maximum values of the deflecting displacement of the five specimens are 0 mm, 0.96 mm, 1.33 mm, 1.50 mm, and 1.71 mm.

(a) S2

(b) S3

(c) S4

(d) S5

Figure 4: Fracture forms of specimens.

Figure 5: Relationships between vertical displacement and deflecting displacement of running cracks.

3.2. Fracture Process

Under impact loading, stress concentration occurs and a caustic spot is generated at the endpoint of the precrack O. With accumulation of energy, the caustic spot becomes larger and the stress intensity factor increases. When the stress intensity factor reaches the fracture toughness of the material under the loading condition, the running crack initiates at endpoint O; the crack initiation time at O is recorded as  μs. Figure 6 shows caustic photographs of S1 and S5 during crack propagation. This figure shows that the fracture of specimens is dominated by mode Ι tensile failure.

(a) S1

(b) S5

Figure 6: Caustic photographs of running cracks during propagation.

For S1, the crack path after initiation is always straight and the fracture process is undisturbed. At  μs, the running crack propagates to the middle of the specimen. The total duration of the running crack propagation is 170 μs.

For S5,  μs, the caustic spot is tangent to the defect boundary (position A, as shown in Figure 7). At  μs, the right part of the caustic spot appears to be absorbed by the defect. Subsequently, the area of caustic spot being “absorbed” increases with further propagation of the running crack. At the same time, the crack path deflects toward the defect. At  μs, the running crack propagates to the middle of the specimen, and absorption of the caustic spot and deflection of the running crack reach a maximum (position B, Figure 7). Then, the running crack propagates away from the defect, the morphology of the caustic spot restores gradually, and the crack path straightens with continuing propagation. At  μs, the caustic spot is tangent to the defect again (position C in Figure 7) and the complete morphology is restored. Then, the running crack propagates toward the impact loading position P. The total duration of the propagation is 170 μs.

Figure 7: Diagram of running crack path near to the defect.

Thus, it can be seen that the intersection of the caustic spot and the defect is a direct cause of the deflection of the running crack, and the degree of deflection is determined by the intersection area. According to the principle of caustics, the caustic spot at the crack tip is actually an area of local stress concentration that is subjected to tensile stress. Caustics are the boundary between this caustic spot and other areas of constant stress. The stress concentration area inside the boundary is relatively balanced and stable with external constant stress area when the running crack propagates without obvious disturbances. However, when the caustic spot intersects a defect, the boundary becomes discontinuous, and the relative equilibrium between the two areas is broken. The strain energy at the crack tip is released into this discontinuous region, resulting in local stress redistribution. Macroscopically, this is exhibited as the running crack being deflected toward the defect. The larger the intersecting area, the greater the local stress redistribution, and the greater the deflection of the running crack toward the defect.

3.3. Stress Distribution

In the present study, extended FEM is adopted for numerical simulation. The maximum circumferential tensile stress is used as the criterion for crack initiation and propagation. The stress distribution around the defect during running crack propagation is analyzed. Figure 8 shows relative position of the running crack path and the defect. Figure 10 shows the change of von Mises stress with calculation time in the observation area. Because the material parameters set in the numerical simulation do not completely correspond to relevant parameters in the experiment, the running crack propagation time in the numerical simulation is not exactly the same as that in the experiment; however, crack path and fracture form are similar to those in the experiment.

Figure 8: Crack path based on the numerical simulation of extended finite element method.

At  μs, stress concentrations occur at the upper and lower boundaries of the defect; stress fields on both sides of the running crack are approximately symmetrical. At  μs, the running crack appears in the observation area, existence of the defect significantly affects stress distribution in the observation area, and stress concentration around the defect is intensified. At the same time, stress fields on both sides of the running crack become nonsymmetrical and deflect to the right of the observation area. At  μs, the running crack propagates to the same horizontal position as the lower boundary of the defect and begins to move toward the defect; stress concentrations at the upper and lower sides of the defect are obviously enhanced. At  μs, the running crack is clearly deflected toward the defect; stress field distributions on both sides of the running crack occur in a clockwise defection and stress concentrations at the upper and lower boundaries of the defect also turn clockwise. This state remains and continues to be strengthened as the running crack propagates further. At  μs, deflection of the running crack to the defect reaches a maximum and the deflection gradually recovers. At  μs, the running crack begins to expand in a straight line (Figure 9).

Figure 9: Change of von Mises stress contour plots with time .

Figure 10: Change of dynamic stress intensity factor with time .

3.4. Fracture Analysis

Dynamic stress intensity factor is an important physical quantity that well characterizes the local stress field at crack tip. The formula of mode Ι dynamic stress intensity factor [21, 22] is

where is the mode Ι dynamic stress intensity factor; is the maximum diameter of the caustic spot; is the distance between the reference plane and the specimen ( m in the experiment); is the effective thickness of the specimen ( m); is the correction factor associated with crack propagation velocity and the surface wave velocity of the specimen, in practice ; is the optical stress constant of PMMA ( m²/N); is the numerical factor (for PMMA, ). For a particular experiment, is determined by .

Figure 10 shows the relationship between the dynamic stress intensity factor and time for the five specimens. The caustic spots in S2, S3, S4, and S5 are absorbed in the vicinity of the defect and the morphologies of these caustic spots are not complete. Therefore, cannot be accurately measured, so the dynamic stress intensity factor cannot be calculated in the period of  μs. According to the trend of the dynamic stress intensity factor with time, the change process is divided into three stages: stage Ι ( μs), stage ΙΙ ( μs), and stage ΙΙΙ ( μs). The dynamic stress intensity factors of running cracks in the five specimens when initiating at endpoint O of precracks are 1.22 MN·m−3/2, 1.21 MN·m−3/2, 1.26 MN·m−3/2, 1.24 MN·m−3/2, and 1.21 MN·m−3/2. Thus, the mean fracture toughness of PMMA material under the experimental loading is about 1.23 MN·m−3/2.

The dynamic stress intensity factor of the running crack in S1 increases continuously after initiation and reaches a peak of 1.66 MN·m−3/2 at  μs. Then, the dynamic stress intensity factor gradually decays until the running crack completely penetrates the specimen. The dynamic stress intensity factors of the other four specimens show significant differences; the existence of the defect has a substantial effect on the dynamic mechanical behavior of nearby running cracks.

In stage Ι, the running crack extends toward the defect. The dynamic stress intensity factor of the running crack clearly increases after initiation under impact loading. However, the dynamic stress intensity factor begins to decrease when the running crack further propagates to the defect. The time at which the dynamic stress intensity factors of S4 and S5 start to decrease is about  μs, but the time for S2 and S3 is a little later, about  μs. In this stage, as the diameter of the defect increases, the effect of the defect on the dynamic stress intensity factor appears earlier.

In stage ΙΙΙ, the running crack propagates away from the defect. When the running crack begins to propagate away from the defect, the dynamic stress intensity factor increases again (except for S2 due to the relatively small diameter that does not produce significant effects). In addition, with increasing defect diameter, the secondary peaks in the dynamic stress intensity factor increase: in S3, S4, and S5 these are 1.58 MN·m−3/2, 1.68 MN·m−3/2, and 1.76 MN·m−3/2, respectively. Subsequently, as propagation continues, the dynamic stress intensity factor gradually decreases until the running crack reaches the edge of the specimen.

Under the experimental conditions, the dynamic stress intensity factor of running crack in stage ΙΙ cannot be obtained by calculation directly. However, it is speculated that the dynamic stress intensity factor should decrease first and then increase (as shown by the dashed line in Figure 10). This is consistent with the trend of deflection of the running crack in this stage.

In general, defects have significant effects on the dynamic mechanical behavior of nearby running cracks. The dynamic stress intensity factor of running cracks is locally inhibited when they propagate toward defects; the repressive effect occurs earlier with increasing defect diameter. Defects locally increase the dynamic stress intensity factor of the running crack when it propagates away from the defect, and the promotive effect is greater for larger defect diameters.

3.5. Error Analysis

For each specimen, more than 4 repeated experiments are conducted. The experimental phenomena are repeatable and the experimental results are reliable. The measurement accuracy of crack tip displacement is mainly restricted by the exposure time of the high-speed camera. However, the moving object continues to move forward in the limited exposure time, which affects the clarity and the identification of effective information, to a certain extent. Therefore, the image inpainting technique of MATLAB, which involves wavelets, denoising, wiener filtering, and image restoration, is used to minimize errors in data processing.

The dynamic stress intensity factor is related to , z0, and according to (2), and , , and are constants. Calculation error of is

where , , and are errors of , , and , respectively.

The incident light is not absolutely perpendicular to the plane of the , which is the main reason for the errors. Assuming the deflection angle of incident light is , then the error is

In the present work, the maximum slant angle of incident light is ; substituting into (4), we get .

is mainly caused by the indistinct pixels on the speckled caustics border. For the caustics, the speckle has  mm, corresponding to 180 pixels; that is, 0.044 mm per pixel. In the experiment, the identification error is no more than one pixel, so is less than 0.044 mm. Accordingly, we have .

is caused by machining errors in the laser cutting; the maximum value of is 0.1 mm for  mm in the present work. That is, . Substituting the above calculations into (3), we obtain the maximum calculation error of , . The calculation error is relatively small, supporting the reliability of our experimental data.

4. Discussion

In the present paper, the effect of defect diameter on the dynamic behavior of running cracks was investigated. The eccentric distance between left edge of the defect and the axis of the specimen was set to 3 mm. Experiments were also conducted with different values of the eccentric distance . When  mm, it was found that the propagation path of the running crack is straight and the defect has little effect on the propagation behavior of the running cracks nearby. However, for  mm, the attractive effect of the defect on the running crack is enhanced when the running crack propagates toward the defect. This causes the running crack to be displaced more strongly toward the defect and the running crack passes through the defect directly. After a period of energy accumulation, the running crack again initiates and propagates from the upper edge of the defect. Thus, the eccentric distance has a strong influence on the dynamic mechanical behavior of nearby running cracks and plays a decisive role in the fracture form of the specimen.

It is concluded that, under similar experimental loading conditions, the eccentric distance determines size of the intersection between the caustic spot and the defect, which in turn affects stress release and redistribution of the local stress field around the running crack tip. When is larger, the caustic spot and the defect never intersect, and the local stress field at the running crack tip is undisturbed and remains relatively stable. When is smaller, the intersection of the caustic spot and the defect is larger. This results in a large displacement of the running crack toward the defect, and the running crack then traverses the defect directly. This is only a preliminary analysis; we intend to conduct rigorous analysis of influence of the eccentric distance on the dynamic behavior of the running crack in depth in future studies.

5. Conclusion

The defect in structure has a significant effect on propagation path of nearby running cracks. When the running crack propagates near to the defect, it deflects toward the defect. The deflecting degree is greater for larger defect diameters. When the running crack propagates away from the defect, the deflection displacement gradually reduces and the original propagation path is restored. The intersection between the caustic spot and the defect is the direct cause of the crack deflection. The intersection area determines the degree of deflection. The larger the interaction area between the caustic spot and the defect, the greater the local stress redistribution, and the stronger the deflection of the running crack toward the defect.

Defects have a significant effect on the dynamic mechanical behavior of nearby running cracks. Defects locally inhibit the dynamic stress intensity factor of the running crack when it propagates toward the defect. This repressive effect occurs earlier with increasing defect diameter. The defect locally promotes the dynamic stress intensity factor of the running crack when it propagates away from the defect, and the promotive effect is stronger with increasing defect diameter. The existence of defects changes the distribution of the stress field. During deflection of the running crack toward the defect, the stress concentration around the defect is strengthened and deflected clockwise.

Finally, the present study also shows that the defect diameter and the eccentric distance are the two variables that directly determine the effect of defects on running cracks. The eccentric distance determines the propagation path shape of the running crack. When is smaller, the running crack passes through the defect directly. When is larger, propagation of the running crack is not affected by the defect. When is moderate, the running crack exhibited is deflected toward the defect. In this case, the defect diameter produces the above-mentioned effects on the propagation path and the dynamic mechanical behavior of nearby running cracks.

Nomenclature

: | Optical stress constant of specimen
: | Maximum diameter of caustic spot
: | Eccentric distance
: | Vertical distance
: | Effective thickness of specimen
: | Dynamic elasticity modulus
: | Correction factor
: | Numerical factor
: | Deflecting displacement of running crack
: | Stress intensity factor of the mode Ι crack
: | Dynamic stress intensity factor of the mode Ι crack
: | Radius of the hole defect
: | Polar radius
: | Vertical displacement of running crack
: | Propagation time of running crack
: | Distance between the specimen plane and the reference plane
: | Deflection angle of incident light
: | Error of
: | Error of
: | Error of
: | Calculation error of
: | Polar angle
: | Dynamic Poisson’s ratio.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This research is supported by the State Key Development Program for Basic Research of China (no. 2016YFC0600903) and the Natural Science Foundation of China (no. 51404273).

References

A. A. Rubinstein, “Macrocrack-microdefect interaction,” Journal of Applied Mechanics, vol. 53, no. 3, pp. 505–510, 1986.

M. F. Kanninen, C. A. Popelar, and H. Saunders, Advanced Fracture Mechanics, Oxford University Press, 1985.

M. F. Ashby and S. D. Hallam, “The failure of brittle solids containing small cracks under compressive stress states,” Acta Metallurgica et Materialia, vol. 34, no. 3, pp. 497–510, 1986.

D. Gross and T. Seelig, “Dynamic fracture mechanics,” in Fracture Mechanics, Mechanical Engineering Series, pp. 207–228, Springer, Berlin, Germany, 2011.

Q. B. Zhang and J. Zhao, “Effect of loading rate on fracture toughness and failure micromechanisms in marble,” Engineering Fracture Mechanics, vol. 102, pp. 288–309, 2013.

P.-P. Cortet, G. Huillard, L. Vanel, and S. Ciliberto, “Attractive and repulsive cracks in a heterogeneous material,” Journal of Statistical Mechanics: Theory and Experiment, vol. 2008, no. 10, Article ID P10022, 2008.

J. Koivisto, M.-J. Dalbe, M. J. Alava, and S. Santucci, “Path (un)predictability of two interacting cracks in polycarbonate sheets using Digital Image Correlation,” Scientific Reports, vol. 6, Article ID 32278, 2016.

M. L. Fender, F. Lechenault, and K. E. Daniels, “Universal shapes formed by two interacting cracks,” Physical Review Letters, vol. 105, no. 12, Article ID 125505, 2010.

R. Yang, P. Xu, Z. Yue, and C. Chen, “Dynamic fracture analysis of crack-defect interaction for mode I running crack using digital dynamic caustics method,” Engineering Fracture Mechanics, vol. 161, pp. 63–75, 2016.

I. V. Simonov and B. L. Karihaloo, “Dislocation model of an asymmetric weak zone for problems of interaction between crack-like defects,” Philosophical Magazine, vol. 85, no. 17, pp. 1847–1864, 2005.

J. Milios and G. Spathis, “Dynamic interaction of a propagating crack with a hole boundary,” Acta Mechanica, vol. 72, no. 3-4, pp. 283–295, 1988.

Z. W. Yue, Y. Guo, and X. Wang, “Influence of empty hole shape on directional fracture controlled blasting of rock,” Rock Soil Mech, vol. 37, pp. 376–382, 2016 (Chinese).

Y. Ju, L. Sudak, and H. Xie, “Study on stress wave propagation in fractured rocks with fractal joint surfaces,” International Journal of Solids and Structures, vol. 44, no. 13, pp. 4256–4271, 2007.

A. K. Soh and C. H. Yang, “Numerical modeling of interactions between a macro-crack and a cluster of micro-defects,” Engineering Fracture Mechanics, vol. 71, no. 2, pp. 193–217, 2004.

Y. Wang, R. Yang, and G. Zhao, “Influence of empty hole on crack running in PMMA plate under dynamic loading,” Polymer Testing, vol. 58, pp. 70–85, 2017.

Z. Wu and L. N. Y. Wong, “Investigating the effects of micro-defects on the dynamic properties of rock using Numerical Manifold method,” Construction and Building Materials, vol. 72, pp. 72–82, 2014.

T. G. Sitharam, J. Sridevi, and N. Shimizu, “Practical equivalent continuum characterization of jointed rock masses,” International Journal of Rock Mechanics and Mining Sciences, vol. 38, no. 3, pp. 437–448, 2001.

T. Kazerani, G. F. Zhao, and J. Zhao, “Dynamic fracturing simulation of brittle material using the distinct lattice spring method with a full rate-dependent cohesive law,” Rock Mechanics and Rock Engineering, vol. 43, no. 6, pp. 717–726, 2010.

R. T. Coates and M. Schoenberg, “Finite-difference modeling of faults and fractures,” Geophysics, vol. 60, no. 5, pp. 1514–1526, 1995.

Y. Zhao, G.-F. Zhao, and Y. Jiang, “Experimental and numerical modelling investigation on fracturing in coal under impact loads,” International Journal of Fracture, vol. 183, no. 1, pp. 63–80, 2013.

A. J. Rosakis, “Analysis of the optical method of caustics for dynamic crack propagation,” Engineering Fracture Mechanics, vol. 13, no. 2, pp. 331–347, 1980.

X. Yao, J. Chen, G. Jin, K. Arakawa, and K. Takahashi, “Caustic analysis of stress singularities in orthotropic composite materials with mode-I crack,” Composites Science and Technology, vol. 64, no. 7-8, pp. 917–924, 2004.