Research Article  Open Access
Shigang Lai, Li Shi, Alex Fok, Haiyan Li, Libin Sun, Zhengming Zhang, "A New Method to Measure Crack Extension in Nuclear Graphite Based on Digital Image Correlation", Science and Technology of Nuclear Installations, vol. 2017, Article ID 2416545, 10 pages, 2017. https://doi.org/10.1155/2017/2416545
A New Method to Measure Crack Extension in Nuclear Graphite Based on Digital Image Correlation
Abstract
Graphite components, used as moderators, reflectors, and coresupport structures in a HighTemperature GasCooled Reactor, play an important role in the safety of the reactor. Specifically, they provide channels for the fuel elements, control rods, and coolant flow. Fracture is the main failure mode for graphite, and breaching of the above channels by crack extension will seriously threaten the safety of a reactor. In this paper, a new method based on digital image correlation (DIC) is introduced for measuring crack extension in brittle materials. Crosscorrelation of the displacements measured by DIC with a step function was employed to identify the advancing crack tip in a graphite beam specimen under threepoint bending. The loadcrack extension curve, which is required for analyzing the Rcurve and tension softening behaviors, was obtained for this material. Furthermore, a sensitivity analysis of the threshold value employed for the crosscorrelation parameter in the crack identification process was conducted. Finally, the results were verified using the finite element method.
1. Introduction
In the HighTemperature GasCooled Reactor (HTGR) design, the fuel and control rods are contained and supported by graphite moderator bricks. Graphite is normally described as a brittle or quasibrittle material with low fracture toughness [1]. Under irradiation, the material properties and dimensions of graphite change dramatically with the increase of dose [2, 3], which will create stresses high enough to cause brittle failure of the graphite components. As a result, the fuel and control rod channels may be breached, jeopardizing the operation and safety of HTGR. Therefore, it is important to consider the fracture behaviors of nuclear graphite at both the design and operation stages.
The fracture behaviors of brittle materials, including ceramic, concrete, rocks, and graphite, have been widely investigated [4–8]. Various failure criteria for nuclear graphite have been proposed and verified by experiments [9–11]. As the nonlinear fracture behavior of nuclear graphite is closely related to its grain size [1, 8], many researchers have focused on the studies of crack initiation, propagation, and failure of graphite materials with different grain sizes. Su et al. [12] carried out threepoint bend tests with singleedgenotched beams (SENB) to measure the tension softening curves and Fracture Process Zone (FPZ) of a finegrain graphite. Similar tests were adopted for a coarsegrain graphite, NBG18 [13]. Important fracture parameters, including the fracture toughness, FPZ, size effect, and Rcurve, which can be used for evaluating crack resistance in quasibrittle materials were obtained. The compact tension [14] and double torsion tests [15], which are able to produce relatively stable cracks, have also been used to investigate the crack extension of graphite. In such studies, one primary requirement was to detect and monitor crack initiation and propagation. However, the failure of graphite is usually associated with fast unstable crack initiation and propagation, making it difficult to capture the process [16]. To date, many techniques have been employed to identify the crack position and measure its extension in situ during loading. Although the loaddisplacement curve can be used to determine the crack extension using the compliance of the cracked specimen [13], the compliance of the equipment and deviation of the material from linearelastic behavior could introduce errors to the measurement. Another method for detecting crack extension is to measure the change in electric resistance between two electrodes placed on the opposite sides of the anticipated crack path [17]. When the crack propagates, the resistivity of the material will change accordingly. However, the results are sensitive to the environment since the resistivity of the sample is also a function of temperature and strain.
Consequently, noncontact imaging techniques have been developed for crack measurement. For example, Electronic Speckle Pattern Interferometry (ESPI) and Xray microcomputed tomography (microCT) have been used to monitor crack growth in nuclear graphite [12, 15, 18]. ESPI can detect very small deformations on an object’s surface but the measurement is affected by external disturbances such as mechanical vibrations. With microCT, the internal microcracks developed around the main crack were observed and the actual shape of the crack front was clearly revealed [18]. However, this method needed specially designed test equipment that could be fitted into the Xray machine, and the scanning after each crack increment took a considerable amount of time to perform.
Digital image correlation (DIC) is an optical method which can also be used to measure the deformation and strain on an object’s surface [19]. By analyzing a series of images recorded through a chargecoupled device (CCD) camera the whole displacement and strain fields can be calculated, from which the cracks on the material’s surface can be identified. Compared to other methods, DIC has many advantages. For example, it is easy to set up, its measurement is in situ, fullfield, and noncontacting, and it can be used to detect large deformation [20]. Moreover, it can be used for measurement at high temperatures, which is often required for studying nuclear graphite. Not surprisingly, DIC has been widely used as an alternative method to monitor in situ the crack growth within structural materials, including nuclear graphite and functionally graded materials [21–24].
Although DIC has been shown to be a powerful tool for measuring crack extension, the accuracy of its results is uncertain. In order to identify the crack tip, a commonly used method is to define a strain threshold, , and the areas with strains higher than are considered to be part of the crack [13–15]. The opening strain, which is an artificial strain caused by the crack opening displacement, was employed by Mostafavi et al. to visualize the cracks of different types of graphite [14]. A strain threshold of 1% was applied to the map of the maximum principal strain in their study. However, the position of the crack tip thus determined is rather sensitive to the value chosen for and is greatly dependent on the resolution of the DIC results [19]. Additionally, methods using strains are not capable of identifying cracks on small specimens. Therefore, it is necessary to develop a less sensitive method for DIC that can be used to accurately identify the crack tip. In this study, the crosscorrelation [25] algorithm was employed for the first time to process the displacement fields obtained from DIC to identify the crack tip position in the specimens. A sensitivity analysis was conducted to test the convergence of this method. The results were verified through linearelastic finite element (FE) simulations and good agreement between theory and experiment was achieved.
2. Experiments
Threepointbend tests had previously been conducted using graphite singleedgenotched beams (SENB) and the fullfield distributions of displacement and strain were obtained through DIC [13]. The results for one group of these specimens were analyzed using the crosscorrelation technique presented in this study.
2.1. Material and Specimens
The test material was a pitchcoke, vibrationally molded, near isotropic nuclear grade graphite NBG18 (SGL Carbon Group, Germany). Its main material properties are listed in Table 1.

