Numerical modelling of pressurized fracture evolution in concrete using zero-thickness interface elements

Numerical modelling of pressurized fracture evolution in concrete using zero-thickness interface elements

Engineering Fracture Mechanics 77 (2010) 1386–1399 Contents lists available at ScienceDirect Engineering Fracture Mechanics journal homepage: www.el...

764KB Sizes 3 Downloads 50 Views

Engineering Fracture Mechanics 77 (2010) 1386–1399

Contents lists available at ScienceDirect

Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech

Numerical modelling of pressurized fracture evolution in concrete using zero-thickness interface elements J.M. Segura, I. Carol * ETSECCPB (School of Civil Engineering), UPC (Technical University of Catalonia), Jordi Girona 1-3, E-08034 Barcelona, Spain

a r t i c l e

i n f o

Article history: Received 2 April 2009 Received in revised form 1 March 2010 Accepted 11 March 2010 Available online 17 March 2010 Keywords: Quasi-brittle materials Plain concrete Finite element analysis Hydro-mechanical coupling Hydraulic fracture

a b s t r a c t The development of cracks due to the effect of fluid pressure is a problem that concerns many areas of engineering, ranging from structural to geotechnical or petroleum. Within the context of the Finite Element Method, the authors have recently proposed a formulation for the coupled hydro-mechanical behaviour of zero-thickness interface elements. This formulation has been verified for pre-existing discontinuities (e.g. natural joints, faults in rock). In this paper, the above formulation, complemented with an appropriate fracture mechanics-based constitutive model, is applied to developing cracks in plain concrete. The numerical results are compared with a series of wedge splitting tests available in the literature, performed on concrete specimens under the influence of fluid pressure at the notch and along the propagating crack. A good agreement is obtained in terms of wedge-splitting force vs. crack mouth opening displacement (CMOD), crack and fluid fronts vs. CMOD, and fluid pressure along the crack vs. time. The influence of splitting rate and input fluid pressure is also systematically analyzed. Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction Coupled hydro-mechanical (HM) fracture processes are important in many fields of engineering (e.g. dams, oil or environmental). In these situations, the pressure of the fluid invading a discontinuity promotes its opening and propagation, which at the same time influences the flow behaviour, through changes in the permeability and the fluid storage capacity of the crack. Hydraulic fracturing in petroleum engineering is a classical example of a HM fracture process. This consists in the initiation and propagation of cracks from the wellbore into the oil reservoir rock, by injecting a highly pressurized fluid that drives the cracks. Such discontinuities are highly permeable paths for the depletion of the reservoir, and they considerably improve the production of the well. Fluid-pressurized cracks are also created in other engineering activities such as the measurement of in situ stresses in rock masses, the creation of impermeable barriers for contaminant containment, the injection of slurry wastes, or in grouting operations. They can also naturally appear in magma-driven fractures or in structures subjected to high water pressure levels like dams. 1.1. Numerical approaches for fracture simulation The numerical approaches, based on the Finite Element Method (FEM), that model the creation and propagation of cracks in quasi-brittle materials, have been mainly classified into the ‘‘smeared crack approach” and the ‘‘discrete crack approach” [1], although more recent developments include other alternatives [2–4]. * Corresponding author. Tel.: +34 93 401 65 09; fax: +34 93 401 72 51. E-mail address: [email protected] (I. Carol). 0013-7944/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.engfracmech.2010.03.014

J.M. Segura, I. Carol / Engineering Fracture Mechanics 77 (2010) 1386–1399

1387

Nomenclature

rm r0m

Wcr

fracture mid-plane total stress vector fracture mid-plane effective stress vector fracture constitutive relationship tensor fracture walls relative displacements vector fracture aperture flux vector along the fracture mid-plane fracture mid-plane fluid pressure fracture longitudinal transmissivity tensor fracture Biot’s modulus fracture Biot’s coefficient fracture mid-plane elevation fluid specific weight fluid dynamic viscosity nodal displacements vector nodal fluid pressure vector interface element length time interface element B-matrix interface element coupling matrix interface element conductivity matrix interface element capacity or storage matrix effect of gravity potential on the flow system nodal forces vector nodal fluid flows vector interface element tensile strength work spent on fracture

GIf

fracture energy in Mode I

GIIa f

fracture energy in Mode IIa

KMn a EM

normal stiffness of the interface element representing the rubber membrane rubber membrane length rubber membrane Young’s modulus

Dm A An jl pm Tl Mj

aj zm

cf lf u p l t B L E S Qz fu fp

v0

