An implicit stress gradient plasticity model for describing mechanical behavior of planar fiber networks on a macroscopic scale

Engineering Fracture Mechanics 77 (2010) 1240–1252

Contents lists available at ScienceDirect

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

An implicit stress gradient plasticity model for describing mechanical behavior of planar fiber networks on a macroscopic scale P. Isaksson * Department of Engineering Physics, Mid Sweden University, SE-851 70 Sundsvall, Sweden

a r t i c l e

i n f o

Article history: Received 10 August 2009 Received in revised form 23 December 2009 Accepted 12 March 2010 Available online 19 March 2010 Keywords: Nonlocal theory Gradient plasticity Network material Crack mechanics

a b s t r a c t The plasticity behavior of fiber networks is governed by complex mechanisms. This study examines the effect of microstructure on the macroscopic plastic behavior of two-dimensional random fiber networks such as strong-bonded paper. Remote load is a pure macroscopic mode I opening field, applied via a boundary layer assuming small scale yielding on the macroscopic scale. It is shown that using a macroscopic classical homogeneous continuum approach to describe plasticity effects due to (macroscopic) singular-dominated strain fields in planar fiber networks leads to erroneous results. The classical continuum description is too simple to capture the essential mechanical behavior of a network material since a structural effect, that alters the macroscopic stress field, becomes pronounced and introduces long-ranging microstructural effects that have to be accounted for. Because of this, it is necessary to include a nonlocal theory that bridges the gap between microscopic and macroscopic scales to describe the material response in homogeneous continuum models. An implicit stress gradient small deformation plasticity model, which is based on a strong nonlocal continuum formulation, is presented here that has the potential to describe the plasticity behavior of fiber networks on a macroscopic scale. The theory is derived by including nonlocal stress terms in the classical associated J2-theory of plasticity. The nonlocal stress tensor is found by scaling the local Cauchy stress tensor by the ratio of nonlocal and local von Mises equivalent stresses. The model is relatively easy to implement in ordinary finite element algorithms for small deformation theory. Fairly good agreements are obtained between discrete micromechanical network models and the derived homogeneous nonlocal continuum model. Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction It is known that classical homogeneous continuum mechanical descriptions of fiber-based network materials, such as paper, textiles or human tissue, cannot fully describe deformations near macroscopic defects, cf. [1]. The fibers introduce longranging microstructural effects that alter the stress field on a macroscopic scale when the inhomogeneous body becomes loaded. The fundamental mechanism behind this phenomenon has received little attention in the literature. There are several candidate mechanisms that may alter the deformation field in a network, as compared to an ideal homogeneous continuum, such as plastic straining of the fibers, fiber and bond breakage or wrinkling and buckling. In tensile experiments on fiber-based materials such as paper it is well known that the stress–strain curve will deviate from linear behavior and upon unloading there will be a permanent deformation, commonly referred to as plasticity. Hence, in the context of elastic–plastic fiber-based materials, a plastic process is defined as the irreversible process when there is an * Tel.: +46 60148882. E-mail address: [email protected] 0013-7944/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.engfracmech.2010.03.023

P. Isaksson / Engineering Fracture Mechanics 77 (2010) 1240–1252

1241

Nomenclature x1, x2; r, h Cartesian coordinates; polar coordinates X; R problem domain; radius to boundary f ðÞ; ðÞ nonlocal counterpart of variable (  ); dimensionless nonlocal variable characteristic internal length; theoretical estimation of l l; l0 r0 ; re ; /; k yield stress; von Mises equivalent stress; yield surface; Lagrange multiplier rij; sij Cauchy stress tensor; deviatoric Cauchy stress tensor eeij ; epij ; w; ui elastic strain; plastic strain; strain energy density; displacements E; m; l; Lijkl macroscopic properties: Young’s modulus; Poisson’s ratio; shear modulus; stiffness tensor L; h; Ef; r0f fiber properties: length; width; Young’s modulus; yield stress c; Ls; N network properties: coverage; average segment length; total numbers of fibers K I ; K Ip ; K I macroscopic stress intensity factor; KI at onset of plastic flow; KI at limit for small scale yielding dij; d(  ) Kronecker delta; dirac delta function G(  ); B0(  ) Green’s function; modified zero-order Bessel function of second kind total potential energy in homogeneous nonlocal model; network model U; U# estimation of macroscopic plastic zone size in homogeneous nonlocal model; network model rp ; r# p macroscopic force; force at which full plastic yielding occur P; Pc

