A non-local model with tension and compression damage mechanisms

Eur. J. Mech. A/Solids 20 (2001) 1–22  2001 Éditions scientifiques et médicales Elsevier SAS. All rights reserved S0997-7538(00)01111-6/FLA

A non-local model with tension and compression damage mechanisms Claudia Comi 1 Department of Structural Engineering, Politecnico di Milano, P. Leonardo da Vinci, 32-20133 Milano, Italy (Received 15 November 1999; revised and accepted 14 June 2000) Abstract – A non-local isotropic damage model is proposed which can be used to predict the behaviour of rock-like materials up to failure. Two isotropic damage variables account for the progressive degradation of mechanical properties under stress states of prevailing tension and compression and two internal lengths, one for tension and the other for compression, are introduced as localisation limiters. A linear bifurcation analysis highlights the regularisation properties of the non-local model. An iterative scheme for the numerical solution of the finite-step problem consisting of a linear global predictor, an averaging phase and a non-linear local corrector is presented. Some illustrative examples of tension and splitting tests show the effectiveness of the proposed model.  2001 Éditions scientifiques et médicales Elsevier SAS damage / non-local model / localisation

1. Introduction The correct modelling of mechanical behaviour of quasi-brittle materials such as concrete and rocks subject to multiaxial tension-compression stress states still represents a challenging theme, especially when the prediction of failure is in point. A key dissipative phenomenon for such materials is certainly microcracking which results in non-symmetric, progressive degradation of mechanical properties under tension and compression conditions. Strength and stiffness degradation can be effectively modelled in the framework of continuous damage mechanics: isotropic damage models with a single damage variable (see, e.g., (Mazars, 1984)) or with a tension damage variable and a compression damage variable (see, e.g., (Mazars, 1986; Mazars and Pijaudier-Cabot, 1989)), anisotropic damage models (see, e.g., (Ramtani et al., 1992; Halm and Dragon, 1996)), damage coupled to plasticity models (see, e.g., (Ortiz, 1985; Simo and Ju, 1987; Lee and Fenves, 1998)) have been proposed, to cite a few. However, as widely studied (see, e.g., (Maier and Hueckel, 1979; Benallal et al., 1988; Desoyer and Cormery, 1994; Comi et al., 1995)) the boundary value problem in the presence of damage-induced softening may become ill-posed and severe difficulties arise in numerical analyses. Various regularisation techniques have been proposed in the literature, leading to enhanced formulations; namely, micropolar (Steinmann and Stein, 1994), non-local integral (Pijaudier-Cabot and Bažant, 1987; Pijaudier-Cabot, 1995), rate dependent (Dubé et al., 1994) and gradient-dependent damage formulations (Frémond and Nedjar, 1996; Peerlings et al., 1996; Comi, 1999). In all these models the explicit or implicit introduction of a material characteristic length fixes the width of the zone in which damage localises, thus preventing strain localisation into a line with consequent zero energy dissipation. The experimental measure of this length is far from being trivial. Under uniaxial tension conditions concretelike materials exhibit a localised failure which is governed by the fracture energy. Comparison of experimental 1 E-mail: [email protected]

2

C. Comi

results concerning tension tests on concrete specimens in which localised failure occurs and tension tests in which a particular experimental apparatus forces the damage to diffuse in the whole specimen allows the definition of a material internal length related to the aggregates size (Bažant and Pijaudier-Cabot, 1989). More recently, experimental studies of Jansen and Shah (1997) concerning simple compression of concrete showed that also the compression failure is localised and that a compression fracture energy and a compression internal length can be measured. This latter length turns out to be very different from the tension internal length. To simulate this behaviour which is non-symmetric both for the homogeneous and localised response, in this paper a two-surface non-local damage model is proposed. In this model the variation in load carrying capacity (hardening followed by softening) and the elastic degradation are described by two scalar damage variables Dt and Dc , ranging from zero to one. The first one measures damage for stress states of prevailing tension and the second one for stress states of prevailing compression. The current elastic domain in stress space is limited by two surfaces Ft and Fc , evolving with damage and non-locality is introduced by means of the weighted averages of the strain invariants which affect the loading functions. The presence in this model of two loading surfaces, besides a correct definition of the elastic domain and of its evolution, gives the possibility of introducing two different internal material lengths, related to the width of the localised band for stress states of prevailing compression or tension. This model represents the non-local enhancement of the fracture-energy based, local model originally proposed in (Comi and Perego, 2000). The paper is organised as follows. In section 2 the underlying local model is briefly presented and the non-local model is formulated. In section 3 a bifurcation analysis for the rate problem is performed and the regularising role of non-locality is highlighted. The numerical aspects related to the use of this model in finite element analyses are discussed in section 4. To solve the time-step problem, a predictor–averaging–corrector scheme is proposed in which the correction can be performed at the local Gauss-point level and only a procedure to evaluate the non-local value of the predicted strain invariants is required in addition to the customary procedure for non-linear local models. In section 5 numerical results obtained with the proposed non-local model are compared with experimental results concerning tension and splitting tests on concrete specimens. As the mesh is refined, convergence of numerical results is observed both in terms of the damage pattern and of the global response. 2. Two-surfaces isotropic damage model 2.1. Local model The local model considered and enhanced in what follows is the bi-dissipative isotropic damage model proposed in (Comi and Perego, 2000). The model is characterised by two scalar damage variables Dt and Dc , ranging from zero to one, which measure damage under stress states of prevailing tension and of prevailing compression. These variables enter in the definition of the free energy ψ: ψ=

1 2µ0 (1 − Dt )(1 − Dc )e : e + K0 (1 − Dt )(tr+ ε)2 + K0 (1 − Dc )(tr− ε)2 , 2

(1)

where ε is the small-strain tensor, e is the strain deviator tensor, µ0 and K0 are the initial, undamaged, shear and bulk elastic moduli, respectively. The first addend represents the deviatoric strain energy which is affected by both damage variables. The volumetric free energy is split into two parts (second and third addends), one ε| depending on the positive part of the volumetric strain tr+ ε (tr+ ε = tr ε+|tr ) and affected by the damage in 2

A non-local damage model

3

ε| tension and the other depending on the negative part of the volumetric strain tr− ε (tr− ε = tr ε−|tr ) and affected 2 by the damage in compression. The assumed expression of the free energy allows to partially model, while preserving isotropy, the unilateral effect due to fracture closure under compression experimentally observed (see, e.g., (Ramtani et al., 1992)). A discussion on the anisotropy induced by a different description of the unilateral effect can be found in (Fichant et al., 1999). The stress–strain relation is obtained from (1) as:

σ=

∂ψ = 2µ e + K+ tr+ ε I + K− tr− ε I ∂ε

(2)

σ being the stress tensor, I the second-order identity tensor, µ ≡ µ0 (1 − Dt )(1 − Dc ) the current damaged value of the shear modulus, K+ ≡ K0 (1 − Dt ) the current damaged value of the bulk modulus when tr ε > 0, and K− ≡ K0 (1 − Dc ) the current damaged value of the bulk modulus when tr ε < 0. The activation and evolution of damage are governed by the following loading–unloading conditions: ft 6 0,

D˙ t > 0,

ft D˙ t = 0,

fc 6 0,

D˙ c > 0,

fc D˙ c = 0

(3)

with ft = J2 − at I12 + bt rt (Dt )I1 − kt rt2 (Dt )(1 − αDc ),

ri (Di ) =

fc = J2 + ac I12 + bc rc (Dc )I1 − kc rc2 (Dc ),
 1 − σσ0e i    (D0i − Di )2 for Di < D0i  2 1 − D0i "   #0.75  Di − D0i ci    for Di > D0i  1− 1 − D0i

(4) (5)

i = t or c.

(6)

The loading functions ft and fc are defined by (4)–(5) in the stress space: J2 = 12 s : s is the second invariant of the stress deviator s, I1 = tr σ is the first invariant of the stress tensor, at , bt , kt , ac , bc , kc and α are non-negative material parameters, rt (Dt ) and rc (Dc ) are the hardening–softening functions defined in equation (6). In these latter, D0t and D0c are the values of damage at the peak σ0t and σ0c of the stress–strain curve in tension and compression, respectively, σet and σec are the elastic limit values, ct > 1 and cc > 1 are parameters defining the slope of the softening branch. The term (1 − αDc ) in ft allows the modelling of the decrease of the load carrying capacity in tension due to compression damage. The tension loading function ft is a hyperboloid around the hydrostatic axis in the Haigh–Westergaard space of principal stresses, with circular deviatoric sections and hyperbolic meridian sections. The compression loading function fc is an ellipsoid (see figure 1). The elastic domain delimited by these functions initially expands due to damage hardening and then shrinks due to damage softening. The boundary of the failure domain, understood as the set of stress states corresponding to the peak of stress–strain curves recorded in tests with radial stress paths, can be easily computed from equations (4), (5) by substituting, in ft and fc , rt (Dt ) = 1 and rc (Dc ) = 1, i.e. the values of the hardening functions corresponding to the peak values of damage D0t and D0c . Under plane stress conditions the resulting rupture domain is represented in figure 2(a) together with the experimental points concerning biaxial tests on concrete specimens taken from (Kupfer et al., 1969). Figure 2(b)

4

C. Comi

Figure 1. Loading surfaces and resulting elastic domain.

Figure 2. Biaxial tests. (a) Failure domain: experimental points and numerical curves; (b) Stress–strain curves under biaxial compression; bold lines: numerical results, light lines: experimental results.

shows a comparison between the experimental stress–strain curves (thin lines) and the curves obtained with the present model (thick lines) for compression tests with different confinements. Using equation (2), the loading functions can be expressed equivalently in the strain space as: ft = ft (ε, Dt , Dc ) = 4µ2 Jε − 9at (K+ tr+ ε + K− tr− ε)2 + 3bt rt (Dt )(K+ tr+ ε + K− tr− ε) − kt rt2 (Dt )(1 − αDc ),

(7)

fc = fc (ε, Dt , Dc ) = 4µ2 Jε + 9ac (K+ tr+ ε + K− tr− ε)2 + 3bc rc (Dc )(K+ tr+ ε + K− tr− ε) − kc rc2 (Dc ),

(8)

A non-local damage model

5

where Jε = 12 e : e is the second invariant of the strain deviator e. These latter expressions of the loading functions will be used in the rest of the paper. In this model the limit values of damage Dt = 1 and Dc = 1 can only be asymptotically reached for infinite strains, however, the specific fracture energy, i.e. the area below the stress–strain curve in one-dimensional conditions, is finite owing to the chosen exponents of the softening function. Remark 2.1: The values D0t , σ0t and D0c , σ0c can in principle be directly measured in homogeneous uniaxial tension and compression tests, respectively. In practice, unless special experimental apparatus are used (Bažant and Pijaudier-Cabot, 1989), the specimen response is not homogeneous due to strain localization. However, the pre-peak experimental behaviour is not significantly affected by strain localisation and therefore the direct identification of D0t , D0c , σ0t and σ0c do make sense. Also the results of localisation analysis for the proposed model (Comi and Rizzi, 2000) under uniaxial conditions confirm that the inception of non-homogeneous responses is around the peak or in the softening regime. Since non-locality affects the overall behaviour only after localisation the meaning of the above parameters also holds true for the non-local model of the next section. 2.2. Non-local model The proposed enhanced version of the above model is based on the concept of non-local damage (PijaudierCabot and Bažant, 1987). As pointed out in (Bažant and Pijaudier-Cabot, 1988), non-locality in damage models can be introduced both through the definition of a weighted average of the damage variable itself or of a strain measure. This second possibility is here exploited in view of a simpler numerical treatment. The basic non-local variables at a point x are here assumed to be the average strain invariants, i.e. the weighted averages over a volume V of the local strain invariants: hJε (x)i = htr+ ε(x)i =

Z V

Z V

W (x − s)Jε (s) ds,

W (x − s) tr+ ε(s) ds,

htr− ε(x)i =

(9)

Z V

W (x − s) tr− ε(s) ds,

(10)

where W (x − s) is the weighting function, adequately defined to normalise the averaging. In equations (9), (10) and in what follows the symbol h•i denotes the weighted average value of the quantity •. In this work W (x − s) is assumed as the normalised Gauss function and the average is extended to the whole body so that V coincides with the body volume: 

W (x − s) =

1 kx − sk2 exp − W0 (x) 2l 2



with W0 (x) =

Z V



exp −

kx − sk2 2l 2



ds.

(11)

The length l is a material parameter which can be related to the width of the zone in which damage phenomena localise. No particular provisions need to be introduced for points near the boundary of the body since W0 (x) in equation (11) already normalises the averaging. For concrete-like materials the localised behaviour is very different in tension and compression (see, e.g., (Bažant and Pijaudier-Cabot, 1989)) for experimental results in tension and (Jansen and Shah, 1997) for experimental results in compression) and it seems reasonable to introduce two different internal lengths lt and lc , one for prevailing tension states and the other for prevailing compression states, to take into account such non-symmetry. Therefore, two different weighting functions are used in (9)–(10) to obtain the mean values of