The smeared crack approach implies a continuum-type representation of the problem with, generally, a fixed FE mesh. It considers an equivalent material where the influence of the opening cracks is incorporated into the constitutive stress/strain law with softening. Macroscopic structural cracks can be identified a posteriori as collections of cracked Gauss points sufficiently aligned [5,6]. In the discrete crack approach each single discontinuity is represented explicitly. Any line of the mesh can potentially become a discontinuity and, in that case, the continuum on each side of the mesh can exhibit different displacement values, leading to the discontinuity opening or sliding. If discontinuities can only develop along mesh lines, the layout of the mesh lines obviously poses an important geometric constraint to a general crack path, unless the mesh is really fine or it can be changed as cracks develop. This has led to two strategies of discrete crack implementation: with evolving mesh and with fixed mesh. 1.2. Discrete crack models in Non-Linear Fracture Mechanics From the mechanics viewpoint, there are two different approaches to analyze fracture processes, one corresponding to Linear Elastic Fracture Mechanics (LEFM), and the other to Non-Linear Fracture Mechanics, normally in the form of Cohesive Crack/Fictitious Crack Model, all expanded next. Models based on LEFM were first developed for general crack propagation problems, by the use of remeshing [7,8]. Such an approach does not consider the description of a fracture process zone (FPZ) ahead of the crack tip, where material behaviour is non-linear and energy is dissipated. In LEFM, the FPZ is actually reduced into a zero-size crack tip. Except for very large structures, where the size of the FPZ is negligible and LEFM is applicable, consideration of a FPZ with dimensions related to the size of micro-structural heterogeneities turns out to be very important in quasi-brittle material analyses [9]. This motivated the development of Non-Linear Fracture Mechanics (NLFM), which considers the existence of the FPZ. Usually, in such an approach, linear elasticity is used for the mechanical response of the continuum material, whereas the non-linearity is included in the formulation along the crack line in the form of a softening law. This softening function relates the cohesive stress (stress acting across the crack faces) to the Cohesive Crack opening (relative displacement between the upper and

1388

J.M. Segura, I. Carol / Engineering Fracture Mechanics 77 (2010) 1386–1399

lower faces of the interface). The integral of the softening function (or graphically, the area enclosed under the softening curve) corresponds to the fracture energy under Mode I ðGIf Þ [10,11]. The NLFM approach, initially conceptualized for Mode I fracture analysis, has also been adapted to mixed-mode fracture with normal/shear coupled interface laws [1,12,13] using the so called ‘‘zero-thickness” interface elements, which were originally proposed to describe the mechanical behaviour of pre-existing joints in rock masses [14]. In a general situation, such a model should be able to simulate both the initiation and propagation of fractures under Mode II or mixed mode loading, and also the influence of joint dilatancy. Most of such ‘‘zero-thickness” interface models are based on the theory of elasto-plasticity incorporating fracture mechanics concepts. Although the Coulomb criterion has been traditionally used for its simplicity, other yield curves have been proposed, as parabolic [15] or hyperbolic [1,13]. Most models consider a non-associated flow or slip rule to describe the plastic dilatancy. The evolution of the yield surface is controlled by some internal variables that vary according to a set of softening rules. Most of the models consider the evolution of the internal variables as a function of the interface fracture work [1,15] or the inelastic displacement [13]. 1.3. Discrete crack models for HM fracture processes HM coupled extensions of the discrete crack approach for the simulation of fluid-pressurized cracks based exclusively in the FEM, are not common in the literature (a review can be found in [16]). Most of the existing formulations are based on partitioned schemes that combine the finite differences (FD) method for the flow problem with the FEM for the mechanics one, either considering some fluid-exchange between the hydraulic fracture and the porous material [17] or no leak-off [18]. Fully-coupled HM formulations based only on the finite differences method have been recently proposed [19]. A fully coupled FEM formulation has been proposed by Simoni and Secchi [20] and Secchi et al. [21], which remeshes the domain according to fracture mechanics principles as the fracture nucleates and propagates. The continuum finite elements are formulated according to poro-elasticity, and the exchange of fluid between the cohesive fracture and the surrounding porous material is considered. Slowik and Saouma [22] combine a Fictitious Crack Model [13] with a transient diffusion formulation, which takes into account the interaction between the penetrating fluid into the crack, and the air in both the initially dry concrete (assumed impermeable) and the propagating fracture. The authors have recently presented a numerical model that describes the coupled HM behaviour of discontinuities using ‘‘zero-thickness” interface elements with double nodes [16]. This model has been obtained by coupling the diffusion formulation of discontinuities proposed earlier [23] with a mechanical formulation of interface elements that can model both preexisting geomechanical discontinuities (e.g. joints or faults) [24,25], and developing cracks [12,26–31], provided that in each case the appropriate constitutive model is used (i.e. rock mechanics based model for pre-existing joints, while fracture mechanics based model for developing cracks). The formulation in [16] considers poro-elasticity in the continuum and allows the exchange of fluid between the discontinuity and the surrounding (porous) material. Each finite element (continuum and interface) is formulated in terms of the displacements (u) and the fluid pressure (p) at the nodes. Two different numerical approaches are presented in [16] to solve the system of coupled equations: staggered and fully coupled. The staggered procedure reorganizes the system of coupled equations so that it can be solved by iterating between a standard mechanical and flow simulators, provided that the appropriate coupling loops between codes are introduced. On the other hand, the fully coupled approach solves the system of equations as a whole, using the Newton–Raphson method, and it is shown that computation of the full Jacobian matrix can considerably improve the convergence of this method [32]. In the present article, the same all-FEM fully-coupled HM formulation is combined with a fracture mechanics-based constitutive model for the interfaces, and the resulting model is used to simulate the propagation of pressurized cracks in concrete specimens. The numerical results are compared with existing well-documented experimental results [22,33,34]. The content of the paper is organized as follows: after this introduction, Section 2 summarizes the main characteristics of the HM coupled model and the interface constitutive model. Section 3 shows the capabilities of the formulation by closely reproducing the experimental results by Brühwiler and Saouma [33,34] and Slowik and Saouma [22]. Section 4 summarizes the concluding remarks. 2. HM model using zero-thickness interface elements The HM coupled formulation of the zero-thickness interface elements was presented in detail in Ref. [16] and it is briefly outlined in this section. Poro-elasticity is considered for the continuum in which the crack is developing, and the non-linear behaviour is restricted to the zero-thickness interface elements representing the discontinuity. After assembly of the continuum and interface elements, a particular form of the well-known ‘‘u–p” FE formulation in geomechanics [35] is obtained, in which the main unknowns are the displacements (u) and fluid pressure (p) at the nodes. 2.1. Governing equations Following standard procedures in the formulation of zero-thickness interface models [24], the governing equations are expressed first in terms of jump variables at the mid-plane of the discontinuity, and then these variables are related to

J.M. Segura, I. Carol / Engineering Fracture Mechanics 77 (2010) 1386–1399

1389

the absolute values of the variable on the discontinuity walls, in order to obtain the discretized FE equations for the interface element. The mechanical formulation of a saturated crack under the influence of an internal fluid pressure distribution is developed in its mid-plane through consideration of the effective stress principle (Eq. (1)), and the constitutive model (Eq. (2)), which links the stresses in its mid-plane to the corresponding relative displacements of the fracture walls:

dr0m ¼ drm þ aj mm dpm 0 m

dr ¼ Dm dA

ð1Þ ð2Þ

where rm and r0m are the mid-plane total and effective stress vectors respectively, aj is the Biot’s coefficient for the joint, mm is a vector introducing the influence of fluid pressure in the normal direction to the discontinuity, pm is the fluid pressure at the mid-plane, Dm is the constitutive relationship tensor and A is the relative displacements vector. The formulation of the flow of fluid along a deforming discontinuity is obtained by combining the fluid mass balance equation (Eq. (3)) with the general form of Darcy’s law (Eq. (4)) along the mid-plane:

