Integral-type regularization of non associated softening plasticity for quasi brittle materials

Computers and Structures 224 (2019) 106120

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

Integral-type regularization of non associated softening plasticity for quasi brittle materials G. Mazzucco a,⇑, G. Xotta a, V.A. Salomoni a,b, C. Majorana a a b

Department of Civil, Environmental and Architectural Engineering (DICEA), University of Padua, Via F. Marzolo 9, 35131 Padua, Italy Department of Management and Engineering, University of Padua, Stradella S. Nicola 3, 36100 Vicenza, Italy

a r t i c l e

i n f o

Article history: Received 27 March 2019 Accepted 2 September 2019 Available online 17 September 2019 Keywords: Concrete behavior Three dimensional modeling Elastoplastic constitutive law Non-local regularization technique

a b s t r a c t Within the framework of infinitesimal elasto-plasticity for describing concrete behaviour, this paper proposes a new non-local constitutive model. Particularly, a non-local integral-type regularization technique is applied to the three-invariant Menétrey-Willam model, here modified by the inclusion of a nonassociated flow rule. The resulting model allows for taking into account the role of the intermediate principal stress and the Lode parameter, as well as to control inelastic dilatancy. A regularization scheme is hence introduced to circumvent strain localization and resolve the incapability of objectively describing localized material failure. One- and three-dimensional Finite Element analyses are developed to prove the soundness of the novel implemented regularized scheme and the upgraded response with respect to a local approach. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction Frictional materials such as concrete are characterized by a complex behavior whose realistic description within the framework of elasto-plasticity requires a non-associated flow rule and a transition from hardening to softening under tension or lowconfined compression. Unfortunately, such requirements make standard constitutive models incapable to objectively describe localized material failure for the outcoming of strain localization [1–6], a phenomenon characterized by the concentration of deformations in infinitely narrow zones. From a mathematical viewpoint, this aspect is related to the loss of ellipticity of the governing differential equations (turning to hyperbolic) and the boundary value problem becomes illposed, i.e. the solution is not unique and depends continuously on the given data. Numerically, this translates into the dependence of the results on the size of the finite elements. In the framework of single-phase solids, several techniques for regularizing the quasi-static (or dynamic), boundary (or initial) value problem with material instability can be found in literature [7–13]. In the numerical simulation of strain localization on multiphase materials, an internal length scale may be naturally intro-

⇑ Corresponding author. E-mail address: [email protected] (G. Mazzucco). https://doi.org/10.1016/j.compstruc.2019.106120 0045-7949/Ó 2019 Elsevier Ltd. All rights reserved.

duced through Darcy’s law and the related constitutive relationship for permeability [14]. When considering multiphase porous media, a series of different techniques should be adopted, such as: gradient plasticity models, including the effect of higher-order deformation gradients within the constitutive equations or the gradient of internal variables [10,15-24], rate dependent constitutive models [25,26], Cosserat or micro-polar models (where the regularizing effect comes from the introduction of additional static and kinematic variables such as couple stresses and micro curvatures [27,7,28,29]) or viscoplasticity [11,30–33]. An alternative solution is the so-called integral or non-local continuum [38,39,34,41,40,9,35–37,43,42], where a certain state variable is replaced by its non-local counterpart or by a linear combination of the local and non-local counterparts, obtained through a space-weighted average which reflects the ability of the microstructure to transmit information to neighboring points. The non-local quantity is computed with an integral format, in which the associated intrinsic length scale influences the weight amplitude in the vicinity of a material point. In the present work the three-invariant Menétrey-Willam yield function [44,45], chosen to consider the significant role of the intermediate principal stress and the Lode parameter, is innovately regularized via a non-local integral-type technique. Octree structures are also originally applied in this context, so to optimize the search algorithm for the integration points. A realistic description of concrete post-peak behaviour, characterized by growth and coalescence of microcracks and openings of

2

G. Mazzucco et al. / Computers and Structures 224 (2019) 106120

Nomenclature

Ag ; Bg c Ce ; Cp e E Er f fe Df e fi fc ft g G J2 ; J3

plastic potential parameters cohesion parameter elastic and plastic constitutive tensors eccentricity elastic modulus elastic domain yield function external forces vector incremental external forces vector internal forces vector compressive strength tensile strength plastic potential function shear modulus second and third invariants of the deviatoric stress tensor local equivalent plastic strain non-local equivalent plastic strain

k ~ k  k

mixed local/non-local equivalent plastic strain increments of the equivalent plastic strains

~ Dk  Dk; Dk; kt K K l m

~ constants in the expression of Dk plastic internal variable hardening and softening functions initial hardening value function that defines the roundness of the failure surface residual vector radius of the search sphere deviatoric stress tensor parameter that controls the slope of the softening function qs displacements vector incremental displacements vector volume position vectors second and fourth order identity tensors non-local weight Gauss-like function Kronecker delta total strain tensor elastic strain tensor plastic strain tensor plastic flow Lamé constants Poisson’s ratio Haigh-Westergaard coordinates stress tensor Helmholtz free-energy

nlA ; nlB

The main symbols used in this work are:

q qh ; qs qh0 r r R s t u Du V x; y; z 1; I