irreversible straining of the fibers themselves, cf. [2,3]. Obviously, this definition relies on the assumption that the fibers can deform irreversible before they eventually break. In particular, for strong-bonded elastic–plastic fiber-based materials (i.e. in situations when the bonds connecting the fibers are strong compared to the fibers) subjected to traction it is reasonable to assume that there is a significant amount of plastic deformation of the fibers before fracture occurs. This is because the strong-bonded fiber-network structure results in a high degree of axial stress in the fibers. Energy dissipation owing to plasticity processes in the network is in this case assumed significant and the macroscopic behavior of the material is strongly non-linear: a behavior that may be captured by an appropriate plasticity theory on the macroscopic scale. Typical examples of strong-bonded elastic–plastic fiber network materials include printing and packaging papers and some textiles. The literature shows a variety of approaches for studying this problem. Even though some important results have been reported, cf. [4,5], the knowledge of deformation processes in e.g. paper materials is still in its infancy and needs further exploration. For example; linear fracture mechanics and non-linear fracture mechanics have commonly been used to characterize macroscopic failure phenomena of paper possessing certain defects, cf. [6–8]. Despite these efforts, it is still not established in the materials physics community whether the current non-linear fracture mechanics formulation of heterogeneous fiber-based materials, which essentially is imported from theories developed for homogeneous materials, can be directly applied to describe fracture and deformation behavior. One should have in mind that the mechanical behavior on a microscopic level in network materials such as paper are significant different from those in, say, ductile metals, which are closely connected to dislocation movements in the material. Not only are the basic phenomena entirely different (permanent deformation of fibers vs. dislocation movements in atomic lattices), the mechanisms themselves are acting on totally different scales (1 mm vs. 10 Å). Addressed in this investigation is a nonlocal plasticity effect on the macroscopic scale that originates from the microstructural properties of the material. Firstly, there is a size-defect dependence in a network material originating from pores and fibers that is not captured by classical homogeneous plasticity continuum theories such as the well-known J2-flow theory. Secondly, since the network consists of relatively stiff fibers connected to other fibers located at remote distances, the fibers themselves introduce nonlocal structural effects in the material as the overlapping fibers apply long-range actions on each other in a complex manner. Thus, as compared to a homogeneous material that often is approximated as a classical continuum, certain inherent lengths influence the mechanical field in a network material. In the last years a number of gradient, or nonlocal, plasticity theories have been presented (among others: Fleck and Hutchinson [9,10]; Fleck et al. [11]; Gao et al. [12]; Huang et al. [13]; Mentzel and Steinmann [14]). The various proposed theories are fairly different with respect to the structure of the equations, however, each aims to capture boundary layer phenomena related to phase and grain boundaries, or slips in single crystals, within small deformation formulations. Gudmundson [15] formulated a small deformation isotropic strain gradient plasticity theory, where the plastic flow direction is governed by a certain microstress and not the deviatoric Cauchy stress that is assumed in the classical J2-flow theory. The microstress is assumed work conjugate to plastic strains. All the abovementioned studies are considered weakly nonlocal in character because the governing equations in a point include gradients explicitly obtained from state variables in an arbitrary small neighborhood of that point. In the present study, however, a strong nonlocal small deformation theory is used, which is derived by implicitly including gradients on a considerable larger distance than in the weakly nonlocal models. This approach is physically motivated since fibers in e.g. cellulosic networks often have the length of millimeters as compared to microns, or even smaller lengths, in the case of dislocation movements in metals, as discussed above.

1242

P. Isaksson / Engineering Fracture Mechanics 77 (2010) 1240–1252

All nonlocal theories provide a constitutive framework that bridges the gap between classical homogeneous continuum theory and micromechanical modeling. A nonlocal theory abandons Saint–Venant’s principle of local action, i.e. the assumption that the mechanical state in a given point in a material is uniquely determined by the state in that point only as it includes interactions with a neighborhood to that point. The concept of strong nonlocal field theory, cf. Pijaudier-Cabot and Bazant [16] or Eringen [17], has emerged as an effective mean for regularizing boundary value problems and capturing size effects. The subject has been reviewed by, among others, Jirásek and Rolshoven [18,19]. Different techniques exist to incorporate strong nonlocality in continuum formulations. Examples of such theories include nonlocal integral formulations (cf. Kröner [20]; Eringen and Edelen [21]; Silling [22]) and implicit gradient theories (cf. Peerlings et al. [23]) that is the foundation of the work in this study. Further, an aim of this work is to develop a constitutive model that acknowledges the underlying physics of fiber-based network materials and is able to quantitatively predict plastic deformations for macroscopic structures containing defects. In addition, the model should be appropriate to incorporate into a finite element framework in a relatively easy manner. 2. Constitutive nonlocal model: an implicit stress gradient plasticity approach 2.1. Foundation  in a point (x1, x2) of an Lets start with the usual assumption for strong nonlocal field theories that a nonlocal counterpart s arbitrary local state variable s over a surrounding two-dimensional domain X is given by

sðx1 ; x2 Þ ¼

1

Z

W

X

uðx01 ; x02 ; x1 ; x2 Þsðx01 ; x02 ÞdX; where W ¼

Z X

uðx01 ; x02 ; x1 ; x2 ÞdX:

ð1Þ

 equals s for a homogenous Here ðx01 ; x02 Þis the position of the infinitesimal area dX and the scaling factor W provides that s state in X. The weight function u is assumed homogenous, isotropic and Gaussian-distributed according to:

u ¼ ð2pl2 Þ1 exp½q2 =ð2l2 Þ;

ð2Þ

where the distance q is given by

q ¼ jðx1 ; x2 Þ  ðx01 ; x02 Þj;

ð3Þ

and a characteristic length l controls the range of nonlocal actions in the two-dimensional continuum. For sufficiently  can be approximated from a Taylor expansion of Eq. (1) around (x1, x2) (cf. Lasry and Belytschko [24]): smooth fields of s, s

sðx1 ; x2 Þ ¼ sðx1 ; x2 Þ þ c2 r2 sðx1 ; x2 Þ þ c4 r4 sðx1 ; x2 Þ þ    ;

ð4Þ

where r2 is the Laplacian operator. Using the symmetric weight function u in Eq. (2), the odd terms in Eq. (4) vanish and 1 2 ½l =2m ; m ¼ 1; 2; . . .. Differentiating Eq. (4) twice and substituting back into (4), a truncated implicit formulation is c2m ¼ m! obtained, cf. [16], according to:

1 2

sðx1 ; x2 Þ  l2 r2 sðx1 ; x2 Þ ¼ sðx1 ; x2 Þ in X:

ð5Þ

 admits at In the differential Eq. (5), derivatives of the fourth order and higher have been neglected and it is assumed that s least a Fréchet derivative. For the special case when l vanishes, the classical (local) continuum solution is recovered. When the domain X has free boundaries, boundary conditions have to be specified when solving Eq. (5). Often, and also in this  is specified equal study, a Neumann type of boundary condition is utilized on the boundary oX, i.e. the normal derivative of s to zero on oX. The differential Eq. (5), first introduced by Peerlings et al. [23] and often referred to as the implicit gradient nonlocal approach, is the foundation of the work in present study. It is a strong nonlocal theory because in each and every point, interactions are accounted for with every point in the whole problem domain X. 2.2. Implicit stress gradient plasticity On the macroscopic scale, the material is considered isotropic elastic perfectly-plastic with the yield stress r0. It is convenient to extend a traditional local plasticity theory and include nonlocal variables to capture a stress gradient behavior. Assuming an associated J2-flow theory of plasticity, the yield surface / can be expressed by

 e  r0 6 0; /¼r