1 @pm @An þ aj ¼0 M j @t @t ! rpm þ rzm jl ¼ Tl

rjl þ

cf

ð3Þ ð4Þ

where t stands for time, An is the crack aperture, jl is the flux vector along the mid-plane, Mj is the joint Biot’s modulus, Tl is the longitudinal transmissivity tensor, zm is the mid-plane elevation, and cf is the fluid specific weight. Parameters aj and Mj in Eqs. (1) and (3) are understood analogously to the case of porous materials, and they depend on joint properties (e.g. roughness, stress state, presence of infill material) [36]. 2.2. FEM formulation of the zero-thickness interface element Extension of the HM equations from the mid-plane to the interface element nodes is detailed in [16], reaching the following system of coupled equations, that describes the HM coupled behaviour of a zero-thickness interface element with double nodes:

Z l

BT r0m dl  Lp ¼ f

Ep þ LT

u

du dp p þS þ Qz ¼ f dt dt

ð5Þ ð6Þ

where the main unknowns are the vectors of absolute displacements (u) and fluid pressure (p) at the joint element nodes, l stands for interface element length, B is the interface element B-matrix, L is the interface element coupling matrix, E is the interface element conductivity matrix, S is the interface element capacity or storage matrix, Qz introduces the effect of gravity potential on the flow system, and fu and fp are respectively nodal forces and fluid flows to cancel with neighbouring elements during finite elements assembly. Eqs. (5) and (6) are discretized in time using the generalized trapezoidal rule [16,37]. Once the continuum finite and the zero-thickness interface elements are assembled, a system of coupled equations is reached which can be solved according to two different procedures [16]: staggered and fully coupled. The fully coupled approach is used in the numerical simulations presented in the following Section 3, due to its robustness and faster convergence to the solution than that of the staggered approach. 2.3. Constitutive models A general model for normal/shear cracking in quasi-brittle materials originally proposed by Carol et al. [1] has been used in this study. From this model, the constitutive stiffness Dm in Eq. (2) is obtained. The model is based on the theory of elastoplasticity. According to this, crack initiation occurs when the stress state reaches the crack surface in the stress space in Fig. 1a. The model is associated in tension and non-associated in compression (Fig. 1a). The evolution of the crack surface (Fig. 1c) is determined by three independent fracture energy softening laws (Fig. 1d and e) that evolve with the work spent on fracture under Mode I and Mode IIa (Fig. 1b). From the point of view of flow, the cubic law [38] is considered as an estimate for the relationship between the fracture transmissivity tensor Tl, in Eq. (4), and the fracture aperture. 3. Combined effect of internal water pressure and external mechanical load for crack propagation in concrete The experimental results used to validate the numerical model presented in Section 2 were performed and reported by Saouma and co-workers at the University of Colorado at Boulder [22,33,34]. These experiments consist of an extensive series

1390

J.M. Segura, I. Carol / Engineering Fracture Mechanics 77 (2010) 1386–1399

(a)

(b)

(c)

(d)

c, χ

(e)

tan φ

c0

tan φ0

=2

1

=2

χ0

0

1

-1

0

=2

1

-2

tan φr

-1

0 -2

-1 -2

Wcr

Wcr

GIf

GIIa f

GIIa f

Fig. 1. (a) Cracking and plastic potential surface and (c) their evolution; (b) modes of fracture; (d–e) softening laws [1].

of wedge splitting tests on concrete specimens under the influence of water pressure at the notch and along the propagating crack. The tests in Brühwiler and Saouma [33,34] were performed under quasi-static conditions, and they analyze the influence of fluid pressure on the degradation of the fracture properties, as well as the evolution of the fluid pressure build-up in the opening of the crack. These experiments were followed by those performed by Slowik and Saouma [22] using the same experimental setup, but this time under dynamic loading conditions. In such tests, the effect of the wedge splitting rate on the HM performance of the specimen is studied, measuring applied forces/CMOD and fluid pressure distribution along the crack. The tests have their application in dam safety analyses, and try to reproduce the effect of the water column height (imposed mechanical load) and hydrostatic fluid pressure (imposed fluid pressure at the notch) acting on the upstream face of concrete dams. Fig. 2a shows the concrete specimens dimensions, and Fig. 2b sketches the experiment setup. The experiments are controlled via crack mouth opening displacement (CMOD) prescribed by the wedge splitting device. This mechanical apparatus deforms the specimen at a given CMOD rate and records the load applied by the loading cell that

Load Cell FLC CMOD

102 mm

FH

FV

FM

σM

83 mm

FLC

pn

288 mm

FW 15O

p FH

154 mm

FR

FH FV

FV

MR 304 mm

FLC

(a)

(b)

(c)

Fig. 2. (a) Concrete specimen dimensions; (b) experimental setup (after [33]); (c) decomposition of the load applied by the cell.

1391

J.M. Segura, I. Carol / Engineering Fracture Mechanics 77 (2010) 1386–1399