6

C. Comi

strain invariants to be introduced in the loading function in tension and compression, respectively: exp − kx−sk 2l 2

2

Wt (x − s) = R V

t

exp − kx−sk 2l 2

2





exp − kx−sk 2l 2

2

, ds

Wc (x − s) = R V

t

c

exp − kx−sk 2l 2

2

c





. ds

(12)

The non-local model is then obtained by replacing in the loading functions the local strain invariants with their averages; in the loading function ft the averages should be performed with the weight Wt , while in loading function fc the averages should be performed with the weight Wc . Denoting by h•it and h•ic the averages of the quantity • performed with weights Wt and Wc , respectively, the resulting non-local loading functions Ft and Fc are: Ft = 4µ2 hJε it − 9at K+ htr+ εit + K− htr− εit


2




+ 3bt rt (Dt ) K+ htr+ εit + K− htr− εit − kt rt2 (Dt )(1 − αDc ), Fc = 4µ2 hJε ic + 9ac K+ htr+ εic + K− htr− εic


2

(13)


+ 3bc rc (Dc ) K+ htr+ εic + K− htr− εic − kc rc2 (Dc ). (14)

The non-local model is completed by the loading–unloading conditions (3) in which ft and fc are replaced by Ft and Fc and by the local stress–strain relation (2). Remarks: 2.2 – In equations (13)–(14) the damage variables are locally defined, but the damage at a point is influenced by the damage of the neighbouring points. In fact if damage tends to grow at a single point, say P , the strains increase at that point; as a result the average values of strain invariants which intervene in the loading functions of points in a proper surrounding of P increase, thus inducing damage also in the neighbouring points of P . Even though strains of all points of the body intervene in the average values defined in equations (9)–(10) the shape of the weight functions allows to practically drop contributions outside a proper region, related to the characteristic lengths lt or lc . 2.3 – Only the inelastic behaviour, involving damage activation, is treated as non-local, while the elastic relations remain local, thus satisfying the requirement already pointed out in (Pijaudier-Cabot and Bažant, 1987). According to this basic assumption only the loading functions are made non-local: the activation of tensile damage is modified by the existence of a tensile internal length and the activation of compression damage is modified by the existence of a compression internal length. For this reason all averages are made using lt in Ft and using lc in Fc . The separation between states with positive volumetric strains and negative volumetric strains concerns the linear elastic behaviour during unloading; this behaviour is local and therefore the internal length does not need to change across the interface tr ε = 0. 2.4 – With the assumed non-local definitions of Ft and Fc , equations (13)–(14), the damage variables, locally defined, remain always limited between zero and one as happens for the local model. On the contrary, the introduction in the loading functions of a non-local definition of the damage variables would require ad hoc provisions to check that locally the damage does not exceed the limit value one. 2.5 – A rigorous thermodynamic formulation of the proposed model can be obtained following the same lines proposed in (Borino et al., 1999) for non-local plasticity models.

A non-local damage model

7

3. Bifurcation analysis 3.1. Rate problem Assuming that at each point of the solid full loading occurs with activation of both modes, i.e. considering the linear comparison solid (Hill, 1959), the rate constitutive equations can be obtained from equation (2) and the consistency conditions. The stress–strain relation in rate form reads: e : ε˙ − AD ˙ t − BD˙ c σ˙ = E

with

(15)

  2µ e E = 2µ I ⊗ I + K+ H(tr ε) + K− H(−tr ε) − I ⊗ I,

(16)

A = 2µ0 (1 − Dc ) e + K0 tr+ ε I,

(17)

3

−

B = 2µ0 (1 − Dt ) e + K0 tr ε I.

(18)

e is the damaged isotropic elastic tensor, I ⊗ I is the fourth order symmetric identity In the above relations E tensor such that its component ijhk is 12 (δih δj k + δik δj h ), ⊗ denotes the tensor product, H( ) is the Heavyside function (H(x) = 1 for x > 0 and H(x) = 0 for x < 0). Relation (16) follows from (2) noting that:

d + (tr ε) = H(tr ε) I : ε˙ , dt 

d − (tr ε) = H(−tr ε) I : ε˙ , dt

(19)



1 e˙ = I ⊗ I − I ⊗ I : ε˙ . 3 The damage increments are then evaluated from the consistency conditions. Consider first the local model. These conditions are: ∂ft ˙ ∂ft ˙ ∂ft f˙t = Dt + Dc + : ε˙ = 0, ∂Dt ∂Dc ∂ε ∂fc ˙ ∂fc ˙ ∂fc f˙c = Dt + Dc + : ε˙ = 0, ∂Dt ∂Dc ∂ε where the partial derivatives are computed from equations (7) and (8): 
 ∂ft = −8µµ0 (1 − Dc )Jε + 3K0 tr+ ε 6at K+ tr+ ε + K− tr− ε − bt rt (Dt ) ∂Dt
drt drt + 3bt K+ tr+ ε + K− tr− ε − 2kt rt (Dt ) (1 − αDc ), dDt dDt 
 ∂ft = −8µµ0 (1 − Dt )Jε + 3K0 tr− ε 6at K+ tr+ ε + K− tr− ε − bt rt (Dt ) + αkt rt2 (Dt ), ∂Dc 
∂fc = −8µµ0 (1 − Dc )Jε − K0 tr+ ε 2ac K+ tr+ ε + K− tr− ε + bc rc (Dc ), ∂Dt 
 ∂fc = −8µµ0 (1 − Dt )Jε − 3K0 tr− ε 6ac K+ tr+ ε + K− tr− ε + bc rc (Dc ) ∂Dc

(20)

(21) (22)

(23)

(24) (25)

8

C. Comi + 3bc


drc drc K+ tr+ ε + K− tr− ε − 2kc rc (Dc ) , dDc dDc

(26)


  ∂ft (27) = 4µ2 e + 3 −6at K+ tr+ ε + K− tr− ε + bt rt (Dt ) K+ H(tr ε) + K− H(−tr ε) I, ∂ε 
  ∂fc (28) = 4µ2 e + 3 6ac K+ tr+ ε + K− tr− ε + bc rc (Dc ) K+ H(tr ε) + K− H(−tr ε) I. ∂ε Solving the system (21)–(22) and substituting the result into (15) one obtains the rate constitutive equations for the local model: e− Et = E