a a1 d

 e p k K; l

equivalent strain at uniaxial concrete strength bulk modulus stiffness matrix characteristic length weighting parameter

m

ðn; q; hÞ

r

W

or in terms of the bulk and shear moduli, K and G:

the micro voids within cement paste and interfacial transition zone, has been attained via a series of numerical analyses, so confirming the potentialities of the chosen regularization technique.

K¼K

2. Theoretical framework

Then, the constitutive relations can be defined as follows:

A convex elastic domain Er in the Cauchy stress space r defined with a yield surface f is considered:

Er ¼ fðr; qÞ 2 S  Rm j f ðr; qÞ 6 0g

ð1Þ

where S is the space of symmetric second order tensors and q is the plastic internal variable. The constitutive relations can be expressed by using the Helmholtz free-energy W, function of the elastic strain e and the plastic variable k (conjugated with q), defined as:

1 1 Wð ; kÞ ¼ W þ W ¼ e : Ce : e þ kCp k 2 2 e

e

p

ð2Þ

where Ce is the fourth order elasticity tensor, defined in terms of the Lamé constants K and l:

Ce ¼ K1  1 þ 2lI:

ð3Þ

while Cp is the plastic constitutive tensor. In Eq. (3), 1 is the second order identity tensor and I is the fourth order identity tensor which is expressed as Iijkl ¼ 12 ðdik djl þ dil djk Þ. The Lamé constants may be expressed in terms of the elastic modulus and Poisson’s ratio, E and m:

K¼

Em ; ð1 þ mÞð1  2mÞ

l¼

E 2ð1 þ mÞ

r¼

2G ; 3

@ Wðe ; kÞ ; @

l ¼ G:

q¼

ð5Þ

@ Wðe ; kÞ @k

ð6Þ

where, if isotropic elasticity is assumed, the constitutive relation for the stress tensor follows the second law of thermodynamics:

r¼

@W ¼ Ce :  e : @

ð7Þ

If small strains are considered, the total infinitesimal strain  can be decomposed into an elastic and a plastic part:

 ¼ e þ p :

ð8Þ

In order to express the rate of the plastic strain _p , it is assumed that a plastic potential function gðr; qÞ exists, generally different from the yield function f in order to define a non-associated flow rule:

_p ¼ k_

@gðr; qÞ @r

ð9Þ

where k is the plastic multiplier. The evolution of the internal variable k is defined as:

@gðr; qÞ k_ ¼ k_ @q

ð10Þ

Moreover the loading conditions need to be satisfied:

ð4Þ

k_ P 0;

f ðr; qÞ 6 0;

_ ¼ 0: kf

ð11Þ

3

G. Mazzucco et al. / Computers and Structures 224 (2019) 106120

2.1. Yield function

g ¼ Ag

When a material reaches its yield limit it undergoes permanent deformation. The onset of yielding for a material is defined by the yield function f ðr; kÞ: the material is in elastic region when the function is strictly less than zero while, if f ðr; kÞ ¼ 0, plastic deformations can trigger. The Menétrey-Willam [44,45] format of the Willam-Warnke yield function [46] is here adopted, characterized by a dependence on the deviatoric stress tensor s through its second and third invariants, J 2 and J 3 , together with a non-associated flow rule defining a plastic potential function dependent on the condition of concrete confinement. For a better description of the plasticity model the HaighWestergaard (H.W.) coordinates of stress are used:

trðrÞ n ¼ pffiffiffi ; 3

q¼

pffiffiffiffiffiffiffi 2J 2 ;

cos3h ¼

pffiffiffi 3 3J 3 2J 23=2

:

ð12Þ

Hence the yield function f can be written as:

f ðn; q; h; kÞ ¼

 2   3 q q c q n þ h pffiffiffi rðh; eÞ þ pffiffiffi  qh qs 6 0 2 fc fc 6 3

ð13Þ

where f c is the compressive strength of the material; qh and qs are the hardening and softening functions, respectively, which depend on the internal variable k. The cohesion parameter of the material c ¼ cðf c ; f t ; eÞ and rðh; eÞ function of the Lode angle h and eccentricity e, defining the roundness of yield surface along the deviatoric plane, are expressed as: 2

cðf c ; f t ; eÞ ¼

rðh; eÞ ¼

2

3e fc  ft ; ðe þ 1Þ f c  f t

ð14Þ

4ð1  e2 Þcos2 ðhÞ þ ð2e  1Þ2

: 1=2 2ð1  e2 ÞcosðhÞ þ ð2e  1Þ½4ð1  e2 Þcos2 ðhÞ þ 5e2  4e ð15Þ

while the variables qh ðkÞ and qs ðkÞ, are defined in agreement with [47]:

qh ðkÞ ¼

8 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2ffi > > k k > if < qh0 þ ð1  qh0 Þ 1  tkt > 1 > > :

if

k < kt k P kt

ð16Þ

where qh0 is the initial value of hardening (initial yield surface), generally assumed equal to 0:2, and kt is the equivalent strain at uniaxial concrete strength. The function qs , controlling the softening branch of the material behaviour, takes the form:

8 1 > > > 2 32 > > < 6 7 qs ðkÞ ¼ 4  1 2 5 > > n1 1 > > 1þ n 1 > 2 :

if

k < kt

if

k P kt

ð17Þ

where n1 ¼ k=kt ; n2 ¼ ðkt þ tÞ=kt and the parameter t controls the slope of the softening curve. The rate of the plastic variable k_ is defined as the rate of the volumetric plastic strain k_ ¼ trð_p Þ. In agreement with Grassl [48], for a better description of concrete material behaviour beyond the elastic limit, a nonassociated flow rule is hence adopted. The plastic potential, defined in terms of the Haigh-Westergaard coordinates n and q, is:

q

!2

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi f c qh ðkÞqs ðkÞ

þ

fc

Bg q þ n pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qh ðkÞqs ðkÞ

ð18Þ

where Ag and Bg are two material parameters, determined by means of the axial strain at maximum stress in uniaxial compression and in a triaxial compressive state (see [48]). 2.2. Non local approach Within the framework of elastoplasticity, the softening branch and the adoption of a non-associated flow rule leads to loss of ellipticity for the governing equations of a continuum model when localization occurs under static loading conditions [7]. In this case the thickness of the localized zone is dependent on the fineness and orientation of the adopted spatial discretization [5]. The localized plastic strain width, related to the heterogeneous material microstructure, can be correctly predicted only by models accounting for an intrinsic parameter with the dimension of a length l [9]. Among the different regularization techniques previously recalled, in the case of a non-local integral-type model the regularization function can be written as:

~ kðxÞ ¼

Z

V

aðy; xÞ kðyÞ dVðyÞ

ð19Þ

~ is the non-local hardening/softening variable and a is a where k weight function of the position vector of the current point x and the position vector of a general neighboring point y inside the volume V of the body. This technique, calledpurely non-local, while regularizing the global response it cannot prevent mesh-dependency [42,9,38,43]. To get around this issue, in the present work a mixed local nonlocal formulation is adopted [9,39,41,40,43,38,42], where the non~ is replaced by a linear combination of the local, k, local variable k ~ counterparts: and non-local, k,

 ~ þ ð1  mÞk ¼ m kðxÞ ¼ mk

Z

V

aðy; xÞ kðyÞ dVðyÞ þ ð1  mÞk

ð20Þ

wheree m is a weighting parameter that controls the size of the plastic zone and the shape of the plastic strain profile [41]. This parameter can be chosen almost arbitrarily but, with m close to 1, the search radius R has to be much larger than the plastic region width, so increasing the computational time needed to perform the procedure. Otherwise, for very large values of m, the search radius would be much smaller than the plastic region width, and a fine mesh discretization is needed to capture the non-local effects. In this work m has been assumed equal to 2, in agreement with [9]. Moreover, it is also worth noting that, for values of m respectively equal to 0 and 1, Eq. (20) reduces to the local and the purely nonlocal formulations. The consistency conditions, reported in Eq. (11), can be rewritten as [43]:

k_ P 0;

 6 0; f ðr; kÞ

 ¼ 0: _ ðr; kÞ kf

ð21Þ

and the weight function a in Eq. (19), so to satisfy the normalizing condition [9]:

Z

V

aðy; xÞdVðyÞ ¼ 1

8x 2 V:

ð22Þ

is set as:

aðy; xÞ ¼ R

a1 ðky  xkÞ a ðkz  xkÞdVðzÞ 1 V

ð23Þ

where a1 ðdÞ is a monotonically decreasing positive function of the distance d ¼ ky  xk.

4

G. Mazzucco et al. / Computers and Structures 224 (2019) 106120

3. Non local numerical procedure Numerically, the mixed local non-local internal variable at the current step n, is evaluated at each integration point through the return mapping procedure [50], given the value of the variable at the previous step ðn  1Þ and its increment for the interval ½ðn  1Þ; n:

 ¼k   k ðnÞ ðn1Þ þ Dk

ð25Þ

 is as follows: From Eq. (20), the increment Dk

DkðnÞ ¼ mDk~ðnÞ þ ð1  mÞDkðnÞ

Fig. 1. Non-local connection scheme in a bidimensional representation.

For the non-local weight a1 , several choices can be found in literature [49,9,33]; in this work a Gauss-like function, in agreement with [33], is assumed:

a1 ¼

8   2  > 2jjyxjj > exp  if > < l > 0 > > :

if

~ðnÞ is the integral part of the non-local procedure: where Dk

Dk~ðnÞ ¼

Z

V

jjy  xjj 6 R jjy  xjj > R

ð26Þ Z

aðy; xÞDkðyÞðnÞ dV ¼

V

aðy; xÞk_

@g dij dV @ rij

ð27Þ

that, during a FE analysis, is replaced by a summation operator [42]

ð24Þ

np X