is necessary to reach a given CMOD value. If a pressurized fluid is introduced in the notch, the fluid tends to invade the opening crack and to pressurize it as well. The resulting fluid pressure distribution acting along the notch and the opening discontinuity helps the external mechanical load to deform the specimen and to propagate the crack. Therefore, the external load applied to reach a certain CMOD value is lower than that in a dry specimen if the fluid pressure is acting within the crack. Also, for higher crack mouth opening rates the fluid pressure distribution cannot build-up so easily within the discontinuity, and higher splitting forces are required. These and other factors are dealt with in Refs. [22,33,34] and numerically reproduced in this paper. A hydraulic loading device maintains the water under the desired pressure inside the notch with the aid of a rubber membrane glued to the surface of the specimen, which prevents leakage of water from the notch and the developing crack. According to the experimental procedure described in Refs. [22,33,34], three different forces contribute to the splitting of the concrete specimen as shown in Fig. 2b: 1. Mechanical load imposed by the wedge: the force FLC, measured by the loading cell, divides into a mechanical splitting (horizontal) FH and a vertical FV component at the point of contact between the wedge and the rolling bearings (Fig. 2c). 2. Hydraulic load FW promoted by the imposed fluid pressure at the notch (pn) and the pressure (p) of the fluid that flows into the developing crack. 3. Reaction force due to stretching of the rubber membrane FM, which opposes to the splitting of the specimen. The influence of these components on the results is explained more in detail for each experiment in the following two subsections. Fig. 3a shows the mechanical boundary conditions applied in the numerical simulations: horizontal and vertical displacements are not allowed along the central bottom line of the specimen (point B in Fig. 3a). The opening displacement is imposed 12 mm below the top specimen face (points A and A0 in Fig. 3), whereas the CMOD is monitored 28.5 mm over the top specimen face. From a flow point of view, the fluid pressure is prescribed along the notch, which in the experiments is controlled by a servovalve. Fig. 3b shows the finite element mesh used for the numerical simulations. The continuum is discretized with 688 nodes and 998 linear triangular elements, with a finer discretization in the vicinities of the interface elements and of points A and A0 , where the splitting force is applied. The mesh is extended 28.5 mm over the top specimen to compute the numerical CMOD as the increment of distance between points C and C0 . Three different types of interface elements are used both along the notch and the expected crack path (i.e. specimen symmetry axis). At the upper notch tip, a small elastic interface element reproduces the effect of the rubber membrane at the mouth of the notch (I1 in Fig. 3b). The notch is also represented via elastic interface elements (I2 in Fig. 3b), although these are assigned a much lower stiffness than the top interface I1 as explained later on. 125 interface elements, governed by the fracture mechanics model [1] cited in Section 2.3, are inserted along the expected crack path (I3 in Fig. 3b). The continuum material is modelled linear elastically and the non-linearity of the system is captured via the interface elements. The coupling between fluid pressure and displacements field is considered at the interface element level only under the main assumption that concrete is almost impermeable in comparison with the opening crack at the time range

28.5 mm C

C’

A

A’

CMOD measurement

A

A’

12 mm I1 I2

I3

B

(a)

(b) Fig. 3. (a) Boundary conditions; (b) mesh used in the analyses.

1392

J.M. Segura, I. Carol / Engineering Fracture Mechanics 77 (2010) 1386–1399

analyzed. From the mechanical point of view, a Poisson ratio of 0.20 is used in all the cases, and the value of the Young modulus is specific for each experiment. The cubic law [38] is used for both interface elements I2 and I3. The model employed [16] incorporates the use of Biot’s modulus for the discontinuity, MJ (see Eq. (3)). In this particular case, the discontinuity is partially saturated at some length and a value for this modulus MJ = 0.0002 MPa/mm, which is lower than the corresponding one for a saturated discontinuity, has been used in all the analyses. This avoids having to impose a vapour pressure at the crack tip as boundary condition for the fluid flow along the fracture. From the mechanical point of view, interface element I1, representing the rubber membrane, is considered linear elastic with a normal stiffness KMn (MPa/mm). Taking into account that the length of the membrane (a) is approximately the thickness of the notch, KMn is related to the rubber material Young’s modulus (EM) through:

r0Mn ¼ K Mn Da ¼ EM Da=a $ EM ¼ K Mn a

ð7Þ

where r0Mn is the normal stress acting on the membrane. This is a simplification of the real behaviour of the rubber membrane in the experiments, since such a membrane gets curved due to the pressure applied by the fluid, and the curvature varies due to the elongation of the membrane as the specimen progressively deforms. Elements I2 representing the notch are linear elastic, and if the specimen were infinitely wide out-of-plane (i.e., the numerical analysis is made under plane strain) these interface elements should exhibit no resistance to deformation, which can be modelled using a nil normal stiffness. However, a certain resistance to crack growth should be incorporated to these interface elements, since the notch is mechanically influenced by the rubber membrane located at the rear and front lateral faces of the specimen. It could be argued that the effect of the membrane should not be significant since in the experiments it was left sufficiently loose so that the fluid pressure would fill in and the membrane would take e a curved approximately semi-cylindrical shape. However, this does not guarantees the absence of tangential reaction at the points (or better lines) of attachment. Since the accurate analysis of the system membrane-fluid is not simple, the first approach has been rather to explore the influence on the results of introducing a membrane stiffness, and seek an inverse calibration of the corresponding value. This procedure has been verified by validating first the model with experiments without pressure or membrane, which has turned out very accurate. Then, the membrane stiffness has been introduced for the calculations with pressure, and, indeed the influence of the membrane on the numerical performance of the system has been found to be important, it cannot be neglected. By varying the thickness and the normal stiffness of the membrane element I1 (with the corresponding normal stiffness for the interface elements I2 specified in each case), very similar numerical results have been obtained. A first group of studies has been performed considering a nil stiffness for the elements I2 along the notch, and therefore joint element I1 concentrates all the membrane resistance to the splitting of the specimen. Three different lengths for the interface element have been considered, which reproduce the membrane thickness. Shorter elements (i.e. thinner membranes) require higher values for the normal stiffness to maintain a certain CMOD. Table 1 shows the values of the interface element normal stiffness KMn adjusted to fit the numerical and the experimental results, and the equivalent range of the membrane Young’s modulus EM according to Eq. (7), considering a notch width of 1 mm or 10 mm. The range obtained for EM in this case gives high values in comparison with the standard Young’s modulus found in the literature for rubber membranes [39], which is of the order of 1–10 MPa at the most, and for rubber material that is usually below 100 MPa. On the other hand, the same normal stiffness could be considered along the entire notch, which would be a too conservative approach. In this case, a normal stiffness of KMn = 0.4 MPa/mm is needed to reproduce the numerical results, which is equivalent to use a rubber membrane with EM = 0.4 MPa in the case of a = 1 mm, or EM = 4 MPa in the case of a = 10 mm. A more reasonable option seems to consider a reduced stiffness along the notch in comparison with the stiffness of the interface element I1. Table 2 shows the normal stiffness for interface element I1 for different interface element lengths. In all the cases, the normal stiffness in I2 is less than 5% the stiffness in I1. If the influence of interface elements I2 is reduced, then

Table 1 Membrane Young modulus considering free deformation of the notch. Membrane thickness (i.e. interface length) (mm)