ð6Þ

 e is a nonlocal equivalent stress given by Eq. (5), i.e. where r

1 2

r e  l2 r2 r e ¼ re ; 3

1=2

ð7Þ

where re ¼ 2 sij sij is the von Mises equivalent stress and the deviatoric Cauchy stress tensor is given by sij = rij  dijrkk/3. Repeated index assumes summation and dij is the Kronecker delta. A nonlocal deviatoric stress tensor is introduced according to:

1243

P. Isaksson / Engineering Fracture Mechanics 77 (2010) 1240–1252

sij ¼ sij

r e : re

ð8Þ

 ij is determined by scaling the local counterpart rij with the scalar valued stress ratio Hence, the nonlocal stress tensor r   r e =re whereupon the relations sij ¼ r ij  dij r kk =3 and r e ¼ 32 sijsij 1=2 holds. The elastic strains eeij are given by Hooke’s generalized law:

rij ¼ Lijkl eekl ;

ð9Þ

where Lijkl = l[(dikdjl + dildjk) + dijdkl2m/(1  2m)] is the macroscopic elastic material stiffness tensor, l = E/(2 + 2m) the shear modulus, E the Young’s modulus and m the Poisson’s ratio. The constitutive relations are written in terms of the incremental strain tensor deij decomposed as follows:

deij ¼ deeij þ depij ;

ð10Þ

where deeij is the linear elastic strain increment and depij is the plastic strain increment. Normality of plastic strain increment requires that

 ij ; depij ¼ k@/=@ r

ð11Þ

 ij ¼ 3sij =ð2r0 Þ and k P 0 is a plastic Lagrange where the outward normal of the associated yield surface is given by @/=@ r multiplier. The mechanical plastic dissipation, given by the well-known Clausius–Duhem inequality, is

r ij depij P 0;

ð12Þ p ij

 ij de P 0 which is similar to the classic plastic dissipation term except that the stress tensor here is nonlocal. The inequality r satisfies that the derived model is admissible from a thermodynamic point of view and is secured by choosing k so that the consistency requirement d/ = 0 together with a loading/unloading condition / = 0 is fulfilled. For an elastic perfectly-plastic material yields

 ij dr  ij ¼ 0: d/ ¼ @/=@ r

ð13Þ

Assuming that the elastic–plastic stiffness remains constant within each stress/strain increment, the strain-driven stress increment yields

drij ¼ Lijkl ðdekl  depkl Þ ¼ Lijkl dekl  kLijmn

@/ :  mn @r

ð14Þ

 ij ¼ drij ðr  e =re Þ and multiply Eq. (14) with @/=@ r  ij , then using relation (13) and solving for k reMake use of the relation dr sults in:

k ¼ sij deij =r0 ;

ð15Þ

 ij Lijkl ¼ 3lskl =r0 have been applied. The final incremental stress–strain relation is written as where the relation @/=@ r

 2e =r40 dekl drij ¼ ½Lijkl  b3lsijskl =r20 dekl ¼ ½Lijkl  b3lsij skl r

ð16Þ

where

 b¼

1 for 0

r e ¼ r0 and dr e P 0

else:

3. The models In order to study structural plasticity effects in a random fiber network and the applicability of the derived implicit stress gradient small deformation plasticity theory to capture this behavior on a larger scale, two different types of models are compared. Both models assume the material is isotropic elastic perfectly-plastic on the macroscopic scale. For simplicity, it is also assumed that the material has a Poisson’s ratio m = 0.3 on the macroscopic scale. The first model is a homogeneous nonlocal continuum model while the second is a micromechanical network model describing each fiber as a structural member and consequently includes components on a considerable smaller scale. The idea is to judge the implicit stress gradient plasticity theory by comparing the mechanical behavior of the model (when utilizing different characteristic lengths l) with the network models (when using different inherent microstructural properties such as fiber length/width/density). The constitutive behavior of the homogeneous nonlocal model follows Section 2 whereas the fibers in the micromechanical network model are described by a classical elastic perfectly-plastic J2-theory with the yield stress r0f (which may be different from the macroscopic yield stress r0). The two different models are used to analyze small deformations in a body containing a macroscopic strain singularity by means of a semi-infinite crack, see Fig. 1.

1244

P. Isaksson / Engineering Fracture Mechanics 77 (2010) 1240–1252

σij x2

θ

r x1

Fig. 1. Body containing a semi-infinite crack.

A Cartesian (x1, x2) and a polar coordinate system (r ¼ ½x21 þ x22 1=2 , h ¼ tan1 ½x2 =x1 ) are introduced with their origins coinciding with the crack tip. The crack occupies the negative part of the x1-axis, i.e. x1 < 0 and x2 = 0. Since the numerical treatment requires a finite body, the body is limited to a disc of radius R (with unit thickness) with its centre located at the crack tip whereupon a boundary layer problem in this manner is formulated. The extent of the body in the out-of-plane direction is small whereby plane stress conditions are assumed to prevail. Hence, the analyzed problem domain X 2 f0 6 r 6 R; jhj < pg. Furthermore, the macroscopic body is loaded so that a pure mode I opening field is acting distant from the crack tip while the shearing mode II field is vanishing. Under the assumption that the Cartesian macroscopic strain tensor eij is linear and infinitesimal, the relation between strains and displacements ui in the macroscopic continuum is given by:

1 2

eij ¼ ðui;j þ uj;i Þ for i ¼ 1; 2 and j ¼ 1; 2:

ð17Þ

At a large distance from the crack tip the macroscopic stress tensor rij is given by:

K

I ffi fij ðhÞ; rij ¼ pffiffiffiffiffiffiffiffi 2pr

ð18Þ

