Theoretical and numerical studies of crack initiation and propagation in rock masses under freezing pressure and far-f i eld stress

Yongshui Kng,Qunsheng Liu,Xioyn Liu,Shiing Hung

aState Key Laboratory of Geomechanics and Geotechnical Engineering,Institute of Rock and Soil Mechanics,Chinese Academy of Sciences,Wuhan,430071, China

bKey Laboratory of Rock Mechanics in Hydraulic Structural Engineering of Ministry of Education,School of Civil Engineering and Architecture,Wuhan University,Wuhan,430072,China

Theoretical and numerical studies of crack initiation and propagation in rock masses under freezing pressure and far-f i eld stress

Yongshui Kanga,*,Quansheng Liua,Xiaoyan Liub,Shibing Huanga

aState Key Laboratory of Geomechanics and Geotechnical Engineering,Institute of Rock and Soil Mechanics,Chinese Academy of Sciences,Wuhan,430071, China

bKey Laboratory of Rock Mechanics in Hydraulic Structural Engineering of Ministry of Education,School of Civil Engineering and Architecture,Wuhan University,Wuhan,430072,China

A R T I C L EI N F O

Article history:

Received 31 December 2013

Received in revised form

12 May 2014

Accepted 29 May 2014

Available online 9 August 2014

Freeze—thaw action

Water-bearing rocks exposed to freezing temperature can be subjected to freeze—thaw cycles leading to crack initiation and propagation,which are the main causes of frost damage to rocks.Based on the Griff i th theory of brittle fracture mechanics,the crack initiation criterion,propagation direction,and crack length under freezing pressure and far-f i eld stress are analyzed.Furthermore,a calculation method is proposed for the stress intensity factor(SIF)of the crack tip under non-uniformly distributed freezing pressure.The formulae for the crack/fracture propagation direction and length of the wing crack under freezing pressure are obtained,and the mechanism for coalescence of adjacent cracks is investigated. In addition,the necessary conditions for different coalescence modes of cracks are studied.Using the topology theory,a new algorithm for frost crack propagation is proposed,which has the capability to def i ne the crack growth path and identify and update the cracked elements.A model that incorporates multiple cracks is built by ANSYS and then imported into FLAC3D.The SIFs are then calculated using a FISH procedure,and the growth path of the freezing cracks after several calculation steps is demonstrated using the new algorithm.The proposed method can be applied to rocks containing f i llings such as detritus and slurry.

?2014 Institute of Rock and Soil Mechanics,Chinese Academy of Sciences.Production and hosting by Elsevier B.V.All rights reserved.

1.Introduction

Low temperatures and harsh weather can be observed in nearly half of the earth’s land area in winter.Frost weathering in cold regions,where a quarter of the land mass is under permafrost or covered by snow,poses a serious threat to the stability of geotechnical engineering endeavors.Frost shattering caused by freeze—thaw cycles in moist rocks,as well as expansion and contraction due to temperature changes,is the major reason for physical weathering of rocks(Yang and Chen,2004;Seto,2010; Kang et al.,2012a,2013).The freeze—thaw deterioration of rock is mainly ref l ected by the propagation of cracks under freezing pressure.Therefore,thewater—ice phase transition and its effect on the evolution of the crack network should be investigated to reveal the mechanism of freeze—thaw damage to rocks.

The freeze—thaw action is caused by temperature f l uctuations around the freezing point of water.Freeze—thaw action,which is also referred to as ice crystal growth or frost shattering,occurs when water in cracked rock is frozen and expands.As shown in Fig.1,when water in cracks is frozen,a volume expansion of approximately 9%occurs,initiating a large freezing pressure that deepens and expands the crack width.When ice in the rock thaws, the water can f l ow further into the rock through these new cracks, facilitating the next freeze—thaw cycles.Thus,the repeated freeze—thaw processes of water in the cracks have long been recognized as a signif i cant incentive for mechanical weathering of rocks(Hilbich, 2010;Su et al.,2010).

Frost damage in rocks is mainly caused by crack propagation under freezing pressure.The heterogeneous features of natural rocks give rise to various types of defects or local stress concentrations that might be the sources of microcracks.Based on the assumption that pre-existing microcracks are considered the most likely sources for crack propagation,Paterson and Wong(2005) analyzed the sliding crack model in detail.

