- Email: [email protected]
Computation of crack extension energy rate for elastic-plastic hardening materials under variable loading
Pergamon
Engineering Fracture Mechanics Vol. 55, No. 5, pp. 813-830, 1996
PII:
COMPUTATION FOR
S0013-7944(96)00055-0
OF CRACK
ELASTIC-PLASTIC UNDER
EXTENSION HARDENING
VARIABLE
Copyright © 1996 Elsevier Science Ltd Printed in Great isritain. All rights reserved 0013-7944/96 $15.00 + 0.00
ENERGY
RATE
MATERIALS
LOADING
M. CHIARELLI, A. FREDIANI Dipartimento di Ingegneria Aerospaziale, Universita' di Pisa, via Diotisalvi 2, Pisa, Italy and M. LUCCHESI Dipartimento di Costruzioni, Universita' di Firenze, piazza Brunelleschi 6, Firenze, Italy Abstract--The Crack Extension Energy Rate (CEER) for an elastic-plastic hardening material is assessed starting from the flow theory of plasticity. Knowledge of the work done by internal forces as a function of the elastic-plastic state variables, permits numerical evaluation of the CEER under general loading conditions. First results show that if monotonic loading conditions are considered, the Rice Jintegral completely describes the fracture behaviour of a ductile material; on the other hand, when dealing with variable loading conditions, the energy available to extend the crack strongly depends on loading history and hardening rules. Copyright © 1996 Elsevier Science Ltd
INTRODUCTION CRACK PROPAGATIONand residual static strength in a cracked body are assumed to depend on the amount of energy released during a unit crack extension; the assessment of this energy is a crucial problem in fracture mechanics. Different expressions of this quantity as a function of load applied (or loading history) and geometry of the cracked body are available for different constitutive models of materials. For a homogeneous, hyperelastic body, subjected to quasi-static deformations in the absence of body forces, the energy required to produce a unit crack extension, named "Energy Release Rate" (ERR), is defined as G(£) = - ~-~ c r d V +
IaBn, • Tn dA,
(1)
where £, a, T, n, u are the crack length, the strain energy density, the Cauchy stress, the outward unit normal to 0B and the displacement field, respectively; lastly u' is the derivative of u with respect to ~. In a bidimensional problem, if ~ is a path around the crack tip, i.e. a smooth non-intersecting path that starts and ends on the crack faces and surrounds the tip (Fig. 1), the quantity
J = e. Iy(aI - VuTT)m ds,
(2)
where e, I and m are the direction of propagation of the crack, the identity tensor and the outward unit normal to y, respectively, is the J-integral corresponding to path y [1]. The energy momentum tensor (aI - VuTT) is solenoidal, i.e. div(crI - VuTT) = 0
(3)
and thus, if the crack is straight, we have G(£) = J
and, furthermore, J is path-independent. 813
(4)
814
M. CHIARELLIet
al.
Fig. 1. Integration path and domains. According to eq. (4), ERR can be simply obtained by calculating the J-integral along any integration path ),; this result is useful for many engineering applications [2]. The definition of J in eq. (2) is meaningful if an elastic potential a exists or, in other words, if the work done on the body during a deformation process depends solely on the initial and final strains. In any cracked body made of a ductile material, like many metals, a region exists at the crack tip where yielding always occurs; the dimension of this region depends on the geometry, the level of loads and the constitutive law of the material. During a loading process in which some unloading occurs in this region, the work done on the body depends on the strain history, and the elastic fracture model is no longer appropriate; this situation occurs during a static monotonic loading where crack growth is stable, as well as under fatigue loads. For elastic-plastic materials with hardening, a definition of the energy available for a unit crack extension can be obtained by substituting the strain energy density a in eq. (1) with the work w done by the internal forces [3],
d lBwdV+ JoBu'.Tnda. G(e)=-~
(5)
As proved in ref. [3], for a straight crack in a bidimensional body we can write G(~) = J(y) - e. f div(wI - Vu]'T) dA, JL
(6)
for every path around the crack tip, with L the region enclosed by y where yielding occurred and .l(y) = e. lv(wI - VuTT)m ds.
(7)
For elastic-plastic materials the symbols J(~,) in eq. (6) and J in eq. (2) denote two completely different quantities; in fact, J is a function of the actual elastic state alone, while J(~,) depends on the whole deformation process. Moreover, Rice's integral J is path-independent while J(y) is not. On the other hand, for elastic materials, relations (6) and (4) coincide and so, the definition of G is an appropriate generalization of Rice's J-integral for elastic-plastic materials. In this paper a method for estimating G is presented that uses the flow theory of plasticity with general hardening laws. In this way, for all loading histories, both local and global unloading and reloading phenomena can be taken into account. The numerical procedure, described in Section 4, uses the results of an elastic-plastic finite element analysis of the cracked body. After each load increment, the stress T, the plastic deformation E p and the Odqvist parameter ~ (i.e. a parameter defined in the next section of the paper which indicates the amount of plastic deformation) are obtained and, by means of eq. (6), the value of G corresponding to each load increment is assessed. A closed-form expression of the work of the internal forces w as a function of T, E p, e, (, allows accurate computation of the surface integral in eq. (6). Examples are considered in which bidimensional cracked bodies are subjected to loads varying over time and the crack does not propagate during the application of the external load. The results obtained show that the difference between the elastic-plastic and elastic G values is very
Computation of crack extension energy rate
815
small when the load applied increases monotonically. On the contrary in the presence of cyclic loads, the difference is quite large. When cyclic loads with a constant amplitude are applied, different situations may occur. If the values of G become negative during unloading, which can be interpreted as a "closure effect" of the crack, the maximum value of G attained in the following load cycle diminishes. Other situations are possible for different material hardening laws and different values of the ratio between minimum and maximum loads applied. In the case of variable amplitude loading, the behaviour of G depends on the loading history and, in particular, on the presence of overloads that produce a decrease of G owing to the residual stresses during compression. The numerical results obtained are in agreement with experimental data when a realistic constitutive behaviour is assumed and this underscores the importance of adopting a suitable model for the hardening laws of the material. THE CONSTITUTIVE M O D E L For the reader's convenience, we shall briefly present certain elements of the flow theory of infinitesimal plasticity; we refer the reader to ref. [4] for a detailed account of the theory dealt with in the following. A deformation process or, more briefly, a history of duration ~ is a mapping E that assigns to each instant r < ~ the value of the infinitesimal strain E(~) reached at a fixed material point; we suppose that E(0) = O, so that all deformation processes are assumed to begin at some fixed initial state. The materials being considered here are elastic-plastic isotropic solids whose mechanical response to deformation processes is described by a rate-independent constitutive functional. We use T(T) to denote the stress at time z associated with history E by the constitutive functional. The kind of constitutive response is further specified by the notion of elastic range and plastic history. Elastic range e(z) at time z corresponding to history E is a set that contains the current value of strain E(z); its points are interpreted as the infinitesimal deformations from the reference configuration to configurations which are elastically accessible from the current one. We suppose that at each instant z the current deformation E(T) can be decomposed into the sum of an elastic part Ee(z) and a plastic part EP(z), E(r) = Ee(r) + EP(r),
(8)
having the following properties. EP('r) is a deviatoric tensor and belongs to e(T). The current stress depends linearly and isotropically on the elastic part of the deformation, T(r) = 2/z[E(r) - EP(r)] + ~.trE(r)I,
(9)
where trE(T) is the trace of E(T), I is the identity tensor and 2, # are the Lame' moduli of the material. Relation (9) reflects the classical hypothesis that the stress response to a purely elastic strain from the unstressed configuration reached after unloading at the current instant ~ is both unaffected by the past deformation process and completely determined by E(~) and EP('c). For each history E, ~'(r) =
IIEP(v')[[ dr'
(10)
is the length of the path described up to instant • by the plastic deformation tensor. ~ is called the Odqvist parameter. In view of the applications we have in mind, we accept the von Mises criterion. That is to say, we suppose that for each history E and for each T, the corresponding elastic range is the cylinder e(r) = {EIJE0 - C(r)ll _< p(U(r))},
(11)
816
M. CHIARELLIet al.
where Eo is the deviatoric part of E and (1) p is a differentiable, non-decreasing real function that depends on the material, is independent of history and describes the isotropic hardening; (2) for each history E, C is a deviatoric history. Moreover, we accept the classic linear kinematic hardening rule. That is to say, we suppose there exists a non-negative material constant ~b, such that for each history E and for each instant z we have C(r) = (1 + ~)EP(r).
(12)
In particular, a material for which p is a constant function and for which we have ~ = 0 is called ideally plastic. The set of constitutive hypotheses is completed by the flow rule which states that when EP('r) is different from zero, it is parallel to N(z), the outward unit normal to the elastic range at E(z), EP(v) = ((r)N(r), N(r) = [p(((v))l-l(E0(v)
- -
C(v)).
(13)
For each history E and for each z, T(z).l~(z) is the stress power and thus the work done by internal forces in the deformation process E up to time z is w(r) =
T(r'). l~(r') dr'.
(14)
In ref. [5] it has been proved that for each history E and for each time z, w ( z ) can be expressed as a function of E(z), EP(z) and ~(z). In fact, we have w(r) -- #ll(E(r) - EP(r))II 2 + lk(trE(r)) 2 + #lPIIEP(r)I[ 2 + 2/zw(((r)),
(15)
where ~o is the primitive of p such that 09(0) = 0.
CRACK EXTENSION ENERGY RATE FOR ELASTIC-PLASTIC MATERIALS
Let B be a regular homogeneous two-dimensional elastic-plastic body. We suppose that B contains an edge crack, represented at every instant z by the image of a smooth non-intersecting curve ~e, parameterized by arc length ct (Fig. 1). The length e = g(z) of curve ~ is a nondecreasing function of time z during the motion of B; ~e(0) and ~(g) represent the intersection of the crack with the boundary aB of B and the tip of the crack, respectively, while e(t~) = dle/do~
(16)
denotes the unit vector field tangent to the crack. In the interval [g0, gl], where the crack advances, we put dl wdA+ J u'.Tnds. G(g) = - ~-~ B ~B
(17)
J] G(~) d~
(18)
For each ~---~1,
0
can be interpreted as the work per unit thickness needed to increase the length of the crack from g0 to e; G is called the crack extension energy rate. For each g ~ [g0, el], let B~(g) = (xEBl((x, e) # 0},
(19)
with ~(x,g) defined by eq. (10), be the set of all points of B in which plastic deformation have
Computation of crack extensionenergyrate
817
taken place. Let us call B~(£) the plastic region; moreover, it should be noted that B~(£), in general, does not coincide with the subset of B in which we have EP(x,£) ~ 0. A curve ? is called a path around the tip if it is a smooth non-intersecting path that starts and ends on the crack, and includes the tip of the crack; we write m for the outward unit vector field normal to y. For each path ~ around the tip, let F be the subset of B enclosed by 7 and let L = bY~B~ be the intersection of F and B~ (Fig. 2). In ref. [3] it was proved that if the crack is straight, for each path 7 around the tip we have
G(e) =
)(y) - e. [ div(wI - VuTT) dA, JL
(20)
where f
J(7) = e. ] ( w I - VuTT)m ds.
(21)
If, for some £ ~ [£0,£1], the distance between B~(£) and aB is positive, there exist paths around the tip that include B~(Q. If 7 is one of these particular paths around the tip, F contains B~ and it therefore follows from eq. (20) that G(e) = J(y)
-
e. [ div(wI - VuTT) dA. JB~(t)
(22)
Moreover, since 7 does not intersect with B~, in view of eq. (17), J(7) can be calculated, rather than from eq. (21), from eq. (2) as in the case of hyperelastic materials, and thus we have
G(e) = J - e. [ JB
div(wI - VuTT) dA. ~(t)
(23)
No analytical assessment of G is available; it is to this purpose that a numerical procedure, presented in the next section, has been developed for bidimensional elastic-plastic bodies.
NUMERICAL EVALUATION OF G Numerical evaluation of G is carried out in two steps: finite element analysis of the cracked body and computation of the line and surface integrals according to eq. (20). The finite element code used in this research is NO.S.A. [6]. This code is particularly suitable for elastic-plastic analyses for both infinitesimal and finite deformations; in fact, it is very versatile in defining the input of the constitutive laws of elastic-plastic materials and, moreover, the values of ( are available at each Gauss point and for each load step. Computation of G is car-
i
fa(~)
i
7, --
-
s
411--
-
t
i
f2(g ) __
i
4"7 s.. 1 i
s'.
'
i 2
~{ (Plasticregion) Fig. 2. Plastic region and its intersectionwith the region F whose boundary is the integrationpath ~,.
Fig. 3. Sketch of the six parabolic functionswhich interpolate Ep and ~. O - - n o d e , * - - G a u s s point.
818
M. CHIARELLI et al.
ried out by GJINT2D, a code developed in-house, for which the finite element results at any load step are the input and the value of G is the output. Both the line and surface integrals in eq. (20) are computed with the help of a Gaussian technique. This is performed by means of an appropriate finite element mesh around the crack tip, defined in such a way that the integration paths pass through a set of Gauss points of the elements. In any Gauss point, the quantities to be integrated are expressed as functions of u, E, E p, T, ( and their derivatives. Assessment of these derivatives plays a central role in the whole numerical procedure implemented in code GJINT2D. The components of tensor E and their derivatives are computed starting from the nodal components of the displacement field by using the shape functions of the elements in the same way as in ref. [2]. Calculation of the derivatives of ( and of the components of E p by means of the shape functions is not reliable, because these quantities are known at the Gauss points only. So, for this purpose we employ a simple technique involving parabolic interpolation along coordinate directions ~ and r/ of the isoparametric space. As an example of this procedure, let us compute the derivatives of ( at the Gauss point labelled as number 5 in Fig. 3. We set f2(O = A~2 + B~ + C,
(24)
g2(0) = Dr/2 + Erl + F,
(25)
where the six coefficients of the quadratic functions (24) and (25) are evaluated by imposing the crossing conditions at points 4, 5 and 6 for f and at points 2, 5 and 8 for g. In this way, the desidered derivatives can be obtained by the chain rule, Og" df2 O~ dg2 Or/ Ox d$ Ox ~ do Ox' 8( dr2 8~ Oy -- d~ Oy
dg2 Or/ do Oy'
(26) (27)
where ax
O0 Or/
(28)
ax
is the inverse of the Jacobian matrix of the isoparametric transformation. Appendices A and B are devoted to solving the main problems that arise in computation of the surface integral; they are: (a) the Gauss points used to calculate the surface integral can be different from the original Gauss points of the element; (b) interpolation of ( and of the components of E p near the boundary of the yield region requires particular care. The general scheme of the whole procedure is shown in Fig. 4. In this paper, the reliability of the numerical procedure in calculating G can not be deduced from comparisons with results available in the literature. Nonetheless, we can point
Computation of crack extension energy rate
I
NO.S.A.F.E. PROGRAM ~ . ~
"<./'//////////////J
819
Solution of constitutive and equilibrium [ equations under the hypothesis of infinitesimal deformations.
Evaluation of ~, E P,E, T, u ] at each load step
i
~GJINT2D
Reading of the geometrical parameters and total displacements from a post file
MODULvIIIIII:~
T I Vectors and matrices containing 1 4, EP, E and T are initialized Loop on elements and their Gauss points
Evaluation of p(~) and c0(t)
I
Evaluation of Line Integral Functions
Evaluation of first derivatives of W, 4, E, EPand second derivatives of displacement components u i
I Y Evaluation of Line Integral
Evaluation of Surface Integral
L
I Evaluation of the Crack Extension Energy Rate
Fig. 4. Flow diagram with the principal steps for numerical evaluation of the Crack Extension Energy Rate.
out two conditions that are necessary (albeit not sufficient) for reliable computation of G and that are verified in all the examples discussed in the next section. They are: (1) the surface integral must be null in the first loading step when the material is elastic, starting from its annealed condition; (2) G must be path-independent, even though the line and the surface integrals are not.
EXAMPLES The examples shown in this section concern computation of G in the case of bidimensional bodies under plane strain conditions, with a crack that does not propagate during the loading process. Different loading histories are considered.
820
M. CHIARELLIet al. 19.65
I
10.9
---Q--" I I
26
Iiiiiiiiiiiiiiiill
82.5
Fig. 5. CT specimen undergoing a monotonic loading history for three differentvalues of kinematic hardening parameter ¢ [dimensionsin mm]. Compact tension specimen subjected to a monotonically increasing load
The CT specimen in Fig. 5 is made of a ductile material whose constitutive uniaxial behaviour is described by a bilinear curve with the following parameter values: Young's modulus
E = 113.000 MPa,
Poisson's ratio
o = 0.32,
Shear yield stress
• o = 461 MPa.
Here we assume the isotropic hardening function
P(O - (vo/2/z).,~(1 + ,80, with fl = 25. For kinematic hardening, we consider the following three values of the material constant if: ~p = 0.1, ¢ = 0.30 and ~p = 0.45. The finite element mesh and the integration paths considered in this analysis are shown in Fig. 6. The specimen is subjected to a monotonic load, increasing up to a value such that large deformations are avoided in the yield region around the crack tip according to the hypotheses of the infinitesimal t ~--ory of plasticity. The results (Fig. ~ show that for all values of ~kconsidered, the surface integral is small and so, in view of eq. (23), G can be well approximated by Rice's J-integral [eq. (2)]; moreover, as the contribution of the surface integral is always negative, Rice's J-integral is an upper bound of G. Experimental assessment of the fracture toughness of materials is carried out on specimens under a monotonically increasing load. According to the previous conclusion, plasticity does not play an important role and, then, when the yield region is not too large, characterization of material can be made with the methods of linear elastic fracture mechanics. The results in Fig. 7 also show that, when ~k increases, the surface integral tends to zero and then, G approaches Rice's J-integral; this agrees with the fact that the material becomes elastic when ~k increases indefinitely.
Computation of crack extensionenergyrate
821
III11111111IIII11111111I tl IIIIIIIIIIIIIII1111111I Imrlll-I I 1 I
II111111tl IIIIII111I]111 II lilllllll IIIIIIIIIit111
11 I!1111111IIIIIIIIt1~11 I
, 4
L~
-rrr~"nt
Illlttli-I-I'll II I II I:Yd-rtll IIIIIIIII I,'1 I I I I IIdllllll IIIIIIIII I,I I I I
7i
a
I
I lllllllll IIIIIIIII1~1 I I
I I rl I I llllllll i I IIIIIIII
III
IITITI1 rl-l'l i i IIIIII
II I I1
..........I I ' 1"',I IIIIIIII
A
crack Fig. 6. Mesh of the CT specimen shown in Fig. 5 and integration paths used for the evaluation of G.
Although it is not possible to compare the present results with analogous data in the literature, we observe that conditions (i) and (ii) pointed out at the end of the previous section are satisfied. In fact, as shown in Table 1, the surface integral in the first load increment is about 10 - 5 times the corresponding line integral; the same holds true for all the integration paths considered. The previous results are valid in the presence of small deformations only. For a CT specimen, when the external load is sufficiently high, large deformations occur in front of the crack tip and the infinitesimal theory of plasticity is no longer appropriate for computing the surface integral. For these problems, although the expression (22) of G is still valid, we need a theory of finite deformation plasticity to assess the work of the internal forces and the stresses; these aspects will be treated in a future paper. Worthy of note is the fact that in all the computations carried out with different loading conditions and constitutive laws, the values of G computed for different paths are the same within an error of about 0.1%, as Table 2 shows.
Center-cracked panel undergoing time variable loads The finite element mesh of a quarter of a CCP and the integration paths are shown in Fig. 8. Two materials are considered. The constitutive law of the first material is modelled by the bilinear curve shown in Fig. 9, having the characteristics of a steel alloy:
Young's modulus
E = 0.206 x 106 MPa,
Poisson's ratio
v -- 0.3,
Shear yield stress
Zo -- 165 MPa.
We suppose that the material is purely isotropically hardening, and apply the same function p (() used in the previous example. The loading history, shown in Fig. 10, consists of three cycles with R = 0. When the load increases during the first cycle, the surface integral is small and the results are similar to those obtained in the previous example. The effect of the loading history becomes very important during subsequent unloading (Fig. 11), when the absolute value of the surface integral increases and, according to eq. (20), the crack extension energy rate is reduced. Then, for 0_~O'~_~0"1, G
822
M. CHIARELLI et al.
Combined Hardeningwith []=25.0 and ¥=0.I 200
Lm~~ o ]
--o-
i
/
[IOhn^2 ] 100
i
l
l
l
l
l
l
}
l
l
l
l
l
l
,
,
,
i
. . . .
10
,
. . . .
i
20
. . . .
t
. . . .
30
40
Load [ KN 1
Combined I-lard~xingwid~~25.0 and ¥=0.3 200
[ KJ/m^2 ]
'"
'
'
i
. . . .
I
. . . .
I
. . . .
I
. . . .
I
. . . .
I--°--
I
. . . .
/
I
. . . .
J"
100
IO
20
30
40
Load [ K N ]
Combined Hardening with [3=25.0and ¥=0.45 200
' ' r ' r l
. . . .
I
. . . .
I
. . . .
i
. . . .
I
r
I
. . . .
.
.
I
. . . .
l
~
t
. . . .
I ' ' '
[I~J~^2 ] lOO
. . . .
i
. . . .
10
.
.
.
L
l
l
20
l
l
l
l
l
30
l
l
h
i
,
,
,
,
40
Load [ I ~ 1 Fig. 7. Results for the CT specimen undergoing a monotonic loading history with three different values of the kinematic hardening parameter ~ (Example 1).
Computation of crack extension energy rate T a b l e 1. O u t p u t o f the G J I N T 2 D
823
r o u t i n e f o r the C T s p e c i m e n u n d e r m o n o t o n i c l o a d i n g c o n d i t i o n . Units: k J / m 2. F i r s t load increment
Circuit number
Line integral
1 2 3 4 5 6 7 8 9 10 11 12
0.89279e0.89189e0.89207e 0.89334e0.89281e 0.89269e 0.89328e 0.89300e0.89286e0.89299e 0.89297e 0.89283e -
Surface integral
C.E.E.R.
0.43184e- 6 0.43159e - 6 0.39360e- 6 0.41509e-6 0.41506e- 6 0.39407e - 6 0.39396e-6 0.39444e- 6 0.39438e - 6 0.39011e-6 0.39488e-6 0.39386e- 6
0.89279e 0.89190e0.89207e 0.89334e0.89281e 0.89269e 0.89328e-1 0.89301e 0.89287e 0.89300e-1 0.89290e0.89283e -
1 1 1 1 1 1 1 I 1 1
becomes negative and, at the end of the first cycle, we have G(e)l~_-0<0
and
aG(e)/atrl,,=o=O.
The condition aG(e)/aal,~ = 0 = 0 at the end of the cycle indicates that during unloading, for small values of tr, we can write G(e)~-O~o+ ~ l a 2 as in the elastic case; then we deduce that during the first unloading the material in the yield region near the crack tip behaves elastically. At the beginning of the second cycle the material continues behaving elastically and, moreover, G(e) remains negative up to tr2. For O
T a b l e 2. O u t p u t o f the G J I N T 2 D
Circuit number 1 2 3 4 5 6 7 8 9 10 11 12
r o u t i n e f o r the C T s p e c i m e n u n d e r m o n o t o n i c l o a d i n g c o n d i t i o n . Units: k J / m 2. Fourth load increment Line integral
Surface integral
C.E.E.R.
0.55879 0.55822 0.55833 0.55913 0.55880 0.55872 0.55909 0.55892 0.55883 0.55892 0.55890 0.55881
- 0.12992e - 0.12992e - 0.12994e- 0.12993e - 0.12993e -0.12994e-0.12994e-2 -0.12994e- 0.12994e- 0.12994e -0.12994e-2 - 0.12994e-
0.55749 0.55692 0.55704 0.55783 0.55750 0.55743 0.55779 0.55762 0.55753 0.55762 0.55760 0.55751
2 2 2 2 2 2 2 2 2 2
824
M. CHIARELLI et al.
III
IIII1
--tll crack
++++++++++o ~,
800
.~
Fig. 8. CCP model used to study two different materials for the loading histories shown in Figs 10 and 13. The finite element mesh and a magnification o f the crack tip zone together with the integration paths are shown [dimensions in ram].
Center-cracked tension specimen undergoing variable amplitude loading In this example, we consider the CCT specimen in Fig. 15 made of a material with the following characteristics: Young's modulus
E = 73000 MPa,
Poisson ratio's
v = 0.32,
Shear yield stress
Zo = 185 MPa;
different hardening laws and load histories are considered. In the first case, the material is supposed to be purely kinematically hardening with = 0.28. Moreover, the load history considered, subdivided into 480 load increments, consists of five cycles followed by three overstresses and, finally, by five cycles with the same amplitude as the first ones (Fig. 16). The maximum values of G reached during the two series of five cycles are constant, while a decreasing trend is seen over the course of the three overload cycles because the overstresses cause an increase in the plastic region near the crack tip (Fig. 17). The shape of the crack faces is associated with the presence of residual strains close to the crack tip; as an example, Fig. 18 shows that the crack is still open when the load is null after the last overstress cycle. The decrease in the crack extension energy rate after an overstress cycle can be interpreted (even though in this paper we are dealing with a stationary crack) as a delay in crack propagation. In the present case of loading history always under tension, no contact between the crack faces ~o = 165 MPa 21.1,= 158.5 GPa T0
tga = 5654 MPa TO P°='4T 2~
= 0.00148
t g a tga) = 25 13: P o ( ~ ,
,
,
,
,
,
,
1 2 3 4 5 6 7
,
im
e(O/O0)
P u r e Isotropic H a r d e n i n g
Fig. 9. Stress-strain curve and numerical parameters used for the first material in Example 2.
Computation of crack extension energy rate
825 Pulsating Load Pure IsotropicHardening 6=25.0
200 Om~, =
108.86 MPa
6min •O
--
• ~aoadt~1 [ ILI/m^2
R=0
]
- - b - - loading 2 , • unloading 2 ---o-- l o a d ~ 3
100
Jd" G2 O 1 ~ P---
170'w~l( c~ If .~f
Gmin' 31 123 183 275 335 418 I
Load Increment N u m b e r Gmax
Fig. 10. First loading history applied to the CCP model shown in Fig. 8.
-100
~
~ 40
80
o [ lVlPa1 Fig. 11. Results for the CCP model undergoing the pulsating loading history shown in Fig. 10.
takes place in the yield region, so this delay can be attributed to the stress and strain fields in the yield region in front of the crack tip, rather than to the kinematics of the crack faces. In the second load history considered here, the maximum stress is constant and three series of three cycles, with R = 0.0, 0.3 and 0.5, are applied (Fig. 19); the third load history is made up of three series with the same amplitude as before and R variable according to Fig. 20. Owing to the assumption of pure kinematic hardening, the decrease in G at the end of each cycle o f the first series is nearly constant in both loading histories, as shown in Figs 21 and 22. In the second series of cycles, the same occurs, even though the decrease in G per cycle is much smaller. In the third series the load amplitude is 50% of the first, the stress variations in the yield region are always inside the elastic range and, so, no decrease in G occurs during these cycles. However, the effects of the loading history on the variation of G in the last series of cycles is very important; for example, in the case of R = 0, the delay effect is so marked that the last cycles produce nearly no crack extension. The importance of the material constitutive law is clearly shown by considering a combined hardening, with: 13 = 25, ~p = 0.28. The response of the material is completely different, as shown in Fig. 23: in the first three cycles the reduction of G is much smaller than in the previous case and the rate of decrease at the end of each cycle is not constant. In the next six cycles, the material behaves elastically owing to the increase in the radius of the elastic range, and the rate of energy available after the first group of cycles is higher than in the case of pure isotropic or pure kinematic hardening. This behaviour seems more realistic.
Young Modulus
73000 MPa
"to = 144.3 MPa Xmaxffi 283 MPa F.max =
0.018
To each of the different ranges of the Odclvist parameter values corresponds a
different isotzopic hardening parameter ~. 10
20
e (o/oo)
Pure lsotropic Hardening
Fig. 12. Stress-strain curve used for the second material treated in Example 2.
120
826
M. CHIARELLI et al.
(~
/'N
~,
~
Omax = 98.4 M P a Omin = 68.8 M P a
~mig
R = 0.7
21
49 77 1 5 133 1 1 Load Increment Number
Fig. 13. Second loading history applied to the CCP model shown in Fig. 8.
Results for the CCP model undergoing variable amplitude loading history Pure Isotropic Hardening 500 .... . . , . . . . . . , • .
[KJlm^2 ] 300
IOO
(followil~ ~ unloatang h
l
i
i
20
l
i
i
l
i
i
i
40
l
i
i
60
i
loadil~"
pnaa~ i
80
l
i
l
100
a [MPa]
Fig. 14. Results for the CCP model, with the stress-strain curve depicted in Fig. 12, undergoing the variable amplitude loading history shown in Fig. 13.
\/ i i
too
i
~
Sl~men
--~ thickness
Ultimate stress ~/ield stress UlIimate elongation
crack
6.a5 mr. 469 MPa
324 MPa 0.20
Fig. 15. CCT specimen used in Examples 3 and 4. The finite element mesh and the crack tip zone with the integration paths used are shown [dimensions in mini.
Computation of crack extension energy rate
827
160 o
[MPal 120
80
40
A ....
i ....
B
L ....
i ....
100
i ....
C
i ....
200
L ....
i ....
aoo
400
Load increment number Fig. 16. First loading history applied to the CCT specimen shown in Fig. 15 (Example 3).
Results for the CCT specimen undergoing loading history with overloads Pure Kinematic Hardening - ¥ = 0.28 15
9 [KJ/m^21 lO
-5
. . . .
I ' ' J l ~
50
. . . .
100
'
. . . .
150
i
. . . .
200
J
. . . .
260
i
. . . .
300
i
. . . .
350
400
Load Increment Number Fig. 17. G vs load increment for the specimen shown in Fig. 15 undergoing the loading history shown in Fig. 16.
1.2 V / rdCL 1 0.8
0.6
.........................
. ........................................................
÷ ............................
i. . . . . . . . . . . . . . . . . . . . . . . . . . .
0.4 ...................................................
i ..............................................................................
0.2 0 0
0.2
0.4
0.6
0.8
(distance from the crack tip) / (half crack length) Fig. 18. Profile of the crack face for the CCT specimen shown in Fig. 15. This profile is obtained unloading the panel completely, after point B of the loading history shown in Fig. 16 is reached.
828
M. C H I A R E L L I et al. 200
200
a [MPa]
(~ [MPa] 160
120
120
81) 4O
40
0
i
0
100
I
I
200
,
I
,
I
300
,
I
0
,
500
400
0
100
200
Load IncrementNumber
300
500
400
Load Increment N u m b e r
Fig. 19. Second loading history applied to the C C T specimen: variable R (Example 3).
Fig. 20. Third loading history applied to the CCT specimen: R = 0 (Example 3).
CONCLUSIONS
On the basis of the flow theory of plasticity, a path independent parameter, the Crack Extension Energy Rate (CEER or G), taking into account local and global unloading of stress state has been defined. A computer code, GJINT2D, has been developed which evaluates G from the results of elastic-plastic finite element analysis. The CEER has been assessed for generic elastic-plastic materials with isotropic, kinematic or combined hardening laws, in the case of plane problems with a stationary crack. The results show that if a monotonic loading condition is considered, the CEER does not greatly differ from the Rice J-integral; moreover, when unloading is considered, a difference, at times large, does appear even in the presence of small scale yielding. In the case of cyclic loading, the influence of the hardening law is so crucial, that the cyclic stress-strain curves of the material must be known very accurately. The occurrence of overloads gives rise to a decrease in CEER and, then, the delay in crack growth can be linked to the G vs load curves.
Results for the loading history with variable R Pure Kinematic [-Iarde~ting- ¥ = 0.28 2O
•
.
i
.
.
,
i
'
'
i
•
'
'
i
.
Results for the loading history with R=0 Pure Kinematic Hardening - ¥ = 0.28 .
j
[ KJImA2 ] 1o
0
40
80 120 a [MPa]
.
.
.
.
2o
i
,
1o
Fig. 21. G vs load for the C C T specimen shown in Fig. 15 undergoing the loading history shown in Fig. 19: pure kinematic hardening (Example 3).
200
-10
'
i
•
'
'
J
'
'
t
j
[ KJlm^2 ]
160
,
i
b
i
I
40
b
i
J
I
i
80
i
i
I
120
*
i
,
I
i
i
160
cr [MPa] Fig. 22. G vs load for the C C T specimen shown in Fig. 15 undergoing the loading history shown in Fig. 20: pure kinematic hardening (Example 3).
*
200
Computation of crack extension energy rate
829
Results for the loading history with variable R Combined Hardening - [~= 25.0, ¥ = 0.28 20
//
[KJlm^2] 10
~
0
t
i
40
i
80 o
i
t
,
i
i
120
,
20o
160
[MPal
Fig. 23. G vs load for the CCT specimen shown in Fig. 15 undergoing the loading history shown in Fig. 19: combined hardening (Example 3).
REFERENCES 1. Rice, J. R., A path independent integral and the approximate analysis of strain concentrations by notches and cracks. 3". appl. Mech., 1968, 35, 379-386. 2. Chiarelli, M. and Frediani, A., A computation of the three-dimensional J-integral for elastic materials with a view to applications in fracture mechanics. Engng Fracture Mech., 1993, 44, 763-788. 3. Chiarelli, M., Frediani, A. and Lucchesi, M., On the crack extension energy rate in elastic-plastic bodies. Defect assessment in components. Fundamentals and applications, ESIS/EGb9, ed. J. G. Blauel and K. H. Schwlbe. Mechanical Engineering Publications, London, 1991, pp. 87-99. 4. Lucchesi, M. and Podio Guidugli, P., Materials with elastic range. Part III: approximate constitutive relations, Arch. Rat. Mech. AnaL, 1992, 117, 53-96. 5. Lucchesi, M., Free-energy functions for elastic-plastic material elements. Q. appl. Math., 1993, 60, 299-318. 6. Guidotti, P., Lucchesi, M., Pagni, A. and Pasquinelli, G., Elastic-plastic behaviour with work-hardening: an appropriate model for structural software. Meccanica, 1984, 19, 43-51.
(Received 6 October 1994) APPENDIX
A
Correct application of Gaussian integration rules in calculation of the surface integral requires all elements crossed by integration paths to be modified and the new Gauss points to be properly placed [2]. For this reason, the interpolation scheme becomes quite complex. As an example, crossing an element in the t/direction, as in Fig. 3, once we get the isoparametric coordinates of new Gauss points, for a generic component E p of plastic strain E p at the point 1' (Fig. A1), we have
e p =A(~I'), ~l, ¢~1, E~ =g,(ol'), ~,' = 0 t ;
tit
I' ?,
i, , l
il
i! ' Y3-
(AI)
integra~onpath
f'
il
il
',
i"i
5"ii
v,,j
I/
I
'I
Fig. AI. Location of Gauss points in a modified element because of a crossing in the positive h direction. • - - generic Gauss point of unmodified element, * - - generic Gauss point of modified element.
Fig. A2. Location of Gauss points in a modified element because of a comer crossing, s - - generic Gauss point of unmodified element, * - - generic intermediate control point, [] - - generic Gauss point of modified element.
M. CHIARELLI et al.
830
V 1
~=-1
2
~i\
interpolating parabola
~=0
!~
3
1 ~=-1
~=1
Fig. BI. If ( = 0 at Gauss point no. 3, with a standard interpolation procedure we would obtain a wrong derivative.
;' E ~
/
parabolawho6e vertex coincides
[ withintermediate pointofthe tern
2 ~,=0
a ~=1
Fig. B2. Interpolation example with two points o f the term for which ( = 0 or E ~ = 0.
the apex indicates that the Gauss point of a modified element must be considered and gl(r/) is the parabolic function interpolating E ~ on points 1', 4' and 7'. So, its derivative with respect to x at 1' is
OE p dfl O~ . dgl Or/ 0x = ~'~xx "t ~ ~x"
(A2)
It is note-worthy that, when the derivatives of E p and ~ are calculated for a modified element, all terms o f the inverse of the Jacobian matrix must be computed in the isoparametric reference frame of the complete element. p For a corner crossing (Fig. A2), E ij and ( and their derivatives are computed by means o f a double interpolation scheme. If we use two apices to indicate the new Gauss points and one apex for intermediate control points (Fig. A2), at 1" we have E~=gj(r/1,,), r/J" # r / l ' = r / l E p =fl(~{'), ~1" = ~ t ' # ~ l
(A3)
and then
OE~jP _ d f l O~ Ox d~ Ox
dgl0r/ dr~ Ox'
(A4)
where, as in the previous case, gl(r/) interpolates E p in l', 4' and 7', whileJ~(O interpolates E p0 in l", 2" and 3".
APPENDIX B Approaching the boundary of the plastic region, E P, ( and their gradients tend to zero. Consequently, the previous interpolation scheme needs to be modified on the elements lying close to this boundary. In fact, supposing ( = 0 at the Gauss points 3, 6 and 9 (Fig. 3), if we simply interpolated ( values in the ~ direction, we would get a parabolic function, resulting in the curve drawn roughly in Fig. B1, whose derivative at Gauss point 3 is different from zero; the same holds for points 6 and 9. To avoid this kind of error, whose numerical influence is greater for a small plastic region size, we must modify the procedure. Cases of interest are:
(1) = 0 o r e 0P = 0 at two Gauss points o f the current interpolation term;
(2) ( = 0 or E ~ = 0 at a single Gauss point of the current interpolation term. In the first case, we select a parabola whose vertex coincides with intermediate integration point: for instance point 2 of terms 1, 2, 3 (Fig. B2). In the second case, we select a parabola as shown in Fig. B3; the possible error introduced at the intermediate point can be accepted. Therefore, before calculating the derivatives, a check is made to locate the Gauss points where ( or E ~ are null.
~
' E~I parabola withvertex , / ~ a t pointforwhich
2
a
{
Fig. B3. Interpolation example with one point of the term for which [ = 0 or E ~ = 0.
Engineering Fracture Mechanics Vol. 55, No. 5, pp. 813-830, 1996
PII:
COMPUTATION FOR
S0013-7944(96)00055-0
OF CRACK
ELASTIC-PLASTIC UNDER
EXTENSION HARDENING
VARIABLE
Copyright © 1996 Elsevier Science Ltd Printed in Great isritain. All rights reserved 0013-7944/96 $15.00 + 0.00
ENERGY
RATE
MATERIALS
LOADING
M. CHIARELLI, A. FREDIANI Dipartimento di Ingegneria Aerospaziale, Universita' di Pisa, via Diotisalvi 2, Pisa, Italy and M. LUCCHESI Dipartimento di Costruzioni, Universita' di Firenze, piazza Brunelleschi 6, Firenze, Italy Abstract--The Crack Extension Energy Rate (CEER) for an elastic-plastic hardening material is assessed starting from the flow theory of plasticity. Knowledge of the work done by internal forces as a function of the elastic-plastic state variables, permits numerical evaluation of the CEER under general loading conditions. First results show that if monotonic loading conditions are considered, the Rice Jintegral completely describes the fracture behaviour of a ductile material; on the other hand, when dealing with variable loading conditions, the energy available to extend the crack strongly depends on loading history and hardening rules. Copyright © 1996 Elsevier Science Ltd
INTRODUCTION CRACK PROPAGATIONand residual static strength in a cracked body are assumed to depend on the amount of energy released during a unit crack extension; the assessment of this energy is a crucial problem in fracture mechanics. Different expressions of this quantity as a function of load applied (or loading history) and geometry of the cracked body are available for different constitutive models of materials. For a homogeneous, hyperelastic body, subjected to quasi-static deformations in the absence of body forces, the energy required to produce a unit crack extension, named "Energy Release Rate" (ERR), is defined as G(£) = - ~-~ c r d V +
IaBn, • Tn dA,
(1)
where £, a, T, n, u are the crack length, the strain energy density, the Cauchy stress, the outward unit normal to 0B and the displacement field, respectively; lastly u' is the derivative of u with respect to ~. In a bidimensional problem, if ~ is a path around the crack tip, i.e. a smooth non-intersecting path that starts and ends on the crack faces and surrounds the tip (Fig. 1), the quantity
J = e. Iy(aI - VuTT)m ds,
(2)
where e, I and m are the direction of propagation of the crack, the identity tensor and the outward unit normal to y, respectively, is the J-integral corresponding to path y [1]. The energy momentum tensor (aI - VuTT) is solenoidal, i.e. div(crI - VuTT) = 0
(3)
and thus, if the crack is straight, we have G(£) = J
and, furthermore, J is path-independent. 813
(4)
814
M. CHIARELLIet
al.
Fig. 1. Integration path and domains. According to eq. (4), ERR can be simply obtained by calculating the J-integral along any integration path ),; this result is useful for many engineering applications [2]. The definition of J in eq. (2) is meaningful if an elastic potential a exists or, in other words, if the work done on the body during a deformation process depends solely on the initial and final strains. In any cracked body made of a ductile material, like many metals, a region exists at the crack tip where yielding always occurs; the dimension of this region depends on the geometry, the level of loads and the constitutive law of the material. During a loading process in which some unloading occurs in this region, the work done on the body depends on the strain history, and the elastic fracture model is no longer appropriate; this situation occurs during a static monotonic loading where crack growth is stable, as well as under fatigue loads. For elastic-plastic materials with hardening, a definition of the energy available for a unit crack extension can be obtained by substituting the strain energy density a in eq. (1) with the work w done by the internal forces [3],
d lBwdV+ JoBu'.Tnda. G(e)=-~
(5)
As proved in ref. [3], for a straight crack in a bidimensional body we can write G(~) = J(y) - e. f div(wI - Vu]'T) dA, JL
(6)
for every path around the crack tip, with L the region enclosed by y where yielding occurred and .l(y) = e. lv(wI - VuTT)m ds.
(7)
For elastic-plastic materials the symbols J(~,) in eq. (6) and J in eq. (2) denote two completely different quantities; in fact, J is a function of the actual elastic state alone, while J(~,) depends on the whole deformation process. Moreover, Rice's integral J is path-independent while J(y) is not. On the other hand, for elastic materials, relations (6) and (4) coincide and so, the definition of G is an appropriate generalization of Rice's J-integral for elastic-plastic materials. In this paper a method for estimating G is presented that uses the flow theory of plasticity with general hardening laws. In this way, for all loading histories, both local and global unloading and reloading phenomena can be taken into account. The numerical procedure, described in Section 4, uses the results of an elastic-plastic finite element analysis of the cracked body. After each load increment, the stress T, the plastic deformation E p and the Odqvist parameter ~ (i.e. a parameter defined in the next section of the paper which indicates the amount of plastic deformation) are obtained and, by means of eq. (6), the value of G corresponding to each load increment is assessed. A closed-form expression of the work of the internal forces w as a function of T, E p, e, (, allows accurate computation of the surface integral in eq. (6). Examples are considered in which bidimensional cracked bodies are subjected to loads varying over time and the crack does not propagate during the application of the external load. The results obtained show that the difference between the elastic-plastic and elastic G values is very
Computation of crack extension energy rate
815
small when the load applied increases monotonically. On the contrary in the presence of cyclic loads, the difference is quite large. When cyclic loads with a constant amplitude are applied, different situations may occur. If the values of G become negative during unloading, which can be interpreted as a "closure effect" of the crack, the maximum value of G attained in the following load cycle diminishes. Other situations are possible for different material hardening laws and different values of the ratio between minimum and maximum loads applied. In the case of variable amplitude loading, the behaviour of G depends on the loading history and, in particular, on the presence of overloads that produce a decrease of G owing to the residual stresses during compression. The numerical results obtained are in agreement with experimental data when a realistic constitutive behaviour is assumed and this underscores the importance of adopting a suitable model for the hardening laws of the material. THE CONSTITUTIVE M O D E L For the reader's convenience, we shall briefly present certain elements of the flow theory of infinitesimal plasticity; we refer the reader to ref. [4] for a detailed account of the theory dealt with in the following. A deformation process or, more briefly, a history of duration ~ is a mapping E that assigns to each instant r < ~ the value of the infinitesimal strain E(~) reached at a fixed material point; we suppose that E(0) = O, so that all deformation processes are assumed to begin at some fixed initial state. The materials being considered here are elastic-plastic isotropic solids whose mechanical response to deformation processes is described by a rate-independent constitutive functional. We use T(T) to denote the stress at time z associated with history E by the constitutive functional. The kind of constitutive response is further specified by the notion of elastic range and plastic history. Elastic range e(z) at time z corresponding to history E is a set that contains the current value of strain E(z); its points are interpreted as the infinitesimal deformations from the reference configuration to configurations which are elastically accessible from the current one. We suppose that at each instant z the current deformation E(T) can be decomposed into the sum of an elastic part Ee(z) and a plastic part EP(z), E(r) = Ee(r) + EP(r),
(8)
having the following properties. EP('r) is a deviatoric tensor and belongs to e(T). The current stress depends linearly and isotropically on the elastic part of the deformation, T(r) = 2/z[E(r) - EP(r)] + ~.trE(r)I,
(9)
where trE(T) is the trace of E(T), I is the identity tensor and 2, # are the Lame' moduli of the material. Relation (9) reflects the classical hypothesis that the stress response to a purely elastic strain from the unstressed configuration reached after unloading at the current instant ~ is both unaffected by the past deformation process and completely determined by E(~) and EP('c). For each history E, ~'(r) =
IIEP(v')[[ dr'
(10)
is the length of the path described up to instant • by the plastic deformation tensor. ~ is called the Odqvist parameter. In view of the applications we have in mind, we accept the von Mises criterion. That is to say, we suppose that for each history E and for each T, the corresponding elastic range is the cylinder e(r) = {EIJE0 - C(r)ll _< p(U(r))},
(11)
816
M. CHIARELLIet al.
where Eo is the deviatoric part of E and (1) p is a differentiable, non-decreasing real function that depends on the material, is independent of history and describes the isotropic hardening; (2) for each history E, C is a deviatoric history. Moreover, we accept the classic linear kinematic hardening rule. That is to say, we suppose there exists a non-negative material constant ~b, such that for each history E and for each instant z we have C(r) = (1 + ~)EP(r).
(12)
In particular, a material for which p is a constant function and for which we have ~ = 0 is called ideally plastic. The set of constitutive hypotheses is completed by the flow rule which states that when EP('r) is different from zero, it is parallel to N(z), the outward unit normal to the elastic range at E(z), EP(v) = ((r)N(r), N(r) = [p(((v))l-l(E0(v)
- -
C(v)).
(13)
For each history E and for each z, T(z).l~(z) is the stress power and thus the work done by internal forces in the deformation process E up to time z is w(r) =
T(r'). l~(r') dr'.
(14)
In ref. [5] it has been proved that for each history E and for each time z, w ( z ) can be expressed as a function of E(z), EP(z) and ~(z). In fact, we have w(r) -- #ll(E(r) - EP(r))II 2 + lk(trE(r)) 2 + #lPIIEP(r)I[ 2 + 2/zw(((r)),
(15)
where ~o is the primitive of p such that 09(0) = 0.
CRACK EXTENSION ENERGY RATE FOR ELASTIC-PLASTIC MATERIALS
Let B be a regular homogeneous two-dimensional elastic-plastic body. We suppose that B contains an edge crack, represented at every instant z by the image of a smooth non-intersecting curve ~e, parameterized by arc length ct (Fig. 1). The length e = g(z) of curve ~ is a nondecreasing function of time z during the motion of B; ~e(0) and ~(g) represent the intersection of the crack with the boundary aB of B and the tip of the crack, respectively, while e(t~) = dle/do~
(16)
denotes the unit vector field tangent to the crack. In the interval [g0, gl], where the crack advances, we put dl wdA+ J u'.Tnds. G(g) = - ~-~ B ~B
(17)
J] G(~) d~
(18)
For each ~---~1,
0
can be interpreted as the work per unit thickness needed to increase the length of the crack from g0 to e; G is called the crack extension energy rate. For each g ~ [g0, el], let B~(g) = (xEBl((x, e) # 0},
(19)
with ~(x,g) defined by eq. (10), be the set of all points of B in which plastic deformation have
Computation of crack extensionenergyrate
817
taken place. Let us call B~(£) the plastic region; moreover, it should be noted that B~(£), in general, does not coincide with the subset of B in which we have EP(x,£) ~ 0. A curve ? is called a path around the tip if it is a smooth non-intersecting path that starts and ends on the crack, and includes the tip of the crack; we write m for the outward unit vector field normal to y. For each path ~ around the tip, let F be the subset of B enclosed by 7 and let L = bY~B~ be the intersection of F and B~ (Fig. 2). In ref. [3] it was proved that if the crack is straight, for each path 7 around the tip we have
G(e) =
)(y) - e. [ div(wI - VuTT) dA, JL
(20)
where f
J(7) = e. ] ( w I - VuTT)m ds.
(21)
If, for some £ ~ [£0,£1], the distance between B~(£) and aB is positive, there exist paths around the tip that include B~(Q. If 7 is one of these particular paths around the tip, F contains B~ and it therefore follows from eq. (20) that G(e) = J(y)
-
e. [ div(wI - VuTT) dA. JB~(t)
(22)
Moreover, since 7 does not intersect with B~, in view of eq. (17), J(7) can be calculated, rather than from eq. (21), from eq. (2) as in the case of hyperelastic materials, and thus we have
G(e) = J - e. [ JB
div(wI - VuTT) dA. ~(t)
(23)
No analytical assessment of G is available; it is to this purpose that a numerical procedure, presented in the next section, has been developed for bidimensional elastic-plastic bodies.
NUMERICAL EVALUATION OF G Numerical evaluation of G is carried out in two steps: finite element analysis of the cracked body and computation of the line and surface integrals according to eq. (20). The finite element code used in this research is NO.S.A. [6]. This code is particularly suitable for elastic-plastic analyses for both infinitesimal and finite deformations; in fact, it is very versatile in defining the input of the constitutive laws of elastic-plastic materials and, moreover, the values of ( are available at each Gauss point and for each load step. Computation of G is car-
i
fa(~)
i
7, --
-
s
411--
-
t
i
f2(g ) __
i
4"7 s.. 1 i
s'.
'
i 2
~{ (Plasticregion) Fig. 2. Plastic region and its intersectionwith the region F whose boundary is the integrationpath ~,.
Fig. 3. Sketch of the six parabolic functionswhich interpolate Ep and ~. O - - n o d e , * - - G a u s s point.
818
M. CHIARELLI et al.
ried out by GJINT2D, a code developed in-house, for which the finite element results at any load step are the input and the value of G is the output. Both the line and surface integrals in eq. (20) are computed with the help of a Gaussian technique. This is performed by means of an appropriate finite element mesh around the crack tip, defined in such a way that the integration paths pass through a set of Gauss points of the elements. In any Gauss point, the quantities to be integrated are expressed as functions of u, E, E p, T, ( and their derivatives. Assessment of these derivatives plays a central role in the whole numerical procedure implemented in code GJINT2D. The components of tensor E and their derivatives are computed starting from the nodal components of the displacement field by using the shape functions of the elements in the same way as in ref. [2]. Calculation of the derivatives of ( and of the components of E p by means of the shape functions is not reliable, because these quantities are known at the Gauss points only. So, for this purpose we employ a simple technique involving parabolic interpolation along coordinate directions ~ and r/ of the isoparametric space. As an example of this procedure, let us compute the derivatives of ( at the Gauss point labelled as number 5 in Fig. 3. We set f2(O = A~2 + B~ + C,
(24)
g2(0) = Dr/2 + Erl + F,
(25)
where the six coefficients of the quadratic functions (24) and (25) are evaluated by imposing the crossing conditions at points 4, 5 and 6 for f and at points 2, 5 and 8 for g. In this way, the desidered derivatives can be obtained by the chain rule, Og" df2 O~ dg2 Or/ Ox d$ Ox ~ do Ox' 8( dr2 8~ Oy -- d~ Oy
dg2 Or/ do Oy'
(26) (27)
where ax
O0 Or/
(28)
ax
is the inverse of the Jacobian matrix of the isoparametric transformation. Appendices A and B are devoted to solving the main problems that arise in computation of the surface integral; they are: (a) the Gauss points used to calculate the surface integral can be different from the original Gauss points of the element; (b) interpolation of ( and of the components of E p near the boundary of the yield region requires particular care. The general scheme of the whole procedure is shown in Fig. 4. In this paper, the reliability of the numerical procedure in calculating G can not be deduced from comparisons with results available in the literature. Nonetheless, we can point
Computation of crack extension energy rate
I
NO.S.A.F.E. PROGRAM ~ . ~
"<./'//////////////J
819
Solution of constitutive and equilibrium [ equations under the hypothesis of infinitesimal deformations.
Evaluation of ~, E P,E, T, u ] at each load step
i
~GJINT2D
Reading of the geometrical parameters and total displacements from a post file
MODULvIIIIII:~
T I Vectors and matrices containing 1 4, EP, E and T are initialized Loop on elements and their Gauss points
Evaluation of p(~) and c0(t)
I
Evaluation of Line Integral Functions
Evaluation of first derivatives of W, 4, E, EPand second derivatives of displacement components u i
I Y Evaluation of Line Integral
Evaluation of Surface Integral
L
I Evaluation of the Crack Extension Energy Rate
Fig. 4. Flow diagram with the principal steps for numerical evaluation of the Crack Extension Energy Rate.
out two conditions that are necessary (albeit not sufficient) for reliable computation of G and that are verified in all the examples discussed in the next section. They are: (1) the surface integral must be null in the first loading step when the material is elastic, starting from its annealed condition; (2) G must be path-independent, even though the line and the surface integrals are not.
EXAMPLES The examples shown in this section concern computation of G in the case of bidimensional bodies under plane strain conditions, with a crack that does not propagate during the loading process. Different loading histories are considered.
820
M. CHIARELLIet al. 19.65
I
10.9
---Q--" I I
26
Iiiiiiiiiiiiiiiill
82.5
Fig. 5. CT specimen undergoing a monotonic loading history for three differentvalues of kinematic hardening parameter ¢ [dimensionsin mm]. Compact tension specimen subjected to a monotonically increasing load
The CT specimen in Fig. 5 is made of a ductile material whose constitutive uniaxial behaviour is described by a bilinear curve with the following parameter values: Young's modulus
E = 113.000 MPa,
Poisson's ratio
o = 0.32,
Shear yield stress
• o = 461 MPa.
Here we assume the isotropic hardening function
P(O - (vo/2/z).,~(1 + ,80, with fl = 25. For kinematic hardening, we consider the following three values of the material constant if: ~p = 0.1, ¢ = 0.30 and ~p = 0.45. The finite element mesh and the integration paths considered in this analysis are shown in Fig. 6. The specimen is subjected to a monotonic load, increasing up to a value such that large deformations are avoided in the yield region around the crack tip according to the hypotheses of the infinitesimal t ~--ory of plasticity. The results (Fig. ~ show that for all values of ~kconsidered, the surface integral is small and so, in view of eq. (23), G can be well approximated by Rice's J-integral [eq. (2)]; moreover, as the contribution of the surface integral is always negative, Rice's J-integral is an upper bound of G. Experimental assessment of the fracture toughness of materials is carried out on specimens under a monotonically increasing load. According to the previous conclusion, plasticity does not play an important role and, then, when the yield region is not too large, characterization of material can be made with the methods of linear elastic fracture mechanics. The results in Fig. 7 also show that, when ~k increases, the surface integral tends to zero and then, G approaches Rice's J-integral; this agrees with the fact that the material becomes elastic when ~k increases indefinitely.
Computation of crack extensionenergyrate
821
III11111111IIII11111111I tl IIIIIIIIIIIIIII1111111I Imrlll-I I 1 I
II111111tl IIIIII111I]111 II lilllllll IIIIIIIIIit111
11 I!1111111IIIIIIIIt1~11 I
, 4
L~
-rrr~"nt
Illlttli-I-I'll II I II I:Yd-rtll IIIIIIIII I,'1 I I I I IIdllllll IIIIIIIII I,I I I I
7i
a
I
I lllllllll IIIIIIIII1~1 I I
I I rl I I llllllll i I IIIIIIII
III
IITITI1 rl-l'l i i IIIIII
II I I1
..........I I ' 1"',I IIIIIIII
A
crack Fig. 6. Mesh of the CT specimen shown in Fig. 5 and integration paths used for the evaluation of G.
Although it is not possible to compare the present results with analogous data in the literature, we observe that conditions (i) and (ii) pointed out at the end of the previous section are satisfied. In fact, as shown in Table 1, the surface integral in the first load increment is about 10 - 5 times the corresponding line integral; the same holds true for all the integration paths considered. The previous results are valid in the presence of small deformations only. For a CT specimen, when the external load is sufficiently high, large deformations occur in front of the crack tip and the infinitesimal theory of plasticity is no longer appropriate for computing the surface integral. For these problems, although the expression (22) of G is still valid, we need a theory of finite deformation plasticity to assess the work of the internal forces and the stresses; these aspects will be treated in a future paper. Worthy of note is the fact that in all the computations carried out with different loading conditions and constitutive laws, the values of G computed for different paths are the same within an error of about 0.1%, as Table 2 shows.
Center-cracked panel undergoing time variable loads The finite element mesh of a quarter of a CCP and the integration paths are shown in Fig. 8. Two materials are considered. The constitutive law of the first material is modelled by the bilinear curve shown in Fig. 9, having the characteristics of a steel alloy:
Young's modulus
E = 0.206 x 106 MPa,
Poisson's ratio
v -- 0.3,
Shear yield stress
Zo -- 165 MPa.
We suppose that the material is purely isotropically hardening, and apply the same function p (() used in the previous example. The loading history, shown in Fig. 10, consists of three cycles with R = 0. When the load increases during the first cycle, the surface integral is small and the results are similar to those obtained in the previous example. The effect of the loading history becomes very important during subsequent unloading (Fig. 11), when the absolute value of the surface integral increases and, according to eq. (20), the crack extension energy rate is reduced. Then, for 0_~O'~_~0"1, G
822
M. CHIARELLI et al.
Combined Hardeningwith []=25.0 and ¥=0.I 200
Lm~~ o ]
--o-
i
/
[IOhn^2 ] 100
i
l
l
l
l
l
l
}
l
l
l
l
l
l
,
,
,
i
. . . .
10
,
. . . .
i
20
. . . .
t
. . . .
30
40
Load [ KN 1
Combined I-lard~xingwid~~25.0 and ¥=0.3 200
[ KJ/m^2 ]
'"
'
'
i
. . . .
I
. . . .
I
. . . .
I
. . . .
I
. . . .
I--°--
I
. . . .
/
I
. . . .
J"
100
IO
20
30
40
Load [ K N ]
Combined Hardening with [3=25.0and ¥=0.45 200
' ' r ' r l
. . . .
I
. . . .
I
. . . .
i
. . . .
I
r
I
. . . .
.
.
I
. . . .
l
~
t
. . . .
I ' ' '
[I~J~^2 ] lOO
. . . .
i
. . . .
10
.
.
.
L
l
l
20
l
l
l
l
l
30
l
l
h
i
,
,
,
,
40
Load [ I ~ 1 Fig. 7. Results for the CT specimen undergoing a monotonic loading history with three different values of the kinematic hardening parameter ~ (Example 1).
Computation of crack extension energy rate T a b l e 1. O u t p u t o f the G J I N T 2 D
823
r o u t i n e f o r the C T s p e c i m e n u n d e r m o n o t o n i c l o a d i n g c o n d i t i o n . Units: k J / m 2. F i r s t load increment
Circuit number
Line integral
1 2 3 4 5 6 7 8 9 10 11 12
0.89279e0.89189e0.89207e 0.89334e0.89281e 0.89269e 0.89328e 0.89300e0.89286e0.89299e 0.89297e 0.89283e -
Surface integral
C.E.E.R.
0.43184e- 6 0.43159e - 6 0.39360e- 6 0.41509e-6 0.41506e- 6 0.39407e - 6 0.39396e-6 0.39444e- 6 0.39438e - 6 0.39011e-6 0.39488e-6 0.39386e- 6
0.89279e 0.89190e0.89207e 0.89334e0.89281e 0.89269e 0.89328e-1 0.89301e 0.89287e 0.89300e-1 0.89290e0.89283e -
1 1 1 1 1 1 1 I 1 1
becomes negative and, at the end of the first cycle, we have G(e)l~_-0<0
and
aG(e)/atrl,,=o=O.
The condition aG(e)/aal,~ = 0 = 0 at the end of the cycle indicates that during unloading, for small values of tr, we can write G(e)~-O~o+ ~ l a 2 as in the elastic case; then we deduce that during the first unloading the material in the yield region near the crack tip behaves elastically. At the beginning of the second cycle the material continues behaving elastically and, moreover, G(e) remains negative up to tr2. For O
T a b l e 2. O u t p u t o f the G J I N T 2 D
Circuit number 1 2 3 4 5 6 7 8 9 10 11 12
r o u t i n e f o r the C T s p e c i m e n u n d e r m o n o t o n i c l o a d i n g c o n d i t i o n . Units: k J / m 2. Fourth load increment Line integral
Surface integral
C.E.E.R.
0.55879 0.55822 0.55833 0.55913 0.55880 0.55872 0.55909 0.55892 0.55883 0.55892 0.55890 0.55881
- 0.12992e - 0.12992e - 0.12994e- 0.12993e - 0.12993e -0.12994e-0.12994e-2 -0.12994e- 0.12994e- 0.12994e -0.12994e-2 - 0.12994e-
0.55749 0.55692 0.55704 0.55783 0.55750 0.55743 0.55779 0.55762 0.55753 0.55762 0.55760 0.55751
2 2 2 2 2 2 2 2 2 2
824
M. CHIARELLI et al.
III
IIII1
--tll crack
++++++++++o ~,
800
.~
Fig. 8. CCP model used to study two different materials for the loading histories shown in Figs 10 and 13. The finite element mesh and a magnification o f the crack tip zone together with the integration paths are shown [dimensions in ram].
Center-cracked tension specimen undergoing variable amplitude loading In this example, we consider the CCT specimen in Fig. 15 made of a material with the following characteristics: Young's modulus
E = 73000 MPa,
Poisson ratio's
v = 0.32,
Shear yield stress
Zo = 185 MPa;
different hardening laws and load histories are considered. In the first case, the material is supposed to be purely kinematically hardening with = 0.28. Moreover, the load history considered, subdivided into 480 load increments, consists of five cycles followed by three overstresses and, finally, by five cycles with the same amplitude as the first ones (Fig. 16). The maximum values of G reached during the two series of five cycles are constant, while a decreasing trend is seen over the course of the three overload cycles because the overstresses cause an increase in the plastic region near the crack tip (Fig. 17). The shape of the crack faces is associated with the presence of residual strains close to the crack tip; as an example, Fig. 18 shows that the crack is still open when the load is null after the last overstress cycle. The decrease in the crack extension energy rate after an overstress cycle can be interpreted (even though in this paper we are dealing with a stationary crack) as a delay in crack propagation. In the present case of loading history always under tension, no contact between the crack faces ~o = 165 MPa 21.1,= 158.5 GPa T0
tga = 5654 MPa TO P°='4T 2~
= 0.00148
t g a tga) = 25 13: P o ( ~ ,
,
,
,
,
,
,
1 2 3 4 5 6 7
,
im
e(O/O0)
P u r e Isotropic H a r d e n i n g
Fig. 9. Stress-strain curve and numerical parameters used for the first material in Example 2.
Computation of crack extension energy rate
825 Pulsating Load Pure IsotropicHardening 6=25.0
200 Om~, =
108.86 MPa
6min •O
--
• ~aoadt~1 [ ILI/m^2
R=0
]
- - b - - loading 2 , • unloading 2 ---o-- l o a d ~ 3
100
Jd" G2 O 1 ~ P---
170'w~l( c~ If .~f
Gmin' 31 123 183 275 335 418 I
Load Increment N u m b e r Gmax
Fig. 10. First loading history applied to the CCP model shown in Fig. 8.
-100
~
~ 40
80
o [ lVlPa1 Fig. 11. Results for the CCP model undergoing the pulsating loading history shown in Fig. 10.
takes place in the yield region, so this delay can be attributed to the stress and strain fields in the yield region in front of the crack tip, rather than to the kinematics of the crack faces. In the second load history considered here, the maximum stress is constant and three series of three cycles, with R = 0.0, 0.3 and 0.5, are applied (Fig. 19); the third load history is made up of three series with the same amplitude as before and R variable according to Fig. 20. Owing to the assumption of pure kinematic hardening, the decrease in G at the end of each cycle o f the first series is nearly constant in both loading histories, as shown in Figs 21 and 22. In the second series of cycles, the same occurs, even though the decrease in G per cycle is much smaller. In the third series the load amplitude is 50% of the first, the stress variations in the yield region are always inside the elastic range and, so, no decrease in G occurs during these cycles. However, the effects of the loading history on the variation of G in the last series of cycles is very important; for example, in the case of R = 0, the delay effect is so marked that the last cycles produce nearly no crack extension. The importance of the material constitutive law is clearly shown by considering a combined hardening, with: 13 = 25, ~p = 0.28. The response of the material is completely different, as shown in Fig. 23: in the first three cycles the reduction of G is much smaller than in the previous case and the rate of decrease at the end of each cycle is not constant. In the next six cycles, the material behaves elastically owing to the increase in the radius of the elastic range, and the rate of energy available after the first group of cycles is higher than in the case of pure isotropic or pure kinematic hardening. This behaviour seems more realistic.
Young Modulus
73000 MPa
"to = 144.3 MPa Xmaxffi 283 MPa F.max =
0.018
To each of the different ranges of the Odclvist parameter values corresponds a
different isotzopic hardening parameter ~. 10
20
e (o/oo)
Pure lsotropic Hardening
Fig. 12. Stress-strain curve used for the second material treated in Example 2.
120
826
M. CHIARELLI et al.
(~
/'N
~,
~
Omax = 98.4 M P a Omin = 68.8 M P a
~mig
R = 0.7
21
49 77 1 5 133 1 1 Load Increment Number
Fig. 13. Second loading history applied to the CCP model shown in Fig. 8.
Results for the CCP model undergoing variable amplitude loading history Pure Isotropic Hardening 500 .... . . , . . . . . . , • .
[KJlm^2 ] 300
IOO
(followil~ ~ unloatang h
l
i
i
20
l
i
i
l
i
i
i
40
l
i
i
60
i
loadil~"
pnaa~ i
80
l
i
l
100
a [MPa]
Fig. 14. Results for the CCP model, with the stress-strain curve depicted in Fig. 12, undergoing the variable amplitude loading history shown in Fig. 13.
\/ i i
too
i
~
Sl~men
--~ thickness
Ultimate stress ~/ield stress UlIimate elongation
crack
6.a5 mr. 469 MPa
324 MPa 0.20
Fig. 15. CCT specimen used in Examples 3 and 4. The finite element mesh and the crack tip zone with the integration paths used are shown [dimensions in mini.
Computation of crack extension energy rate
827
160 o
[MPal 120
80
40
A ....
i ....
B
L ....
i ....
100
i ....
C
i ....
200
L ....
i ....
aoo
400
Load increment number Fig. 16. First loading history applied to the CCT specimen shown in Fig. 15 (Example 3).
Results for the CCT specimen undergoing loading history with overloads Pure Kinematic Hardening - ¥ = 0.28 15
9 [KJ/m^21 lO
-5
. . . .
I ' ' J l ~
50
. . . .
100
'
. . . .
150
i
. . . .
200
J
. . . .
260
i
. . . .
300
i
. . . .
350
400
Load Increment Number Fig. 17. G vs load increment for the specimen shown in Fig. 15 undergoing the loading history shown in Fig. 16.
1.2 V / rdCL 1 0.8
0.6
.........................
. ........................................................
÷ ............................
i. . . . . . . . . . . . . . . . . . . . . . . . . . .
0.4 ...................................................
i ..............................................................................
0.2 0 0
0.2
0.4
0.6
0.8
(distance from the crack tip) / (half crack length) Fig. 18. Profile of the crack face for the CCT specimen shown in Fig. 15. This profile is obtained unloading the panel completely, after point B of the loading history shown in Fig. 16 is reached.
828
M. C H I A R E L L I et al. 200
200
a [MPa]
(~ [MPa] 160
120
120
81) 4O
40
0
i
0
100
I
I
200
,
I
,
I
300
,
I
0
,
500
400
0
100
200
Load IncrementNumber
300
500
400
Load Increment N u m b e r
Fig. 19. Second loading history applied to the C C T specimen: variable R (Example 3).
Fig. 20. Third loading history applied to the CCT specimen: R = 0 (Example 3).
CONCLUSIONS
On the basis of the flow theory of plasticity, a path independent parameter, the Crack Extension Energy Rate (CEER or G), taking into account local and global unloading of stress state has been defined. A computer code, GJINT2D, has been developed which evaluates G from the results of elastic-plastic finite element analysis. The CEER has been assessed for generic elastic-plastic materials with isotropic, kinematic or combined hardening laws, in the case of plane problems with a stationary crack. The results show that if a monotonic loading condition is considered, the CEER does not greatly differ from the Rice J-integral; moreover, when unloading is considered, a difference, at times large, does appear even in the presence of small scale yielding. In the case of cyclic loading, the influence of the hardening law is so crucial, that the cyclic stress-strain curves of the material must be known very accurately. The occurrence of overloads gives rise to a decrease in CEER and, then, the delay in crack growth can be linked to the G vs load curves.
Results for the loading history with variable R Pure Kinematic [-Iarde~ting- ¥ = 0.28 2O
•
.
i
.
.
,
i
'
'
i
•
'
'
i
.
Results for the loading history with R=0 Pure Kinematic Hardening - ¥ = 0.28 .
j
[ KJImA2 ] 1o
0
40
80 120 a [MPa]
.
.
.
.
2o
i
,
1o
Fig. 21. G vs load for the C C T specimen shown in Fig. 15 undergoing the loading history shown in Fig. 19: pure kinematic hardening (Example 3).
200
-10
'
i
•
'
'
J
'
'
t
j
[ KJlm^2 ]
160
,
i
b
i
I
40
b
i
J
I
i
80
i
i
I
120
*
i
,
I
i
i
160
cr [MPa] Fig. 22. G vs load for the C C T specimen shown in Fig. 15 undergoing the loading history shown in Fig. 20: pure kinematic hardening (Example 3).
*
200
Computation of crack extension energy rate
829
Results for the loading history with variable R Combined Hardening - [~= 25.0, ¥ = 0.28 20
//
[KJlm^2] 10
~
0
t
i
40
i
80 o
i
t
,
i
i
120
,
20o
160
[MPal
Fig. 23. G vs load for the CCT specimen shown in Fig. 15 undergoing the loading history shown in Fig. 19: combined hardening (Example 3).
REFERENCES 1. Rice, J. R., A path independent integral and the approximate analysis of strain concentrations by notches and cracks. 3". appl. Mech., 1968, 35, 379-386. 2. Chiarelli, M. and Frediani, A., A computation of the three-dimensional J-integral for elastic materials with a view to applications in fracture mechanics. Engng Fracture Mech., 1993, 44, 763-788. 3. Chiarelli, M., Frediani, A. and Lucchesi, M., On the crack extension energy rate in elastic-plastic bodies. Defect assessment in components. Fundamentals and applications, ESIS/EGb9, ed. J. G. Blauel and K. H. Schwlbe. Mechanical Engineering Publications, London, 1991, pp. 87-99. 4. Lucchesi, M. and Podio Guidugli, P., Materials with elastic range. Part III: approximate constitutive relations, Arch. Rat. Mech. AnaL, 1992, 117, 53-96. 5. Lucchesi, M., Free-energy functions for elastic-plastic material elements. Q. appl. Math., 1993, 60, 299-318. 6. Guidotti, P., Lucchesi, M., Pagni, A. and Pasquinelli, G., Elastic-plastic behaviour with work-hardening: an appropriate model for structural software. Meccanica, 1984, 19, 43-51.
(Received 6 October 1994) APPENDIX
A
Correct application of Gaussian integration rules in calculation of the surface integral requires all elements crossed by integration paths to be modified and the new Gauss points to be properly placed [2]. For this reason, the interpolation scheme becomes quite complex. As an example, crossing an element in the t/direction, as in Fig. 3, once we get the isoparametric coordinates of new Gauss points, for a generic component E p of plastic strain E p at the point 1' (Fig. A1), we have
e p =A(~I'), ~l, ¢~1, E~ =g,(ol'), ~,' = 0 t ;
tit
I' ?,
i, , l
il
i! ' Y3-
(AI)
integra~onpath
f'
il
il
',
i"i
5"ii
v,,j
I/
I
'I
Fig. AI. Location of Gauss points in a modified element because of a crossing in the positive h direction. • - - generic Gauss point of unmodified element, * - - generic Gauss point of modified element.
Fig. A2. Location of Gauss points in a modified element because of a comer crossing, s - - generic Gauss point of unmodified element, * - - generic intermediate control point, [] - - generic Gauss point of modified element.
M. CHIARELLI et al.
830
V 1
~=-1
2
~i\
interpolating parabola
~=0
!~
3
1 ~=-1
~=1
Fig. BI. If ( = 0 at Gauss point no. 3, with a standard interpolation procedure we would obtain a wrong derivative.
;' E ~
/
parabolawho6e vertex coincides
[ withintermediate pointofthe tern
2 ~,=0
a ~=1
Fig. B2. Interpolation example with two points o f the term for which ( = 0 or E ~ = 0.
the apex indicates that the Gauss point of a modified element must be considered and gl(r/) is the parabolic function interpolating E ~ on points 1', 4' and 7'. So, its derivative with respect to x at 1' is
OE p dfl O~ . dgl Or/ 0x = ~'~xx "t ~ ~x"
(A2)
It is note-worthy that, when the derivatives of E p and ~ are calculated for a modified element, all terms o f the inverse of the Jacobian matrix must be computed in the isoparametric reference frame of the complete element. p For a corner crossing (Fig. A2), E ij and ( and their derivatives are computed by means o f a double interpolation scheme. If we use two apices to indicate the new Gauss points and one apex for intermediate control points (Fig. A2), at 1" we have E~=gj(r/1,,), r/J" # r / l ' = r / l E p =fl(~{'), ~1" = ~ t ' # ~ l
(A3)
and then
OE~jP _ d f l O~ Ox d~ Ox
dgl0r/ dr~ Ox'
(A4)
where, as in the previous case, gl(r/) interpolates E p in l', 4' and 7', whileJ~(O interpolates E p0 in l", 2" and 3".
APPENDIX B Approaching the boundary of the plastic region, E P, ( and their gradients tend to zero. Consequently, the previous interpolation scheme needs to be modified on the elements lying close to this boundary. In fact, supposing ( = 0 at the Gauss points 3, 6 and 9 (Fig. 3), if we simply interpolated ( values in the ~ direction, we would get a parabolic function, resulting in the curve drawn roughly in Fig. B1, whose derivative at Gauss point 3 is different from zero; the same holds for points 6 and 9. To avoid this kind of error, whose numerical influence is greater for a small plastic region size, we must modify the procedure. Cases of interest are:
(1) = 0 o r e 0P = 0 at two Gauss points o f the current interpolation term;
(2) ( = 0 or E ~ = 0 at a single Gauss point of the current interpolation term. In the first case, we select a parabola whose vertex coincides with intermediate integration point: for instance point 2 of terms 1, 2, 3 (Fig. B2). In the second case, we select a parabola as shown in Fig. B3; the possible error introduced at the intermediate point can be accepted. Therefore, before calculating the derivatives, a check is made to locate the Gauss points where ( or E ~ are null.
~
' E~I parabola withvertex , / ~ a t pointforwhich
2
a
{
Fig. B3. Interpolation example with one point of the term for which [ = 0 or E ~ = 0.