according to LEFM (cf. Williams [25]) where KI is the mode I stress intensity factor and fij(h) are known functions. The boundary conditions, acting at jhj 6 p, are given as prescribed displacements given by Eq. (18):

ui ¼

KI 2l

rffiffiffiffiffiffiffi r g ðhÞ; 8p i

ð19Þ

where gi(h) are known functions. 3.1. Random fiber network The fibers in the discrete networks are highly idealized. This means that each network consists of totally N fibers having length L, width h and unit thickness. In a large network area A, the network’s material coverage c and the average distance Ls, or segment length, between adjacent crossings along a fiber is given by, cf. [2]:

c ¼ NLh=A;

Ls ¼ ph=2c:

ð20Þ

Either the material coverage c or the average segment length Ls specifies completely a two-dimensional RFN when its fiber properties are constant. The Young’s modulus of each fiber (Ef) and their Poisson’s ratio (mf = 0.3) are the same for all fibers in all network models. The fibers are rigidly coupled at each intersection. In all considered networks the fibers length-to-width ratio is in the range L/h 2 [50, 200], which is representative of paper materials. The distance R from the crack tip to the problem domain’s boundary is made relatively large in comparison to the length of the individual fibers and in all analyzed networks L = R/4. 4. Finite element implementation The two models are solved using well-established finite element algorithms for small deformation theory. The numerical schemes are implemented in the Matlab [26] code. A state of plane stress is considered throughout in the computations. The constitutive laws described in Section 2 have been implemented for a two-dimensional solid material element. Four-node isoparametric elements have been used in the homogeneous stress gradient plasticity model while three-node isoparametric elements have been used in the network models. The reason for using three-node elements in the network models is that it is

P. Isaksson / Engineering Fracture Mechanics 77 (2010) 1240–1252

1245

Fig. 2. Illustration of fiber modification to ensure that fiber and bond densities are similar over the network and have a periodic arrangement. The quadratic domain has side-length 2R. The fiber is cut at the boundaries and the deserted part is moved one side-length and placed within the domain. This procedure is repeated until all fiber parts are located within the domain.

a

b

c

Fig. 3. (a) Example of a network mesh consisting of around 106 triangular elements. (b) Close-up of near-tip region. One may observe the crack (here illustrated with broken lines) extending to the left from the centre (tip) of the picture. (c) Mesh constructed of 8000 four-node elements and used in the nonlocal homogeneous model.

significantly easier to design algorithms filling a heterogeneous media with triangular elements than with four-node elements. All elements have two degrees of freedom, translation in the x1- and x2-directions. An iterative technique has been employed to solve the discretized equilibrium equations for each load step. The result obtained – after each iteration – then corresponds to estimates of the incremental displacements from which the currents stress and strain states can be calculated at the integration points of each finite element, as well as estimations of the nonlocal stress variables. Since deformations are assumed to remain small, all derivatives and integrals are evaluated with respect to the initial topology of the considered body. A modified Newton–Raphson iteration scheme for the solution of non-linear finite-element equations has been employed, [27]. The governing finite-element equation is

KA ¼ R  F;

ð21Þ

where K is the global stiffness matrix, A is the global vector consisting of nodal displacement variables, R is the vector of externally applied nodal point loads and F represents the vector of nodal point forces equivalent to the element stresses. The finite element implementations follows mainly the one discussed in [23] for the implicit gradient model. 4.1. Meshes To ensure that the fiber density is evenly distributed over the network, some precautions are introduced when generating the networks. With reference to Fig. 2, the fiber sections that extend outside the quadratic domain with side-length 2R in a mesh generation procedure are moved one side-length in the preferred direction and so located inside the problem domain. When all fibers have been positioned, the fibers are cut at the radius R whereupon a circular domain in this manner is obtained (the origin is located at the centre of the quadratic domain). Cutting the fibers that are crossing the crack-line into two separate parts then produces the macroscopic crack. Finally, the fiber-domains are ‘‘filled” by triangular elements having approximately the same size over the whole network. The crack surfaces in the networks are thus formed by the elementsides located along x1 < 0 and |x2| ? 0. Fig. 3 shows an example of a sparse network mesh and the mesh used for the homogeneous nonlocal model. 5. Analysis 5.1. Asymptotic equivalent stress at tip Eq. (7) is recognized as a modified inhomogeneous Helmholtz equation and it is well known that a solution in a two-dimensional infinite domain can be obtained using a Green’s function, i.e.

1246

P. Isaksson / Engineering Fracture Mechanics 77 (2010) 1240–1252

r e ðx1 ; x2 Þ ¼

Z

Z

1 1

1

1

re ðx01 ; x02 ÞGðx1 ; x2 ; x01 ; x02 Þdx01 dx02 ;

ð22Þ

where GðÞ is the Green’s function, i.e. the solution to Eq. (7) with the source term equaling the Dirac delta function dðÞ, so that G satisfies

1 2 Gðx1 ; x2 ; x01 ; x02 Þ  l r2 Gðx1 ; x2 ; x01 ; x02 Þ ¼ dðx1  x01 Þdðx2  x02 Þ: 2

ð23Þ

Its fundamental solution Gp=ffiffiffiG0 in the two-dimensional full-space can be found in textbooks on mathematical physics and 2 reads G0 ðx1 ; x2 ; x01 ; x02 Þ ¼ B0 ð 2q=lÞ=½pl  where B0 ðÞ is the modified zero-order Bessel function of the second kind. With use of G0 the general solution to Eq. (7) is given by

r e ðx1 ; x2 Þ ¼

Z

1

pl

2

1

Z

1

1

1

pffiffiffi

re ðx01 ; x02 ÞB0 ð 2q=lÞdx01 dx02 :

ð24Þ

Under the present loading situation, the von Mises equivalent stress is given by re ¼ K I ½7 þ 4 cos h  3 cos 2h=ð16prÞ1=2 and  e0 :  e captures the finite value r at the crack tip r

r e0 ¼

Z

KI 4p

3=2 l2

1

