CEMENT and CONCRETE RESEARCH. Vol. 20, pp. 92-I02, 1990. Printed in the USA. 0008-8846/90. $3.00+00. Copyright (c) 1990 Pergamon Press plc.
A S I M P L I F I E D M O D E L TO STUDY F R A C T U R E B E H A V I O U R IN COHESIVE M A T E R I A L S
J. LLorca and M. Elices Department of Materials Science Polytechnical University of Madrid E.T.S. Ingenieros de Caminos Ciudad Universitaria. 28040 - Madrid
(Communicated by F.H. Wittmann) (Received March 22, 1989)
ABSTRACT A simplified model to study crack growth resistance in cohesive materials has been developed. The model relies on the assumption that the crack displacements are known a priori for different crack lengths. The load-displacement curves are obtained for three point bend specimens and compared with experimental results and numerical calculations performed using the finite element method. T h e range of validity of the model is discussed. Introduction Fracture of cementitious materials has been modelled successfully by using the cohesive crack concept, developed by Hillerborg and coworkers (1, 2), and calculations with this model are intimately related with numerical methods; mainly because of the lack of analytical solutions for strain-softening relations which make numerical methods of fitting strain-softening relations very appealing. A typical strain-softening curve is sketched in figure 1, where stress a is a function of crack displacement w. In this model the area under the curve is considered a material property, named the fracture energy GF, and several methods have been devised for its measurement (see, for example, 3, 4). Unfortunately, measurements using direct tensile tests of the shape of the strain-softening curve have not been so successful, (5), and indirect methods to evaluate ~ - w relations are nowadays in use. A simple indirect m e t h o d consists in making an educated guess of the strain-softening curve and check numerical results with experimental ones. As it is well known, experimental results for concrete show some scatter and in order to have confidence with the choosen strain-softening curve it is recommended to check against different sizes. Moreover, if one also has different geometries reliability will be improved. This trial and error method to estimate the ~ - w curve involves a huge effort of numerical calculations due to the many sizes and geometries to solve and because every solution needs iterative procedures. This is why a simplified indirect method, if accurate enough and fast, will be welcome. 92
Vol.
20,
No.
!
93
FRACTURE, MODEL, COHESIVE MATERIALS
The aim of this paper is to develop one of such simplified models, based on findings of Mai and coworkers (6), although the equations of the model are derived following a different procedure. The model is applied to a notched beam loaded at three points and numerical results for different sizes are compared with experimental ones and with numerical ones obtained by a more accurate procedure based on a finite element calculation. The Model This simplified model belongs to the family of cohesive crack models (4), where material properties outside of the crack are elastic and the crack behaves according a strain-softening relation ~ - w , with the restrictions shown in figure 1.
O"
a
il!!!i:::::,
o
[
/
~.
w _-',v (c~., a o , ~
Oo
t
Wc
W
",
FIG. i. Strain-softening curve.
Q - -
-r'
FIG. 2. Crack and displacement notation.
The central hypothesis of the present model is to assume that during crack propagation, displacements of crack faces are a k n o w n function of a, a0 and x as sketched in figure 2.
w=~(.,.o,~)
O)
This assumption allows a straightfoward procedure to compute the external load P in equilibrium with the cohesive crack. This is done by assuming that the stress intensity factor, due to external load P and to cohesive forces, vanishes at the tip of the crack, i.e.,
Ke = Kao
(2)
where Kp and K,o are respectively the stress intensity factors due to P and to cohesive forces. Computation of P is done according the following steps: It begins with the knowledge of w(z), according equation (1). From this value, ~c(x) is known, from the assumed strainsoftening relation. Once the cohesive forces are known, the corresponding stress intensity factor Koo may be computed. Finally, from equation (2) P can be deduced. Obviously, the accuracy of this model depends on the choice of the crack displacement function w = w(a, ao,x). If this function is close to actual displacement this method will provide a fast and accurate solution of the problem. Mai and coworkers (6, 7) make the simplest assumption, a lineal relation between w and the distance z, i.e.: = 2~(. - ~)
(3)
94
Vol. 20, No. l J. Llorca and M. Elices
The slope c~ is known once we (the critical opening displacement, and ape (critical size of the process zone) are known (see figure 3). Computation of %c is done following an iterative procedure and wc is associated with the chosen strain-softening behaviour. From knowledge ofw(x,a) one can compute, for every crack length a, values of the load P as well as corresponding CTOD and displacement.
W¢/2
-
-
X
-
-
a"
-
-
(3. O
--
-
-
+6
QO
FIG. 3. Crack and displacement notation.
Using this linear relation, Foote et al., (6) obtained very accurate results for double cantilever beam (DCB) geometries with a wide range of strain-softening curves (from elasto-plastic to elastic-brittle materials). Three point notched bend beams (NBB) were not explored in that paper and it is unknown if equation (3) is suitable for this geometry, very common in testing cementitious materials. Numerical Procedure This section is devoted to the numerical procedure needed for implementing the model. Detailed calculations are included in the Appendix, and for simplicity a notched beam loaded in three points, NBB, (Figure 4) was selected. These results can be easily extended to other geometries by using the appropriate stress intensity factors.
u! P 1
d
oll
tL
i i cMoo
b-)FIG. 4. Three point bend specimen dimensions.
For a NBB specimen loaded with P = 1 (per unit thickness), the stress intensity factor is (8), (4) I(_p(a) = " ~ vl'~F(a/d)
Vol. 20, No. 1
95 FRACTURE, MODEL, COHESIVE MATERIALS
where S, d and a are shown in figure 4 and F(a/d) appears in the quoted reference as a function of d / S . The stress intensity factor for c~,(z) = 1, as depicted in figure 4, can be also found in the same reference for a beam of infinite length, as
(8)
Ko(a,x) = - - ~ H ( a / d , x / a )
The first step is to assume a reliable displacement function w(x). To start with, a linear relation, like equation (3), will be used.
~(,)= yT...~&~_~) i
\
(6)
e
ape
The critical displacement w, is known and the critical size of the process zone ape is guessed. The actual value will be found by an iterative process, as shown. (Let us suppose we are in the i-th iteration and the critical process zone is ape). If w(z) is known, the cohesive forces are also known via the strain-softening relation and the stress intensity factor due to cohesive forces (according to equation (5)) will be given by aO+at'e
'
.,no
.
aic
~
ae]
(7)
where a,, i a0 and ape i are shown in figure 3. From this result, the external load P~ that balances the stress intensity factor, equation (7), may be computed. Once known, the value of the w(ao) can be evaluated from w(ao) = C T O D i = CPao(aie,ao)P~b -
t
{l
C o o , ( ai
(8~)
o
where
Cpao(ae,ao) = "~
Kp(a). Ka(a,ao)da
(sb)
and 2 [°:
C~o~(aie,ao,x) = ~ y~ K~,(a,x)Ka(a,ao)da These results are derived in the Appendix, jointly with Cp~ o and C~o,. (Expressions A.8, A.9a and A.9b, respectively). The computed C T O D ~ has to be compared with the critical crack tip opening displacement we. If C T O D i is higher than we, one needs to increase a~c (in order to close the cohesive crack) and repeat the whole procedure. If C T O D i is lower than we, the reverse should be done. The actual value of ape is found by iteration following this procedure. Once %~ is known, it is possible to compute the load-displacement curves (P - u) and load-crack mouth opening displacement ( P - C M O D ) . Such results are also measured during experimentation and are needed to check this model;
96
Vol. 20, No. I J. Llorca and M. Elices
Displacement u, is given in the Appendix by: u
Cpp(a)Pb -
=
ja a Cp~(a,x)~r(x)dx
(9)
o
where
2 /o~
Cpp(a) = --~
S3
(10a)
K~(a)da + .4d3E
and Cp~(a, x) = "~
(10b)
Kp(a)Kq(a, x)da
and C M O D is given by the same expressions as for CTOD --equations (8a), (8b) and (8c)-where a0 has to be changed to 0. All these integrals have to be computed numerically and care should be taken with some singularities from the stress intensity factors. All this procedure is independent from the choosen displacement function w(x). The linear relation, shown in equation (3), used by Mai and coworkers is the simplest one, although other relations like: w__.sc
(11)
w -~ a2pe(a -- x)(a-- ao)
were also proven useful, as shown in the next section. Numerical and Exoerimental Results In order to check the accuracy of this simplified model, numerical results were compared with experimental ones, already published by one of the authors (9). Three sets of notched beams were tested, with the same geometry sketched in figure 4. Dimensions, as well as age, and number of samples for each set are summarized in table 1. TABLE 1. Dimensions, age and number of samples for each set tested.
Set
b
d
ao
S
Age(days)
No. o f Tests
FT1 FT2 FT3
100 i00 100
100 200 300
50 I00 150
800 1131 1368
79 45 43
12 6 7
Dimmensions ram.
The average value of direct tensile tests performed with cylindrical specimens, casted at the same time when moulding beams, was 3.04 MPa and the Young's modulus, infered from beams compliance, was 20 GPa. Fracture energy GF, measured according the RILEM proposal, was fairly constant, for all sizes, with an average value of 125 J/rn 2. More details are given in the quoted reference. Records of load P versus displacement u are summarized in figures 5, 6 and 7, for sets
Vol. 20, No. I
97 FRACTURE, MODEL, COHESIVEMATERIALS
FT1, FT2 and FT3, respectively. In fact, figures show a scatterband, and all experimental records lie inside. t25
=
o
0.75
O.SO
SliT
FT-I
' 61UNEAR SOFT;'NING ~ ~ -- ~(I~0NENTIAL SOFTENING
(10
OI
02
03
04
0S
00
r
I
I
I
~1
Q2
03
~,
OI~FLECT~N .{n~n)
~.
z w
:
0E~ECTIGN.Imm)
~.~
2 Experimental and calculated (FEM) results. FIG. 5. Set FTI. FIG. 6. Set FT2. FIG. 7. Set FT3. O0
~1
0.2
03
04
05
OEFLECTION , (ram)
Before performing numerical calculations with the simplified method a strain-softening relation was needed. Because it was not measured during the experimental research, the chosen relation was one that gave the best fit with the experimental results of P versus u. To this end we used a finite element procedure, known as the influence matrix (2), with 40 rectangular elements of constant strain along the beam depth. Two strain-softening relations were found to fit well the experimental results. The first one was a bilinear curve, following Petersson's suggestions, and the second one a exponential curve, according the experimental work carried out by Reinhardt and coworkers (5). Both relationships have been drawn in figure 8. The values of ~, and GF were the experimental measured ones. The correspondent numerical values of P versus u, are also drawn for each size in figures 5, 6 and 7. J Computations done with the simplified method, using the strain-softening relations depicted in figure 8 and expression ~ for the P versus u curve are shown in figures 9, 10 and 11, for the three experimental and for two displacement functions w(x), given by equations (3) and (11). In these figures the experimental scatterband is also shown. Also, the P versus CMOD curves were calculated for the three sizes, according to equation (8) with a0 = 0. Such results are plotted in figures 12, 13 and 14, for the two displacement functions w(x) and the exponential softening curve, already mentioned. For comparison purposes the P - C M O D curves computed using the influence matrix finite element method, are also included. Similar results are obtained with the bilinear softening. Results using eq. (3) are always below those computed using finite element analysis or eq. (11). As
98
Vol. 20, No. l J. Llorca and M. Elices
0"
O"I .
~ -- - - ~
~l
BILINEAR SOFTENING (FROM [2] ) EXPONENTIAL SOFTENING (FROM [B) ) Gt = 3.04 MPO G F = 125
~
Jim 2
FIG. 8. Bilinear and exponential softening curves.
t
0"8~,
3.6N
125 I
o
/
t/ 00
Q1
02
f
no
03 04 DEFLECTION ,(ram)
0.5
0.2
0.3 {]& DEFLECTION .(mm)
Q1
0.2
EO(3)
(11)
0.3 04 OEFLECTION .(ram)
05
Load-Displacement curves. Experimental and calculated results. FIG. 9. Set FTI. FIG. I0. Set FT2. FIG. Ii. Set FT3.
I - - - - - E O (~1
01
0,0
I ....
EXPONENTIALSOFTENING I . ~ . E O
05
already mentioned, accuracy of this m e t h o d depends on the chosen displacement function w; when eq. (3) is used to analyze a DCB good results are obtained, but for NBB eq. (11) gives a better description of the crack profile in the cohesive zone and consequently better values of load-displacement curves.
Vol.
20,
No.
l
99
FRACTURE, MODEL, COHESIVE MATERIALS 125
100
o
075 050I
05 10~1# "
o:sL a
COO0
oo2s
I EO (3) THIS MODEL . . . . EO (11) o~
oor~s
o~o
o12s
o'~o
r
0O0
o~s
CM00 ,(ram)
SET FT-2 I .... FINITEELEMENT51 ~ 005
i 010
i OlS
--~(11)4, EQ 020
..J
0,25
CM.OD .[mm )
z
Load-CMOD curves. Exponential softening. FIG. FIG. FIG.
o
12. Set FTI. 13. Set FT2. 14. Set FT3.
000
o~s
o,~
.... THJ$ MODEL .... o;s o~o o~s
EO (3} EO (?1)
CMO.O
o~o
o.3s
.(ram}
Final Comments A simple model for concrete cracking, based on work of Mai and coworkers, has been developed. It allows a relative simple way to compute the values of P versus u, or CMOD, during crack propagation. 1.- This model relies on the assumption of a simplified crack displacement function: w = w(z, a, ao) and is very sensitive to this choice. Values of P versus u are shown in figures 9, 10 and 11 for the two choices already discussed in this paper, i.e. equations (3) and (11): Values of P versus CMOD are also shown in figures 12, 13 and 14. In both cases (P versus u, and P versus MOD) the curves computed using equation (3) can be shifted inside the experimental band at the cost of using an artificially high tensile strength for the concrete. 2.- With the aim of checking the simplified model in different conditions, various sets of computations were carried out. In order to have a reference value, the same computations were also done using more accurate and lengthy, finite element method already described. The depth of the beams d was changed between 100 and 1000 mm, usual values for engineering purposes. The span to depth ratio was equal to 4. The initial notch to depth ratio ao/d was 0.3, 0.5 and 0.7. The displacement function choosen for the crack was w from equation (11). The maximum load values with the simplified model and the finite element calculations are shown in figures 15 and 16 for the hilinear and the exponential strain-softening curves, respectively. The depth of the beams is presented dimensionless by using the material characteristic length l~h EGF (12)
100
Vol.
20, No. 1
J. L1orca and H. E l i c e s
--~
FINITE ELEMENTS ~ THiS MODEL
"ss
BtLLNEAR SOFTENING
p~..~,
J
s ~ "" "~ ° ° :0.3 d
/ ~
2!f
EXPONENTLAL
1.~ w =.
FINITE ELEMENTS m THIS MOOEL
~m
s"
SOFTENING
,.I jt
**s///p
s
~'
O0
Q5
tO
1.5
20
25
30
35
L.O
C0
0.5
"~:O3
/ / / /
~a 07
tO
1-~
2.0
2.S
30
lS
/- 0
ditch
dilci~
FIG. 15. Maximum load value vs. beam depth. Bilinear softening.
FIG. 16. Maximum load values yes. beam depth. Exponential softening.
As is shown in these figures, the simplified model using equation (11) gives accurate results especially for beam depths between 200 and 500 mm. With smaller sizes, the simplified model underestimates the maximum load and with bigger sizes, overestimates. Accuracy of the results is very similar with both softening curves. Taking into account that, for an elastic-perfectly plastic material, the results of the simplified model are equal to the exact solution, it is reasonable to suppose that the simplified model gives the same accurate results for different strain-softening relationships. Acknowledgements The authors wish to thank Prof. Y.W. Mai for his useful comments and suggestions on this paper. Financial support for this research provided by the CICYT, Spain, under grant PB86-0494 is also acknowledged. References .- A. Hillerborg, M. Modeer and P.E. Peterson (1976). Analysis of Crack Formation and Crack Growth in Concrete by Means of Fracture Mechanics and Finite Elements. Cement and Concrete Research, 6,773-782. .- P.E. Petersson (1981). Crack Growth and Development of Fracture Zones in Plain Concrete and Similar Materials. Report TVBM-1006. Division of Building Materials. Lund Institute of Technology. Sweden. .- F.H. Wittmann, Ed. (1986). Fracture Toughness and Fracture Energy of Concrete. Chapter 7. Elsevier Sc. Publ. The Netherlands. .- J. Planas and M. Elices (1988). Conceptual and Experimental Problems in the Determination of the Fracture Energy of Concrete.Fracture Toughness and Fracture Energy. Test Methods for Concrete and Rock. Tohoku University. Sendai. Japan. .- H.A.W. Cornelissen, D.A. Hordijk and H.W. Reinhardt (1986). Experimental Determination of Crack Softening Characteristics of Normalweight and Lightweight Concrete. Heron, 31, 45-56. .- R.M.L. Foote, Y.W. Mai and B. Cotterell (1986). Crack Growth Resistance Curves in Strain-Softening Materials. J. Mech. Phys. Sohds, 34, 593-607. 7.- B. Cotterell and Y.W. Mai (1987). Crack Growth Resistance Curve and Size Effect in
Vol. 20, No. 1
I01
:RACTURE, MODEL, COHESIVE MATERIALS
the Fracture of Cement Paste. J. Materials Science, 22, 2734-2738. 8.- H. Tada, P. Paris and G. Irwin (1985). The Stress Analysis of Cracks Handbook. Del Research Corporation. 9.- J. Planas and M. Elices (1986). Towards a Measure of GF: An Analysis of Experimental Results. Fracture Toughness and Fracture Energy of Concrete. (F.H. Wittmann, Ed.), Elsevier Sc. Publ., 381-390. ADDendix
Figure A.1 represents a symmetrical solid of a linear elastic material, loaded symmetrically by two forces/'1 and P~. The stored elastic energy n is n
1
= ~[Plul + P2u2]
(A.1)
where Ul and u2 are, respectively, the displacements of P1 and P2. The energy release rate G, will be given by:
(OA~
1 (8u1~ 1p [bu2~ P = ~P1 +
(A.2)
or
a = 1p21C~l(a ) + [P.CL(a) 1 ,
+ P~P.Ch(a)
(A.3)
where ill : Cll(a)P1 -t- C12(a)P2
(A.4a)
II 2 --- C21(a)P1 -~- 622(a)P2
(A.4b)
C12(a) = Cm(a) because the material is linear. Compliances are depending on crack length a as is well known, and C' stands for dC/da. When the energy release rate is written using the stress intensity factors, we have:
P Up O"
(Y
,\
u(x)
u2 FIG.
AI.
Elastic in m o d e
NBB I.
specimen
loaded
FIG.
A2.
S t r e s s and d i s p I a c e m e n t s . NomencIature.
I02
Vol. 20, No. I J. Llorca and M. Elices
G -- (K1 + K2) 2
E"
K~
K~
-E*+~
2K1K2
-+
(A.5)
E-----:--
where E* is the generalized elasticity modulus (depending if we consider plane strain or generalized plane stress). From equations (A.3) and (A.5) one obtains, after integration,
2//
Cij(a) = ~
KiKjda
(A.6)
where Ki are the stress intensity factors corresponding to P~ = 1, per unit thickness. Once known such stress intensity factors, displacements u~ can be computed using equations (A.4a) and (A.4b). When there is a distributed load along the crack faces, as in cohesive cracks (Figure A.2), displacements are given by
up = Cpp(a)Pb + ~a a Cp~(a, x)a(x)dx o
u(x') = C,,p(a, x')P +
~aa C,,,(a, x, x')o'(x)dz
(A.7)
(A.8)
0
where Cpx(a, x) means the displacement of point P when an unitary a(x) is applied at point x. These compliances are given by expressions similar to equation (A.6), i.e.:
2//
c>~(a,x)= ~
c.,.(~,~,~')= ~
KpK.da
o~(.,*')K~, K.da
(A9a)
(A.9b)
The compliance due to the applied load Cpp(a) has to be evaluated by adding to the compliance of the uncracked body the contribution --like equation (A.6) - - of the crack of length a.