Copyright @ 2013 ACMECohesive zone model is a well known concept in nonlinear fracture mechanics of elasto-plastic materials. In contrast to that, we discuss a development of the cohesive zone model to linear, but time and history dependent, materials. The stress distribution over the cohesive zone satisﬁes a history dependent rupture criterion for the normalised equivalent stress, represented by a nonlinear Abel-type integral operator. The cohesive zone length at each time step is determined from the condition of zero stress intensity factor at the cohesive zone tip. It appeared that the crack starts propagating after some delay time elapses since a constant load is applied to the body. This happens when the crack tip opening displacement reaches a prescribed critical value. A numerical algorithm to compute the cohesive zone and crack length with respect to time is discussed and graphs showing the results are give

Hakim, L

Mikhailov, SE

English

Brunel University Research Archive

International Conference on Computational Mechanics (CM13) 25-27 March 2013, Durham, UKCohesive Zone Models in History Dependent MaterialsLayal Hakim∗, Sergey MikhailovDept. of Mathematics, Brunel University London, UKAbstractCohesive zone model is a well known concept in nonlinear fracture mechanics of elasto-plastic materials. In contrast tothat, we discuss a development of the cohesive zone model to linear, but time and history dependent, materials. The stressdistribution over the cohesive zone satisfies a history dependent rupture criterion for the normalised equivalent stress,represented by a nonlinear Abel-type integral operator. The cohesive zone length at each time step is determined from thecondition of zero stress intensity factor at the cohesive zone tip. It appeared that the crack starts propagating after somedelay time elapses since a constant load is applied to the body. This happens when the crack tip opening displacementreaches a prescribed critical value. A numerical algorithm to compute the cohesive zone and crack length with respect totime will be discussed and graphs showing the results will be given.Keywords: Cohesive zone, Time dependent fracture, Abel integral equation1. IntroductionThe cohesive zone in a material is the region between 2 surfaces ahead of the crack tip that are separating due toexternal loading but are pulled together by cohesive stresses σ, while the crack faces are traction-free, see Figure 1.Eventually, the crack will propagate through the cohesive zone points, where the criterion of the cohesive zone breakageis reached.Figure 1: Cohesive zoneThe Dugdale-Leonov-Panasyuk (1959-1960) (DLP) model was the first model introducing perfectly plastic cohesivezones, with constant cohesive stresses, σ = σy. This model has been modified and used in many applications in fracturemechanics. Another popular cohesive model is the Barenblatt (1962) model. The 3 main components needed to studycohesive zone models of the DLP type are: (i) the constitutive equations for the bulk of the material; (ii) the criterion onthe stress (history) on the cohesive zone for the cohesive zone to appear and propagate; (iii) the criterion for the cohesivezone to break and the crack to propagate.∗Corresponding authorEmail address: layal.hakim@brunel.ac.uk (Layal Hakim)Draft extended abstract submitted to CM13 January 29, 20132. Problem FormulationThe following problem is an extension of the DLP model to linear, but time and history dependent, materials. To thisend, we will replace the DLP cohesive zone stress condition, σ = σy, with the conditionΛ(σ; t) = 1, (1)whereΛ(σ; t) = βbσβ0∫ t0|σ(τ)|β(t − τ) βb−1dτ1β(2)is the normalised equivalent stress, |σ| is the maximum of the principal stresses, and t denotes time. The parameters σ0and b are material constants in the assumed power-type relation t∞(σ) = (σ/σ0)−b between the rupture time t∞(σ) and theconstant uniaxial tensile stress applied to a body without cracks. The parameter β is a material constant in the nonlinearaccumulation rule for durability under variable load, see [1].Note that relations (1)-(2) were implemented in [2] and [3] to solve the corresponding crack propagation problem butwithout a cohesive zone, i.e. it was assumed that when condition (1) is reached at a point, this point becomes part ofthe crack. However, such an approach appeared to be inapplicable for b ≥ 2. In this paper, a cohesive zone approach isdeveloped instead, in order to cover the larger range of b values relevant to structural materials. In this approach, whencondition (1) is reached at a point, this point becomes part of the cohesive zone.Let the problem geometry be as in Figure 1, i.e, the crack occupies the interval [−a(t), a(t)] and the cohesive zoneoccupies the intervals [−c(t),−a(t)] and [a(t), c(t)] in an infinite plane loaded at infinity by traction q in the directionnormal to the crack, which is constant in x, applied at the time t = 0 and kept constant in time thereafter.The cohesive zone condition (1)-(2) at a point x on the cohesive zone can be rewritten as∫ ttc(x)σβ(x, τ)(t − τ) βb−1dτ = bσβ0β−∫ tc(x)0σβ(x, τ)(t − τ) βb−1dτ, t ≥ tc(x), (3)where a(t) ≤ |x| ≤ c(t). Here, tc(x) denotes the time when the point x becomes part of the cohesive zone. Equation (3) isan inhomogeneous nonlinear Volterra integral equation of the Abel type with unknown function σ(x, t), t ≥ tc(x).Assuming that the bulk of the material is linearly elastic and applying the results by Muskhelishvili (see [4], Section120), we have for the stresses ahead of the crack,σ(x, t) =1√x − c(t)x√x + c(t)q − 2pi∫ c(t)a(t)√c2(t) − ξ2x2 − ξ2 σ(ξ, t)dξ , t ≤ tc(x), (4)where |x| > c(t). A sufficient condition for the normalised equivalent stress, Λ, to be bounded at the cohesive zone tipis that the stress intensity factor, denoted by K is zero at the cohesive zone tip. Multiplying the stress in equation (4) by√x − c(t) and taking the limit as x tend to c(t) yieldsK(t) =√c(t)2q − 2pi∫ c(t)a(t)1√c2(t) − ξ2σ(ξ, t)dξ .We will normalise σ, t and x, a and c using the following transformations but keeping the same notations:σ→ σq, t → tt∞(q)= t(qσ0)b, x→ xa(0), a→ aa(0), c→ ca(0).Thus, after the normalisation, we state the following principle equations for the considered problem:(a) the cohesive zone condition on stresses using equation (1):∫ ttc(x)σ(x, τ)β(t − τ) βb−1dτ = bβ−∫ tc(x)0σβ(x, τ)(t − τ) βb−1dτ a(t) ≤ |x| ≤ c(t); (5)(b) the expression for the stress ahead of the cohesive zone:σ(x, t) =x√x2 − c2(t)1 − 2pi∫ c(t)a(t)√c2(t) − ξ2x2 − ξ2 σ(ξ, t)dξ |x| > c(t); (6)2(c) the zero stress intensity factor at the cohesive zone tip:K(t) = −√2c(t)pi∫ c(t)a(t)1√c2(t) − ξ2σ(ξ, t)dξ +√c(t)√2= 0. (7)3. Cohesive Zone Growth for a Stationary CrackFirst, we will consider the stage, when the crack is stationary, a(t) = 1, and thus only the cohesive zone grows withtime. This is related to the times before the crack starts propagating. Our aim is to find the cohesive zone tip position c(t)and the crack opening at the crack tip a(t) = 1.3.1. Numerical MethodFirst, we will introduce a time mesh, with nodes ti, for i = 0, 1, 2, 3, ...n. At each time step ti, we use the secant methodto find the roots, c(ti), of the equation K(c(ti), ti) = 0. For each approximation within the secant iterations, we calculatethe stress intensity factor using equation (7). To compute the integral∫ c(ti)a01√c2(ti) − ξ2σ(ξ, ti)dξ, (8)we linearly interpolate σ(ξ, ti) on the cohesive zone between ξ = c(tk) and ξ = c(tk+1), where k = 0, 1, 2, ...i − 1. On theother hand, to find σ(c(tk), ti), at each c(tk) for ti > tk, we use the Abel integral equation (5). After linearly interpolatingthe function σβ(c(tk), τ) between τ = t j and τ = t j+1, where j = 0, 1, 2, ...k − 1, in the integral∫ tk0σβ(c(tk), τ)(t − τ) βb−1dτin the right hand side of equation (5), we use an analytical formula to solve the equation. To this end, in turn, we need tofind σβ(c(tk), t j) for t j < tk. This is obtained using equation (6) (since c(tk) > c(t j)), where the integral∫ c(t j)a0√c2(t) − ξ2c(tk)2 − ξ2 σ(ξ, t j)dξwill be calculated similar to integral (8). This means we will linearly interpolate σ(ξ, t j) between ξ = c(tm) and ξ = c(tm+1)for m = 0, 1, ... j − 1.All programming was implemented in MATLAB.3.2. Cohezive Zone Growth Numerical ResultsFigure 2 shows the results obtained for the case b = 4, using various mesh sizes.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1100101102tc(t)  n=25n=50n=100n=200n=400n=8000.665 0.67 0.675 0.68 0.6857.58587.76257.94338.12838.3176tc(t)  n=25n=50n=100n=200n=400n=800Figure 2: CZ tip position vs time for b = 4 and different meshes (non-propagating crack)3The graphs in Figure 3 show the results obtained when considering 3 different cases for the parameter β using n = 500time steps on the interval 0 ≤ t ≤ 1. The graph on the right hand side of Figure 3 is a zoomed part of the graph on the lefthand side in the vicinity of the point t = 0.56.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1100101102103104tc(t)  β=b/2β=b/4β=b/80.55 0.555 0.56 0.565 0.574.89785.01195.12865.2481tc(t)  β=b/2β=b/4β=b/8Figure 3: CZ tip position vs time for b = 4 and different β (non-propagating crack)Figure 4 shows the stress behaviour with respect to time at the point x = c(0.6), i.e., tc(x) = 0.6.0 0.2 0.4 0.6 0.8 10.80.911.11.21.31.41.5tσ(c 0.6,t)  β=b/2β=b/4β=b/8Figure 4: σ(c(0.6), t) vs time for b = 4, n = 5003.3. Crack OpeningNow, we want to find the displacement jump at the crack shores, which is also referred to as the crack opening. Wedenote the crack opening by δ(x, t), see Figure 1. Using the representations by Muskhelishvili (see [4], Section 120), itcan be deduced that the crack opening is given byδ(x, t) = [u2(x; c(t))] = [u(q)2 (x; c(t))] + [u(σ)2 (x; c(t), t)]where[u(q)2 (x; c)] =4q(1 − ν2)E√c2 − x2and[u(σ)2 (x; c, t)] =4(1 − ν2)piE(∫ caσ(ξ, t)Γˆ(x, ξ; c)dξ),whereΓˆ(x, ξ; c) = ln2c2 − ξ2 − x2 − 2 √(c2 − x2)(c2 − ξ2)2c2 − ξ2 − x2 + 2 √(c2 − x2)(c2 − ξ2).4In the above expressions, E and ν denote Young’s modulus of elasticity and Poisson’s ratio respectively. We will furthernormalise δ using the transformationδ→ δEaq(1 − ν2) . (9)Therefore, the normalised crack opening δ(x, t) at the crack tip, x = a(t), is given asδ(a(t), t) =4pi(pi√c2(t) − a2(t) +∫ c(t)a(t)σ(ξ, t)Γˆ(a(t), ξ; c(t))dξ).Linearly interpolating σ(ξ, ti) between c(tk) and c(tk+1) for k = 0, 1, ...i − 1, we can calculate the crack opening, δ(1, ti), atthe tip a(t) = 1 of the stationary crack. The graphs in Figure 5 show the results obtained.0 0.2 0.4 0.6 0.8 10510152025tδ  β=2β=1β=1/2Figure 5: Crack tip opening vs time in the stationary crack for b = 4, n = 5004. Crack PropagationWe have, so far, assumed that the crack is stationary and only the cohesive zone is growing ahead of the crack.However, at some time instant, td, that will be referred to as the delay time, the crack opening reaches a critical value,δ = δc, and the crack starts to propagate while still having the cohesive zone ahead of it. The crack and cohesive zonewill not necessarily grow at the same rate. The parameter δc is known as the critical crack tip opening. It depends on thematerial and can be experimentally measured along with other material parameters. For example, for PMMA (also knownas plexiglas) one can take from literature δc = 0.0016mm, the Poisson ratio ν = 0.35 and the Young modulus E = 3.1GPa.Thus, normalising δc according to (9), we will have δc = 0.113. After determining td, the aim is to find the crack lengthand the corresponding cohesive zone length for times t > td.4.1. Numerical MethodConsider a time mesh with uniformly spaced time nodes, ti. During crack growth, at each step, the crack openingδ(a(ti), ti) equals to the critical value δc, which leads to the following equation4pi(pi√c(ti)2 − a(ti)2 +∫ c(ti)a(ti)σ(ξ, ti)Γˆ(a(ti), ξ; c(ti))dξ)− δc = 0, ti ≥ td. (10)For each ti, we apply the following algorithm to obtain a(ti) and c(ti). We use the secant method to find a(ti) by settingthe crack opening equal to the critical crack opening value, i.e., by solving equation (10) for a(ti). To do this, we need toknow c(ti) at each iteration. It is obtained using the secant method to set the stress intensity factor to 0; i.e., to solve theequation−√2c(ti)pi∫ c(ti)a(ti)1√c2(ti) − ξ2σ(ξ, ti)dξ +√c(ti)√2= 0, (11)5for c(ti). Note that we choose previous cohesive zone tip positions, c(tm), as initial approximations for a(ti) within thesecant algorithm. The advantage of doing this is that we already know the stress history at these previous points. Forinstance, for the first step after crack growth begins, it is reasonable to take the corresponding a(ti) approximations equalto c(t1) and c(t2). At some point, we will come across the step where a(ti) will exceed c(ti−1); this means that we can onlyuse c(ti−1) as one of the initial approximations for a(ti). For the steps when this occurs, we fix a(ti) to be equal to c(ti−1)and compute the corresponding ti and c(ti) by solving equation (10) and K(ti) = 0 respectively.4.2. Numerical resultsTo check the algorithm, we first used in our calculations the value δc = 1.13 for the normalised critical crack tipopening, which is 10 times larger than the corresponding value estimated for PMMA. For the case of b = 4, and β = 1,we obtained td = 0.0102 and present in Figure 6 the graphs showing coordinates of the crack tip and the cohesive zone tipas well as the cohesive zone length against time.0 0.005 0.01 0.015 0.02100101TimeLength (log scale)Length (log scale) vs Time for b=4, β=1  Crack LengthCZ tip coordinate0 0.005 0.01 0.015 0.0200.020.040.060.080.10.120.14Timecz lengthCZ Length vs Time for b=4, β=1  Figure 6: Length vs. time for b = 4, β = 15. ConclusionsThe solutions for c(t), a(t) and σ(x, t) converge as the mesh becomes finer. For the stationary crack stage, t < td, therate of the cohesive zone length growth decreases when β decreases. As expected, we have an increase, with time, of thecrack opening. Moreover, as β becomes smaller, the crack opening increases more slowly with time. For the growingcrack stage, t > td, we can see from Figure 6, that the crack growth rate increases, while the cohesive zone length decreaseswith time. The time, when the cohesive zone time becomes 0 seems to coincide with the time when the crack becomesinfinite and can be associated with the complete fracture of the body.References[1] Mikhailov S. E., Namestnikova I. V. History-sensitive accumulation rules for life-time prediction under variableloading, Archive of Applied Mechanics, Vol. 81, 2011, 1679-1696.[2] Mikhailov S.E., Namestnikova I.V., Local and non-local approaches to creep crack initiation and propagation, in:Proceedings of the 9th International Conference on the Mechanical Behaviour of Materials, Geneva, Switzerland,2003.[3] Hakim L, Mikhailov S.E., Nonlinear Abel type integral equation in modelling creep crack propagation, in: IntegralMethods in Science and Engineering: Computational and Analytic Aspects, editors: C. Constanda C. and P. Harris,Springer, 2011, 191-201.[4] Muskhelishvili N.I., Some Basic Problems of the Mathematical Theory of Elasticity, Noordhoff International Publish-ing, The Netherlands, 1954.6

History-sensitive accumulation rules for life-time prediction under variable loading,

Local and non-local approaches to creep crack initiation and propagation, in:

Nonlinear Abel type integral equation in modelling creep crack propagation, in:

Some Basic Problems of the Mathematical Theory of Elasticity,

Cohesive zone models in history dependent materials

http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.426.6013

Cohesive zone models in history dependent materials

Abstract

Similar works

Full text

Available Versions

Brunel University Research Archive