σ˙ = Et : ε˙ ;

A⊗

∂ft ∂fc ∂Dc ∂ε

h

−

∂fc ∂ft ∂Dc ∂ε




−

B⊗

∂fc ∂ft ∂Dt ∂ε

∂ft − ∂D t

∂fc
∂ε

h

(29)

with: ∂ft ∂fc ∂ft ∂fc − . (30) ∂Dt ∂Dc ∂Dc ∂Dt Since two modes have been considered, the tangent tensor Et is a two rank one update of the damaged elastic tensor. If only one loading function is active the second or the third term in equation (29b) vanishes and the classical form of the tangent operator of single-surface models is recovered. Consider now the non-local model. The damage increments are again evaluated from the consistency conditions but using in this case the non-local loading functions (13)–(14). These conditions are expressed by the following integral–differential equations: h=

∂Ft ˙ ∂Ft ˙ F˙t = Dt + Dc + Gt (ε˙ ) = 0, ∂Dt ∂Dc ∂Fc ˙ ∂Fc ˙ F˙c = Dt + Dc + Gc (ε˙ ) = 0, ∂Dt ∂Dc

(31) (32)

where: 
 ∂Ft = −8µµ0 (1 − Dc )hJε it + 3K0 htr+ εit 6at K+ htr+ εit + K− htr− εit − bt rt (Dt ) ∂Dt
drt drt + 3bt K+ htr+ εit + K− htr− εit − 2kt rt (Dt ) (1 − αDc ), dDt dDt

(33)


 ∂Ft = −8µµ0 (1 − Dt )hJε it + 3K0 htr− εit 6at K+ htr+ εit + K− htr− εit − bt rt (Dt ) + αkt rt2 (Dt ), (34) ∂Dc 
 ∂Fc = −8µµ0 (1 − Dc )hJε ic − 3K0 htr+ εic 6ac K+ htr+ εic + K− htr− εic + bc rc (Dc ) , (35) ∂Dt 
 ∂Fc = −8µµ0 (1 − Dt )hJε ic − 3K0 htr− εic 6ac K+ htr+ εic + K− htr− εic + bc rc (Dc ) ∂Dc
drc drc + 3bc K+ htr+ εic + K− htr− εic − 2kc rc (Dc ) , dDc dDc

Gt (ε˙ ) = 4µ

Z

2 V






Wt e : ε˙ dV − 3 6at K+ htr+ εit + K− htr− εit − bt rt (Dt )



(36)

A non-local damage model 

× K+ Z

Gc (ε˙ ) = 4µ2 

V

Z V

Wt H(tr ε) I : ε˙ dV + K−



Z



9

V

Wt H(−tr ε) I : ε˙ dV ,


(37)

Wc e : ε˙ dV + 3 6ac K+ htr+ εic + K− htr− εic + bc rc (Dc )

× K+

Z V

Wc H(tr ε) I : ε˙ dV + K−

Z V





Wc H(−tr ε) I : ε˙ dV .

(38)

Substituting the damage increments into equation (15) leads to the rate constitutive equation: e : ε˙ − σ˙ = E

 





∂Ft ∂Fc 1 ∂Fc ∂Ft A Gc (ε˙ ) − Gt (ε˙ ) + B Gt (ε˙ ) − Gc (ε˙ ) H ∂Dc ∂Dc ∂Dt ∂Dt



(39)

with: ∂Ft ∂Fc ∂Ft ∂Fc − . (40) ∂Dt ∂Dc ∂Dc ∂Dt Note that, due to the non-local nature of the model, a local tangent operator Et such that σ˙ = Et : ε˙ can no longer be defined. The rate problem describing the quasi-static evolution of an elasto-damageable solid is given by the rate constitutive equations (equation (29) for the local model and equation (39) for the non-local one), the field equations of equilibrium and compatibility: H=

div σ˙ = 0, (41)
1 ε˙ = grad u˙ + gradT u˙ , (42) 2 where u is the displacement vector and small strains and zero rates of body forces have been assumed, and boundary conditions on displacements and traction rates. 3.2. Bifurcation conditions Starting from a homogeneous equilibrium state, possible bifurcated solutions of the rate problem with the non-local model of section 2.2 are sought in the form: u˙ = v0 e−iq n·x ,

(43)

where q is the wave number and n denotes the normal direction. The corresponding strain rate is: 1 ε˙ = − iq(v0 ⊗ n + n ⊗ v0 ) e−iq n·x . 2

(44)

This perturbation is admissible if it satisfies the rate equilibrium equations, with σ˙ given by equation (39). As shown in the Appendix, this admissibility condition can be expressed in the form: (n · L(q) · n) · v0 = 0,

(45)

where L(q) is a kind of non-local tangent operator defined as: 







∂ft ∂fc ∂fc ∂ft ∂fc ∂ft ∂ft ∂fc A ⊗ W c (q) + B ⊗ W t (q) − W t (q) − W c (q) H ∂Dc ∂ε ∂Dc ∂ ε ∂Dt ∂ε ∂Dt ∂ε

e− 1 L(q) = E



. (46)

10

C. Comi

The functions W t (q) and W c (q) for an infinite body are the Fourier transforms of the weight functions Wt and Wc , respectively: W t (q) = W c (q) =

Z V

Z

V

Wt (ξ ) e−iq n·ξ dξ = e− Wc (ξ ) e−iq n·ξ dξ = e−

lt2 q 2 2

,

(47)

lc2 q 2 2

.

(48)

Thus, a perturbed solution of the form (43) exists if, and only if, the coefficient matrix of system (45) is singular, i.e.: det(n · L(q) · n) = 0.

(49)

If one considers the local model, standard arguments lead to the following critical condition: det(n · Et · n) = 0.

(50)

The tensor n · L(q) · n plays, for the non-local model here considered, a role similar to the one of the acoustic tensor n · Et · n for the local model, but the wave number q enters in the definition (46) through the Fourier transforms W t (q) and W c (q) of the weighting functions. Therefore unlike the local medium, in the non-local one the wavelength is not arbitrary at a bifurcation point. Localisation into a band of zero width, which occurs in the local medium as soon as condition (50) is fulfilled (Rice, 1976; Borré and Maier, 1989), cannot occur in the non-local medium. In fact for zero wavelength (q → ∞) the Fourier transforms W t (q) and W c (q) of e which is positive definite the weighting functions go to zero and L(q) tends to the elastic damaged tensor E for 0 6 Dt < 1 and 0 6 Dc < 1. Therefore the ill-posedness of the boundary value problem arising when using local damage models is avoided by the proposed non-local enhancement. The tensor n · L(q) · n reduces to the local acoustic tensor when W t (q) = W c (q) = 1, i.e. when the weights are Dirac distributions. To compare the bifurcation conditions of the local and non-local model assume that the internal length under compression lc be equal to the internal length under tension lt , then W t (q) = W c (q) = W (q) for any q and one has: 