KMn (MPa/mm)

EM (MPa) a = 1 mm

EM (MPa) a = 10 mm

0.2 1 2

90 20 9

90 20 9

900 200 90

Table 2 Membrane Young modulus considering a certain constraint deformation of the notch. Membrane thickness (i.e. interface length) (mm)

KMn (MPa/mm)

EM (MPa) a = 1 mm

EM (MPa) a = 10 mm

0.2 1 2

10 8 5

10 8 5

100 80 50

1393

J.M. Segura, I. Carol / Engineering Fracture Mechanics 77 (2010) 1386–1399

the stiffness for I1 would tend to the values in Table 1, whereas if the influence of interfaces I2 is increased, the normal stiffness for I1 and I2 would be more similar and would tend to the value of KMn = 0.4 MPa/mm previously specified. The constitutive model of Carol et al. [1] introduced in Section 2.3 is considered in interface elements I3 shown in Fig. 3b, which really capture the propagation of the crack. Although the model can reproduce both I and II fracture Modes, the problem under analysis is a pure fracture Mode I case (Fig. 1b). Then, at a given Gauss point of an interface element the effective normal stress will increase until it reaches the interface element strength, v0 . This increase is made along the rn axis in Fig. 1c via an elastic normal stiffness of the interface, with very high value, which has the meaning of a penalty coefficient. Once the normal stress reaches the initial tensile strength v00 , the initial cracking surface denoted as 0 in Fig. 1c starts softening towards the final cracking surface, denoted as 1 in Fig. 1c. At this final stage under pure tension, the interface element does not transfer any stress between the two joint faces. The evolution between states 0 and 1 is governed by the softening law shown in Fig. 1d, which relates the normal effective stress (initially equal to v00 ) to the work spent on fracture Wcr. This analysis is under tension and the effective stresses are always positive according to the sign convention used in this work. The characteristics of the fluid (water) are specific weight (cf) equal to 1  105 MPa/mm and dynamic viscosity (lf) equal to 1  109 MPa s. The effect of gravity is neglected under the assumption that the effect of water weight on the water flow is negligible compared to the reservoir pressure head [22]. As a matter of fact, the tests try to reproduce the effect of water loads on concrete dams. 3.1. Effect of the mechanical loading rate: comparison with Slowik and Saouma experiments Two groups of wedge splitting tests are considered: one set applies a slow CMOD rate of 2 lm/s, and the other a fast rate of 200 lm/s. As explained by the authors [22], this experimental procedure allows them to analyze the influence of the crack opening rate on the pressure distribution in the crack, and the velocity of the pressure build-up within the discontinuity in comparison with the crack propagation. The fluid pressure distribution within the developing discontinuity will clearly condition the external mechanical load that has to be applied to reach a given CMOD. The mechanical and flow parameters for the continuum and the zero-thickness interface elements have the following values. The Poisson ratio is 0.2 for both slow and fast loading cases, whereas Young’s modulus are 20,000 MPa and 29,000 MPa for the slow and fast opening rate cases, respectively. This difference in the elastic moduli is to capture the higher pre-peak slope that the fast loading case shows in the load-CMOD curves of Fig. 5. The reason for the difference in initial slope could correspond to the effect of the loading rate, which if increased usually produces an increase of the specimen strength and elasticity modulus. Zero-thickness interface elements that capture the developing crack (I3 in Fig. 3b) are assigned a fracture energy of GIf ¼ 0:125 MPa mm for the slow loading case and GIf ¼ 0:1535 MPa mm for the fast one. The initial interface tensile strength is v00 = 3.05 MPa. The resulting cohesive functions for the joint elements are compared in Fig. 4 with those used by Slowik and

3.5 coh. function used

3

sigma'_n (MPa)

sigma'_n (MPa)

3.5

coh. function [5]

2.5 2 1.5 1 0.5 0

coh. function used

3

coh. function [5]

2.5 2 1.5 1 0.5 0

0

0.05

0.1

0.15

0.2

0.25

0.3

0

0.05

crack opening (mm) 0

0.05

0.1

0.15

0.2

0.15

0.2

0.25

0.3

0.25

0.3

crack opening (mm) 0.25

0.3

0

0.05

0.1

0.15

0.2

1.0E+05

1.0E+00 1.0E-05 1.0E-10 conductivity used

1.0E-15

conductivity [5]

1.0E-20

(a)

conductivity (mm2/s)

1.0E+05

conductivity (mm2/s)

0.1

1.0E+00 1.0E-05 1.0E-10 conductivity used

1.0E-15

conductivity [5]

1.0E-20

(b)

Fig. 4. Cohesive function and hydraulic conductivity curves for (a) slow and (b) fast splitting tests.

1394

J.M. Segura, I. Carol / Engineering Fracture Mechanics 77 (2010) 1386–1399

Saouma [22] in their simulations. Fig. 4 also compares the hydraulic conductivity used in our model (cubic law based) with that used by Slowik and Saouma in their simulations. The results shown in Figs. 5–7 have been obtained with a normal stiffness KMn = 20 MPa/mm and a length of 1 mm for interface element I1 representing the membrane, and a nil stiffness for interfaces I2, although almost identical results were obtained with the other values in Tables 1 and 2. A backward scheme is used for the time discretization. Fig. 5 shows the applied cell load vs. CMOD curves for both slow (a) and fast (b) crack opening rates. From an experimental point of view, two main distinctions appear when comparing both (a) and (b) curves: the pre-peak/peak behaviour and the post-peak response. The pre-peak response takes place in the interval 0 6 CMOD 6 0.15 mm (Fig. 5), and, as shown by the experimental curves in Fig. 6, for these values of CMOD the fluid pressure has not almost developed along the discontinuity, with a gap

experimental fast

4

numerical fast

3.5

experimental slow

Cell load (kN)

3

numerical slow

2.5 2 1.5 1 0.5 0 0

0.2

0.4

0.6

0.8

1

CMOD (mm)

140 120

location in ligament (mm)

numerical crack front crack front [5] numerical fluid front fluid front [5]

160

100 80 60 40 20 0 0

0.2

0.4

0.6