Z p

pffiffiffi 0 1=2 ½r 0 ð7 þ 4 cosðh0 Þ  3 cosð2h0 ÞÞ1=2 B0 ð 2r0 =lÞdr dh0 ¼ aK I l

ð25Þ

p

0

pffiffiffi pffiffiffi 02 1=2 where r 0 ¼ ½x02 , h0 ¼ tan1 ½x02 =x01  and a ¼ Cð34 Þ2 21=4 ½6 þ 3 logð2 þ 3Þ=ð6p3=2 Þ is a constant. Eq. (25) differs slightly 1 þ x2  from another published result, Tovo and Livieri [28]. However, this discrepancy is explained by their utilization of a principal stress instead of the von Mises equivalent stress as a local stress. The relation (25) will be used later to estimate the material’s inherent gradient length l.  e =re0 , ~r ¼ r=l, ~ ~e ¼ r ~ ¼ q=l, which when applied to (24) leads to a xi ¼ xi =l and q It is convenient to introduce the scaling: r dimensionless integral formula for nonlocal equivalent stress:

r~ e ð~x1 ; ~x2 Þ ¼

3 pffiffiffi pffiffiffi  25=4 Cð34 Þ2 ½6 þ 3 logð2 þ 3Þ

Z

1

1

Z

1

1

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi pffiffiffi 7 þ 4 cosðh0 Þ  3 cosð2h0 Þ ~ Þd~x01 d~x02 pffiffiffiffi B0 ð 2q ~r0

ð26Þ

5.2. Linking micro/macro structural properties A convenient way to link a network’s fiber Young’s modulus Ef to the material’s macroscopic Young’s modulus E in a homogeneous material description, is to relate the two systems’ potential energies to each other: naturally the energies have to be equal since the two models consider the very same body loaded in equal manner. In the present loading case, the local macroscopic strain energy density w ¼ rij eeij =2 yields

w¼

K 2I cos2 ðh=2Þ½3  m  ð1 þ mÞ cosðhÞ: 4Epr

ð27Þ

The total potential energy U of a circular disc with radius R and thickness t is readily obtained using Eq. (27), i.e.:

U¼t

Z

R

wrdr

Z p

0

dh ¼

p

K 2I Rt ð10  6mÞ: 16E

ð28Þ

The potential energy U# for a network loaded with boundary displacements according to (19) is directly obtained from the finite element analysis, U # ¼ 12 AT KA. Then, setting U# = U, the macroscopic E is directly obtained by

E¼

K 2I Rt ð10  6mÞ: 16U #

ð29Þ

Thus, for a sufficiently low remote load KI, i.e. so low so that plastic flow not has initiated in the network, Eq. (29) relates E to every specific network solution with inherent fiber Young’s modulus Ef. Subsequently, making use of the scaling ratio E/Ef, the macroscopic yield stress r0 is assumed linked to the fiber’s yield stress r0f through the rough estimation:

r0  r0f ½E=Ef :

ð30Þ

The estimation in (30) is motivated by that the network’s microscopic load carrying area is assumed to affect its macroscopic quantities E and r0 in an equal manner.

P. Isaksson / Engineering Fracture Mechanics 77 (2010) 1240–1252

1247

6. Results 6.1. Nonlocal equivalent stress field ~ e , prior to yielding, in the near-tip region. The posiFig. 4 shows the field of the dimensionless nonlocal equivalent stress r ~ e is numer~ e g  1:486 are also marked, which are located at ~r  0:910 and jhj  69:538 degrees. The stress r tions of maxfr ically estimated at discrete locations throughout the problem domain using the integral formula (26) and assuming a very large radius of the problem domain, in comparison to any other length, so that the kernel in (26) vanishes far away from the dominant macroscopic strain singularity located at the origin. Essentially, the routine dblquad is utilized in Matlab [26] to estimate (26) numerically. 6.2. Estimation of inherent gradient length l LEFM is based on the concept of small scale yielding, i.e. in situations when the plastic zone at the crack tip is sufficiently small such that the macroscopic K-dominated field still gives a good approximation of the actual field in an annular region surrounding the tip. As the remote load increases, the asymptotic condition is increasingly violated. Several researchers, among others Hutchinson [29], reports that small scale yielding appears to be a reasonable assumption as long as the applied load is below one half of the load at which full plastic yielding occurs (i.e., the limit load for an elastic perfectly-plastic body). Fig. 5a shows a typical calculated, by the finite element method, macroscopic load–displacement curve of a network subjected to a remote pure KI-field. Here, P is the load that is work conjugate to the applied KI, given by displacements according to (19). The applied KI is normalized with K I , the stress intensity factor KI corresponding to the limit load for the macroscopic small scale yielding assumption, i.e. the load of one half of the load Pc at which full plastic yielding occurs. Fig. 5b shows the calculated, by the finite element method, macroscopic load at corresponding remote KI when using the homogeneous nonlocal plasticity model utilizing a gradient length l = l0 contrasted with the classical homogeneous J2-plasticity model (i.e. l = 0). The classical J2-plasticity model deviates noticeable from the network solution while, for a certain chosen gradient

a

b

~ e in the crack tip region, obtained by evaluating Eq. (26). The locations of maxfr ~ e g  1:486 are also marked (+). (b) r ~ e along crack Fig. 4. (a) Contours of r ~ e is finite as r ! 0 in contrast to its local counterpart re that have a square-root singularity. plane (h = 0). Notice that r

a

b

Fig. 5. (a) Typical network solution. For this particular network, the segment length ratio Ls/L = 0.09 and the length-to-width ratio L/h = 100. (b) Homogeneous nonlocal continuum plasticity model with gradient length l = l0, contrasted with a classical solution (i.e. l = 0) and the network solution in (a). The macroscopic load P is normalized with Pc, the load at which full plastic yielding occurs in the network.

1248

P. Isaksson / Engineering Fracture Mechanics 77 (2010) 1240–1252