W (q) ∂ft ∂fc ∂fc ∂ft ∂fc ∂ft ∂ft ∂fc A⊗ +B⊗ . (51) − − H ∂Dc ∂ε ∂Dc ∂ε ∂Dt ∂ε ∂Dt ∂ε Compare this tensor with the local tangent tensor of equation (29b): the two tensors are identical except for the presence of the term H /W (q) instead of h. Therefore, if perturbated solutions exist, for the local and the non-local models the direction n is the same. Denoting by Hcr the maximum value of H , at varying n, which fulfils condition (49) and by hcr the maximum value of h, at varying n, which fulfils condition (50), being W (q) 6 1 for any q, one has: e− L(q) = E

hcr > Hcr .

(52)

Note that, due to the presence of two loading surfaces, h and H are not directly related to a single scalar hardening–softening parameter monotonically decreasing during loading. Therefore from equation (52) one cannot conclude that the criterion of bifurcation for the local continuum gives a lower bound for the corresponding non-local continuum in general. However, this is the case when bifurcation is sought starting from a state in which only one mode is active and the results obtained in (Benallal and Pijaudier-Cabot, 1993) are recovered.

A non-local damage model

11

4. Algorithmic aspects 4.1. Finite-step problem Consider the evolution of a damageable body of volume V and boundary S subject to body forces F(t), surface forces f(t) on the free part of the boundary Sf , endowed with an outward normal m, and imposed displacements U(t) on the constrained part of the boundary Su and let the material behaviour be described by the non-local damage model of section 2.2. Numerical solution of this evolutive problem requires time and space discretizations. Time dicretization is performed by dividing the time interval of interest into finite steps. At the beginning, say tn , of a time step 1t all quantities are known and the solution must be computed at the end of the step, say tn+1 , for given increments of external actions 1F, 1f and 1U. The finite step problem is then obtained by enforcing equilibrium and compatibility at the end of the step and integrating the constitutive law according to a backward difference scheme (subscript n + 1 marks quantities computed at tn+1 and 1 indicates the increment of a variable over the step): div σ n+1 + Fn+1 = 0 in V ; m σ n+1 = fn+1 on Sf ,
1 grad un+1 + gradT un+1 = εn+1 in V ; un+1 = Un+1 on Su , 2 σ n+1 = 2µ0 (1 − Dtn+1 )(1 − Dcn+1 ) en+1 + K0 (1 − Dtn+1 ) tr+ εn+1 I + K0 (1 − Dcn+1 ) tr− ε n+1 I

(53) (54) in V , (55)

(Ft )n+1 = Ft (ε n+1 , Dtn+1 , Dcn+1 ) 6 0;

1Dt > 0;

(Ft )n+1 1Dt = 0 in V ,

(56)

(Fc )n+1 = Fc (εn+1 , Dtn+1 , Dcn+1 ) 6 0;

1Dc > 0;

(Fc )n+1 1Dc = 0

(57)

in V .

Space discretization is performed following a standard finite element displacement approach in which only ¯ and strains are computed through displacements are modelled in terms of their nodal values u¯ (u(x) = Nu (x)u) compatibility. Note that even though the model is non-local, no additional variables should be modelled nor additional continuity requirements arise as it happens in gradient enhanced damage models (see, e.g., (Comi, 1999)). 4.2. Iterative procedure The finite-step problem concerning the time interval 1t = tn+1 − tn is non-linear and can be solved by an iterative scheme. The scheme here adopted is similar to the linear predictor–non-linear corrector scheme used for local models except that an intermediate averaging phase should be added between predictor and corrector. In what follows the index n + 1 of the time step will be omitted and the superscripts i and i + 1 will refer to the iteration. • Predictor phase – At the iteration i + 1 for fixed damage variables values, compute the estimate of the nodal displacement vector u¯ i+1 from the linearized equations of the finite element model. Denoting by Ri+1 the vector of residual forces and by S i+1 the prediction stiffness matrix, the equations to be solved read:


Si+1 u¯ i+1 − u¯ i = Ri+1 .

(58)

The prediction matrix can be chosen in different ways. If in this phase the damage variables are set equal to the values reached at the end of the previous step, Dt = Dtn , Dc = Dcn , the procedure outlined above corresponds to a secant prediction matrix.

12

C. Comi • Averaging phase – At each Gauss point xg , compute the weighted mean values of the strain invariants by adding contributions of strains at all the NG Gauss points of the mesh: hJε (xg )i =

NG X

ωi W (xg − xi )Jε (xi ),

(59)

i=1 +

htr ε(xg )i =

NG X

+

ωi W (xg − xi ) tr ε(xi ),

−

htr ε(xg )i =

i=1

NG X

ωi W (xg − xi ) tr− ε(xi ).

(60)

i=1

In equations (59)–(60) ωi is the Gauss weight of the ith Gauss point and W is the weighting function of the non-local model which is equal to Wt (equation (12a)) when computing average quantities to be introduced in the tension loading function and equal to Wc (equation (12b)) when computing average quantities to be introduced in the compression loading function. Note that the terms ωi Wt (xg − xi ) and ωi Wc (xg − xi ) can be computed and stored once for all at the beginning of the analysis. • Corrector phase – At the Gauss point xg , with the current (of iteration i + 1), predicted values of the mean strain invariants and the values of damage of the previous time step, compute the trial loading functions i+1 i+1 Ft trial and Fctrial : Fttrial

i+1




= Ft ε i+1 , Dtn , Dcn ,

Fctrial

i+1

= Fc εi+1 , Dtn , Dcn




(61)