In recent years,many studies have concentrated on the frost damage induced by freeze—thaw action in rock masses.However, as a result of the complicated freeze—thaw process and mathematical diff i culties,most of the previous research concentrated onspecial cases,and the evolution of the crack network caused by freeze—thaw action has rarely been reported(Matsuoka,1990, 2008;Neaupane et al.,1999;Matsuoka and Murton,2008).

In this regard,the present paper tries to present a theoretical method to predict crack propagation under freezing pressures and far-f i eld stresses.The direction and length of crack propagation are derived based on the Griff i th theory of brittle fracture mechanics. The conditions necessary for different coalescence modes of cracks are also studied.

Fig.1.Sketch of freezing-induced crack.

2.Stress intensity factors of freezing-induced rock cracks

2.1.Stress intensity factors

Rock is basically regarded as a brittle material and as such it cannot resist large tensile stresses.In general,cracks in rocks under freezing pressures and far-f i eld stresses can be regarded as a mixed mode of tensile crack(Mode I)and shear crack(Mode II).The freezing pressure caused by volume expansion of ice can be considered as a normal pressure imposed on the crack plane,while the far-f i eld stress might induce both normal and shear stresses on the crack plane.The sliding crack model considers the source of tensile stress concentrations located at the tips of inclined preexisting cracks.The far-f i eld stress would induce shear traction on the crack plane.If this shear traction is suff i ciently high to overcome the frictional resistance along the interface between the crack and the ice,frictional slip occurs and a tensile stress concentration is induced at the tips of the sliding crack,with the appearance of wing cracks(Matsuoka,2008).The driving force for nucleation is characterized by the stress intensity factor(SIF),KI,at the position of wing crack initiation.With the increasing load or freezing pressure,KIwill increase to the critical value KIC(the fracture toughness)at which point a wing crack is nucleated.

Fig.2.Freezing pressure imposed on the crack.(a)Non-uniform stress.(b)Crack with concentrated force.

Numerically,the freezing pressure on the crack might not be uniform;therefore attention should be focused on the SIF induced by non-uniform stress.As shown in Fig.2a,the model is conceived as a plane stress problem.Normal and shear stresses are nonuniformly distributed on the faces of the crack.Therefore,the integration method should be used for this problem.The nonuniform stress is roughly transformed into numerous point forces as shown in Fig.2b,which has been studied in previous research.

The SIFs at a crack tip caused by the point force on the crack plane can be expressed as follows(Yin,1992;China Aviation Academy,1993):

where P(x)andτ(x)are the normal and shear forces at point x, respectively;a is the half length of the crack;dKIand dKIIare the Mode I and Mode II SIFs at the crack tip caused by the concentrated force.For the plane stress problem,κ=（3-ν）/（1+ν）,whereνis the Poisson’s ratio of the rock.Thus,the total SIFs caused by the non-uniform stress can be expressed as

It is necessary to def i ne the effective SIFs when the results are applied to numerical simulations.As shown in Fig.3,the shear and normal stresses are considered to act on elements.For the plane stress problem,the effective SIFs can be calculated as follows:

Fig.3.Elements of the crack.(a)Sketch of a freezing-induced crack.(b)Elements of the crack interface.

2.2.Freezing pressure imposed on cracks

Crack initiation and propagation of rocks induced by freezing pressure are quite complex processes.If the crack is initially saturated,unfrozen water might migrate out of crack as the ice crystal grows and expands.On the other hand,if the crack is unsaturated, the water could migrate to the freezing edge of the crystal.Theoretically,freezing pressure is signif i cant only when the crack is conf i ned or the freezing rate is high enough to form an ice plug that conf i nes the unfrozenwater,and consequentlya signif i cant amount of water can leak out of the crack tip(Fig.4).

The volume expansion of conf i ned water in crack tip can be written as

where t is the time,k is the f l uid permeability coeff i cient,?P is the freezing pressure gradient,φTdenotes the segregation potential,?T is the temperature gradient,ρwis the density of water,g is the gravity acceleration vector,and n is the unit direction vector of f l uid water in the crack.

The maximum value of the freezing pressure in the crack tip can be written as follows:

3.Stress distribution in the crack tip region