length l = l0, the homogeneous nonlocal plasticity model produces almost a perfect agreement with the load curve obtained from the network model. However, how is an optimal value of l0 determined? Let KIp denote at which remote macroscopic load KI onset of plastic  e0 g is estimated to 1.486, Eq. (25) can be directly applied to estimate the mate e =r flow occur in the body. As the ratio maxfr  e0 ¼ 1:486 resulting in: rial’s inherent gradient length l0 by setting r0 =r

 2  2 K Ip K Ip l0 ¼ 1:486a  0:433

r0

ð31Þ

r0

To judge Eq. (31), a numerical investigation is performed. The gradient length l0 is estimated in an iterative procedure for several different discrete networks and is determined as the l, which when used in the gradient plasticity finite element model, result in a macroscopic load P at K I that captures the limit value for macroscopic small scale yielding in the network solution (i.e. P ¼ 12 P c at K I ). The iterative procedure is performed on the nonlocal homogeneous model at each network configuration as follows: a too low l0 will, when used in the homogeneous gradient plasticity model, produce a too low load P at K I compared to the network solution, as is shown in Fig. 5b (dotted line), whereas utilizing a too high l0 produces a higher load P at K I than the network solution. The length l0 is successively evaluated until a perfectly match (i.e., less than 0.1% relative error) is achieved between the two models. Fig. 6a shows these numerically estimated gradient lengths l0 together with the theoretical value given by Eq. (31) for 30 different networks where N and L/h have been varied. The ratio (KIp/r0)2 and length l0 are both normalized by the fiber length L = R/4 used in all the networks. The theoretical value of l0 represents a characteristic gradient length for a nonlocal homogeneous continuum having certain macroscopic properties KIp and r0. Having in mind that each discrete network is randomly generated, Fig. 6a reveals that Eq. (31) reproduces the numerically estimated l0 remarkably well and thus lends confidence to the theory. It should be mentioned that the numerically determined l0 resulted not only in a load P ¼ 12 P c at K I , but also resulted in a macroscopic load P at onset of plastic flow (Fig. 5) that captures the corresponded load in the network solution remarkably well (within 5% in all analyzed cases). In the classical local continuum solution, this yielding load is considerable lower as a consequence of the square-root singularity in stress field (see Fig. 5). Hence, a properly chosen l0 produces agreement throughout the whole loading cycle up to the classical limit value for macroscopic small scale yielding. Fig. 6b shows the computed macroscopic ratio (KIp/r0)2 as function of material coverage c (Eq. (20)). Since, according to Eq. (31), the length l0 scales linearly to the ratio (KIp/r0)2 one realizes, with support from Fig. 6b, that the inherent gradient length controlling the stress gradient plasticity actions depends roughly on the level of coverage in the network rather than the average fiber segment length or equivalent quantity. In the limit c ? 1, i.e. for very dense networks, one may expect that l0 will vanish. A future analysis will be performed to fully examine this hypothesis. 6.3. Estimation of plastic zone size The size of the plastic zone in a homogeneous nonlocal continuum can be approximated by a first order estimation. The  e g is given by ~r  0:910, according to the findings in Section 6.1. To a first order approximation this disdistance to maxfr tance is the size ~r p of the macroscopic plastic zone at onset of plastic flow, i.e. when applied KI = KIp. When the remote load KI further increases, the plastic zone increases in size. Utilizing the nondimensional equivalent stress field, integral formula (26), and the finite tip-stress, Eq. (25), apfirst order estimation of plastic zone size ~rp is given by estimating at which ffi  e0 ¼ r0 l=ðaK I Þ prevails at increasing remote load K I P K Ip . Fig. 7 shows in which way ~r P 0:910 the condition r ~ e ¼ r0 =r ~rp to a first order approximation depends on the applied KI/KIp. As illustrated, the approximated plastic zone size in a nonlocal homogeneous continuum is larger for larger gradient lengths l at onset of plastic flow since ~r p ¼ 0:910, i.e. r p ¼ 0:910l, when remote load KI = KIp. The size of the macroscopic plastic zone increases relatively faster for larger gradient lengths as the remote load factor KI/KIp further increases. Fig. 7 further implies that ~rp with good accuracy can be approximated with a

a

b

Fig. 6. (a) Numerically estimated gradient lengths l0 for different networks. (b) The computed length ratio (KIp/r0)2/L at different coverage c. An in a leastsquare sense fitted line is also drawn to indicate a linear behavior.

P. Isaksson / Engineering Fracture Mechanics 77 (2010) 1240–1252

1249

6 4 2 0

1

1.2

1.4

1.6

1.8

2

Fig. 7. First order estimation of nondimensional plastic zone size ~rp . The broken thin line is a fitted straight line capturing the exact ~r p at KI = KIp and KI = 2KIp, Eq. (32).

straight line for low ratios KI/KIp. A simple linear function is easily obtained by using an ansatz with a polynomial of the form ~r p ¼ 0:910 þ a0 ðK I =K Ip  1Þ and seeking a0 that produces an exact agreement with the numerical solution of ~rp at the points KI = KIp and KI = 2KIp. This result in a0 = 6.811, i.e.

~rp  0:091 þ 6:811½K I =K Ip  1;

for K I P K Ip :

ð32Þ