~ Þ ¼ Dkðx i ðnÞ where l is the non-local characteristic length. In Eq. (24) the variable R is the radius of the ‘‘search sphere” (SS) where the non-local interaction is evaluated; only the Gauss points inside SS are assumed connected with the current integration point and a1 is consequently found. An example of the SS is reported in Fig. 1; the current integration point xi is located at the center of the search sphere and only the connected Gauss points are used to evaluate the non-local weight. The SS procedure has been implemented to reduce computational costs related to the non-linear solver. As specifIed in [9] the boundary may affect the weight function, reported in Eq. (23). Indeed, for a weight function with a bounded support, this influence is present if the distance of the plastic zone from the boundary (db ) is smaller than 2R; the search sphere then needs to be scaled to a value equal to Rs so to obtain a db greater than 2Rs , in agreement with [9]. Otherwise, if no boundary scaling is performed or if db > 2R, the solution is the same as for an infinite domain.





a1 ð yp  xi ÞÞDV p Dkðyp ÞðnÞ

p¼1

ð28Þ

np X

a1 ð yp  xi ÞDV p

p¼1

where np is the total number of Gauss points inside a search sphere of radius R centered in xi (see Fig. 1), while yp and DV p are respectively the position vector and the volume related to each integration point inside the SS. Expression (28) can be rewritten and compacted, by expressing the weight function in agreement with Eq. (23) and evaluating the local increments Dkðyp Þ at the previous step ðn  1Þ except when yp = xi :

~ Þ ¼ Dkðx i ðnÞ

np X

!

aðyp ; xi ÞDV p Dkðyp Þðn1Þ

p¼1

þ yp –xi

ð29Þ

þ aðxi ; xi ÞDV i Dkðxi ÞðnÞ but, for the Gauss point with a position vector equal to xi , the weight function is aðxi ; xi Þ ¼ 1; so obtaining:

Fig. 2. Speed test comparison in terms of solution time: classic loop vs. octree structure.

5

G. Mazzucco et al. / Computers and Structures 224 (2019) 106120

~ Þ ¼ Dkðx i ðnÞ

np X

!

aðyp ; xi ÞDV p Dkðyp Þðn1Þ

p¼1

þ yp –xi

ð30Þ

þ Dkðxi ÞðnÞ DV i

as done in the present work, only the variables referred to the current point must pass in the material response function, according to the classic finite element implementation. The scheme of the non-local numerical procedure is reported in Box 1.

In Eq. (30) it can be noted that the unique variable part is the local increment of the equivalent strain Dkðxi ÞðnÞ :

~ Þ ¼ nlA þ nlB Dkðx Þ Dkðx i ðnÞ i ðnÞ

ð31Þ

while nlA and nlB are constant and can be defined as:

nlA ¼

np X p¼1

!

aðyp ; xi ÞDkðyp Þðn1Þ DV p yp –xi

nlB ¼ DV i

ð32Þ

The terms in Eq. (32) can then be computed outside the integration ~ point, thus reducing the computational time to evaluate DkðxÞ ðnÞ . Indeed, if the non-local integral is executed at the integration point, it means that the global array, where all the local variables are stored (in this case k), needs to pass into the function where the material response is estimated (stress state and tangent constitutive matrix), but normally at this level of the code the global vector, where the state variables are located, is not accessible and further code modifications are needed to evaluate the non-local procedure. Othewise, if nlA and nlB are computed outside the integration point,

Fig. 5. Punch test geometry and boundary conditions.

Table 1 Constitutive parameters for the elasto-plastic material. E ½MPa

m

Fig. 3. Integration points’ positions.

f c ½MPa f t ½MPa e

30000 0:15 32 3 0:55

t kt K0 Ag Bg

Fig. 4. Normalized equivalent strain vs integration points’ position with: R = 50 mm (a); R = 300 mm (b) and R = 600 mm (c).

0:02 0:0091 0:6 7:46 10:384

6

G. Mazzucco et al. / Computers and Structures 224 (2019) 106120

Box 1. Non-local Newton-Raphson scheme. function nonLinearSolver() . . . P x //global array of the Gauss points’ position vectors ½totP  dimSpace P dV //global array of the volumes related to each Gauss point ½totP ... ... for n ¼ 1 to totSteps (i) compute global external force vector f e f eðnÞ ¼ f eðn1Þ þ Df e (ii) set initial guess for displacements u and initial residual r uðnÞ ¼ uðn1Þ jjrjj ¼ 1 (iii) compute the parameters, P nlA and P nlB , for calculating ~ (see Eqs. (31, 32)) the non-local variables Dk ðnÞ

[P nlA ; P nlB ] = getNLparameters(P x; P dV; P Dkðn1Þ ) while jjrjj > tolerance (iv) assemble global stiffness matrix K KðnÞ ¼

@rðuðnÞ ;P nlA ;P nlB ;...Þ @uðnÞ

(v) compute global internal force vector f i and update residual r ¼ f eðnÞ  f iðnÞ (vi) compute incremental displacements Du and update total displacements u