A crack,as shown in Fig.5a,is considered,whereσ1andσ2are the normal stresses imposed on the rock,andβis the angle between the crack plane and the direction ofσ1.The model is also considered as a plane stress problem,and the crack is initially water saturated.As the temperaturedecreases tothe freezing point,a thin ice layer will appear in the crack,which would induce freezing pressure on the crack face.In addition,a shear stress and a normal stress would be induced byσ1andσ2on the crack plane,as shown in Fig.5b.

The shear and normal stresses on the crack plane can be written as

Fig.4.Ice plug effect in freezing-induced crack.

Fig.5.Composite crack under compression and shearing conditions.

According to theory of superposition for Mode I and Mode II cracks,the stress components in Fig.6 can be expressed as follows (Yin,1992):

Fig.6.The stress f i eld in the crack tip region.

4.Growth path of the wing crack

4.1.Direction of crack initiation

A wing crack grows out of the initial plane of the sliding crack when nucleated.The prediction of crack propagation behavior requires the computation of the SIF at each stage as well as the wing crack trajectory.The formulation and solution of this fracture mechanics problem can be very complicated(Matsuoka,2008).This paper adopts the maximum circumferential stress criterion to analyze the freezing crack initiation criterion.It is assumed that crack initiation takes place in the direction vertical tothe maximum circumferential stress.Therefore,the partial derivative ofσθtoθcan be calculated as follows(Zhuang,2006;Yang,2010):

Whenσθθ,θ≤0 andσθθ,θθ＜0,σθθreaches the maximum value. Since0≤θ＜π,weobtaincos(θ/2)≠0.Therefore, KIsinθ+2KII［cosθ-sin2（θ/2）］=0.Then we have

Dividing both sides of Eq.(12)by 2 cos2（θ/2）,we get the following equation:

The solution to Eqs.(13)and(14)can be written as

Generally,the symbol“±”means that we get two values ofθ0whenσθθ,θ=0.Substitutingθ=θ0into Eq.(16),we can get the value ofσθθ,θθ.

Ifσθθ,θθ＜0,θ0is supposed to be the direction of crack initiation. It is evident that one of the values ofθ0is positive and the other is negative when KII≠0.According to the expression ofσθθ,the sign ofθ0should be opposite to that of KIIin order to get the maximum value ofσθθ.

Specif i cally,when KI=0 and KII＞0,the case becomes the problem of a Mode II crack;from the present analysis,we can obtain the crack initiation angle:

According tothe conclusion of Erdogan and Sih(1963),the crack initiation criterion can be expressed as follows:

where KICis a critical parameter referred to as the crack toughness, ref l ecting rock’s ability to resist crack propagation.When the calculated KIexceeds KIC,crack initiation commences.In the case of a pure tensile crack,θ0=0 and KI=KIC=σ(πa)1/2when KII=0; whilstθ0=-70.53°,KII=0.866KICwhen KI=0.

4.2.Length of the wing crack

When the direction of crack initiation is determined,the crack length should be specif i ed to demonstrate the path of crack growth. As shown in Fig.7,the angle between the wing crack and main crack is designated asθ0.Crack propagation is regarded as unstable crack growth and each step is a transient process.

According to the theory proposed by Erdogan and Sih(1963),a wing crack can be considered as a pure tensile crack(Mode I).As shown in Fig.7,Ftis the traction force normal to the face of wing crack,which is caused by stresses imposed on the face of main crack at the intersection:

where b is the width of the crack.Considering the plane stress condition,we have b=1.

The wing crack is considered to be a pure tensile crack and the SIF at the tip of the wing crack can be expressed as

The circumferential stress will gradually be reduced along the path of crack propagation;propagation ceases when the SIF is reduced to the critical value KIC:

Fig.7.Sketch of a wing crack.

Then,L can be deduced from Eq.(21):

Numerically,L can be calculated using the following formula:

5.Coalescence mechanism of neighboring cracks

Neighboring cracks will affect each other under freezing pressures and far-f i eld stresses.From a micromechanical point of view, the f i nal stage of crack propagation is crack coalescence.The complex interaction of stress f i elds from numerous neighboring cracks plays an important role in the evolution of the crack network.The individual cracks may coalesce with each other as the wing cracks propagate in the rock.Three coalescence patterns are considered in the freezing cracked rock:(1)coalescence of tension cracks,which is considered as tensile coalescence,(2)coalescence of shear cracks,considered as shear coalescence,and(3)a mixed pattern of tensile and shear coalescence.