numerical crack front crack front [5] numerical fluid front fluid front [5]

160 140 120 100

0.8

80 60 40 20 0 0

0.2

0.4

CMOD (mm)

CMOD (mm)

(a)

(b)

Fig. 6. Crack and water front vs. CMOD for (a) slow loading and (b) fast loading.

0.25 0.2

340 s 290 s

p (MPa)

location in ligament (mm)

Fig. 5. Cell load vs. CMOD for slow and fast loading.

240 s

0.15

140 s

190 s

0.1 0.05

90 s

0 0

50

100

150

location from notch tip (mm) Fig. 7. Fluid pressure along the ligament for slow loading.

0.6

0.8

J.M. Segura, I. Carol / Engineering Fracture Mechanics 77 (2010) 1386–1399

1395

of about 40 mm between the crack and the fluid fronts in the slow loading case (Fig. 6a), and much higher for the fast loading case (Fig. 6b). Therefore, the effect of fluid pressure in the pre-peak curves is minimal and the different behaviour may be only attributed to mechanical effects of the loading rate. As a matter of fact, the loading rate effect in dry tests affects the concrete strength, showing lower peak loads for slow loading than for fast loading, but giving similar post-peak responses [40]. The effect of fluid pressure along the ligament appears in the post-peak response. Fig. 5 shows that a higher splitting force is needed in the fast loading case than in the slow one to reach the same CMOD value. This is explained by the higher buildup water pressure along the crack in the slow loading specimen that contributes to the deformation of the specimen. Then, if the propagation of the crack is fast enough, the fluid pressure does not have time to develop, and the opposite for a slow crack propagation. This can be clearly observed in Fig. 6, which shows how the gap between the fluid and crack fronts remains approximately constant throughout the slow experiment, but it increases in the fast loading analysis. As observed in Fig. 5, the numerical results reproduce quite accurately the experimental load vs. CMOD response of the specimens for both slow and fast opening rates. The numerical simulations also capture the trend of all the curves in Fig. 6, which show the evolution of the fluid and crack fronts with the CMOD. A reasonable agreement between the numerical and experimental fluid pressure profiles along the discontinuity has been obtained. As an example, the computed and measured fluid pressure profiles for different values of time are compared in Fig. 7 for the slow loading case. 3.2. Effect of the fluid pressure at the crack mouth: comparison with Brühwiler and Saouma experiments In contrast to the previous case, these wedge splitting tests are run under quasi-static conditions. The objective is now to capture the effect of the fluid pressure build-up within the crack on the FPZ. The loading cell and specimen dimensions are the same as in the previous section experiments (Fig. 2). Therefore, the same mesh and problem setup (Fig. 3) is considered for these numerical simulations. In the experiments described in Section 3.1, the load cell vs. CMOD curves incorporated the effect of fluid pressure and membrane stretching. In this case, the specimen is repeatedly pressurized and unpressurized, giving a saw tooth-like curve (Fig. 9), in which the peak points (unpressurization) represent values of load cell and CMOD without fluid pressure in the crack and with no membrane stretching, but tacking into account the indirect effect of the fluid pressure having predamaged the specimen during the pressurization stages. This is the basic objective of these experiments: to evaluate the degradation of fracture properties due to the action of fluid pressure. The tests are run according to different input fluid pressure values at the notch. Comparison between the experimental and numerical results is shown in Figs. 10 and 11. First, adjustment of the mechanical properties for both concrete and interface elements is made according to the experimental curve in Fig. 8, which shows the wedge splitting test for a dry specimen (i.e. no input fluid pressure at the notch). According to the experimental procedure presented in the article, the stretching of the membrane only takes place when the notch is pressurized. The numerical simulations are performed according to the mesh presented in Fig. 3, in which the stretching of the membrane takes place also in dry conditions, which is a simplification of the problem. A proper modelization of the problem would require ‘‘ad/remove” element capabilities in the coupled code, not available at the time. Then, the reference curve is adjusted under the effect of the membrane represented by interface I1 in Fig. 3, and the objective of the numerical simulations is to prove that, effectively, when the specimens are pressurized and depressurized, the trend of the splitting force vs. CMOD curves adjusts to those predicted experimentally, giving progressively lower values of peak force and CMOD, and also progressively more vertical curves after the peak load as the fluid pressure at the notch is increased. Fig. 8 shows how for

experimental numer. with memb. numer. without memb.

5 4.5 4

Fs (kN)

3.5 3 2.5 2 1.5 1 0.5 0 0

0.1

0.2

0.3

0.4

0.5

CMOD (mm) Fig. 8. Splitting force vs. CMOD curve for the test with no water pressure (reference test).

1396

J.M. Segura, I. Carol / Engineering Fracture Mechanics 77 (2010) 1386–1399

Fig. 9. (a) Experimental procedure and (b) experimental results for the test with pressure at the notch 0.5 MPa [33].

5

5

4.5

4.5 4

4 3 2.5

pn = 0.0

2

pn = 0.5

1.5 1

pn = 0.9

0.5

pn = 0.3

pn = 0.1

Fs (kN)

Fs (kN)

3.5

3.5 3 pn = 0.0

2.5 2 1.5

pn = 0.3

1 0.5

pn = 0.7

0

pn = 0.5

pn = 0.9

pn = 0.1

pn = 0.7

0 0

0.1

0.2

0.3

0.4

0.5

0

0.1

0.2

0.3

CMOD (mm)

CMOD (mm)

(a)

(b)

0.4

0.5

Fig. 10. Splitting force vs. CMOD curve for different fluid pressure values at the notch (a) experimental (b) numerical.