The dimensions of the SENB specimen analyzed in this study were 45 mm length, 10 mm width, and 5 mm thickness. A prenotch was cut into the specimen with a diamond blade (Buehler diamond wafering blade, 102 mm × 0.3 mm). The specimen had an initial crack lengthtowidth ratio of 0.4, and the load span was 40 mm.
2.2. ThreePointBend Test
The specimen was tested in threepoint bending using a universal MTS test machine (858 Mini Bionix II, MTS, US). The test setup is shown in Figure 1. A vertical displacement was applied to the specimen until it fractured completely. During loading, the notch mouth opening displacement was monitored continuously using an extensometer for feedback control purposes so that stable postpeak loaddisplacement data were obtained.
During loading, a CCD camera (Point Grey Grasshopper GRAS20S4CC, Point Grey Research, Inc.) was used to take a series of pictures of the front surface of the specimen at a frequency of 7.5 Hz (Figure 1). To ensure accuracy of the DIC results, the surface of the specimen facing the camera was sprayed with a white paint to produce more irregularly shaped speckles, which facilitated the tracking of the surface displacements. The camera started to record at the same time as the loading process began.
2.3. Digital Image Correlation
The images recorded during the loading process were then analyzed by proprietary image correlation software (DAVIS, LaVision). Wholefield displacement and strain distributions of the front surface of the specimens were obtained. Multiple interrogations were used with a window size of 64 × 64 pixels at an overlap of 50% with two passes, followed by the use of 32 × 32pixel windows at an overlap of 75% with two passes.
Figure 2 shows the distributions of the horizontal displacement and strain under two different loads from the DIC analysis. The images were taken at the postpeak stage when the crack was extending. Since the crack propagated mainly along the vertical direction, only the horizontal displacement and strain were used for crack identification.
(a) Horizontal displacement (mm), load = 18.128 N
(b) Horizontal displacement (mm), load = 3.357 N
(c) Horizontal strain, load = 18.128 N
(d) Horizontal strain, load = 3.357 N
Figure 3 shows the horizontal displacements along different horizontal paths. It can be seen that the displacements on the left side of the beam’s centerline are negative, and those on the right side are positive. Within the upper uncracked area (Path m3), the horizontal displacement is close to zero, whereas within the lower cracked area (Paths m1 and m2), there is a step change in the displacement curves at the position of the crack.
3. CrossCorrelation Analysis
In signal processing, crosscorrelation is a measure of similarity between two waveforms as a function of a varying timelag applied to one of them. It is widely used to detect the time delay between two similar signals [25]. Crosscorrelation is similar in nature to the convolution of two functions. As an example, consider two real valued functions and which are different only by an unknown shift along the axis. Crosscorrelation can be used to find how much must be shifted along the axis to make it identical to . The procedure essentially slides the function along the axis, calculating the summation of their product at each position. When the two functions match, the value of their crosscorrelation is maximized. This is because when the peaks and troughs are aligned, they both make a large positive contribution to the summation [22].
For discrete functions and , the crosscorrelation between them is defined aswhere is the discrete timelag and and are the th element of function and the th element of function , respectively.
The crosscorrelation method was used here to trace the crack in the graphite specimen under threepoint bending based on the displacement distributions from the DIC analysis. The scientific software Matlab [26] was used for programming the method.
3.1. Data Format
Each displacement field data obtained through DIC is stored as a twodimension matrix, as shown in Figures 2(a) and 2(b); the horizontal displacement at each sampling point was extracted and stored in the form of by matrix:or where component is the horizontal displacement in the th row and th column and are the total number of rows and columns, respectively, which are dependent on the resolution of the DIC analysis.
3.2. CrossCorrelation Process
Using each matrix extracted from the DIC results, a crosscorrelation procedure was conducted to accurately map out the crack and to measure its extension. The process included two steps:
Step 1 (establishing the crosscorrelation function and crack position array ). A discrete sign function is defined in (4). As shown in Figure 4, the shape of the Gfunction is similar to the distribution of horizontal displacement along a horizontal path crossed by the crack (Figure 3(b)). For the rows of horizontal displacement, the crosscorrelation function between the Gfunction and the matrix was constructed as where is the crosscorrelation function, is the row number, is the column number, and is an integer variable representing the shift. Therefore, by shifting , the closer it was in matching its step to that of a particular row of entries in , the larger the corresponding value in the crosscorrelation function is. In this way, the dissimilarity between the cracked and uncracked areas was actually enlarged with the crosscorrelation analysis. Note that the crosscorrelation function was also an by matrix.
Along the th row, reached its maximum value where the shift in matched its step with that in the horizontal displacement curve best. This coincided with the crack position. So the crack position was identified by finding the maximum value of along each row. An array was thus constructed to store the crack position (column number) in each row: where is the column number with the maximum in the th row; that is, . Figure 5 shows the maps of the matrix and the array (marked by red crosses) under the same Load 1 and Load 2 as in Figure 2.
(a)
(b)
Step 2 (identification of the crack tip). After obtaining the crack position array , the next task was to accurately identify the crack tip. To do so, an array was introduced to store the variance of the crosscorrelation function on each row:The variance function was defined by the following equation:where and is the mean value of .
With the variance array , the dissimilarity between the cracked and uncracked areas was further enlarged. The crack tip was then positioned by defining the cracked area as the area where , being a threshold of (see later). And the complete crack was fully captured by combining the crack position array and the variance array that defined the crack tip.
Figures 6(a) and 6(c) show the distributions of along the vertical direction for images under Load 1 (18.128 N) and Load 2 (3.357 N), respectively. In Figure 6(a), the value of tends towards 9 × 10^{−4} as it enters the uncracked region. Figures 6(b) and 6(d) show the delineated cracks using a threshold for of .
(a)
(b)
(c)
(d)
3.3. Crack Extension versus Load
The above analysis was repeated for all the images and the crack extensions under different load levels were obtained, as shown in Figure 7. Each point on the curve represents the results from one DIC image. The curve shows that crack propagation started just before the peak load and the crack propagated steadily thereafter with a decrease in load.
3.4. Sensitivity Analysis
In the crack identification process, a threshold, , for the variance was required to locate the crack tip. In order to study how much the results were dependent on , a sensitivity analysis was conducted by measuring the crack extension with different , as shown in Figure 8(a). It can be seen that as decreased, the calculated crack length increased. But the results converged quickly, and when was smaller than 1.5 × 10^{−3}, the curves were almost the same.
(a)
(b)
(c)
For comparison, the crack extension was also calculated using two other methods that were based on (1) the displacement field without the crosscorrelation analysis and (2) the strain fields. These two methods also required a threshold value for their respective parameters, displacement and strain, to identify the crack tip. The results are shown in Figures 8(b) and 8(c). It can be seen that their speeds of convergence were much slower than that based on the crosscorrelation of displacement functions, as shown in Figure 8(a), indicating that their results were more sensitive to the threshold.
4. Verification with the Finite Element Method
In order to verify the results obtained from the crosscorrelation process, the finite element (FE) method was used to predict the crack extension versus load for the same specimen. The dimensions of the specimen were given in Section 2.1.
Cracks with different lengths, from 0 mm to 5 mm, were considered in the FE analysis using ABAQUS. Figure 9(a) shows one of the FE models with a crack of 1 mm long ahead the notch. The mesh around the crack tip was refined. Loading and supporting components were considered in the model and hard contacts were assumed between these components and the graphite specimens. The FE model contained 6904 8node, planestrain elements (CPE8) and 20939 nodes. The material was assumed to be linearelastic, with the same parameters as those in Table 1.
(a)
(b)
(c)
The FEpredicted loaddisplacement curves (green) for different crack lengths were plotted together with one of the test results (blue) in Figure 9(b). It can be seen that with the increase of crack extension, the slopes of the green lines decreased, indicating a decrease in stiffness of the specimen. The loads corresponding to the different crack extensions were obtained from the intersection points between the experimental loaddisplacement curve and the FEsimulated loaddisplacement lines, as shown in Figure 9(b). The crack length as a function of load in the test can thus be estimated. The intersection points in the figure are listed in Table 2. The procedure is similar to that used by Hodgkins et al. to observe the microstructural damage and crack morphology in graphite [27]. In their study, the specimen was unloaded and reloaded during crack propagation and the crack length was calculated from the unloading compliance. In the present study, FE simulations were used instead to obtain the compliance of the cracked specimen.