5.1.Tensile coalescence

As shown in Fig.8,the tensile coalescence pattern is caused by tensile wing cracks.When the wing cracks coalesce with each other,or the tip of awing crack intersects with another pre-existing crack,tensile coalescence occurs.The growth path of wing cracks can be derived from the analyses in Section 4.

5.2.Shear coalescence

As shown in Fig.9,the rock bridge is broken due to shear stress. This case is considered as a plane stress problem.According to the equilibrium principle,the resultant forces in two directions equal zero as expressed by the following equation:

Fig.8.Sketch of tensile coalescence.

Fig.9.Sketch of shear coalescence.

As shown in Fig.9,the following relations can be derived (Vásárhelyi and Bobet,2000;Zhuang,2006;Yang,2010;Ning et al., 2011;Stakhovsky,2011):

Then,we can get

where

According to the Mohr—Coulomb criterion,whenσbandτbsatisfy the following relation,shear coalescence occurs:

where cband φ are the cohesion and internal friction angle of the rock,respectively.

5.3.Mixed pattern of tensile and shear coalescence

AsshowninFig.10,wingcrackswillnotintersectwitheachother. However,asthewingcrackgrows,therockbridgebetweenthecracks becomes thinner.The ability to resist shear damage declines.

Similar to the analysis of shear coalescence,we obtain the following equations:

Then,we can get

Fig.10.Sketch of mixed pattern of tensile and shear coalescence.

The solution of Eq.(31)can be written as

Whenσbandτbsatisfy Eq.(29),shear coalescence occurs.

6.A new algorithm for frost crack propagation

The ice pressure is the primary driving force of crack initiation and propagation.As shown in Fig.11a,the ice-f i lled crack is def i ned as a double-faced structure with a certain open width,while the tip is simplif i ed as an equilateral triangle.As shown in Fig.11b,the elements intersected by the new crack can be classif i ed into two categories:fully cracked elements(e.g.A1,A2,A6)and partially cracked elements(e.g.A3).The update of elements in the cracked region is the most important step for the evolution of the rock crack network.The algorithm for frost crack propagation should have the capability to update the grids in the cracked regions and to automatically distinguish the crack coalescence.

In order to achieve this,the following steps should be taken:

(1)The growing paths of each crack tip are calculated and def i ned according to Eqs.(15)and(23).

(2)Crack coalescence is judged by the topology method.

(3)The cracked elements are identif i ed by the element ergodic loop method using a computing procedure.

(4)The cracked elements are updated and the parameter values of the new elements are def i ned for the next calculating cycle.

6.1.Basic concepts of the algorithm

In order to illustrate the algorithm explicitly,the following concepts should be introduced as shown in Fig.12.

(1)Base points B1—B5:the key points on the crack face that control the shape of the new crack.

(2)Edge points E1—E3:the outer nodes of the new elements.

(3)Base lines lB1and lB2:the connecting lines of the base points in sequence.

Fig.11.An example to illustrate cracked elements.(a)The ice-f i lled crack.(b)Cracked elements.

Fig.12.Basic concepts of the proposed algorithm.

(4)Edge lines lE1and lE2:the connecting lines of the edge points in sequence.

(5)Mapping points M1and M2:the auxiliary points generated on the base line according to the distribution of edge points.The mapping points would be taken as the inner nodes of the new elements.

(6)Auxiliary base line:the line def i ned near the crack tip to optimize the elements in the crack tip region.

(7)Scanning plane:the plane consists of base lines and auxiliary lines.

The growth path of the crack is predicted after a certain number of steps.The path contains information on the initiation angle and propagation length.As shown in Fig.13,the tip point of the initial crack is def i ned as one of the base points of the new crack when θ＜30°orθ≥30°.

6.2.Update of cracked elements

According to the distribution position,the edge points are divided into f i ve zones,as shown as Fig.14.

The methods used toupdateelements in each zone are different. The following conditions should be discussed:

(1)In zones I and V,B1(B5)and the edge points are taken as the nodes of the new elements,as shown in Fig.15a.