and calculate damage using the loading–unloading conditions (56)–(57). i+1 i+1 (i) if Fttrial 6 0 and Fctrial 6 0, unloading occurs, the damage increments are zero, then Dti+1 = Dtn and Dci+1 = Dcn ; i+1 i+1 (ii) if Fttrial > 0 and/or Fctrial > 0, compute the damage increments; the following three possibilities at most should be considered: tension activation: 1Dt > 0 and 1Dc = 0, with 1Dt computed from Ft (εi+1 , Dtn + 1Dt , Dcn ) = 0; this solution is valid if Fc (εi+1 , Dtn + 1Dt , Dcn ) 6 0; compression activation: 1Dt = 0 and 1Dc > 0, with 1Dc computed from Fc (εi+1 , Dtn , Dcn + 1Dc ) = 0; this solution is valid if Ft (εi+1 , Dtn , Dcn + 1Dc ) 6 0; vertex activation: 1Dt > 0 and 1Dc > 0, with 1Dt and 1Dc computed from Ft (εi+1 , Dtn + 1Dt , Dcn + 1Dc ) = 0 and Fc (εi+1 , Dtn + 1Dt , Dcn + 1Dc ) = 0. With the calculated values of damage variables, update the stresses. • Convergence test – Compute residual force vectors Ri+1 with the current value of stress and check convergence comparing a norm of the residuals with a fixed tolerance. If the convergence test is not passed, perform a new iteration. The scheme of the iteration is summarized in box 1. It is worth noting that, despite the non-local nature of the material model, the corrector phase can be performed locally, at the Gauss point level. This is due to the particular choice to introduce non-local values of strain measures keeping damage variables locally defined. Thus in a strain driven process the non-local quantities are the input of the local corrector phase which furnishes updated estimates of the locally defined damage variables. 5. Numerical tests In order to investigate the predictive capability and the regularisation properties of the proposed model and the effectiveness of the numerical strategy, two simple preliminary numerical tests have been performed.

A non-local damage model

13

Box 1. Flow chart of the finite-step iterative scheme.

In both examples, within a given time step, the secant prediction matrix computed at the end of the previous step has been used for all iterations. This choice results in a large number of iterations, but is quite efficient since it does not require the updating of the stiffness at each iteration. The proper definition of a consistent tangent matrix as proposed in (Jirásek, 1999) is rather cumbersome for this model and further work is required in this direction. 5.1. Tension test on a four edge notched specimen The first example concerns the simulation of the experimental tension test performed by Hassanzadeh (1991) on four edge notched concrete specimen (figure 3(a)). This test has already been considered in the literature as a benchmark (di Prisco et al., 2000). The material experimental compressive strength given by the author is 50 MPa, which corresponds to a concrete grade C40 in the CEB-FIP model code; therefore, the other data used to identify the model parameters have been taken from this code. Note that during the test only damage in tension is activated, therefore only the material parameters defining the loading function in tension are of interest. Both material and resulting model parameters are collected in figure 3. Even though the presence of the notch on the four edges of the specimen makes the response typically threedimensional, the finite element analyses have been performed under plane strain conditions; to account for the three dimensionality of the problem, the total reaction R has been obtained assuming ‘axial-symmetric-like’ conditions as follows. From the computed nodal forces Fi (see figure 4(a)), forces per unit area qi have been computed dividing Fi by half of the support bi of node i (see figure 4(b)); these load densities have then been multiplied by suitable influence areas (the square frames of width bi /2 in figure 4(c)) instead of the areas

14

C. Comi

Figure 3. Tension test: (a) geometrical data; (b) 472-elements mesh; (c) 908-elements mesh.

Figure 4. Tension test: evaluation of the total reaction: (a) nodal forces; (b) forces per unit area; (c) influence areas adopted; (b) influence areas in plane strain conditions.

A non-local damage model

15

Figure 5. Tension test: experimental (dotted line) and numerical global response.

of 70 mm long strips (see figure 4(d)). The specimen has been discretized by two meshes of constant strain triangular elements (figure 3(b)–(c)). The computed total reactions R versus the imposed displacement of the top of the specimen u are plotted in figure 5 together with the experimental curve (dotted line) taken from (Hassanzadeh, 1991). One can observe that the two meshes give similar responses which are in quite good agreement with the experimental results. Note that the bump experimentally observed in the softening branch is due to rotational instability, occurring when the test apparatus has a low rotational stiffness (Hillerborg, 1989). The two parts of the specimen, on both sides of the damaged zone rotate with respect to each other. This instability effect is not reproduced by the numerical analysis in which equal displacements are imposed to all points on the top face of the specimen. Also the final pattern of damage is almost mesh-independent, see figure 6, the non-local model allows for a spread of the damage over all the finite elements belonging to a zone whose dimensions are fixed by the

(a)

(b)

Figure 6. Tension test: final damage distribution (a) 472-elements mesh (b) 908-elements mesh.

16

C. Comi

(a)

(b)

(c)

(d)

Figure 7. Tension test: damage pattern evolution during the analysis (908-elements mesh).

material length and not by the element size. The evolution of the damaged zone, simulating the formation of macro-cracks leading to the failure of the specimen, is shown in figure 7 for the refined mesh. 5.2. Splitting test on concrete prism The second example concerns the splitting test on prisms of plain micro-concrete, a test which is carried out under monotonically increasing imposed displacement. This example is chosen because both damage mechanisms in tension and compression are activated during this simple test; therefore its simulation requires the correct modelling of not only damage due to tension, but also damage due to compression. A similar test has been simulated with a two-surfaces plasticity model by Feenstra and de Borst (1996). The geometrical data, taken from (Rocco et al., 1995), are shown in figure 8 together with the model parameters employed in the numerical analyses; these latter are identified on the basis of the material data given in that reference (E = 31 000 MPa, fc = 38 MPa, Gf = 0.072 N/mm) and of the data taken from the CEB-FIP

A non-local damage model

17

Figure 8. Splitting test: geometrical and material data.