Du ¼ K1 ðnÞ r uðnÞ ¼ uðnÞ þ Du end while end for ... end function function [P nlA ; P nlB ] = getNLparameters(P x; P dV; P Dkðn1Þ ) for ip ¼ 1 to totP P nlA ðipÞ // (see Eq. (32)) P nlB ðipÞ // (see Eq. (32)) end for end function function [KðnÞ ] = assembly(uðnÞ ; P nlA ; P nlB , . . .) for e ¼ 1 to totElements ... for P ¼ 1 to totP [ e P r; e P C ] = getConstMat(P ; P D; P nlA ; P nlB , . . .) ... end for ... end for end function ~ i Þ , it is necessary to evaluate the distance When defining Dkðx ðnÞ of the current integration point (xi ) with respect to the other ones belonging to the search sphere; to perform this procedure n2 checks are required, so considerably increasing the computational time in case of more refined mesh discretizations. The use of an optimized procedure for searching the integration points is therefore fundamental. To this end, a more complex structure named octree, normally used to store data, can be adopted so to increase the computational performance; through this method, the three-dimensional space is recursively subdivided into eight octants where the integration

points are located. This method is frequently and successfully used in the Discrete Element Formulation to increase the collision detection algorithm efficiency [51]. In the present work, the adoption of octree structures allows to reduce the number of checks, for the search of the integration points, from n2 to n  logðnÞ. A speed test is reported in Fig. 2, displaying the solution times, classic loop vs. octree structure, necessary to define the integration for a search sphere with a radius R = 30 mm in a cubic space with different numbers of Gauss points randomly located (no parallel computing has been taken into account). Particularly, a 100 mm side cube is investigated, for different thickening levels. Each octree can contain a maximum of 64 points; if this limit is exceeded, the subspace of the ‘‘parent” cube is divided into 8 ‘‘children” cubes and so on until a maximum number of 64 points is achieved within each cube. 4. Numerical examples 4.1. One dimensional non-local behaviour A strip with 26 integration points (see Fig. 3) has been considered, characterized by a constant spacing equal to 100 mm (the geometric position of the integration points ranges from 0 to 2500 mm). An initial strain field is hypothetically assumed so that a localization is imposed in the middle of the sample (Fig. 4); subsequently three different radius values and characteristic lengths are used to evaluate their effects on the non-local procedure (Eq. (20)). By first adopting a radius R = 50 mm, less than the distance between two consecutive points, clearly non-localization can not be performed, independently on the value of the characteristic length (all the curves are superimposed, Fig. 4a; on the contrary, when R exceeds the distance between two integration points (e.g. R = 100 mm), the non-local procedure is activated and the decrease in the strain peak depends on the characteristic length (Fig. 4b and c). Through Eq. (24), the maximum and minimum a1 can be defined. The maximum value occurs when the distance jjy  xjj ¼ 0 and in this case the weight a1 ¼ 1. On the other hand, a non-zero minimum value occurs when jjy  xjj ¼ R; therefore, once assumed R = 300 mm, a1 for different non-local lengths is respectively equal to

a1 ðl ¼ 100Þ ’ 0; a1 ðl ¼ 300Þ ¼ 1:8  101 ; a1 ðl ¼ 600Þ ¼ 3:7  101 . This means that, if R is large enough to contain the strain band, by increasing the non-local length also the weight of the farthest integration points (characterized by lower strain values) increases, and consequently a decrease in the maximum strain occurs when y is related to the maximum value of k. Differently, if y is outside the strain band but R intercepts some integration points within the strain band (e.g. when l = 300 or 600 mm as in Fig. 4b), the strain value estimated in y increases. More importantly, and as evidenced by comparing the results of Fig. 4b and c, the non-local solution becomes substantially independent on R once the radius is (even slightly) larger than the strain band, whereas the dependence on the characteristic length is clearly maintained (in fact in non-local models such a term becomes a constitutive parameter). This has a relevant consequence in terms of efficiency of the numerical procedure, being the computational costs higher for a higher R, but hence being possible to fix an ‘‘optimal” SS to achieve adequate computational times. Hence the role of R is only that of optimizing the speedup of the procedure. 4.2. Punch test A punch length L = 100 mm is now assumed, together with two different materials that characterize a 3D domain (indicated with 1

G. Mazzucco et al. / Computers and Structures 224 (2019) 106120

and 2 in Fig. 5), to evaluate the response of the model. A very rigid linear elastic material (1) is used to define the plate on which the load is applied, with a Young’s modulus equal to E = 2000000 MPa and a Poisson’s ratio m = 0.15. Material 2 is instead assumed to nonlinearly behave, in agreement with the constitutive law described in the previous Sections (data are reported in Table 1). Parameters Ag and Bg , used for the definition of the plastic potential, are in agreement with [48] and a low confinement condition is considered. Three different discretizations are adopted, characterized by hexahedral finite elements with linear shape functions and full integration. The corresponding element size chosen for the three meshes is reported in Table 2. A constant domain thickness equal to 25 mm is assumed, together with symmetry conditions along the vertical axis containing the load force P. Fig. 6 depicts the resulting plastic strain localizations with k P 5:5  104 , with the peak of plastic strain located within a single finite element. To stabilize the solution, the non-local procedure previously described needs to be activated. Normally the characteristic length

7

is taken as a material characteristics, but in this case this value is chosen to stabilize the strain band starting from the coarser mesh (MESH02). The radius of the search sphere is assumed to be equal to R = 30 mm (see Fig. 7). When reducing the size of the adopted finite element, the number of integration points contained into a single search sphere increases; hence for model MESH02 a sphere intercepts about 244 points, for MESH01 about 896 and for MESH00, with an element size equal to 25 mm, the number of integration points contained into a single SS is about 3584. Thus, by refining the mesh (the characteristic length is kept equal to 20 mm), the procedure becomes computationally heavier. The results in terms of non-local equivalent strain are shown in Fig. 8 where it is possible to notice a stabilization of the solution. Fig. 9 reports the spatial evolution of the equivalent strain (along X direction) for the three different meshes and the two different approaches at the end of the analysis, for a reference Z level. If the non-local procedure is not activated (continuous lines in Fig. 9), at a decrease of the element size clearly corresponds a decrease of the strain bandwidth, while the peak increases. Differently, when taking into account a non-local formulation, the results for all models are mesh-size independent, as illustrated by dotted lines of Fig. 9.

Table 2 Element sizes adopted for the three models. Model MESH02 MESH01 MESH00

el size ½mm 100 50 25

4.3. Tensile test A tensile test is reconstructed now adopting linear tetrahedral elements with reduced integration and different levels of mesh discretization, similarly as in the previous example.

Fig. 6. Strain localization for MESH02 (a); MESH01 (b) and MESH00 (c).

Fig. 7. Radius of the search sphere for MESH02 (a); MESH01 (b) and MESH00 (c).

8

G. Mazzucco et al. / Computers and Structures 224 (2019) 106120

Fig. 8. Non-local equivalent strain for MESH02 (a); MESH01 (b) and MESH00 (c).

Fig. 9. Spatial evolution of equivalent strain, with and without non-local procedure, for different mesh discretizations.

Fig. 10. Sample geometry (a), and search sphere (b).

G. Mazzucco et al. / Computers and Structures 224 (2019) 106120

Fig. 11. MESH02 (a); MESH01 (b); MESH00 (c).

Fig. 12. Strain localization for MESH02 (a); MESH01 (b) and MESH00 (c).

Fig. 13. Spatial evolution of the equivalent strain, with and without non-local procedure, for different mesh discretizations.

9

10

G. Mazzucco et al. / Computers and Structures 224 (2019) 106120

Fig. 14. Normalized load vs. normalized displacement in non-local models for different mesh discretizations.

Fig. 15. Split test geometry (a); example of mesh discretization (b).

G. Mazzucco et al. / Computers and Structures 224 (2019) 106120

Taking into account only one integration point for each element, the time needed to define the number of points included within a SS is reduced if compared to hexahedral finite elements having the same mesh density. A rectangular wall, having dimension L = 1000 mm and thickness t = 100 mm, is considered (Fig. 10); a notch is present in the middle of its longitudinal edge to fix the location of plastic strain triggering, characterized by a size of 2t and h = 20 mm. The sample is constrained to avoid rigid motion and a constant vertical displacement u = 0.25 mm is applied to top and bottom edges. The same constitutive parameters as of Section 4.2 are adopted, apart from Poisson’s ratio which is assumed equal to 0 in order to remove plastic localization phenomena that can occur close to the restrained nodes. Model MESH02 now consists in 1200 linear tetrahedral finite elements, MESH01 in 9600 elements and MESH00 in 76800 ele-

11

ments (Fig. 11), with one, two and four elements rows, respectively, along the wall thickness. When the non-local procedure is inactive, as evidenced in the previous test the plastic strain localizes within a bandwidth affecting one or two elements rows (Fig. 12). A search radius having almost the same size of the notch length (i.e. R = 150 mm; see Fig. 10b) and an intrinsic length l = 100 mm are assumed. With such a choice, MESH02 presents a maximum of 18 integration points, with an average distance among them equal to 50 mm and consequently an efficient regularization is not obtained. This phenomenon is visible in Fig. 13 (reference x level of 100 mm), where the equivalent strain curves referred to the local/non-local approach for MESH02 are nearly superimposed; differently, for MESH01 and MESH00 the non-local regularization is correctly obtained. The results are reported at 78% of the total applied load.

Fig. 16. Splitting test meshes and equivalent strain contour maps for local (L) and non-local (NL) approaches.

Fig. 17. Load vs. vertical displacements comparisons.

12

G. Mazzucco et al. / Computers and Structures 224 (2019) 106120

Table 3 Mesh details: el size is the average element size; d.o.f is the total number of degrees of freedom; tot IP is the total number of integration points. model MESH01 (case a) MESH02 (case b) MESH03 (case c)

el size ½mm

d.o.f.

tot IP

5.0 3.0 2.5

4914 20880 34809

9928 47328 80640

Similar considerations arise when analysing the normalized loads F=F max (Fig. 14): the non-local solution from MESH02 is unable to catch the real sample behavior, whereas MESH01 and MESH00 show a stable solution during the softening phase. In Fig. 14 a contour of the non-local equivalent strain distribution at the 82% of the analysis is additionally reported.

4.4. Splitting tensile strength test A concrete splitting test, normally used to assess the tensile strength of the material, is here numerically reproduced so to confirm the soundness of the proposed non-local procedure. The test consists in a cylindrical concrete sample, with a diameter D ¼ 100 mm and a length l ¼ 2D, which splits along the vertical diameter, so as to indirectly obtain its tensile strength f t ; details of the specimen geometry are given in Fig. 15 a), while an example of the adopted mesh in the numerical analyses is depicted in Fig. 15 b).Particularly, different mesh discretizations are adopted for the models (see details in Table 3) and the analyses are performed by activating or not the non-local algorithm. The same material characteristics as those of the previous example are assumed; in this case the characteristic length is taken equal to 30 mm, in agreement with the aggregates dimension [49], while the search sphere radius R ¼ 50 mm. When regularization is not activated, the equivalent strain contours appear strongly mesh dependent, as it can be appreciated in Fig. 16 (cases ‘‘L”), comparing the results obtained for the various discretizations. Differently, this can be overcome if a non-local approach is adopted (Fig. 16 (cases ‘‘NL”)). Load vs. displacement trends are reported in Fig. 17, showing a first peak related to the first cracks triggering near the zone where the load is applied, but more importantly a uniform response when non-localization is activated.