(2)In zones II and IV,B1(B5)and B2(B4)as well as mapping points are taken as the nodes of the new elements,as shown in Fig.15b.The newelements will be generated in sequence by the computing procedure.The method of generating mapping points will be discussed in the following section.

(3)In zone III,the points B2,B3,B4and the edge points in zone III are taken as the nodes of the new elements,as shown in Fig.15c.

(4)In the intermediate zone,the adjacent edge points of neighboring zones,as well as the closest base point,are taken as the nodes of the new elements,as shown in Fig.15d.

In addition,the following two cases should be considered to avoid malformed elements:

(1)The edge point is too close to the base line,in which case the original element is not actually cracked.Therefore,the new element is quite diff i cult to generate,or malformed elements might be generated under this condition.In order to optimize the elements,the original element should be treated as a cracked one.

(2)The base point near the crack tip is too close to the edge line,as shown in Fig.16a.Again,it is diff i cult to generate newelements. In this case the auxiliary line should be def i ned to enlarge the cracked zone near the tip.

Fig.13.Base line of the crack.(a)θ0＜30°,(b)θ0≥30°.

Fig.14.Five zones of the edge points.

Fig.15.The methods used to generate new elements in different zones.(a)Zones I and V,(b)Zones II and IV,(c)Zone III,(d)Intermediate zone.

Fig.16.The base point is too close to the edge line.

6.3.Selecting cracked elements

When the original element is far away from the scanning plane, the element is not cracked,as shown in Fig.17a.On the other hand, the element is considered to be cracked when the scanning plane crossesovertheelementplane,asshowninFig.17b.Thecoordinates of the nodes would be recorded and taken as edge points.Specifically,whenthenodeslieonthescanningplane,orquiteclosetothe scanning plane,they should be deleted,as shown in Fig.17c and d.

Both the scanning plane and the original element plane consist of several straightline segments.Therefore,the relationship between the scanning plane and the original element plane can be determined by the line segments.

As shown in Fig.18,the endpoint coordinates of the line segments are(x1,y1),(x2,y2),(x3,y3)and(x4,y4).The following equation set can be solved to get the intersection points(Shi,1988):

where r and R are both constants.

If 0≤r≤1 and 0≤R≤1,the line segments intersect with each other,as shown in Fig.18a.Otherwise,there is no intersection,as shown in Fig.18b.

6.4.Generating mapping points

Fig.17.Four types of elements.

Fig.18.Relationship between two line segments.(a)Real intersection,0≤r≤1 and 0≤R≤1.(b)False interaction,r＜0,or r＞1,or R＜0,or R＞1.

Fig.19.Generating mapping point.

The direction angle of the vector in global coordinate can be calculated as follows:

6.5.Coalescence of cracks

As shown in Fig.20,three cases of coalescence should be considered,i.e.coalescence at the trunks,tips and intersections, respectively.The predicted position should be adapted to the new one.The intersection points can be obtained according to Eqs.(35) and(36).

7.Case study

The most essential feature of frost damage is the evolution of a crack network,which is a complex process involving coupling of thermal,hydraulic and mechanical f i elds.As a result,proper simplif i cations are necessary.In this study,the following assumptions are made:

(1)Both the ice and the rock material satisfy the features of a Mohr—Coulomb material and their thermal conductivities are isotropic.

(2)Leakage of water during the freezing process can be neglected. The frozen ratio u(T),which represents the percentage of water that has been transformed into ice,satisf i es the following relationship with temperature T(°C):

(3)The freezing point of water is assumed to be 0°C and the supercooling stage of phase transition is neglected,suggesting that the freezing process starts as soon as the temperature drops to the freezing point.The volume expansion of water is 9%when the temperature drops from 0°C to-20°C.Using the method of equivalent thermal expansion coeff i cient to simulate the expansion of the ice,we obtain a value for the thermal expansion coeff i cient of ice of-0.0015°C-1.The negative sign stands for expansion rather than shrinking of the ice volume when the temperature drops(Kang et al.,2012b). The effective SIFs of crack tips are calculated using a computing procedure.The direction and length of crack propagation are def i ned by Eqs.(15)and(23),respectively.

7.1.Modeling and calculating procedure

This paper uses ANSYS for modeling and meshing,and then the model is imported into FLAC3Dby a conversion procedure.The modeling and calculating processes are shown in Fig.21.