Figure 9(c) shows the load versus crack extension behavior predicted by the FE analysis, together with those from the crosscorrelation analysis. It can be seen that the two sets of results agreed well with each other. The results given by the crosscorrelation method were slightly lower than those from the FE simulation.
5. Discussion and Conclusion
Digital image correlation (DIC) has been widely used to monitor the failure of engineering materials because of its capability to provide in situ, fullfield, and noncontacting deformation measurement.
In a DIC analysis, the displacement field is the primary result obtained directly from the correlation analysis of the raw images, while the strain field is a secondary result which is obtained through the differentiation of the displacement field. Therefore, the displacement field provides more detailed information on the crack profile than the strain field. However, unlike strain, the values of displacement could not be directly related to the crack. Therefore, a crosscorrelation method was developed in this study to more accurately identify the tip of an advancing crack and measure its extension at different stages of loading from the displacement field.
The crack length identified is dependent on the threshold set for the parameter used in the identification, which may in turn depend on the material properties, specimen size, and loading conditions. This is reflected in the different strain thresholds, , employed in the literature. For the threepoint bend test [13], as fracture occurred rapidly, the detected crack length was found to be very sensitive to . To find the proper , crack extension curves for different threshold values were compared and the one which could describe the experimental results best was adopted. Strain thresholds of 0.4% and 0.7% were chosen for large and smallsize specimens, respectively. Compared to the SENB specimen, compact tension and double torsion specimens were able to produce more stable crack growth and the measured crack length based on the strain map was less sensitive to . However, the selection of the threshold strain is still based on experience or trial and error, depending on the application. A threshold value of 1% was used to calculate crack length and 0.5% was chosen to identify microcracking in [14]. In a later study using the double torsion technique [15], the threshold strain for identifying microcracking was 0.35%, which was about twice the maximum predicted strain of the material assuming linear elasticity.
Similar to methods that are based on strain maps, a threshold for the variance is also required in the crosscorrelation method in order to locate the crack tip. However, the sensitivity analysis showed that the results with the crosscorrelation analysis were much less dependent on the threshold adopted. As shown in Figure 8, the results based on the crosscorrelation method converged fastest. In addition, the crack profile could be more accurately defined with this method, as shown in Figure 5.
The crack extensionload curve obtained with the crosscorrelation method was verified through FE analysis. As shown in Figure 9(c), the crack extensions given by the FE simulations were slightly longer than those from the crosscorrelation method. This was reasonable because the linearelastic material behavior assumed in the FE model could result in a higher stiffness of the specimen, causing an overestimation of the load for a certain crack extension. In reality, the stressstrain curve of graphite is slightly nonlinear and there is a Fracture Process Zone ahead of the crack tip which leads to inelastic behavior.
Based on the characteristics of the displacement distribution across a crack, a Gfunction (see (1)) was defined to enlarge the difference between the cracked and uncracked areas through crosscorrelation. It should be noted that the Gfunction in (1) is only suitable for mode I cracks, such as those produced under tensile, threepoint, and fourpoint bending tests. For specimens with other fracture modes, the aforementioned Gfunction should be revised according to the characteristics of the displacement distribution.
Although the crosscorrelation method presented in this paper was developed to trace the single crack, it can be extended to identify multiple cracks in specimens under complex loading. For example, crosscorrelation analyses along the  and directions can be conducted jointly to map out crack trajectories in different directions.
In conclusion, the crosscorrelation method developed in this study is a robust way to identify the path of an advancing crack based on DIC results. The main advantages of this method include the following: (1) the crack extension can be measured accurately; (2) the profile of the crack can be clearly identified; and (3) the results are less sensitive to the calculation parameters. The developed method will improve the current measurement of crack propagation of nuclear graphite and lead to a better understanding of its fracture behaviors, which is essential for the design and safe operation of HTGR.
Nomenclature
CCD:  Chargecoupled device 
DIC:  Digital image correlation 
ESPI:  Electronic Speckle Pattern Interferometry 
FE:  Finite element 
FPZ:  Fracture Process Zone 
HTGR:  HighTemperature GasCooled Reactor 
SENB:  Singleedgenotched beams 
:  Displacement field 
:  Variance array 
:  Strain 
:  Strain threshold. 
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was performed using funding received from the DOE Office of Nuclear Energy’s Nuclear Energy University Programs. The authors would like to acknowledge Minnesota Dental Research Center for Biomaterials and Biomechanics (MDRCBB) for providing the testing devices and the Minnesota Supercomputing Institute (MSI) for providing the computing services.
References
 T. D. Burchell, Ed., Carbon Materials for Advanced Technologies, Elsevier, Amsterdam, Netherlands, 1999.
 H. Matsuo, T. Saito, and Y. Sasaki, “Effects of high temperature neutron irradiation on dimension and thermal expansion coefficient of nuclear grade graphites,” Journal of Nuclear Science and Technology, vol. 18, no. 11, pp. 863–869, 1981. View at: Publisher Site  Google Scholar
 H. Matsuo and T. Saito, “Irradiation behaviors of nuclear grade graphite in commercial reactor, (I) dimensional change and thermal expansion,” Journal of Nuclear Science and Technology, vol. 22, no. 2, pp. 139–146, 1985. View at: Publisher Site  Google Scholar
 S. P. Shah, S. E. Swartz, and C. Ouyang, Fracture Mechanics of Concrete: Applications of Fracture Mechanics to Concrete, Rock and Other QuasiBrittle Materials, John Wiley & Sons, 1995.
 A. Bornhauser, K. Kromp, and R. F. Pabst, “Rcurve evaluation with ceramic materials at elevated temperatures by an energy approach using direct observation and compliance calculation of the crack length,” Journal of Materials Science, vol. 20, no. 7, pp. 2586–2596, 1985. View at: Publisher Site  Google Scholar
 M. R. Ayatollahi, M. R. M. Aliha, and M. M. Hassani, “Mixed mode brittle fracture in PMMA—an experimental study using SCB specimens,” Materials Science and Engineering A, vol. 417, no. 12, pp. 348–356, 2006. View at: Publisher Site  Google Scholar
 F. Berto and P. Lazzarin, “Recent developments in brittle and quasibrittle failure assessment of engineering materials by means of local approaches,” Materials Science and Engineering R: Reports, vol. 75, no. 1, pp. 1–48, 2014. View at: Publisher Site  Google Scholar
 A. Hodgkins, T. J. Marrow, M. R. Wootton, R. Moskovic, and P. E. J. Flewitt, “Fracture behaviour of radiolytically oxidised reactor coregraphites: a view,” Materials Science and Technology, vol. 26, no. 8, pp. 899–907, 2010. View at: Google Scholar
 S. Sato, H. Awaji, K. Kawamata, A. Kurumada, and T. Oku, “Fracture criteria of reactor graphite under multiaxial stesses,” Nuclear Engineering and Design, vol. 103, no. 3, pp. 291–300, 1987. View at: Publisher Site  Google Scholar
 M. R. Ayatollahi, F. Berto, and P. Lazzarin, “Mixed mode brittle fracture of sharp and blunt Vnotches in polycrystalline graphite,” Carbon, vol. 49, no. 7, pp. 2465–2474, 2011. View at: Publisher Site  Google Scholar
 P. Lazzarin, F. Berto, and M. R. Ayatollahi, “Brittle failure of inclined keyhole notches in isostatic graphite under inplane mixed mode loading,” Fatigue and Fracture of Engineering Materials and Structures, vol. 36, no. 9, pp. 942–955, 2013. View at: Publisher Site  Google Scholar
 R. K. L. Su, H. H. Chen, S. L. Fok et al., “Determination of the tension softening curve of nuclear graphites using the incremental displacement collocation method,” Carbon, vol. 57, pp. 65–78, 2013. View at: Publisher Site  Google Scholar
 H. Li, J. Li, G. Singh, and A. Fok, “Fracture behavior of nuclear graphite NBG18,” Carbon, vol. 60, pp. 46–56, 2013. View at: Publisher Site  Google Scholar
 M. Mostafavi, M. J. J. Schmidt, B. J. Marsden, and T. J. Marrow, “Fracture behaviour of an anisotropic polygranular graphite (PGA),” Materials Science and Engineering A, vol. 558, pp. 265–277, 2012. View at: Publisher Site  Google Scholar
 T. H. Becker, T. J. Marrow, and R. B. Tait, “Damage, crack growth and fracture characteristics of nuclear grade graphite using the Double Torsion technique,” Journal of Nuclear Materials, vol. 414, no. 1, pp. 32–43, 2011. View at: Publisher Site  Google Scholar
 R. Moskovic, P. J. Heard, P. E. J. Flewitt, and M. R. Wootton, “Overview of strength, crack propagation and fracture of nuclear reactor moderator graphite,” Nuclear Engineering and Design, vol. 263, pp. 431–442, 2013. View at: Publisher Site  Google Scholar
 B. Stalder, P. Béguelin, and H. H. Kausch, “A simple velocity gauge for measuring crack growth,” International Journal of Fracture, vol. 22, no. 2, pp. R47–R50, 1983. View at: Publisher Site  Google Scholar
 A. Hodgkins, T. J. Marrow, P. Mummery, B. Marsden, and A. Fok, “Xray tomography observation of crack propagation in nuclear graphite,” Materials Science and Technology, vol. 22, no. 9, pp. 1045–1051, 2006. View at: Publisher Site  Google Scholar
 DANTEC, “Measurement Principles—Digital Image Correlation (DIC),” 2012. View at: Google Scholar
 R. Lichtenberger and H. Schreier, “Noncontacting measurement technology for component safety assessement,” in Proceedings of the 6th International Symposium and Exhibition on Sophisticated Car Occupant Safety Systems (AIRBAG '02), Karlsruhe, Germany, 2002. View at: Google Scholar
 J. AbantoBueno and J. Lambros, “Investigation of crack growth in functionally graded materials using digital image correlation,” Engineering Fracture Mechanics, vol. 69, no. 14–16, pp. 1695–1711, 2002. View at: Publisher Site  Google Scholar
 S. J. Orfanidis, Optimum Signal Processing, McGrawHill, New York, NY, USA, 2007.
 M. R. Ayatollahi and M. Moazzami, “Digital image correlation method for calculating coefficients of Williams expansion in compact tension specimen,” Optics and Lasers in Engineering, vol. 90, pp. 26–33, 2017. View at: Publisher Site  Google Scholar
 R. Markgraaff, Relationship between the natural frequencies and fatigue life of NBG18 graphite [M.S. thesis], NorthWest University, Potchefstroom, South Africa, 2010.
 Wikipedia, “Crosscorrelation,” 2013. View at: Google Scholar
 MathWorks, “Matlab the language of Technical Computing,” 2013. View at: Google Scholar
 A. Hodgkins, J. Marrow, P. Mummery, A. Fok, and B. J. Marsden, “Crack propagation resistance and damage mechanisms in nuclear graphite,” in Fracture of Nano and Engineering Materials and Structures: Proceedings of the 16th European Conference of Fracture, Alexandroupolis, Greece, July 3–7, 2006, E. E. Gdoutos, Ed., pp. 891–892, Springer, Dordrecht, The Netherlands, 2006. View at: Google Scholar
Copyright
Copyright © 2017 Shigang Lai 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.