Available online at www.sciencedirect.com Chinese

ScienceDirect Sautics

Chinese Journal of Aeronautics 23(2010) 46-53 www.elsevier.com/locate/cja

Mesh-dependence of Material with Softening Behavior

Fu Li*, Du Xingwen

Center for Composite Materials, Harbin Institute of Technology, Harbin 150001, China Received 1 December 2008; accepted 16 February 2009

Abstract

The issue of mesh-dependence emerges when the conventional continuum damage model is applied to handling the softening behavior. In order to circumvent the mesh-dependence, the non-local theory is introduced into the conventional damage model and the finite element formulas are derived for two-dimensional gradient-enhanced damage model. A new element is proposed in which the basic unknown quantities are displacement, non-local equivalent strain and the gradient of non-local equivalent strain. The element and constitutive equation proposed in this article are added to the finite element software ABAQUS through user subroutine UEL. Numerical results show that the gradient-enhanced damage model can eliminate the mesh-dependence and is effective for dealing with the issue of softening behavior.

Keywords: non-local theory; gradient dependence; softening behavior; mesh-dependence; Hermitian element; finite element method

1. Introduction

Softening behavior is a common phenomenon in the field of aerospace structure design. For example, the laminated composites damage due to bird impacts on airplane falls into this case. The development of softening is controlled by the variation and interaction of microstructure. Consequently, the scale of microstructure plays an important role in forming the scale of damage zone[1]. In the conventional continuum damage model, softening behavior can cause local loss of ellipticity of the differential equations when the scale of observation is close to the scale of microstructure. The mathematical description becomes ill-posed and the numerical solutions are dependent on mesh size. T. Drabek, et al.[2] and N. Germain, et al.[3] analyzed composites damage by coupling non-local theory with damage mechanics, which overcame the difficulty of mesh-dependence. The success of this research laid a theoretical foundation for the structural design of astronautical and aeronautical fields. Therefore, it is necessary to conduct in-depth studies.

Aiming at the issue of mesh-dependence, several

Corresponding author. Tel.: +86-451-86402366. E-mail address: bearfl333@163.com

1000-9361/$-see front matter © 2010 Elsevier Ltd. All rights reserved. doi: 10.1016/S1000-9361(09)60186-2

regularization methods have been proposed, such as Cosserat model[4], integral model[5-7] and gradient model[8-10] to remedy the constitutive equation. From the standpoint of computation, the gradient model has been proved to be simpler and more convenient. Therefore, in this article the gradient model is adopted to solve the problem of mesh-dependence for materials with softening behavior. The gradient-dependent plastic model proposed by E. C. Ai-fantis[8] is a simpler gradient model. The Laplacian of non-local plastic strain ep describes the dislocation in the finite and small neighborhood of considered point and the parameter c denotes localization character. The model is widely applied to practice because of its practicality and effectiveness. R. de Borst, et al.[9] coupled the gradient-dependent plastic model proposed by Aifantis with conventional damage constitutive model and formulated gradient-dependent plastic damage model and applied it to analyzing the plastic damage of concrete. P. H. J. Peerlings, et al.[10] developed the gradient-dependent damage model. Aiming at quasi-brittle material, the implicit gradient-dependent damage model was derived further based on integral model and the finite element formulations were derived using weighted residual method.

The analytic solutions of governing equations are difficult to be obtained although the gradient-dependent damage model is more convenient for

computation. Consequently, it is necessary to investigate the numerical methods. J. Petera, et al. [11] proposed Ci continuous Hermitian element and proved that the second-order mixed derivative of the degree of freedom is invariable when it is transformed from local coordinate system to global coordinate system. Due to the implementation of the Laplacian term of plastic strain in the gradient-dependent plastic damage model, J. D. Robert[12] chose Hermitian element to discretize the continuous body and deduced the finite element formulations of governing equations. He compiled the user subroutine UEL and incorporated UEL into ABAQUS. Numerical results showed that the gradient-dependent plastic damage model overcame the difficulty of mesh-dependence.