In the network solutions, the size of the macroscopic plastic zone is truly diffuse as shown in Fig. 8a. At the tip there are some fiber segments that remains unloaded and does consequently not exhibit yielding even when fibers at remote distances from the tip have yielded. It is therefore difficult to evaluate the plastic zone size in network materials. However, it is possible to obtain an estimation of the diffuse macroscopic plastic zone size in a network by defining r # p as the maximum distance to all yielded locations (i.e. integration points) computed at each applied load factor KI/KIp. Fig. 8b shows a typical result for r# p in a network solution together with the numerically estimated ~rp . To obtain comparable quantities, r # p is normalized with l0, given by Eq. (31) for the particular network. The reason for the nearly horizontal sets for r # p , in Fig. 8b, is that the plastic zone in the network does not grow homogeneously under continued load (as in a classical homogeneous material). The plastic flow occurs in the ‘‘weakest part at the moment”, – a phenomenon that is entirely due to the heterogeneous structure. This means that, depending on the macroscopic load; plastic flow sometimes takes place in several regions simultaneously while plastic flow sometime occurs only in a few places. Fig. 8 reveals that the first order estimation of plastic zone size ~r p is a fair approximation of the diffuse plastic zone size in a network. Of course, since the yielded locations in the network are somewhat randomly distributed in the near-tip region, one should have in mind that ~rp is only a very rough approximation. Nevertheless, ~rp may be of interest in order to judge the extension of a non-linear macroscopic process zone in networks at specific applied loads. Also shown in Fig. 8a is the wellknown first order estimation of plastic zone size according to Irwin [30], i.e. r p ¼ ½K I =r0 2 =ð2pÞ. Obviously, the size of the macroscopic diffuse plastic zone in a fiber network is considerable larger as compared to classic homogeneous continuums, which is a result of the long-ranging actions that limit the size of gradients and thereby significantly alters the stress field in the near-tip region. 6.4. Evaluation of the J-integral The well-known J-integral, introduced by Rice [31], plays an important role in non-linear fracture mechanics. Consider an arbitrary counter-clockwise path S enclosing the tip of a crack, as illustrated in Fig. 9.

a

b

Fig. 8. (a) Diffuse macroscopic plastic zone in a network computed at remote load KI/KIp = 2. The yielded locations in the network are shown as black dots. For this particular network, the fiber segment length ratio Ls/L = 0.01 and the length-to-width ratio L/h = 50. The first order estimation of plastic zone size according to Irwin [30] is illustrated as a circle in front of the tip. (b) Comparison between approximations of (normalized) diffuse macroscopic plastic zone # ~ sizes r# p =l0 in a network and the first order estimation r p at increased remote load. Here, r p is defined as the maximum distance to all yielded locations computed at each applied load factor KI/KIp and l0 is given by Eq. (31) for the particular network.

1250

P. Isaksson / Engineering Fracture Mechanics 77 (2010) 1240–1252

2

1

Fig. 9. Arbitrary contour S around the tip of a crack.

Fig. 10. The nonlocal J – integral, according to (34), computed at several contours for some different gradient lengths l. The results are normalized with J0 ¼ K 2I =E, which is the theoretical value for a linear-elastic material according to classical local theory. The value of J is computed at a remote load KI < KIp, i.e. before onset of plastic flow has occurred in the body, and l0 is here equal to 0.15R.

In small scale yielding, the contour S can be chosen to fall within the annular region in which the K-fields holds, i.e. outside the plastic zone if such is present, cf. [29]. The J-integral, on classical form, is given by:

J¼

Z 

wdx2  rij

S

@ui nj dS ; @x1

ð33Þ

where nj are components of the unit vector normal to S. Applying the theory derived in Section 2, a slightly modified J-integral in the present model is given by:

J ¼

Z 

 ij  2r wdx

S

@ui nj dS ; @x1

ð34Þ

where the nonlocal total stress work density (per unit volume) is given by

 ¼ w

Z

eij

0

r ij deeij :

ð35Þ

Thus, it is here assumed that the homogeneous nonlocal continuum can be characterized by a nonlocal strain energy func ij ¼ @ w=@  eeij . The derivation of integral (34) is straightforward and follows exactly the derivation of (33),  eeij Þ such that r tion wð   ij and w. which can be found in any textbook on fracture mechanics, with the exception of utilizing the nonlocal quantities r Rice [31] showed that the value of J-integral is independent of the integration path around the crack in a local classical theory. To illustrate that path independence does not prevail for strong nonlocal material models, J is computed along several contours enclosing the tip utilizing different gradient lengths l when remote load KI < KIp (i.e. the material is still linearlyelastic). The computations are straightforward and follows established routines, cf. [32]. Fig. 10 illustrates that J is path dependent when l > 0 while the classical solution is captured when l = 0. Thus, when l > 0 the result is in contradiction with the classical theory. However, this is not surprising since stress gradients are included when calculating the integral (34) in the body, which completely alter the field around the tip (Fig. 4). One may observe in Fig. 10 that for small magnitudes of the gradient length l, the position of maxfJg is computed along a contour positioned at a radius r approximately equal to l. 7. Conclusions The plasticity behavior of fiber networks is governed by complex mechanisms. This study examines the effect of microstructure on the macroscopic plastic behavior of two-dimensional random fiber networks such as strong-bonded elastic– plastic paper. Remote load is a pure macroscopic mode I opening field, applied via a boundary layer.

P. Isaksson / Engineering Fracture Mechanics 77 (2010) 1240–1252

1251