5. Concluding remarks A non-local integral-type regularization technique has been described and successfully implemented within a 3D F.E. research code, for resolving localization issues when modelling non-linear concrete behavior accounting for non-associated plasticity and softening. Particularly, the three-invariant Menétrey-Willam model has been chosen, together with a non-associated flow rule so to control inelastic dilatancy, and innovately regularized via a non-local algorithm whose efficiency is enhanced by an ad hoc octree-type algorithm and that is based on the definition of a Search-Sphere (SS) iterative procedure depending on two parameters: the sphere radius R and the characteristic length l, the latter having the meaning of a constitutive parameter. Finite element analyses have revealed a prominent dependence of the procedure on the internal length only, being the solution substantially independent on R once it is appropriately chosen to be larger than the plastic strain band; such a choice for R improves the speedup of the technique so reducing the required computational times and supporting the efficiency of the strategy, even

demonstrated by the attained stabilization of the numerical results.

Acknowledgements Financial support from the University of Padua, Italy, in the framework of DOR 2018 Projects - DOR1844491/18 and DOR1814885/18, is gratefully acknowledged.

References [1] De Borst R, Sluys LJ, Mühlhaus HB, Pamin J. Fundamental issues in finite element analyses of localization of deformation. Eng Computat 1993;10:99–121. [2] Borja RI, Regueiro RA. Strain localization in frictional materials exhibiting displacement jumps. Comput Methods Appl Mech Engrg 2001;190:2555–80. [3] Bazˇant ZP, Jirásek M. Nonlocal integral formulations of plasticity and damage: survey of progress. J Eng Mech 2002;128(11):1119–49. [4] Cervera M, Chiumenti M. Size effect and localization in J2 plasticity. Int J Solids Struct 2009;46:3301–12. [5] Tejchman J, Bobin´ski J. Continuous and Discontinuous Modelling of Fracture in Concrete Using FEM. Berlin-Heidelberg: Springer; 2013. [6] Shi J, Ling D, Hu C, Gong S. An improved nonlocal strain-softening model. Electron J Geotech Eng 2017;22:1447–58. [7] De Borst R. Simulation of strain localization: a greenreappraisal of the Cosserat continuum. Eng Computat 1991;8:317–32. [8] di Prisco C, Imposimato S, Aifantis EC. A visco-plastic constitutive model for granular soils modified according to non-local and gradient approaches. Int J Numer Anal Meth Geomech 2002;26:121–38. [9] Jirásek M, Rolshoven S. Comparison of integral-type nonlocal plasticity models for strain-softening materials. Int J Eng Sci 2003;41(13–14):1553–602. [10] Mühlhaus HB, Aifantis EC. A variational principle for gradient plasticity. Int J Solids Struct 1991;1:417–38. [11] Needleman A. Material rate dependence and mesh sensitivity on localisation problems. Comput Methods Appl Mech Engng 1988;67:69–86. [12] Oka F, Adachi T, Yashima A. Instability of an elasto-viscoplastic constitutive model for clay and strain localization. Mech Mater 1994;18:119–29. [13] Sluys L.J., ‘‘Wave Propagation, Localization and Dispersion in Softening Solids, Ph.D. Dissertation Delft University of Technology: Delft, Netherlands, (1992). [14] Schrefler BA, Sanavia L, Majorana CE. A multiphase medium model for localization and post localization simulation in geomaterials. Mech CohesiveFrictional Mater 1996;1:95–114. [15] De Borst R, Mülhaus HB. Gradient-dependent plasticity: formulation and algorithmic aspects. Int J Numer Meth Eng 1992;35(3):521–39. [16] Borino G, Fuschi P, Polizzotto C. A thermodynamic approach to nonlocal plasticity and related variational principles. J Appl Mech 1999;66(4):952–63. [17] Fleck NA, Hutchinson JW. A reformulation of strain gradient plasticity. J Mech Phys Solids 2001;49(10):2245–71. [18] Engelen RAB, Geers MGD, Baaijens FPT. Nonlocal implicit gradient-enhanced elasto-plasticity for the modelling of softening behaviour. Int J Plast 2003;19:403–33. [19] 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. [20] 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. [21] Marotti de Sciarra F. Nonlocal and gradient rate plasticity. Int J Solids Struct 2004;41(26):7329–49. [22] Marotti de Sciarra F. Novel variational formulations for nonlocal plasticity. Int J Plast 2009;25(2):302–31. [23] Polizzotto C. A nonlocal strain gradient plasticity theory for finite deformations. Int J Plast 2009;25:1280–300. [24] Poh LH, Swaddiwudhipong S. Gradient-enhanced softening material models. Int J Plast 2009;25:2094–121. [25] Ehlers W, Graf T, Ammann M. Deformation and localization analysis of partially saturated soil. Comput. Methods Appl. Mech. Eng. 2004;193:2885–910. [26] Loret B, Prevost JH. Dynamic strain localization in fluid-saturated porous media. J Eng Mech 1991;117(4):907–22. [27] Bardet JP, Proubet JA. A numerical investigation of the structure of persistent shear bands. Géotechnique 1991;41(4):599–613. [28] Alsaleh MI, Voyiadjis GZ, Alshibli KA. Modelling strain localization in granular materials using micropolar theory: mathematical formulations. Int J Numer Anal Meth Geomech 2006;30:1501–24. [29] Xotta G, Beizaee S, Willam KJ. Bifurcation investigations of coupled damageplasticity models for concrete materials. Comput Methods Appl Mech Eng 2016;298:428–52. [30] Loret B, Prevost JH. Dynamic strain localization in elasto-(visco-)plastic solids, Part 1. General formulation and one-dimensional examples. Comput Method Appl Mech Eng 1990;83(3):247–73.