In this article, two-dimensional gradient damage model is derived for the softening behavior of materials. The finite element formulations of the model are derived based on weighted residual method. A new element M. 2D. 4 is proposed to make the continuous body be discrete and compiled by UEL. The softening problem is analyzed by using the method proposed in this article. Numerical results show that the mesh-dependence is removed. Additionally, the effects of different characteristic lengths on the distribution of strain localization are discussed.

2. Method of Conventional Continuum Damage Mechanics

2.1.Two-dimensional damage constitutive equation

The viscous, thermal or other non-mechanical effects are not taken into account and strain and rotation are assumed to be small. The damage is considered to be isotropic, thus a scalar quantity, the damage variable 0< rn <1, suffices to describe the damage process. The undamaged material is characterized by ^=0, while ^=1 corresponds to the complete loss of material coherence. Utilizing the concept of effective stress and the hypothesis of strain equivalence, leads to the stress-strain relation

a = (1 -d) Ds

where m is damage variable and denotes the area density of micro-defects of material, 8 second-order strain tensor, o second-order Cauchy stress tensor, and D fourth-order elasticity tensor.

2. 2. Definition of equivalent strain

The scalar quantity of equivalent strain £ can substitute for the strain tensor. The strain components multiplied by different weighted coefficients which weight the different functions of these components

during damage evolution.

There are a lot of definitions of equivalent strain[13-14] and the adopted equivalent strain expression in this article is[15]

£ = 72^ -s2f + /3 (2)

where s1 and s2 are strains along x and y directions.

2. 3. Model of damage evolution

The damage state is governed by the scalar history parameter k, which represents the most severe deformation the material experienced. The loading function f is used to judge whether the damage is extended.

f (£, k) =£- k (3)

where f < 0 denotes that the damage is not extended, and f = 0 represents that the damage is extended. The

relationship between f and the rate of k suffices to Kuhn-Tucker loading-unloading condition:

k > 0, f < 0, fk = 0 (4)

The damage occurs if the equivalent strain is larger than the damage threshold. The increment of damage is described by bilinear damage evolution model[16] which is shown in Fig.1. The relationship between damage variable œ and accumulated damage strain k is

œ( k ) =

kc k - ki

k k^ ki

k < k < k

k > kc

where ki is damage threshold, kc the critical value of complete loss of material coherence, and k is equal to £ when damage occurs.

Fig.1 Bilinear damage evolution relationship.

2. 4. Example

The experiment shown in Ref. [17] is simulated and the relationship between stress and strain shown in Fig.1 is adopted to simplify the simulation. The geometry model of the specimen is shown in Fig.2. The length of the specimen is 100 mm and the width is 50 mm. The specimen is subjected to a pure compressive loading with prescribed displacement

AZ=0.05 mm at both ends. The constitutive behavior is described by the elastic modulus E=!8 GPa, Poisson ratio d=0.209, and the initial threshold value of the material ^=0.004 35. The value of ki is set being equal to 0.004 for the middle shadow area in order to trigger the localization.

y o jt 50 mm

100 mm

Fig.2 Geometry model of specimen.

In Fig.3, the equivalent strain distributions have been plotted as computed with 50*100 and 25*50 element meshes by ABAQUS. The results show that the material failure is directly correlated with the element size. The main reason for this phenomenon is that the conventional damage constitutive model neglects the effects of higher-order terms on stress and strain. Actually, there is larger strain gradient in localization zone. It will lead to mesh-dependence of the simulation if the higher-order strains are neglected.

-60 -40 -20 0 20 40 60 x/mm

Fig.3 Distribution of strain of conventional continuum damage model.

3. Method of Non-local Theory

3. 1. Non-local damage model

From the study of softening behavior, it is clear that the stress or strain state of a point is influenced by its neighbourhood, that is, deformations are generally non-local. As suggested by Z. P. Bazant, et al.[5], the value of non-local equivalent strain £ of x for a certain material is a weighted average of the local strain over the surrounding area:

I(x) = ^ *: x) La(*': *)£(*(6)

where a(x': x) is weight function and x' denotes the relative position vector of the infinitesimal area

Based on the method of gradient damage model proposed by P. H. J. Peerlings, et al.[10], the equivalent strain is expanded into Taylor series and substituted into Eq.(6). Then, the relationship between non-local equivalent strain and equivalent strain is derived as

s = £- cV 2£ (7)

where c is gradient parameter, and V 2 denotes Laplacian operator.

The additional boundary conditions are required because the model includes the Laplacian term of non-local equivalent strain. From the mathematical point of view, it is necessary to specify the W for every boundary point of the considered configuration. In order to avoid the singularity of stiffness matrix, the following boundary conditions are adopted[12]

E,r = 0 (8)

I,rm = 0 (9)

where r and m are normal and tangent unit vectors.

3. 2. Numerical solution of gradient damage model

It is difficult to obtain the analytical solution for gradient model because of the existence of the Lapla-cian term of non-local equivalent strain. Moreover, it cannot be directly solved by commercial finite element software. The element and constitutive model proposed in this article are incorporated into ABAQUS by compiling UEL[18]. The ABAQUS is used to solve the nonlinear equations.

(1) Degree of freedom and interpolating function of element

The Laplacian term of non-local equivalent strain is introduced into the gradient damage model. In order to ensure the continuity of non-local equivalent strain gradient, a new element M.2D.4 (see Fig. 4) is proposed, where 'M' indicates a multifield displacement/non-local equivalent strain element, '2D' indicates a two-dimensional element, and the last number indicates the number of nodes. The element basically consists of two separate elements and the basic unknown quantities are displacement, non-local equivalent strain and non-local equivalent strain gradient. The C0 shape functions are used for the interpolation of displacements and coordinates, and the C1 shape functions are used for the interpolation of non-local equivalent strain.

Mesh 50X100 Mesh 25X50

1 • U. V. ¿\.v, C.y. ¿\.v,

(a) Global coordinate system

• 11, V, S^, S^JJ

(b) Local coordinate system

Fig.4 Two-dimensional element M.2D.4.

The degrees of freedom for node n include displacements un and vn, non-local equivalent strain Wn and non-local equivalent strain gradients £nx,

s_ny and ^nxy . The degrees of freedom for

four-node element with respect to the global coordinates are defined as

C = [ x1 y1 x2 y2 x3 y3 X4

a = [u1 v1 u2 v2 u3 v3 u4

e = [£1 £1,x £1,y £1,xy

74lT V4lT

(10) (11)

£4 £4, x £4, y

£4, xy ]

where C is a vector of nodal coordinates, a a vector of nodal displacements, and e a vector of nodal

non-local equivalent strains.

The same shape functions are used to interpolate the coordinates of nodes and displacements. The relationships are given as follows:

where N is the matrix of shape function, and

N 0 n2

0 n 0 N2

0 N3 0 N4

0 N 0 N

£ - [ H H 2 H 3 H4 ]

where h is the shape function matrix of non-local equivalent strain and the sub-matrix Hn is written as

Hn = [hn hn,x hn,y hn,xy] (17)

In order to derive the Hermitian shape functions in Eq.(17), a local coordinate system as shown in Fig.4(b) is used. The vector of reference nodal coordinates can be defined for the £ coordinates and the ц coordinates as follows:

5n = [ -1 1 1 -!]T (18)

rjn = [-1 -1 1 1]T (19)

Two-dimensional regular element M.2D.4 is adopted to discretize the continuous body and the degrees of freedom include displacement, non-local equivalent strain and non-local equivalent strain gradient. In considering an individual element, the derivative of degrees of freedom can be the derivative of the local coordinate system. However, for assuring the continuity between the elements, the derivatives must be the derivative of the global coordinate system.

The generalized non-local equivalent strain can be expressed as

1 = he (20)

The matrix of Hermitian shape function and the set of degree of freedom are respectively defined as

h = [ H1 H2 H3 H4] (21)

e =[ e1 £2 e3 !u]T (22)

where the sub-matrixes and the sub-vectors represent the shape functions and degrees of freedom at node n with respect to the local coordinate system and are

written as

Hn = [ hin h

In = [!n In

en,f7]T

To transform the shape functions to the global coordinate system, we rewrite the local coordinate system degrees of freedom, Eq.(22), in terms of the global coordinate system degrees of freedom, Eq. (16), through the use of the chain rule. For a given node n, we have

e = Ke (25)

where the transformation matrix K is defined as

Additionally, Hermitian shape functions are used to interpolate non-local equivalent strain. The relationship is written as

The shape function matrix with respect to the global coordinate system is derived by Eq.(16), Eq. (20) and Eq.(25) . It can be written as

Hn = Hn K (27)

The components of the shape functions of displacement and Hermitian with respect to the local configuration coordinate system are given as[10]

Nn = rn ,i(#)rn .„W (28)

h- = fn,i(#) fnv(V) (29)

Kv = fn,i(#) gn,,W (31)

hn,?V = gn,4& gn,n(ri) (32)

where, by introducing a variable 5 and taking 5 as either £ or r\, these shape functions are then obtained through the following polynomials

V = 1(1 + sns) (33)

fn,, (5) = (35 - 53) + 2] (34)

(5) = -4[53 + 5n52 - 5 - 5n ] (35)

(2) Matrixes of strain and non-local equivalent strain gradient

The derivation of strain matrix is not discussed. You can refer to Ref.[19] if you want to know it. In this article, we focus on the derivation of the matrix of non-local equivalent strain. For the two-dimensional model, the gradient reduces to the derivatives of the non-local equivalent strain with respect to the x direction and to the y direction, such that

V £ =

where the matrix of non-local equivalent strain p is expressed as

= № P2 P3 P4]

and the components of the matrix of non-local equivalent strain Pn are written as

Using the chain rule, the gradient of non-local equivalent strain with respect to the local coordinate system is written as

V £ =

^ _ d^

dx ~'v dx d£ _ dn dy

The inverse matrix of Jacobian matrix is written as

dx d77

According to Eq.(16), Eq.(39) and Eq.(40), we yield the gradient of non-local equivalent strain dh

V £ = G

The derivative of Hermitian shape functions with respect to the local coordinate system is given as

= [ Pi P2

where the sub-matrix at node n is written as

~ dH„

= [ P n

Pn,Çri ]

By transforming the shape function in Eq.(41) according to Eq.(27). The shape function in Eq.(42) can be defined with respect to the global coordinate system as

dÇ d77

d£, d77

By taking the appropriate derivatives of the Her-mitian shape functions in Eq.(29) to Eq.(32) with respect to the local coordinates £ and r\, the components of the matrix of Hermitian gradient shape functions are written as

Pn,i =

gi&ffi fn,V(V) 8n4(^) fn,vv(n)_

"fn,ii(#) gn,V(V) /n,i(#) gn,,,W.

gi&ffi gn,V(V)

where the gradient shape functions are written in terms of the following functions

fn,ss (s) = Sn (3 - 3s2)]

fn,sss (s) = 4(-6SnS)

gn,ss (S) = -4(3s2 + 2SnS - 1)

gn.sss (S) = 4(6S + 2Sn )

(3) Finite element equation

According to the weighted residual method, we obtain the following discretized equilibrium equation and the relationship between equivalent strain and non- local equivalent strain

Jß BT (1 - )DBdn^a, -\nBT m1_1Ds1_1 HdßSIi = -f BV-,dß+f NTbdS

in ' 1 JS.

TBdn5at HTH + pTcp)d^5£. = Jß H Te,_1dß + Jß (HT H + pT cp)d^5fi_1 (54)

The sub-matrixes in Eq.(53) and Eq.(54) are defined as

Kuu = BT (1 - ai_1) DBd-ß

Kue=-f BTm, 1 De", Hdn

Jo 11 11

TT-uu Ki-1 TT-uS K i-1 fu J ext f u J int, i—1

v-su _ Ki-1 K £ _ .SIi _ _ fU _ _ K-I h _

Kee = H T H +

K EE1i_1 = ¡n (H T H + / i

Keu = i HT Bdn Ja

feXt = { N TbdS

fUu 1 =LB ^-1d^ fU =lH v^

Consequently, the finite element formula of gradient damage model is written as

' (55)

where i denotes load step and b is nodal load.

(4) Example

We adopt the example of Section 2.4 again to verify that the gradient damage model can overcome the difficulty of mesh-dependence. All the data used for this example, except the value of c, are unchanged. To examine the effect of gradient parameter variations on the response, the 3*11 elements simulation has been conducted for c = 0. 05 mm2, c = 0. 14 mm2 and c = 1.00 mm2. The damage evolution Eq.(5) is still taken into account but the equivalent strain s_ is replaced by the non-local equivalent strain £. The element stiffness matrix (AMATRX) and the element residual (RHS) of the governing equations derived in Section 3.2 are incorporated into ABAQUS by the user subroutine UEL. The flowchart of the simulation is shown in Fig.5.

Characteristic length plays an important role in the simulation and will affect the computing results[20-21]. Characteristic length is related to gradient parameter. Therefore, the effects of different gradient parameters on the distribution of strain localization are analyzed first in this article. The distributions of non-local equivalent strain with different gradient parameters (c = 0.05 mm2, c = 0.14 mm2 and c = 1.00 mm2) are plotted in Fig.6.

In order to highlight the distribution of strain in the localization area, only the strains within the range of [-25, 25] are plotted. It is observed that the large value of characteristic length (c =1.00 mm2) results in wider localization zone and the lower peak of non-local equivalent strain. It can be explained by the fact that more material points within the wider localization zone contribute their residual load- carrying capacities to resist the further deformation. As revealed by Fig.6 that within the frame of gradient damage, the strain distribution is controlled by the value of characteristic length during analyzing the softening behavior.

c = — r 2

Fig.5 Computation flowchart of gradient damage model.

•7 4.5

-m- c = 0.05 mm2

—c = 0.14 mm2

c= 1.00 mm2

-15 -5 5 15 25

Fig. 6 Distributions of non-local equivalent strain for different gradient parameters c.

Therefore, the characteristic length should be determined first before simulation. Characteristic length is close to the maximum size of the heterogeneous body within material. It can be obtained by combining experiment with mesomechanics or by inverse method according to the observed localization width[22]. The latter method is adopted to obtain the characteristic length in this article.

A formula[23] is adopted to estimate the localization width:

w = 2nr (56)

where r denotes the characteristic length related to gradient parameter c[24], namely

Y. S. Pan[17] determined the localization width (w=4.687 5 mm) of the specimen analyzed in Section 2.4 by white light digital speckle correlation method. Then, based on Eqs.(56)-(57) , the gradient parameter c is equal to 0.14 mm2 when w is 4.687 5 mm. Therefore, based on the non-local damage model, this gradient parameter is used to analyze the damage of specimen under compression load.

Fig.7 illustrates the widths of localization and the strain distributions along x axis for different element meshes (50*100 and 25*50). Based on the conventional damage model, the distribution of strain obtained by the finite element simulation is shown in Fig.7(a) whereas Fig.7(b) is the distribution of strain obtained by the non-local damage model simulation for the gradient parameter c of 0. 14 mm2. In order to highlight the distribution of strain in localization area, only the strains within the range of [-15, 15] are plotted.

Fig.7 Contrast figure of distribution of strain localization.

The localization widths are w1=2 mm for 50*100 elements and w2= 4 mm for 25*50 elements in Fig.7(a). The error of the maximum equivalent strain is Ero=13.79% for two cases with different element meshes. Numerical results reveal that based on the conventional continuum damage mo-del, a strong mesh-dependence of strain distribution is obtained through the finite element simulation.

In this article, coupling non-local theory with conventional damage model leads to gradient damage model. The distribution of non-local equivalent strain obtained through the gradient damage model is shown in Fig.7(b). The error of the maximum non-local equivalent strain for two cases of element meshes is 0.78%. Numerical results in Fig.7(b) display that by using gradient damage model the localization width is independent of the finite element meshes, which coincides with Ref.[12].

Fig.7 reveals that the gradient damage model eliminates the mesh-dependence in simulating softening behavior. The strain and localization width within the weaken zone are unchanged when different finite element meshes are chosen.

4. Conclusions

(1) The constitutive equation and the element proposed in this article eliminate the mesh-dependence in simulating softening behavior. Stable localization width and strain can be obtained and provide reliable guarantee for structural design.

(2) Based on the gradient damage model, the localization width obtained by simulation is close to that determined by the experiment described in reference, which verifies the correctness of the finite element method proposed in this article. Consequently, the method can be used to analyze other softening behavior of brittle material.

(3) Based on the non-local damage model, the characteristic length will affect the distribution of strain localization. In order to enhance the precision, it is necessary to determine the characteristic length by experiment.

References

[1] Bazant Z P, Planas J. Fracture and size effect in concrete and other quasi-brittle materials. Boca Raton: CRC Press, 1998.

[2] Drabek T, Bohm H J. Damage models for studying ductile matrix failure in composites. Computational Materials Science 2005; 32(3-4): 329-336.

[3] Germain N, Besson J, Feyel F. Simulation of laminate composites degradation using mesoscopic non-local damage model and non-local layered shell element. Modelling Simulation Material Sciences and Engineering 2007; 15(4): 425-434.

[4] Cosserat E, Cosserat F. Théorie des corps déformable. Paris: A Herman et Fils, 1909.

[5] Bazant Z P, Pijaudier-Cabot G. Nonlocal continuum damage, localization instability and convergence. Journal of Applied Mechanics 1988; 55(2): 287-293.

[6] Jirâsek M, Marfia S. Non-local damage model based on displacement averaging. International Journal for Numerical Methods in Engineering 2005; 63(1): 77-102.

[7] Korsunsky A, Nguyen G, Houlsby G, et al. Analysis of essential work of rupture using non-local damage-plasticity modelling. International Journal of Fracture 2005; 135(1-4): 19-26.

[8] Aifantis E C. On the role of gradient in the localization of deformation and fracture. International Jour-

nal of Engineering Science 1992; 30(10): 1279-1299.

[9] de Borst R, Pamin J, Marc G. On coupled gradient-dependent plasticity and damage theories with a view to localization analysis. European Journal of Mechanics: A/Solids 1999; 18(6): 939-962.

[10] Peerlings P H J, de Borst R, Brekelmans W A M, et al. Gradient enhanced damage for quasi-brittle materials. International Journal for Numerical Methods in Engineering 1996; 39(19): 3391-3403.

[11] Petera J, Pittman F T. Isoparameteric Hermite elements. International Journal for Numerical Methods in Engineering 1994; 37(20): 3489-3519.

[12] Robert J D. A nonlocal model for coupled damage-plasticity incorporating gradients of internal state variables at multiscales. PhD thesis, Baton Rouge, LA: Louisiana State University, 2006.

[13] Lemaitre J, Chaboche J L. Mechanics of solid materials. Cambridge: Cambridge University Press, 1990.

[14] de Borst R, Geers M G D, Peerings R H J. Computational fracture mechanics in concrete technology. Southampton, UK: WIT Press, 1999.

[15] Xue S Y. Theory of elasticity and plasticity. 1st ed. Beijing:China Building Material Industry Publishing House, 2006. [in Chinese]

[16] Kachanov L M. Introduction to continuum damage mechanics. Leiden, the Netherlands: Martinus Ni-jhoff Publishers, 1986.

[17] Pan Y S. Study on rockburst initiation and failure propagation. PhD thesis, Beijing: Tsinghua University, 1999. [in Chinese]

[18] Hibbitt K, Sorensen I. ABAQUS manuals version 6.3. Pawtucket, RI, USA: Hibbitt Karlsson and Sorensen Inc. , 2003.

[19] Wang X C, Shao M. Fundamentals and numerical methods for finite element method. 2nd ed. Beijing: Tsinghua University Press, 2003. [in Chinese]

[20] Bazant Z P, Pijaudier-Cabot G. Measurement of the characteristic length of nonlocal continuum. Journal of Engineering Mechanics 1989; 115(4): 755-767.

[21] le Bellego C, Francois D J, Pijaudier-Cabot G, et al. Calibration of nonlocal damage model from size effect tests. European Journal of Mechanics:A/Solids 2003; 22(1) : 33-46.

[22] Wang X P, Meng G T. Nonlocal elastoplastic theory and its application. Chinese Journal of Rock Mechanics and Engineering 2007; 26 (S1): 2964-2967. [in Chinese]

[23] de Borst R, Sluys L J, Muhlhaus H B, et al. Fundamental issues in finite element analyses of localization of deformation. Engineering Computations 1993; 10(2): 99-121.

[24] Tovo R, Livieri P, Benvenuti E. An implicit gradient type of static failure criterion for mixed-mode loading. International Journal of Fracture 2006; 141(3-4): 497-511.