It is shown that using a classical homogeneous continuum approach to describe plasticity effects due to macroscopic singular-dominated strain fields in fiber networks can lead to erroneous results. The classical continuum description is too simple to capture the essential mechanical behavior of network materials since a structural effect, that alters the macroscopic stress field, becomes pronounced and introduces long-ranging microstructural effects that have to be accounted for. Because of this, it is necessary to include a nonlocal theory that bridges the gap between microscopic and macroscopic scales to describe the material response when using homogeneous continuum models. An implicit stress gradient small deformation plasticity model, which is based on a strong nonlocal continuum formulation, is presented here that has the potential to describe the plasticity behavior of fiber networks on a macroscopic scale. The theory is derived by including nonlocal stress terms in the classical associated J2-theory of plasticity. The nonlocal stress tensor is found by scaling the local Cauchy stress tensor by the ratio between calculated nonlocal and local von Mises equivalent stresses. The model is relatively easy to implement in ordinary finite element algorithms for small deformation theory. Plastic flow in the neighborhood of a macroscopic crack in a homogeneous nonlocal continuum are calculated in finite element models that, for a certain characteristic inherent length, produces plastic zones which sizes agree well with the diffuse macroscopic plastic zones calculated in micromechanical network models. Hence, fairly good agreements are obtained between discrete micromechanical models and the homogeneous nonlocal continuum model. The size of the diffuse macroscopic plastic zone in a fiber network is considerable larger, as compared to classic homogeneous continuums loaded in an equal manner, which is a result of long-ranging actions that limit the size of gradients and thereby significantly alters the macroscopic stress field in the body. A properly chosen characteristic length in the nonlocal plasticity model produces agreement with the micromechanical network model throughout the whole loading cycle up to the limit value for macroscopic small scale yielding. An appropriate characteristic length can be determined by a simple relation in which one need to measure, in experiments, the material’s macroscopic yield stress and the remote load at which plastic flow initiates in a cracked body. In the classical J2-plasticity continuum model, as a consequence of the square-root singularity in stress field, the remote load at which plastic flow initiates in the body is substantially lower than in the nonlocal continuum and network models. Finally, it is demonstrated that path independence of the well-known J-integral does not prevail for this class of material models since long-ranging actions alters the stress field in the cracked body. Only for the special situation of a vanishing inherent length, i.e. in an ideal homogeneous continuum that can be approximated by Saint–Venant’s principle of local action, can the plastic flow be described by the classical theories of plasticity. Acknowledgment The Swedish Research Council is acknowledged for funding this work. References [1] Isaksson P, Hägglund R. Structural effects on deformation and fracture of random fiber networks and consequences on continuum models. Int J Solids Struct 2009;46:2320–9. [2] Niskanen KJ. Paper physics (book 16 of series on papermaking science and technology). Helsinki, Finland: Fapet OY; 1998. [3] Jentzen CA. The effect of stress applied during drying on some of the properties of individual pulp fibers. Tappi 1964;47:412–8. [4] Niskanen K, Kettunen H, Yu Y. Damage width: a measure of the size of fracture process zone. In: 12th Fundamental research symposium, Oxford UK; 2001. [5] Yu Y, Kettunen H, Niskanen KJ. Connection between the damage profile and cohesive stress-widening curve of paper. J Pulp Paper Sci 2002;28:J176–81. [6] Swineheart D, Broek D. Tenacity and fracture toughness of paper and board. J Pulp Paper Sci 1995;21:J389–97. [7] Seth RS. Measurement of fracture resistance of paper. Tappi 1979;62:92–5. [8] Wellmar P, Fellers C, Nilsson F, Delhage L. Crack-tip characterization in paper. J Pulp Paper Sci 1997;23:J269–76. [9] Fleck NA, Hutchinson JW. A phenomenological theory for strain gradient effects in plasticity. J Mech Phys Solids 1993;41:1825–57. [10] Fleck NA, Hutchinson JW. Strain gradient plasticity. Adv Appl Mech 1997;33:295–361. [11] Fleck NA, Muller GM, Ashby MF, Hutchinson JW. Strain gradient plasticity: theory and experiment. Acta Metall Mater 1994;42:475–87. [12] Gao H, Huang Y, Nix WD, Hutchinson JW. Mechanism-based strain gradient plasticity – I. Theory. J Mech Phys Solids 1999;47:1239–63. [13] Huang Y, Gao H, Nix WD, Hutchinson JW. Mechanism-based strain gradient plasticity – II. Analysis. J Mech Phys Solids 2000;48:99–128. [14] Mentzel A, Steinmann P. On the continuum formulation of higher order gradient plasticity for single and polycrystals. J Mech Phys Solids 2000;48:1777–96. [15] Gudmundson P. A unified treatment of strain gradient plasticity. J Mech Phys Solids 2004;52:1379–406. [16] Pijaudier-Cabot G, Bazant ZP. Nonlocal damage theory. J Engng Mech 1987;113:1512–33. [17] Eringen AC. Nonlocal continuum field theories. NY: Springer-Verlag; 2002. [18] Jirásek M, Rolshoven S. Localization properties of strain-softening gradient plasticity models. Part I: strain-gradient theories. Int J Solids Struct 2009;46:2225–38. [19] Jirásek M, Rolshoven S. Localization properties of strain-softening gradient plasticity models. Part II: theories with gradients of internal variables. Int J Solids Struct 2009;46:2239–54. [20] Kröner E. Elasticity theory of materials with long range cohesive forces. Int J Solids Struct 1967;3:731–42. [21] Eringen AC, Edelen DGB. On nonlocal elasticity. Int J Engng Sci 1972;10:233–48. [22] Silling SA. Reformulation of elasticity theory for discontinuities and long-range forces. J Mech Phys Solids 2000;48:175–209. [23] Peerlings RHJ, de Borst R, Brekelmans WAM, de Vree JHP. Gradient enhanced damage for quasi-brittle materials. Int J Num Meth Engng 1996;39:3391–403. [24] Lasry D, Belytschko T. Localization limiters in transient problems. Int J Solids Struct 1988;24:581–97. [25] Williams ML. On the stress distribution at the base of a stationary crack. J Appl Mech 1957;24:109–14. [26] Matlab. Version 7.8. Natick, MA USA: The MathWorks Inc.; 2009.

1252

P. Isaksson / Engineering Fracture Mechanics 77 (2010) 1240–1252

[27] [28] [29] [30] [31] [32]

Bathe KJ. Finite element procedures in engineering analysis. USA: Prentice-Hall; 1982. Tovo R, Livieri P. An implicit gradient application to fatigue of sharp notches and weldments. Engng Fract Mech 2007;74:515–26. Hutchinson JW. Nonlinear fracture mechanics. Dept. solid mechanics, Technical university of Denmark; 1979. Irwin GR. Analysis of stresses and strains near the end of a crack traversing a plate. J Appl Mech 1957;24:361–4. Rice JR. A path independent integral and the approximate analysis of strain concentration by notches and cracks. J Appl Mech 1968;35:379–86. Adina. Theory and modelling guide. Watertown, MA USA: ADINA R&D Inc.; 1995.