Using the above-mentioned method,the freezing model in FLAC3Dis generated as shown in Fig.22.The model contains three cracks,two of which are parallel to each other and perpendicular to the third one.The size of the block is 2.0 m×2.0 m×0.1 m.The length of the cracks is 0.5 m,and the open width is 2.0 mm.The cracks are initially water saturated.The main physico-mechanical and thermal parameters of the model are listed in Table 1.

Interface elements are generated in the gaps between the ice layer and crack faces.Parameters of the ice,such as cohesion and internal friction angle,are modif i ed by a FISH procedure during calculation.Such parameters are roughly assumed to satisfy the following equation:

where F0and Feare the values of the parameters when the temperature T equals 0°C and-20°C,respectively.

Fig.20.Crack coalescence.(a)Trunk coalescence,(b)Tip coalescence and(c)Intersection coalescence.

Fig.21.Flowchart for modeling and calculation.

7.2.Calculation results

A vertical stress of 1.0 MPa and a horizontal pressure of 0.5 MPa are applied to the model.According to the primary assumptions, the crack initiation angle and propagation length are def i ned by the theoretical analyses in Sections 2 and 3.KICis assumed to be 4.0×105Pa m1/2.A FISH procedure is modif i ed to calculate the effective SIFs of each crack tip and to judge whether crack initiation starts.Furthermore,the angle and length of crack propagation are calculated.According to the new algorithm proposed in Section 6, the coalescence is demonstrated as shown in Fig.23.

Another crack tip then propagates without coalescence;details of the element update are shown in Fig.24.The separated new elements in the cracked region are shown in Fig.25.

The f i nal model after the f i rst freezing process is shown in Fig.26.

Fig.22.The pre-existing cracks.

Table 1The main physico-mechanical and thermal parameters of the model.

Fig.23.Coalescence of cracks.(a)Before coalescence,(b)After deleting cracked elements.

Fig.24.Contrast of original and new elements.(a)Original elements,(b)New elements.

Fig.25.Detail of new elements in the cracked zone.

Fig.26.The cracks after element updates.

8.Discussion and conclusions

Frost damage is an important issue in rock weathering.The most essential feature of frost damage is crack propagation caused by temperature f l uctuations around the freezing point of water. Considering the freezing-induced process in a rock mass,this paper analyzes the freezing crack propagation under compression-shear conditions.The freezing-induced cracking under compressionshear is considered to be a mixed mode of tensile cracks(Mode I) and shear cracks(Mode II).The expressions of Mode I and Mode II SIFs are given to predict the growth path of the freezing crack. Three coalescence patterns of adjacent frost cracks are considered, i.e.tensile coalescence,shear coalescence and a mixed pattern of tensile and shear coalescence.

A new algorithm for frost crack evolution is proposed,which could def i ne the crack growth path and identify and update the cracked elements.The algorithm can cope with the evolution of multiple two-dimensional cracks.Further efforts should be made to generalize the method for three-dimensional cracks.

Finally,a numerical example is given to show the growth path of three freezing cracks under compression-shear conditions.The model is generated by ANSYS and is then imported into FLAC3Dbya conversion procedure.The method utilizing the equivalent thermal expansion coeff i cient is applied to simulate the expansion of ice when frozen.A FISH procedure is modif i ed to calculate the effective SIFs of each crack tip and to assess whether crack initiation starts, and both the angle and length of crack propagation are def i ned. Using this new algorithm,the cracked elements are continually updated.Theworkis completedautomaticallyonce the cycle step is specif i ed.The model contains only three cracks and the result proves that the new crack-evolving algorithm is prof i cient in updating the cracked elements.However,when hundreds or thousands of cracks are planted in the model,the computing process could be more complicated,even though the problem could still be solved by the algorithm.Of course,the procedure would have to be modif i ed and a larger matrix built to store details of each crack tip.

This research makes some simplif i cations and assumptions to study the crack propagation mechanism under freezing conditions. Real freezing cracks in rock masses can be much more complicated and,since it is hard to achieve a model that ref l ects the real scenario,a further study is needed.

Conf l ict of interest

We wish to con fi rm that there are no known con fl icts of interest associated with this publication and there has been no signi fi cant fi nancial support for this work that could have in fl uenced its outcome.

Acknowledgments