code. The analyses have been performed discretizing only a quarter of the specimen (shown in figure 8) by two meshes of 1 092 and 1 548 triangular plane strain elements. The only experimental result given in (Rocco et al., 1995) is the value of the reaction at the peak (26.7 kN) which is in good agreement with the computed value (27.2 kN). At the peak, both damage mechanisms are activated: figure 9 shows the corresponding contour plots of the damage variables Dt and Dc obtained with the coarse (left hand side of the contour plots) and refined mesh (right hand side of the contour plots). One can see that damage in tension develops at the centre of the specimen, where a crack forms, while damage in compression is activated below the bearing strips, giving rise to a punch. Again no spurious mesh dependence is visible. A map of the active modes during the analysis is displayed in figures 10(a)–(e). Only the right, upper quarter of the specimen is represented. In black are represented regions where compression is active, in grey are represented regions where tension is active; light grey corresponds to elastic behaviour. Following the sequence, one can see that damage in tension begins at the middle of the specimen (figure 10(a)), then it spreads in both vertical and horizontal directions (figures 10(b)–(c)) and finally it tends to localize in a narrow vertical zone, thus simulating a macrocrack, and in the regions outside this zone elastic unloading occurs (figures 10(d)– (e)). Damage in compression develops below the bearing strips ( figures 10(a) and (b)) and it forms a punch (figures 10(c)–(e)). Figure 11 shows at each Gauss point, the active modes at the same step of the analysis of figure 10(d) (u = 0.0259 mm); the different post-treatment of the data allows the visualisation of the Gauss points where the vertex is active, marked by black circles. It can be observed, and this is true also for all the other steps, that the vertex is active in a very small number of points, differently to what would happen using multi-surfaces plasticity models. This is due to the fact that in the damage models there is no constraint of normality to be fulfilled, thus the line of vertex (or the points of vertex in two dimensions) is not a privileged set of points where the constraint is relaxed, but it is only a single line on a surface of potentially active points.

18

C. Comi

(a)

(b)

Figure 9. Splitting test: (a) tension damage distribution at the peak, left 1 092-elements mesh, right 1 548-elements mesh; (b) compression damage distribution at the peak, left 1 092-elements mesh, right 1 548-elements mesh.

Figure 10. Splitting test: active modes at different steps of the analysis.

A non-local damage model

19

Figure 11. Splitting test: activation of vertex, compression or tension at the peak.

6. Conclusions A non-local isotropic damage model with two damage mechanisms has been developed. Two internal lengths have been introduced as localisation limiters of damage in tension and compression. For the proposed non-local model a bifurcation analysis has shown that the criterion for bifurcation is very similar to that of the underlying local continuum, but with the non-local model the boundary value problem remains well-posed at the onset of bifurcation. In fact while a perturbation with zero wavelength is an admissible solution to the rate problem in the local case, in the non-local case the wavelength is fixed and related to the material internal lengths. For the numerical solution of structural problems in the presence of this enhanced damage model, a simple iterative solution scheme of the finite time-step has been proposed. The effectiveness of the non-local enhancement, which allows the avoidance of spurious mesh dependence of numerical results, has been shown by some numerical tests. The correct identification of the internal, material lengths from experimental tests still remains an open question which, maybe, can be tackled by properly defined inverse problems. The simplicity of the underlying local model, in particular its isotropy, and the simplicity of the non-local enhancement, in which the damage variables rest locally defined, make the proposed model apt to be used in structural problems involving large computations as in dam engineering. The damage induced anisotropy linked to the formation of aligned micro-cracks is not introduced in the material model and is obtained only as a structural effect at the macro-scale. Even though a provision to describe the unilateral effect has been already included, the model is not able to correctly predict the behaviour in situations where loading reversal occurs, due to the lack of permanent strains in its present version. The coupling of the proposed damage model with plasticity and thus the inclusion of permanent strains is now in progress.

20

C. Comi Acknowledgements

This work has been carried out in the framework of the project: “Structural safety assessment of large dams” supported by the Italian Ministry for Universities and Scientific Research. Appendix Consider the rate equilibrium equations for the non-local medium: 

     ∂Ft ∂Fc 1 ∂Fc ∂Ft e div E : ε˙ − A Gc (ε˙ ) − Gt (ε˙ ) + B Gt (ε˙ ) − Gc (ε˙ ) = 0.

H

∂Dc

∂Dc

∂Dt

∂Dt

(A.1)

Note first that, if bifurcation of the strain rates in the form (44) is sought starting from a homogeneous state, the averages of variables at this state do coincide with the corresponding local values, furthermore, only rate quantities give rise to non-zero terms in the divergence. Therefore equation (A.1) can be rewritten as:      ∂Ft ∂Fc 1 ∂Fc ∂Ft e div E : ε˙ − A div Gc (ε˙ ) − div Gt (ε˙ ) + B div Gt (ε˙ ) − div Gc (ε˙ ) = 0.

H

∂Dc

∂Dc

∂Dt

∂Dt

(A.2)

To compute the term div Gt (ε˙ ), taking into account the definitions (37) and (44), one has to evaluate the following terms: term 1 = div

Z 

Z 

V



1 − iqWt (x − s)e : (v0 ⊗ n + n ⊗ v0 ) e−iqn·s ds, 2

(A.3)



1 term 2 = div − iqWt (x − s)H(±tr ε) I : (v0 ⊗ n + n ⊗ v0 ) e−iqn·s ds. (A.4) 2 V Analogous terms, with Wc instead of Wt , intervene in the evaluation of div Gc (ε˙ ). Perform the change of variables x − s = −ξ and note that, being Wt symmetric, Wt (−ξ ) = Wt (ξ), from equations (A.3) and (A.4) one has: term 1 = −iq e : (v0 ⊗ n)

Z V

Wt (ξ ) e−iqn·ξ dξ div e−iqn·x = −q 2 W t (q)(n · e · v0 )n e−iqn·x ,

term 2 = −iqH(±tr ε) I : (v0 ⊗ n)

Z

(A.5)

Wt (ξ ) e−iqn·ξ dξ div e−iqn·x

V

= −q 2 W t (q)H(±tr ε)(n · I · v0 )n e−iqn·x ,

(A.6)

where W t (q) is the Fourier transform of Wt . Using equations (A.5), (A.6) and (37), remembering that the starting state is homogeneous (and, hence, htr+ εit ≡ htr+ εic ≡ tr+ ε, etc.), one can compute the divergence of Gt (ε˙ ) as: 






div Gt (ε˙ ) = −q 2 W t (q) 4µ2 (n · e · v0 ) − 2at K+ tr+ ε + K− tr− ε − bt rt (Dt ) 





× K+ H(tr ε) + K− H(−tr ε) (n · I · v0 ) n e−iqn·x . Using the definition (27) of

∂ft ∂ε



(A.7)

, one can rewrite (A.7) as: 

div Gt (ε˙ ) = −q W t 2



∂ft n· · v0 n e−iqn·x . ∂ε

(A.8)

A non-local damage model

21

Fully analogous passages allow to compute the divergence of Gc (ε˙ ) as: 

div Gc (ε˙ ) = −q 2 W c n ·



∂fc · v0 n e−iqn·x . ∂ε

(A.9)

Substituting (A.8) and (A.9) into (A.2), condition (45) for the admissibility of the perturbation is arrived at.

