## Abstract

Providing a quantitative description of the whole-field stress evolution in complex structures subjected to continuous loading processes using traditional photoelastic approaches is a significant challenge because of the difficulties with fabricating complex structures, identifying the stress distribution and evolution, and unwrapping isochromatic phase maps. To overcome the challenges, we proposed a novel method to quantify the continuous whole-field stress evolution in a complex porous structure that was fabricated with 3D printing technology. The stress fringes were identified by analysing a series of continuous frames extracted from a video recording of the fringe changes and determining the valleys of the light intensity change curve over the entire loading process. The experimental data were compared with the numerical results of the complex model with identical pore geometries, physical properties, and loading conditions to evaluate the accuracy and effectiveness of the method. In principle, the applicability of the reported method for identifying and unwrapping the continuous whole-field stress is not affected by the complexity of a structure.

© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Quantifying and visualising the continuous distribution and evolution of the whole-field stress in complex pore structures are of vital significance to solving a variety of engineering problems where the stress distribution and evolution govern the deformation and failure of a structure. For instance, when exploiting underground resources, the efficient and safe exploration of oil, natural gas, coal methane, and mineral resources greatly relies on knowledge of the distribution and evolution of stress fields in the complex structures of natural rock masses [1–4]. In geological and environmental engineering, the efficient prediction of earthquake occurrence, waste geological disposal, CO_{2} geological storage, water resource conservation, and contaminant transportation and diffusion are extremely dependent on recognising the continuous distribution of the whole-field stress in the complicated pore fissure structures inside rock [5–8]. In civil and underground engineering, the stability of underground tunnels and caverns is closely associated with the dynamic evolution of a continual stress field [9, 10]. For these applications, the quantification and visualisation of the continuous distribution and dynamic evolutionary processes related to the whole-field stress in a complex pore structure are critical to understanding the mesoscopic damage, deformation, and failure mechanisms in underground engineering problems.

Various techniques have been proposed to quantitatively determine the stress distribution and evolution in structures; these can be roughly categorised into theoretical models [11], numerical simulations [12, 13], in situ monitoring [14, 15], and laboratory tests [16, 17]. However, no appropriate approaches are currently capable of accurately determining the continuous whole-field stress evolution in a complex structure that is characterised by numerous and geometrically irregular pores [18]. Numerical simulation methods can be used to determine and visualise the continuous stress evolution in complex structures. However, the calculation results of these methods strongly depend on the input of physical and mechanical properties, as well as constitutive models, thus increasing the difficulty of experimental verification [19, 20]. In situ monitoring methods such as hydraulic fracturing, the overcoming method in boreholes, and the Kaiser method can measure the in situ stresses of rock strata. However, the distribution and evolution of whole-field stress can hardly be obtained with a limited number of measurement points; such methods are not capable of determining the continuous distribution of the stress field induced by resource exploration or tunnel construction activities [21–23]. Moreover, there is controversy in the data accuracy of in situ monitoring methods because of the limitations with installing measurement devices [24, 25].

Laboratory testing is one of the primary methods that are employed to evaluate stress fields. In traditional laboratory approaches, however, the stress field is usually deduced from the deformation and failure properties of samples rather than observed directly [26, 27]. The obtained results are not a representation of the continuous whole-field stress but only represent the state of the local area [28, 29]. With regards to methods for quantifying the stress field that essentially governs the deformation and failure of a complex system, only a few experimental methods are capable of determining the field stress of solids, but it is challenging for them to quantify the distribution and evolution of the whole-field stress in complex structures [30–33]. Photoelasticity is an optical method that uses fringe contours to indicate the principal stress difference and direction. It is widely used for two-dimensional (2D) and three-dimensional (3D) stress analyses in mechanical engineering [34, 35], medical science [36, 37], rock mechanics and rock engineering [38–40], etc. Especially in recent years, the combination of 3D printing technology with computed tomography (CT) imaging technology has simplified the process of preparing complex models. Studies have shown that photosensitive resin is a feasible choice as a printing material [39]; this exhibits transient birefringence phenomena, which is a prerequisite for photoelastic testing during loading [38–40]. However, various difficulties with existing photoelastic methods need to be overcome in order to properly unwrap the photoelastic fringes and quantify the continuous stress field evolution in complex structures. For instance, for phase shift methods, a certain number of images are required, which leads to difficulty in determining the dynamic continuous evolution of stress field [41–43]. The determination accuracy is greatly dependent on a high unwrapping quality, which is extremely difficult to achieve for complex structures [44, 45]. Because the red–green–blue (RGB) photoelastic method requires one colour image to determine whole-field fringe orders [46], the determination precision is vulnerable to the impact of the measuring facilities and requires special equipment for the observable maximum fringe order to be detected [47–49]. The step-loading method and Fourier transform technique also cannot be used for the long duration of a continuous loading process [50, 51].