All authors gratefully acknowledge the f i nancial support from the National Natural Science Foundation of China(Grant Nos.41302237 and 41130742),and the State Key Development Program for Basic Research of China(Grant No.2014CB046900).

China Aviation Academy.The manual of stress intensity factors.Beijing:Science Press;1993(in Chinese).

Erdogan F,Sih GC.On the crack extension in plates under plane loading and transverse shear.Journal of Fluids Engineering 1963;85(4):519—25.

Hilbich C.Time-lapse refraction seismic tomography for the detection of ground ice degradation.The Cryosphere 2010;4:243—59.

Kang YS,Liu QS,Zhao J,Zhang F.Research on frost deformation characteristics of rock and simulation of tunnel frost deformation in cold region.Chinese Journal of Rock Mechanics and Engineering 2012a;31(12):2518—26(in Chinese).

Kang YS,Liu QS,Shi K,Liu XY.Research on modeling method for freezing tunnel with fractured surrounding rock.Advanced Materials Research 2012b;455—456:1591—5.

Kang YS,Liu QS,Huang SB.A fully coupled thermo-hydro-mechanical model for rock mass under freezing/thawing condition.Cold Regions Science and Technology 2013;95:19—26.

Matsuoka N.Mechanisms of rock breakdown by frost action:an experimental approach.Cold Regions Science and Technology 1990;17(3):253—70.

Matsuoka N.Frost weathering and rockwall erosion in the southeastern Swiss Alps: long-term(1994—2006)observations.Geomorphology 2008;99(1—4):353—68.

Matsuoka N,Murton J.Frost weathering:recent advances and future directions. Permafrost and Periglacial Processes 2008;19(2):195—210.

Neaupane KM,Yamabe T,Yoshinaka R.Simulation of a fully coupled thermo-hydromechanical system in freezing and thawing rock.International Journal of Rock Mechanics and Mining Sciences 1999;36(5):563—80.

Ning YJ,An XM,Ma GW.Footwall slope stability analysis with the numerical manifold method.International Journal of Rock Mechanics and Mining Sciences 2011;48(6):964—75.

Paterson MS,Wong TF.Experimental rock deformation—the brittle f i eld.2nd ed. Berlin-Heidelberg:Springer-Verlag;2005.

Seto M.Freeze-thaw cycles on rock surfaces below the timberline in a montane zone:f i eld measurements in Kobugahara,Northern Ashio Mountains,Central Japan.CATENA 2010;82(3):218—26.

Shi GH.Discontinuous deformation analysis—a new numerical model for the statics and dynamics of block systems.PhD Thesis.Berkeley:Department of Civil Engineering,University of California;1988.

Stakhovsky IR.A multifractal model of crack coalescence in rocks.Physics of the Solid Earth 2011;47:371—8.

Su XT,Yang ZJ,Liu GH.Finite element modelling of complex 3D static and dynamic crack extension by embedding cohesive elements in abaqus.Acta Mechanica Solida Sinica 2010;23(3):271—82.

Vásárhelyi B,Bobet A.Modeling of crack initiation,propagation and coalescence in uniaxial compression.Rock Mechanics and Rock Engineering 2000;33(2):119—39.

Yang H.Study on fracture mechanism of multi-crack in rock mass under water-rock interaction.PhD Thesis.Changsha:Central South University;2010(in Chinese).

Yang ZJ,Chen JF.Fully automatic modelling of cohesive discrete crack propagation in concrete beams using local arc-length methods.International Journal of Solids and Structures 2004;41(3—4):801—26.

Yin SZ.Fracture and damage theories and their application.Beijing:Tsinghua University Press;1992(in Chinese).

Zhuang N.Propagation rule of crack rock mass and its mathematics model research under the seepage-stress coupling.PhD Thesis.Shanghai:Tongji University; 2006(in Chinese).

*Corresponding author.Tel.:+86 27 87199185.

E-mail address:kys85@163.com(Y.Kang).

Peer review under responsibility of Institute of Rock and Soil Mechanics,Chinese Academy of Sciences.

1674-7755?2014 Institute of Rock and Soil Mechanics,Chinese Academy of Sciences.Production and hosting by Elsevier B.V.All rights reserved.

http://dx.doi.org/10.1016/j.jrmge.2014.05.004

Freezing pressure

Stress intensity factor(SIF)

Crack propagation