a specimen without the stretching of the membrane, the splitting force vs. CMOD would show lower values, and therefore a higher fracture energy should be included in the fracture mechanics model to accurately predict the splitting curve. The same hydraulic and membrane parameters as in previous numerical simulations are used for the seek of coherence. Poisson’s ratio is considered 0.2 also, and Young’s modulus has been adjusted to 19,000 MPa. A fracture energy of GIf ¼ 0:11 MPa mm and an initial tensile strength of v00 ¼ 1:9 MPa are used for the zero-thickness interface model that capture the fracture. Fig. 9 [33] describes the experimental procedure for the case of a water pressure input of 0.5 MPa. The test commences with no fluid pressure at the notch and under CMOD control. At point A, the CMOD is locked and the notch is fluid pressurized. Since water pressure contributes to the deformation of the specimen, the external load necessary to maintain the CMOD is reduced to point B (Fig. 9b). Then the experiment continues deforming the specimen under CMOD control until point C. At this point the CMOD is again locked and the fluid pressure is released. In order to maintain the CMOD, the external load must now increase till point D, since it does not have the help of water pressure. Once depressurized, the notch is again pressurized and the external applied load reduces now to point E. This procedure is repeated throughout the test giving a saw toothlike curve in which the upper points represent the mechanical load necessary to deform a specimen predamaged by water pressure, and the lower points represent the external mechanical load necessary to deform a fluid pressurized specimen. Five series of these tests are made, by considering five different input fluid pressure values: 0.10, 0.30, 0.50, 0.70 and 0.90 MPa. In the numerical simulations, splitting of the specimen is performed at the slow rate introduced in previous section (2 lm/s) and the fluid pressure increase and decrease at the notch is imposed at a rate of 0.005 MPa/s. Also, in the numerical simulations the pressurization and depressurization is made always at the same CMOD values. Such parameters (loading rate, pressurization/depressurization rate and CMOD values at which pressurization/depressurization takes place) consider-

1397

J.M. Segura, I. Carol / Engineering Fracture Mechanics 77 (2010) 1386–1399

ably affect the shape of the splitting force vs. CMOD, since they control the fluid pressure distribution within the crack during pressurization, depressurization and crack propagation. For each input fluid pressure at the notch, four different specimens were experimentally tested in [33], and then, the resulting fluid pressure release points for each test (upper points in Fig. 9) were curve fitted to obtain the splitting force against CMOD curves represented in Fig. 10a. Since at the beginning of the analysis no fracture process zone exists, water pressure cannot penetrate into it and the initial slope of the experimental curves is the same in all the cases examined. Once the fracture process zone has been created, these curves show that as the input fluid pressure increases, both the peak force

5

experimental

5

numerical

4

pn = 0.0 MPa

4

pn = 0.1 MPa

3.5

Fs (kN)

3.5

Fs (kN)

experimental numerical

4.5

4.5

3 2.5 2

3 2.5 2 1.5

1.5 1

1

0.5

0.5 0

0 0

0.1

0.2

0.3

0.4

0

0.5

0.1

experimental numerical

4.5 4

0.4

0.5

experimental

4.5

numerical

4

pn = 0.3 MPa

3.5

pn = 0.5 MPa

3.5

3

Fs (kN)

Fs (kN)

0.3

5

5

2.5 2

3 2.5 2 1.5

1.5 1

1

0.5

0.5

0

0

0

0.1

0.2

0.3

0.4

0

0.5

0.1

5

5

4.5

4.5

experimental

4

0.2

0.3

0.5

experimental

4

numerical

numerical

3.5

3.5

Fs (kN)

pn = 0.7 MPa

3

0.4

CMOD (mm)

CMOD (mm)

Fs (kN)

0.2

CMOD (mm)

CMOD (mm)

2.5 2

2.5 2

1.5

1.5

1

1

0.5

0.5

0

pn = 0.9 MPa

3

0

0

0.1

0.2

0.3

CMOD (mm)

0.4

0.5

0

0.1

0.2

0.3

0.4

0.5

CMOD (mm)

Fig. 11. Numerical simulation of the complete experimental procedure. The experimental curves represent the experimental pressure release points (upper points in the saw-tooth curves).

1398

J.M. Segura, I. Carol / Engineering Fracture Mechanics 77 (2010) 1386–1399