G. Mazzucco et al. / Computers and Structures 224 (2019) 106120 [31] Wang WM, Sluys LJ, de Borst R. Viscoplasticity for instabilities due to strain softening and strain-rate softening. Int J Numer Meth Eng 1997;40 (20):3839–64. [32] Carosio A, Willam K, Etse G. On the consistency of viscoplastic formulations. Int J Solids Struct 2000;37:7349–69. [33] Lazari M, Sanavia L, Schrefler BA. Local and non-local elasto-viscoplasticity in strain localization analysis of multiphase geomaterials. Int J Numer Anal Meth Geomech 2015;39:1570–92. [34] Bazˇant Z, Feng-Bao L. Non-local yield limit degradation. Int J Numer Meth Eng 1988;26(8):1805–23. [35] Grassl P, Jirasek M. Nonlocal plastic models for cohesive-frictional materials. In: Herrmann J, Ehlers W, Vermeer PA, Ramm E, editors, Modelling of cohesivefrictional materials: proceedings of second international symposium on continuous and discontinuous modelling of cohesive-frictional materials (CDM 2004), 2004. pp. 323–337. [36] Marotti de Sciarra F. A general theory for nonlocal softening plasticity of integral-type. Int J Plast 2008;24:1411–39. [37] Lu X, Bardet J-P, Huang M. Numerical solutions of strain localization with nonlocal softening plasticity. Comput Methods Appl Mech Engrg 2009;198:3702–11. [38] Brinkgreve R.B.J. ”Geomaterial models and numerical analysis of softening”, Ph.D. thesis, Delft Univ. of Technology, (1994). [39] Strömberg L, Ristinmaa M. FE-formulation of a nonlocal plasticity theory. Comput Methods Appl Mech Engrg 1996;136:127–44. [40] Jirásek M. Objective modeling of strain localization. Revue francaise de genie civil 2002;6(6):1119–32. [41] Rolshoven S. Nonlocal plasticity models for localized failure. EPFL: Lausanne; 2003.