In this study, we proposed a novel photoelastic method to determine and visualise the continuous evolution of the whole-field stress of the complex pore structures that were fabricated with 3D printing technology. The stress contours of the complex structure models that comprise geometrically irregular pores were quantitatively determined by analysing a series of continuous frames extracted from a video recording of the fringe changes and counting the valleys of the light intensity change curve over the entire loading process. The novel algorithm was introduced by the analytic solution of the diametrical compression disk. Based on the proposed method, the continuous evolution of the whole-field stress in pore structures subjected to various levels of continual compressive loads was distinctly determined. The experimental results were compared with the numerical results of the complex model with identical pore geometries, physical properties, and loading conditions to evaluate the accuracy and effectiveness of the suggested method. The comparison indicates that the proposed method can well describe the whole-field stress evolution of complex porous structures, and its applicability for identifying and unwrapping the whole-field stress is not affected by the complexity of a structure.

## 2. Models and methods

#### 2.1 Model preparation

The pore structure tested in this study was taken from the CT results of real rocks provided by Imperial College London (http://www.imperial.ac.uk/earth-science/research/researchgroups/perm/research/pore-scale-modelling/micro-ct-images-and-networks/berea-sandstone/). Part of one image from the 400 Micro-CT scan images of Berea sandstone was selected to reconstruct our model as shown in Fig. 1. The local part had dimensions of 120 × 120 pixels and a resolution of 5.345 *µ*m. A binarization process was employed to segment the selected image for easier reconstruction of the digital model. Note that the distribution and evolution of the stress field around irregular pores were the major objectives in our work. Because of the limited camera resolution, only the local part was used for the photoelastic test. Choosing other parts would have no influence on the proposed method; however, the computing time involved in model reconstruction and stress field identification may increase. To facilitate the model reconstruction and experimental observation, the local part was magnified to 50 mm × 50 mm and the size of the image was magnified to 500 × 500 pixels. In addition, to ensure pressures uniformly loading on the top and bottom surfaces, a mask was added to shape the edges. Then, 80 replicas of the local image were input into the software Mimics^{®} to create the 3D model with a certain thickness. This was done because it is unsuitable for the software to reconstruct a thickened model using only a single image. The reconstructed models were output to the stereolithographic (STL) format, which could be directly read by the Object Studio program to produce the 3D models. STL is a file format that combines 3D systems to describe the entire geometry or structure of a 3D object without colour or texture; it has been extensively used in rapid prototyping and 3D printing systems. The resolution of the STL file is 100 *μ*m × 100 *μ*m × 100 *μ*m which is calculated by the thickness 8 mm of the model divide the image number 80 and the length 50 mm of the model divide the pixel number 500. With the help of 3D printing technology, the transparent photopolymer material VeroClear was employed to print the model, as shown in Fig. 2. Table 1 presents the physical/mechanical parameters of the printing material.

#### 2.2 Photoelastic testing setup

Figure 3 shows the photoelastic setup for capturing and identifying the fringes of the printed models. The arrangement of the dark-field circular polariscopes used in the experiments is shown in Fig. 3(a).

The experimental devices, which included polarisers, two quarter-wave plates, and an analyser, are shown in Fig. 3(b). The angles of the optical axes of the polariser and analyser and the fast axes of the two quarter-wave plates to the vertical direction are shown in Fig. 3(b). A digital servo-control universal testing machine with a loading capacity of 100 kN and loading rate of 0.5 mm/min was used to exert continuous compressive pressure on the model. A digital camera (SONY HDR-SR12E) with a dynamic resolution of 1440 × 1080 pixels and a frame rate of 25 f/s was employed to record the fringe contours during the loading process, and a monochromatic light filter, whose central wavelength is 532 nm and light transmittance is 73%, was installed in the lens of the digital camera for a clearer observation.

#### 2.3 Quantification method for isochromatic fringe maps

The method used in the present study was developed by employing the analytic solution of the diametric compression of a circular disk in the dark field of circular polarized light. Based on the theoretical resolution of the diametric compression of a circular disk [52], the equation describing the stress components can be written as

*r*is the disk radius,

*t*is the thickness of the circular disk, and

*P*is the loading pressure. Theoretically, the information shown in the photoelasticity fringe patterns is the principal stress difference, so Eq. (1) can be rewritten as

*σ*and

_{1}*σ*are the principal stresses. Based on the stress-optical law, the relationship between principal stress difference and fringe orders can be written aswhere the

_{2}*f*is the optical fringe constant of the photoelastic material and

*N*is the fringe order. Then, the retardation

*δ*caused by the principal stress difference isAccording to the matrix theory of photoelasticity [53], the relationship between light intensity and retardation can be written asSubstituting Eqs. (4) and (2) into Eq. (5), the light intensity can be rewritten as

To determine the whole-field fringe orders under a specific loading stage, the light intensity maps of the continuous loading process have to be obtained. Figure 4 shows the light intensity maps (the photoelastic fringe patterns calculated by Eq. (6)) where the photoelastic fringe constant of the material and radius are identical to those of the physical disc, and the intensity of light source is 255 (this is just an example of the light intensity; other values of light intensity are also suitable). The light intensity curve of the centre point of the disk is taken as an example to introduce the method to determine the fringe orders shown in Fig. 5. Under the circular polarized light of the dark field, the darkest fringes represent the integral fringe orders. First, the light intensity curve is filtered according to the following conditions: the light intensity will take the value *I _{p}*, which is larger than the medium line value, if it is above the medium line value

*I*= 127.5; otherwise, the light intensity will take the value

_{m}*I*, which is smaller than the medium line value. Second, the fringe orders of the integral parts can be determined from the flow chart shown in Fig. 6. Finally, the integer parts

_{v}*N*of the fringe orders can be written as

_{i}For the ending parts after the last valley, there are three cases that must be evaluated to calculate the corresponding phase: C_{1}, C_{2}, and C_{3} [see Fig. 5]. C_{1} is located in the part after the last peak, C_{2} represents the point before the last peak but higher than the value of the medium line (red line), and C_{3} refers to the point before the last peak but lower than the value of the medium line (red line).

Based on Eq. (5) and by taking the background light intensity into consideration, the cosine function can also be described by the peak value *I _{p}* and valley value

*I*:

_{v}Both the equations that make up Eq. (9) are used to calculate the retardation of the three cases shown in Fig. 5. The first equation is applied to C_{1}, and the second equation is applied to C_{2} and C_{3}. Because of the non-uniform wave length of the light source in the experiment, the incoherent superposition in the interference of the lights with different wave lengths causes the peak and valley value to be close to the medium line. Figure 7 shows a simple example of the interference of two different waves, and the blue line represents the interference caused by the black- and red-dashed lines. Thus, the peak and valley before the calculated point and close to it are selected as the values in Eqs. (9) and (10). The calculated equations are listed in Table 2.

As shown in Fig. 5, *I′ _{p}* is the light intensity of the last peak, and

*I′*is the light intensity of the last valley. Then,

_{v}*I′*and

_{p}*I′*can be substituted into Eq. (5), and the retardation on the scale of [0, π] for the three simulations shown in Fig. 5 can be calculated with the equations listed in Table 2.

_{v}As aforementioned, the integer and fractional parts of the fringe orders at every pixel can be determined using a raster scanning approach. Hence, the fringe orders *N* can be calculated with

*N*and

_{i}*N*are the integer and fractional parts, respectively, of the fringe orders.

_{f}#### 2.4 Illustration of method effectiveness

The calculated results of the experimental and theoretical principal stress differences in the horizontal and vertical directions were counted according to the curves shown in Figs. 8(a) and 8(b), where the red line indicates the experimental results and the blue line shows the theoretical results. The experimental and theoretical lines are clearly close to each other, which indicates that this method can accurately determine the whole-field principal stress difference in a homogeneous and continuous model. Note that the fringe orders at the centre point were calculated with the proposed method. The photoelastic fringe value of the printing material was determined by averaging the values of the three stages, resulting in a measured value of 38.76 ± 2 N/mm. Table 3 presents a comparison of the mechanical parameters between printing materials and conventional photoelastic materials, thus indicating that the printing material used in the experiment is suitable for photoelastic tests.

To evaluate the accuracy of the method, especially in the zones of stress concentration and near the hole boundaries, the comparison of the photoelastic patterns directly observed from the experiment and the fringe order contours calculated using the proposed method of two typical models, i.e., a circular disk and a cubic plate with a circular hole, that underwent diametrical compression loads, is shown in Fig. 9. It is evident that the results are very close to each other. In addition, the analytical errors of the two typical models, using the proposed unwrapping method and traditional photoelastic method, were compared. Under the circular polarized light of the dark field, the darkest fringes in the model represent the integral fringe orders and the brightest fringes indicate the half integer fringe orders. As shown in Fig. 10, the integral and half integer fringe orders of some pixels on the routes of the typical model obtained using the conventional technique were calculated and regarded as the standard results. Then, the fringe orders of the same pixels on the same routes of the models using the proposed method were evaluated and compared with the standard results. The accuracy of the proposed method was assessed using the difference between the results of the two different methods, formulated by the Eq. (12)

where*n*

_{1}represents the fringe orders calculated using the conventional photoelastic method and

*n*

_{2}represents those calculated using the proposed method. Figure 10 shows the comparison results. Even when the fringe order increases to 12, the bias of the results between the two methods is less than 5%, as shown in Figs. 2(a) and 2(b). This means that the proposed method exhibits good accuracy when extracting and characterizing the high-density fringes that occur around the sub-structures of media.

## 3. Photoelastic testing results and discussion

#### 3.1 Stress fringe patterns

The photoelastic technique showed its superiority when used to quantitatively visualise and determine the stress distribution in complex heterogeneous pore structures. Figures 11(a)–11(d) depict the isochromatic patterns of the test model subjected to four different stages of continual compressive stress by using circular polariscopes in the dark-field configuration. These fringe patterns represent the principal stress difference contours. The darkest and brightest fringes, which reveal the equivalent zones of σ_{1} - σ_{2}, were easily observed in monochromatic cases. The distance between the adjacent dark or bright fringes is indicative of the variable gradient. The zones with small gradients indicate that the principal stress difference was small in these areas. Zones with high-density fringes represent the formation of a stress concentration. The local stress concentration clearly appeared near regional boundaries of pores where the radius of curvature was relatively small. The stresses concentrated in these local boundaries increased rapidly with added load, and even fringes in some areas could not be distinguished before the end of the elastic stage on the stress–strain curve because the fringe density was too high or plastic zones formed. Note that the printing material used to fabricate the specimen exhibited good brittleness accompanied by apparent deformability; thus, it could fail by shear stress. Detailed information can refer to the relevant work [38, 57]. Overall, the changes in the photoelastic patterns indicate the distribution of the whole-field principal stress difference. Denser fringes indicate a greater stress concentration.

#### 3.2 Identification and quantification of fringe orders

The stage with a pressure of 7.5 MPa is taken as an example to illustrate the calculated results of whole-field fringe orders. The calculated results of the other stages will be discussed in the next section. Figure 12(a) is the monochromatic fringe pattern directly captured from the photoelastic test, and Fig. 12(b) is the calculated result of the wrapped fringe orders that are equivalent to the wrapped phase pattern. In the dark field configuration of circular polariscopes, the darkest fringe represents the integer portion of the fringe orders, and the lightest fringe denotes one half of the fringe orders. Hence, in the contour of wrapped fringe orders, a fringe order equal to zero or one represents the integer portion; a fringe order of 0.5 denotes one half of the fringe orders. The uniformity of the distributions for the integer portion and one half of the fringe orders in Figs. 12(a) and 12(b) indicates that the method worked for evaluating the wrapped fringe orders. Figure 12(c) shows the contour of the unwrapped fringe orders, which is equivalent to the contour of the principal stress difference. It clearly illustrates the concentration and discontinuity of the principal stress difference in porous structures. For instance, the principal stress difference concentrates in the local parts of pore boundaries where the radius of curvature is relatively small, and a discontinuity forms because of the pores.

#### 3.3 Principal stress difference

The principal stress difference can be obtained by multiplying the isochromatic fringe orders with the photoelastic fringe constant of the printing material. The distributions of the principal stress difference as calculated with the proposed method on horizontal and vertical lines of the test model are shown in Figs. 13(a) and 13(b). The existence of a complex pore structure clearly had a great influence on the stress distribution. For instance, a high concentration of the principal stress difference formed on the line in the horizontal direction. Furthermore, the degree of concentration increased rapidly with pressure, and the stress concentrations had a very serious effect. This indicates that the plastic failure caused by shear stress easily occurs on the boundary of pores where a high stress concentration forms.

## 4. Numerical simulation and comparison

#### 4.1 Finite element models and simulation

Photoelastic tests are extremely advantageous in visualising and measuring the whole-field stress of complex structures. Finite element models (FEMs) with identical pore geometries, material properties, and loading conditions were computed to validate the accuracy of the results of the photoelastic tests. The commercial software Mimics® (http://www.materialise.com/en/medical/mimics-innovation-suite) was employed as a tool to create an FEM from the porous structure (see Fig. 14). The model used tetrahedron elements with six degrees of freedom (DOF) at each node, thus resulting in 1,282,602 DOFs in the model. The elasticity modulus (*E*) and Poisson’s ratio (ν) of the elements were assigned in accordance with the mechanical properties of the printing materials, as listed in Table 1. In the numerical approach, pressure comparable to that exerted on the tested model was applied to the top elements of the numerical models, and the bottom elements were constrained in all directions. Therefore, the numerical simulation parameters and boundary conditions were the same as those used in the experiments. The meshed model was subsequently analysed with the finite element analysis program ANSYS^{©} 15.0, and the elastic constitution relationship was used in the calculation.

#### 4.2 Simulation Results

Figures 15(a)–15(d) plot the contours of the principal stress difference for the four points of the elastic stage of the stress–strain curve from the numerical simulation. The principal stress difference (σ_{1} − σ_{2}) is calculated based on the stresses at nodes of the simulation model. The calculation equations can be written as

Figures 15(e)–15(f) present the principal stress differences calculated by the proposed method at the same four stages in the photoelastic test for comparison with the numerical simulation. The whole-field stress evolution for the four stages can be explicitly displayed by choosing a uniform colour bar. The distributions of the principal stress difference in the experiment and simulation at each stage were highly similar, which indicates that the proposed method worked well at calculating the principal stress difference. However, the stress field evolutions for the experiment and simulation were also close to each other. Obviously, the whole-field principal stress difference increased when the model was exposed to a higher compression stress. Especially, a greater stress concentration was easily observed around pores, and even potential plastic zones where the fringes could not be distinguished formed.

To accurately verify the effectiveness of the experimental results, the calculated results of the experimental and numerical principal stress differences in the horizontal and vertical directions of the four different stages (see Figs. 15 (a)–15(h)) were determined by the curve in two directions (see Figs. 16(a) and 16(b)). The red lines indicate the experimental results, and the blue lines indicate the numerical results. The experimental and theoretical lines for the same stage are clearly close to each other, which indicates that this method can accurately determine the whole-field principal stress difference. However, high-density fringe zones were caused by the high stress concentration, and plastic zones formed at the irregularly shaped pore boundaries where the radius of curvature was smaller. This resulted in fringes that could not be clearly distinguished and captured. In summary, the experimental results were similar to the numerical results, which could distinctly display the distribution and evolution of the whole-field stress except in some high-stress regions.

#### 4.3 Comparison between the fringe contours of the experiments and numerical simulations

Allowing the verification of numerical simulations is an outstanding feature of photoelasticity. For another distinct comparison, it is critical to transfer the contours of the principal stress difference obtained from the FEM results to fringe patterns that are uniform to the experimental results calculated by the proposed method. Based on the optic-stress law, the fringe order can be written as

where*F*is the fringe order of the entire field,

_{N}*d*is the thickness of the model, and

*f*is the photoelastic material fringe constant at room temperature. The experimental fringe contours calculated by the proposed method (see Fig. 17(a)) that correspond to the fringe patterns directly observed from experiments represent the fractional part of the fringe orders. The fractional fringe order values (

*δ*) of the numerical simulation can be evaluated by using the following equation:where INT(

_{N}*F*) is the integer part of

_{N}*F*. Thus,

_{N}*δ*lies in the range of 0 <

_{N}*δ*< 1. The contour of

_{N}*δ*provided a wrapped isochromatic fringe pattern (see Fig. 17(b)), which represents an identical mechanical meaning with the experimental result as shown in Fig. 17(a). Thus, these fringe contours facilitated a comparison between the experimental observations and numerical solutions.

_{N}The stage with a pressure of 7.5 MPa was chosen as an illustration to compare the fringe contours, as shown in Figs. 17(a) and 17(b). Based on the comparison between the experimental calculation and numerical simulation results, the fringe contours calculated with the two methods clearly matched very well with each other in most regions. However, there were some local areas close to the model boundaries where the fringe contour calculated by the numerical simulation was different from the experimental results. In the experiments, the influences of the optic effect of the boundary, the evenness of the top and bottom surfaces of models, and the friction between the contact surfaces were the main reasons for the difference in fringe contours for the model boundaries with the numerical results.

## 5. Conclusions

Based on the natural porous structure of Berea sandstone, a complex porous model was prepared with a 3D printer using photopolymer materials. The continuous distribution and evolution of the whole-field stress in the complex porous model subjected to continual compressive loads were calculated with the proposed method. The applicability and accuracy of the method at evaluating the stress field in a complex pore structure were verified by comparing the experimental results with numerical simulation data. The results indicated that the proposed method worked well at quantifying and visualising the whole-field stress distribution and evolution in complex pore structures. The main findings are as follows:

- (1) The contours of whole-field stress fringes and wrapped fringe orders in the complex pore structures were directly obtained by computing the photoelastic patterns observed from continuous compression loading. Therefore, all of the fringe orders and wrapped fringe orders in the complex pore structures with different pressures in the elastic loading stage can be determined. It was shown that the applicability of the reported method to identify and unwrap the continuous whole-field stress is not affected by the complexity of a structure.
- (2) The experimental and numerical results for the stress distribution and evolution of the complex porous models with identical geometries, physical properties, and loading conditions indicated that the proposed method is applicable to determining the whole-field stress evolution in complex pore structures and has a good precision.
- (3) The experimental and numerical results indicated that high-stress-concentration areas will appear around porous boundaries with a small radius of curvature, and the degree of concentration will rapidly increase with the load. Thus, the existence of irregular pores has an enormous influence on the distribution of the entire stress field.

It is noteworthy that, in this study, we concentrated on the distribution and evolution of the whole-field stress in complex porous structures where the fringes can be clearly observed and the photoelastic optic-stress law is applicable. Comparison indicates that the proposed method is capable of extracting and unwrapping the highly dense fringes even when the fringe order increases to 12, which means the domain of interest may enter the regime of plastic deformation, which is hard to achieve using traditional photoelastic methods. Nevertheless, the prerequisite is that the dense fringes can be identified and distinguished using high- resolution recording systems. If the fringes of the domain are too dense to be identified, such as the domains where large plastic deformation occurs, the proposed method may not be applicable for this problem. Moreover, the optic-stress law of photoelasticity will not be applicable in these domains and a new law may need to be proposed to translate the fringe orders into the state of stress. This is a significant challenge in the study of photoelasticity and photoplasticity.

Although the preliminary results in this study indicate that the proposed method is capable of quantitatively characterising the whole-field stress and its continuous evolution inside complex pore structures, some challenges remain to be solved. For instance, in the areas where the fringe density was too high to distinguish the individual fringes, the potential plastic zones need to be clarified. Moreover, for the purpose of solving engineering problems, better similitude not only in geometric structures but also physical properties between the printed models and natural objects must be achieved [38, 57]. Solutions to these problems rely on further development of the mechanics of photoplasticity, improving the resolution and film speed of digital cameras and the progress of three-dimensional printing technology. Methods to distinguish the high-density fringe areas and quantify the strain in photoplastic zones will be discussed in future work.

## Funding

National Natural Science Foundation of China (NSFC) (51727807, 51674251, 51374213 and 51125017); State Key Research Development Program of China (2016YFC0600705); National Major Project for Science and Technology (2017ZX05049003-006); Fund for Creative Research and Development Group Program of Jiangsu Province (2014-27); Priority Academic Program Development of Jiangsu Higher Education Institutions (Grant No. PAPD2014).

## References and links

**1. **K. Sun, J. Tan, and D. Wu, “The research on dynamic rules of crack extension during hydraulic fracturing for oil shale in-situ exploitation,” Procedia Environ. Sci. **12**, 736–743 (2012). [CrossRef]

**2. **L. Yang, J. Zhao, W. Liu, Y. Li, M. Yang, and Y. Song, “Microstructure observations of natural gas hydrate occurrence in porous media using microfocus X-ray computed tomography,” Energy Fuels **29**(8), 4835–4841 (2015). [CrossRef]

**3. **M. Zhang, H. Shimada, T. Sasaoka, K. Matsui, and L. Dou, “Evolution and effect of the stress concentration and rock failure in the deep multi-seam coal mining,” Environ. Earth Sci. **72**(3), 629–643 (2014). [CrossRef]

**4. **R. Zhang, T. Ai, H. W. Zhou, Y. Ju, and Z. T. Zhang, “Fractal and volume characteristics of 3D mining-induced fractures under typical mining layouts,” Environ. Earth Sci. **73**(10), 6069–6080 (2015). [CrossRef]

**5. **L. Valoroso, L. Chiaraluce, and C. Collettini, “Earthquakes and fault zone structure,” Geology **42**(4), 343–346 (2014). [CrossRef]

**6. **G. W. Bird and W. S. Fyfe, “The nuclear waste disposal problem—An overview from a geological and geochemical perspective,” Chem. Geol. **36**(1–2), 1–13 (1982). [CrossRef]

**7. **P. Bénézeth, B. Ménez, and C. Noiriel, “CO_{2} geological storage: Integrating geochemical, hydro-dynamical, mechanical and biological processes from the pore to the reservoir scale,” Chem. Geol. **265**(1–2), 1–2 (2009). [CrossRef]

**8. **A. J. Desbarats, D. R. Boyle, M. Stapinsky, and M. J. Robin, “A dual‐porosity model for water level response to atmospheric loading in wells tapping fractured rock aquifers,” Water Resour. Res. **35**(5), 1495–1505 (1999). [CrossRef]

**9. **W. R. Judd, “Underground excavations in rock,” Eng. Geol. **19**(3), 244–246 (1983). [CrossRef]

**10. **H. Ma, C. Yang, Y. Li, X. Shi, J. Liu, and T. Wang, “Stability evaluation of the underground gas storage in rock salts based on new partitions of the surrounding rock,” Environ. Earth Sci. **73**(11), 6911–6925 (2015). [CrossRef]

**11. **J. G. Ramsay, “The crack-seal mechanism of rock deformation,” Nature **284**(5752), 135–139 (1980). [CrossRef]

**12. **L. S. Dimas, T. Giesa, and M. J. Buehler, “Coupled continuum and discrete analysis of random heterogeneous materials: Elasticity and fracture,” J. Mech. Phys. Solids **63**(63), 481–490 (2014). [CrossRef]

**13. **C. Farhat, P. Avery, T. Chapman, and J. Cortial, “Dimensional reduction of nonlinear finite element dynamic models with finite rotations and energy-based mesh sampling and weighting for computational efficiency,” Int. J. Numer. Methods Eng. **98**(9), 625–662 (2014). [CrossRef]

**14. **C. Lempp, K. M. Shams, and N. Jahr, “Approaches to stress monitoring in deep boreholes for future CCS projects,” Environ. Earth Sci. **67**(2), 435–445 (2012). [CrossRef]

**15. **A. Mazaira and P. Konicek, “Intense rockburst impacts in deep underground construction and their prevention,” Can. Geotech. J. **52**(10), 1426–1439 (2015). [CrossRef]

**16. **M. Bai, J. C. Roegiers, and D. Elsworth, “Poromechanical response of fractured-porous rock masses,” J. Petrol. Sci. Eng. **13**(3–4), 155–168 (1995). [CrossRef]

**17. **D. Li and L. N. Y. Wong, “The Brazilian disc test for rock mechanics applications: review and new insights,” Rock Mech. Rock Eng. **46**(2), 269–287 (2013). [CrossRef]

**18. **S. P. Timoshenko, *Theory of Elastic Stability* (McGraw-Hill, 1961).

**19. **P. Kulatilake, J. Liang, and H. Gao, “Experimental and numerical simulations of jointed rock block strength under uniaxial loading,” J. Eng. Mech. **127**(12), 1240–1247 (2001). [CrossRef]

**20. **R. Katsman, E. Aharonov, and H. Scher, “Numerical simulation of compaction bands in high-porosity sedimentary rock,” Mech. Mater. **37**(1), 143–162 (2005). [CrossRef]

**21. **B. C. Haimson, L. W. Tunbridge, M. Y. Lee, and C. M. Cooling, “Measurement of rock stress using the hydraulic fracturing method in Cornwall, U.K.-Part II. Data reduction and stress calculation,” Int. J. Rock Mech. Min. Sci. Geomech. Abstr. **26**(5), 361–372 (1989). [CrossRef]

**22. **R. Ikeda, Y. Iio, and K. Omura, “In situ stress measurements in NIED boreholes in and around the fault zone near the 1995 Hyogo-ken Nanbu earthquake Japan,” Isl. Arc **10**(3–4), 252–260 (2010).

**23. **A. Lehtonen, J. W. Cosgrove, J. A. Hudson, and E. Johansson, “An examination of in situ rock stress estimation using the Kaiser effect,” Eng. Geol. **124**(1), 24–37 (2012). [CrossRef]

**24. **K. Sakaguchi, Y. Obara, T. Nakayama, and K. Sugawara, “Accuracy of rock stress measurement by means of conical-ended borehole technique,” Shigen-to-Sozai **108**(6), 455–460 (1992). [CrossRef]

**25. **D. A. Demer, M. A. Soule, and R. P. Hewitt, “A multiple-frequency method for potentially improving the accuracy and precision of in situ target strength measurements,” J. Acoust. Soc. Am. **105**(4), 2359–2376 (1999). [CrossRef]

**26. **V. Vajdova, P. Baud, L. Wu, and T. F. Wong, “Micromechanics of inelastic compaction in two allochemical limestones,” J. Struct. Geol. **43**(7), 100–117 (2012). [CrossRef]

**27. **V. Vajdova, W. Zhu, T. M. N. Chen, and T. F. Wong, “Micromechanics of brittle faulting and cataclastic flow in Tavel limestone,” J. Struct. Geol. **32**(8), 1158–1169 (2010). [CrossRef]

**28. **Y. Liu, F. J. H. G. Kessels, W. D. V. Driel, J. A. S. V. Driel, F. L. Sun, and G. Q. Zhang, “Comparing drop impact test method using strain gauge measurements,” Microelectron. Reliab. **49**(9–11), 1299–1303 (2009). [CrossRef]

**29. **R. Y. Makhnenko, J. Harvieux, and J. F. Labuz, “Paul-Mohr-Coulomb failure surface of rock in the brittle regime,” Geophys. Res. Lett. **42**(17), 6975–6981 (2015). [CrossRef]

**30. **A. S. Voloshin and T. Leng-Tsun, “Investigation of the stress singularities by enhanced moiré interferometry,” Eng. Fract. Mech. **43**(4), 477–486 (1992). [CrossRef]

**31. **S. Takao, S. Yoneyama, and M. Takashi, “Minute displacement and strain analysis using lensless Fourier transformed holographic interferometry,” Opt. Lasers Eng. **38**(5), 233–244 (2002). [CrossRef]

**32. **Z. Lei, M. Fu, and H. Yun, “Experimental study on interfacial shear transfer in partially-debonded aluminum/ epoxy joint,” Int. J. Adhes. **31**(2), 104–111 (2011). [CrossRef]

**33. **J. T. Pinto, F. Touchard, S. Castagnet, C. Nadot-Martin, and D. Mellier, “DIC strain measurements at the micro-scale in a semi-crystalline polymer,” Exp. Mech. **53**(8), 1311–1321 (2013). [CrossRef]

**34. **P. Forte, A. Paoli, and A. V. Razionale, “A CAE approach for the stress analysis of gear models by 3D digital photoelasticity,” Int. J. Interact. Des. Manuf. **9**(1), 31–43 (2015). [CrossRef]

**35. **A. J. Rosakis, O. Samudrala, and D. Coker, “Cracks faster than the shear wave speed,” Science **284**(5418), 1337–1340 (1999). [CrossRef] [PubMed]

**36. **I. A. Takacs, A. I. Botean, M. Hardau, and S. Chindris, “Displacement-stress distribution in a femoral bone by optical methods,” Procedia Technol. **19**, 901–908 (2015). [CrossRef]

**37. **K. Ramesh, M. P. Hariprasad, and S. Bhuvanewari, “Digital photoelastic analysis applied to implant dentistry,” Opt. Lasers Eng. **87**, 204–213 (2016). [CrossRef]

**38. **Y. Ju, L. Wang, H. Xie, G. Ma, Z. Zheng, and L. Mao, “Visualization and transparentization of the structure and stress field of aggregated geomaterials through 3D printing and photoelastic techniques,” Rock Mech. Rock Eng. **50**(6), 1383–1407 (2017). [CrossRef]

**39. **Y. Ju, H. Xie, Z. Zheng, J. Lu, L. Mao, F. Gao, and R. Peng, “Visualization of the complex structure and stress field inside rock by means of 3D printing technology,” Chin. Sci. Bull. **59**(36), 5354–5365 (2014). [CrossRef]

**40. **Y. Ju, Z. Zheng, H. Xie, J. Lu, L. Wang, and K. He, “Experimental visualisation methods for three-dimensional stress fields of porous solids,” Exp. Tech. **41**(4), 331–344 (2017). [CrossRef]

**41. **M. Ramji and K. Ramesh, “Whole field evaluation of stress components in digital photoelasticity-Issues, implementation and application,” Opt. Lasers Eng. **46**(3), 257–271 (2008). [CrossRef]

**42. **M. J. Huang and P. C. Sung, “Regional phase unwrapping algorithm for photoelastic phase map,” Opt. Express **18**(2), 1419–1429 (2010). [CrossRef] [PubMed]

**43. **P. A. Júnior, F. G. Vieira, C. A. Magalhães, J. S. Ribeiro, and I. G. Rios, “Numerical method to digital photoelasticity using plane polariscope,” Opt. Express **24**(12), 12617–12624 (2016). [CrossRef] [PubMed]

**44. **Z. K. Lei, R. X. Bai, W. Qiu, L. B. Deng, and B. C. Huang, “Flexural effects of sandwich beam with a plate insert under in-plane bending,” Opt. Laser Technol. **44**(5), 1223–1231 (2012). [CrossRef]

**45. **V. S. Prasad, K. R. Madhu, and K. Ramesh, “Towards effective phase unwrapping in digital photo-elasticity,” Opt. Lasers Eng. **42**(4), 421–436 (2004). [CrossRef]

**46. **A. Ajovalasit, G. Petrucci, and M. Scafidi, “Review of RGB photoelasticity,” Opt. Lasers Eng. **68**, 58–73 (2015). [CrossRef]

**47. **J. A. Quiroga, A. García-Botella, and J. A. Gómez-Pedrero, “Improved method for isochromatic demodulation by RGB calibration,” Appl. Opt. **41**(17), 3461–3468 (2002). [CrossRef] [PubMed]

**48. **A. Ajovalasit, G. Petrucci, and M. Scafidi, “RGB photoelasticity: review and improvements,” Strain **46**(2), 137–147 (2010). [CrossRef]

**49. **A. Ajovalasit and G. Petrucci, “Developments in RGB photoelasticity,” Appl. Mech. Mater. **3–4**, 205–210 (2005). [CrossRef]

**50. **T. W. Ng, “Photoelastic stress analysis using an object step-loading method,” Exp. Mech. **37**(2), 137–141 (1997). [CrossRef]

**51. **B. Zuccarello and G. Tripoli, “Photoelastic stress pattern analysis using Fourier transform with carrier fringes: influence of quarter-wave plate error,” Opt. Lasers Eng. **37**(4), 401–416 (2002). [CrossRef]

**52. **J. H. Yang, F. Q. Wu, and J. Z. Sun, “Estimation of the tensile elastic modulus using Brazilian disc by applying diametrically opposed concentrated loads,” Int. J. Rock Mech. Min. Sci. **46**(3), 568–576 (2009). [CrossRef]

**53. **B. P. S. Theocaris and E. E. Gdoutos, *Matrix Theory of Photoelasticity* (Springer-Verlag, 1979).

**54. **A. Pawlak and A. Galeski, “Determination of stresses around beads in stressed epoxy resin by photoelasticity,” J. Appl. Polym. Sci. **86**(7), 1436–1444 (2002). [CrossRef]

**55. **R. P. Singh, J. Lambros, A. Shukla, and A. J. Rosakis, “Investigation of the mechanics of intersonic crack propagation along a bimaterial interface using coherent gradient sensing and photoelasticity,” P. Roy Soc. A-Math. Phys. **453** (1967), 2649–2667 (1997). [CrossRef]

**56. **A. Kobayashi, *Handbook on Experimental Mechanics* (Prentice Hall, 1987).

**57. **L. Wang, Y. Ju, H. Xie, G. Ma, L. Mao, and K. He, “The mechanical and photoelastic properties of 3D printable stress-visualized materials,” Sci. Rep. **7**(1), 10918 (2017). [CrossRef] [PubMed]