and the corresponding CMOD decreases, giving also a steeper descending branch, the last one related to a reduction of the fracture process zone size. These characteristic behaviours are captured by the numerical simulations of the experiments represented in Fig. 10b and in Fig. 11. Fig. 11 reproduces the complete numerical simulations corresponding to each input pressure value, obtaining the characteristic saw-tooth curves, and which upper points must compare with the experimental graphs. 4. Concluding remarks A numerical coupled model for the HM behaviour of discontinuities recently proposed by the authors, and successfully applied to pre-existing discontinuities in rock mechanics problems, is applied in this paper for the simulation of HM processes in developing cracks in concrete. Validation of the model is accomplished by numerically reproducing a series of wedge splitting tests under the influence of fluid pressure at the notch performed by the group of Prof. Saouma at the University of Colorado at Boulder, and which are available in literature. An elasto-plastic constitutive model, with softening laws that evolve with the work spent on fracture, is assumed for the developing crack. The numerical simulation considers also the stretching of the membrane when the specimen is deformed, by using interface elements with an elastic constitutive model at the notch. These have found to be necessary to accurately reproduce the experimental results. The experimental results, which include splitting force vs. CMOD, crack and fluid fronts vs. CMOD, and fluid pressure distribution along the crack at different times, are efficiently represented by the numerical model. Acknowledgements The first author wishes to acknowledge the FI doctoral Fellowship from AGAUR-DURSI (Generalitat de Catalunya, Barcelona), and support from ETSECCPB and UPC for completing his doctoral thesis. Partial support for this research was also obtained from research projects BIA2006-12717 and BIA2009-10491 funded by MEC (Madrid). The valuable comments and information from Prof. Saouma and Prof. Slowik are greatly appreciated. References [1] Carol I, Prat P, López CM. A normal/shear cracking model. Application to discrete crack analysis. ASCE J Engng Mech 1997;123(8):765–73. [2] Oliver J, Huespe AE, Pulido MDG, Chaves EWV. From continuum mechanics to fracture mechanics: the strong discontinuity approach. Engng Fract Mech 2001;69:113–36. [3] Jirásek M, Zimmerman T. Embedded crack model: I. Basic formulation. Int J Numer Methods Engng 2001;50(6):1269–90. [4] Oñate E, Idelsohn S, Zienkiewicz OC, Taylor RL, Sacco S. Stabilized finite point method for analysis of fluid mechanics problems. Comput Methods Appl Mech Engng 1996;139:315–46. [5] Rashid Y. Analysis of prestressed concrete pressure vessels. Nucl Engng Des 1968;7:334–44. ˇ ervenka V. Inelastic finite element analysis of reinforced concrete panels under in-plane loads. PhD thesis, University of Colorado, Boulder (CO, USA); [6] C 1970. [7] Ingraffea A. Discrete fracture propagation in rock: laboratory tests and finite element formulation. PhD thesis, University of Colorado at Boulder, Boulder (CO, USA); 1977. [8] Saouma VE. Automated non-linear finite element analysis of reinforced concrete: a fracture mechanics approach. PhD thesis, University of Colorado at Boulder, Boulder (CO, USA); 1980. [9] Bazˇant Z. Mechanics of distributed cracking. Appl Mech Rev 1986;39(5):675–705. [10] Hillerborg A, Modéer M, Petersson PE. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cem Concr Res 1976;6:773–82. [11] Elices M, Guinea GV, Gómez J, Planas J. The cohesive zone model: advantages, limitations and challenges. Engng Fract Mech 2002;69:137–63. [12] Carol I, López CM, Roa O. Micromechanical analysis of quasi-brittle materials using fracture-based interface elements. Int J Numer Methods Engng 2001;52:193–215. ˇ ervenka J, Chandra JM, Saouma V. Mixed mode fracture of cementitious bimaterial interfaces. Part II: numerical simulation. Engng Fract Mech [13] C 1998;60:95–107. [14] Goodman RE, Taylor RL, Brekke TL. A model for the mechanics of jointed rock. ASCE J Soil Mech Found Div 1968;94(SM3):637–59. [15] Stankowski T. Numerical simulation of progressive failure in particle composites. PhD thesis, University of Colorado at Boulder, Boulder, CO, USA; 1990. [16] Segura JM, Carol I. Coupled HM analysis using zero-thickness interface elements with double nodes. Part I: theoretical model. Int J Numer Anal Methods Geomech 2008;32(18):2083–101. [17] Boone TJ, Ingraffea AR. A numerical procedure for simulation of hydraulic-driven fracture propagation in poroelastic media. Int J Numer Anal Methods Geomech 1990;14:27–47. [18] Papanastasiou P. The influence of plasticity in hydraulic fracturing. Int J Fract 1997;84:61–79. [19] Murdoch LC, Germanovich LN. Analysis of a deformable fracture in permeable material. Int J Numer Anal Methods Geomech 2006;30:529–61. [20] Simoni L, Secchi S. Cohesive fracture mechanics for a multi-phase porous medium. Engng Comput: Int J Comp-Aid Engng 2003;20(5):675–98. [21] Secchi S, Simoni L, Schrefler BA. Mesh adaptation and transfer schemes for discrete fracture propagation in porous materials. Int J Numer Anal Methods Geomech 2007;31:331–45. [22] Slowik V, Saouma VE. Water pressure in propagating concrete cracks. ASCE J Struct Engng 2000;126(2):235–42. [23] Segura JM, Carol I. On zero-thickness interface elements for diffusion problems. Int J Numer Anal Methods Geomech 2004;28(9):947–62. [24] Carol I, Gens A, Alonso E. A three dimensional elastoplastic joint element. In: Stephansson, editor. Fundamentals of rock joints. Lulea: Centek Publishers; 1985. p. 441–51. [25] Gens A, Carol I, Alonso EE. A constitutive model for rock joints: formulation and numerical implementation. Comput Geotech 1990;9:3–20. [26] Carol I, Idiart A, López CM, Caballero A. Multiaxial behavior of concrete, a mesomechanical approach. Rev Eur Génie Civil 2007;11:907–26. [27] López CM, Carol I, Aguado A. Mesostructural study of concrete fracture using interface elements. I: numerical model and tensile behavior, and II: Brazilian test, compression and HS concrete. Mater Struct RILEM 2008;41:583–620. [28] Caballero A, Willam K, Carol I. Consistent tangent formulation for 3D interface modeling of cracking/fracture in quasi-brittle materials. Comput Methods Appl Mech Engng 2008;197:2804–22.

J.M. Segura, I. Carol / Engineering Fracture Mechanics 77 (2010) 1386–1399

1399

[29] Caballero A, López CM, Carol I. 3D meso-structural analysis of concrete specimens under uniaxial tension. Comput Methods Appl Mech Engng 2006;195(52):7182–95. [30] Caballero A, Carol I, López CM. A meso-level approach for the 3D numerical analysis of cracking and fracture of concrete materials. Fract Fatigue Engng Mater Struct 2006;29:979–91. [31] Caballero A, Carol I, López CM. 3D mesomechanical analysis of concrete specimens under biaxial loading. Fract Fatigue Engng Mater Struct 2007;30:877–86. [32] Segura JM, Carol I. Coupled HM analysis using zero-thickness interface elements with double nodes. Part II: verification and application. Int J Numer Anal Methods Geomech 2008;32(18):2103–23. [33] Brühwiler E, Saouma V. Water fracture interaction in concrete – part I: fracture properties. ACI Mater J 1995;92(3):296–303. [34] Brühwiler E, Saouma V. Water fracture interaction in concrete – part II: hydrostatic pressure in cracks. ACI Mater J 1995;92(3):383–90. [35] Lewis RW, Schrefler BA. The finite element method in the static and dynamic deformation and consolidation of porous media. Chichester: John Wiley & Sons; 1998. [36] Noorishad J, Ayatollahi MS, Witherspoon PA. A finite-element method for coupled stress and fluid flow analysis in fractured rock masses. Int J Rock Mech Min Sci Geomech Abstr 1982;19:185–93. [37] Zienkiewicz OC. El método de los elementos finitos. Barcelona: Editiorial Reverté; 1980. [38] Snow D. A parallel plate model of fractured permeable media. PhD Dissertation. University of California, Berkeley; 1965. [39] Lancelot L, Shahrour I, Al Mahmoud M. Failure and dilatancy properties of sand at relatively low stresses. ASCE J Engng Mech 2006;132(12):1396–9. [40] Slowik V, Plizzari G, Saouma VE. Fracture of concrete under variable amplitude fatigue loading. ACI Mater J 1996;93(3):272–83.