13

[42] Bobin´ski J, Tejchman J. Numerical simulations of localization of deformation in quasi-brittle materials within non-local softening plasticity. Comput Concr 2004:1(4). [43] Kenawy M, Kunnath S, Kolwankar S, Kanvinde A. Fiber-based nonlocal formulation for simulating softening in reinforced concrete beam-columns. J Struct Eng 2018;144(12):04018217. [44] Menétrey P, Willam KJ. Triaxial failure criterion for concrete and its generalization. ACI Struct J 1995;92(3):311–8. [45] Mazzucco G, Xotta G, Pomaro B, Salomoni VA, Faleschini F. Elastoplasticdamaged meso-scale modelling of concrete with recycled aggregates. Compos B: Eng 2018;140:145–56. [46] Willam K, Warnke EP. Constitutive Models for the Triaxial Behavior of Concrete. Proceedings of the International Assoc. for Bridge and Structural Engineering 1975;19:1–30. ˇ ervenka J, Papanikolaou VK. Three dimensional combined fracture-plastic [47] C material model for concrete. Int J Plast 2008;24(12):2192–220. [48] Grassl P, Karin Ln, Kent G. Concrete in compression: a plasticity theory with a novel hardening law. Int J Solids Struct 2002;39(20):5205–23. [49] Pijaudier-Cabot G, Bazˇant Z. Nonlocal damage theory. J Eng Mech 1987;113 (10):1512–33. [50] Zeng LF, Horrigmoe G, Andersen R. Numerical Implementation of constitutive integration for rate-independent elastoplasticity. Comput Mech 1996;18:387–96. [51] Schäfer BC, Quigley SF, Chan AH. Acceleration of the discrete element method (DEM) on a reconfigurable co-processor. Comput Struct 2004;82(20– 21):1707–18.