References Bažant Z.P., Pijaudier-Cabot G., 1988. Non-local continuum damage, localisation instability and convergence. J. Appl. Mech. 55, 287–293. Bažant Z.P., Pijaudier-Cabot G., 1989. Measurement of characteristic length of non-local continuum. J. Eng. Mech. ASCE 115, 755–767. Benallal A., Billardon R., Geymonat G., 1988. Some mathematical aspects of the damage softening problem. In: Mazars J., Bažant Z.P. (Eds.), Cracking and Damage, Elsevier, Amsterdam, pp. 247–258. Benallal A., Pijaudier-Cabot G., 1993. Strain localisation and bifurcation in a non-local continuum. Int. J. Solids Structures 30 (13), 1761–1775. Borino G., Fuschi P., Polizzotto C., 1999. A thermodynamic approach to non-local plasticity and related variational principles. J. Appl. Mech. 66 (4), 952–963. Borré G., Maier G., 1989. On linear versus nonlinear flow rules in strain localisation analysis. Meccanica 24, 36–41. CEB-FIP Model code 1990, Bullettin d’information, CEB, Lausanne. Comi C., Berthaud Y., Billardon R., 1995. On localisation in ductile-brittle materials under compressive loadings. Eur. J. Mech., A/Solids 14 (1), 19–43. Comi C., 1999. Computational modelling of gradient-enhanced damage in quasi-brittle materials. Mech. Cohes.-Frict. Mater. 4 (1), 17–36. Comi C., Perego U., 2000. A bi-dissipative damage model for concrete with applications to dam engineering. In: Proc. ECCOMAS 2000, Barcelona, Spain. Comi C., Rizzi E., 2000. On bifurcation in local and non-local materials with tension and compression damage. In: Proc. ECCOMAS 2000, Barcelona, Spain. Desoyer T., Cormery F., 1994. On uniqueness and localisation in elastic-damage materials. Int. J. Solids Structures 31, 733–744. di Prisco M., Ferrara L., Meftah F., Pamin J., de Borst R., Mazars J., Reynouard J.M., 2000. Mixed mode fracture in plain and reinforced concrete: some results on benchmark tests. Int. J. Fracture 103 (2), 127–148. Dubé J.-F., Pijaudier-Cabot G., La Borderie C., Glynn J., 1994. Rate dependent damage model for concrete- Wave propagation and localisation. In: Mang H., Bicanic N. and de Borst R. (Eds.), Computer Modelling of Concrete Structures, pp. 313–322. Feenstra P.H., de Borst R., 1996. A composite plasticity model for concrete. Int. J. Solids Structures 33, 707–730. Fichant S., La Borderie C., Pijaudier-Cabot G., 1999. Isotropic and anisotropic descriptions of damage in concrete structures. Mech. Cohes.-Frict. Mater. 4, 339–359. Frémond M., Nedjar B., 1996. Damage, gradient of damage and principle of virtual power. Int. J. Solids Structures 33, 1803–1103. Halm D., Dragon A., 1996. A model of anisotropic damage by mesocrack growth; unilateral effect. Int. J. Damage Mech. 5, 384–402. Hassanzadeh M., 1991. Behaviour of fracture process zones in concrete influenced by simultaneously applied normal and shear displacements. PhD Thesis, Lund Institute of Technology, Sweden. Hill R., 1959. Some basic principles in the mechanics of solids without a natural time. J. Mech. Phys. Solids 7, 209–225. Hillerborg A., 1989. Stability problems in fracture mechanics testing. In: Shah S.P., Swartz S.E., Barr B. (Eds.), Fracture of concrete and rock, Elsevier Applied Science, pp. 369–378. Jansen D.C., Shah S.P., 1997. Effect of length on compressive strain softening of concrete. J. Eng. Mech. ASCE 123, 25–35. Jirásek M., 1999. Computational aspects of non-local models. In: Proc. European Conf. Comp. Mech., München. Kupfer H., Hilsdorf H.K., Rush H., 1969. Behavior of concrete under biaxial stress. ACI J. 66, 656–666. Lee J., Fenves G.L., 1998. Plastic-damage model for cyclic loading of concrete structures. J. Eng. Mech. ASCE 124, 892–900. Maier G., Hueckel T., 1979. Nonassociated and coupled flow rules of elastoplasticity for rock-like materials. Int. J. Rock Mech. Min. Sci. & Geomech. Abstr. 16, 77–92. Mazars J., 1984. Application de la mécanique de l’endommagement au comportement non linéaire et à la rupture du béton de structure. Dr of Science thesis University of Paris 6. Mazars J., 1986. A model of a unilateral elastic damageable material and its application to concrete. In: Fracture Toughness and Fracture Energy of Concrete, Elsevier, New York. Mazars J., Pijaudier-Cabot G., 1989. Continuum damage theory – application to concrete. J. Eng. Mech. ASCE 115, 345–365. Ortiz M., 1985. A constitutive theory for the inelastic behavior of concrete. Mech. of Materials 4, 67–93. Peerlings R.H.J., de Borst R., Brekelmans W.A.M., de Vree J.H.P., 1996. Gradient-enhanced damage for quasi-brittle materials. Int. J. Num. Meth. Eng. 39, 3391–3403. Pijaudier-Cabot G., 1995. Non local damage. In: Mülhaus H.-B. (Ed.), Continuum Models for Materials with Microstructure, Wiley, pp. 105–143. Pijaudier-Cabot G., Bažant Z.P., 1987. Non-local damage theory. J. Eng. Mech. ASCE 113, 1512–1533.

22

C. Comi

Ramtani S., Berthaud Y., Mazars J., 1992. Orthotropic behavior of concrete with directional aspects: modelling and experiments. Nuclear Eng. Design 133, 97–111. Rice J.R., 1976. The localisation of plastic deformation. In: Koiter W.T. (Ed.), Theoretical and Applied Mechanics, North-Holland, Amsterdam, 207–220. Rocco C., Guinea G.V., Planas J., Elices M., 1995. The effect of the boundary conditions on the cylinder splitting strength. In: Wittman F.H. (Ed.), Fracture Mechanics of Concrete Structures, AEDIFICATIO Publishers, Freiburg, pp. 75–84. Simo J.C., Ju J.W., 1987. Strain- and stress-based continuum damage models – I. Formulation. Int. J. Solids Structures 23, 821–840. Steinmann P., Stein E., 1994. Finite element localisation analysis of micropolar strength degrading materials. In: Mang H., Bicanic N., de Borst R. (Eds.), Computer Modelling of Concrete Structures, pp. 435–444.