A computational framework for constitutive modelling

Computers and Structures 187 (2017) 1–23

Contents lists available at ScienceDirect

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

A computational framework for constitutive modelling Lapo Gori ⇑, Samuel Silva Penna, Roque Luiz da Silva Pitangueira Department of Structural Engineering, Engineering School, Federal University of Minas Gerais (UFMG), Avenida Antônio Carlos, 6627, Pampulha, 31270-901 Belo Horizonte, MG, Brazil

a r t i c l e

i n f o

Article history: Received 23 September 2016 Accepted 19 January 2017

Keywords: Finite element method Object-oriented programming Elastic-degradation Elasto-plasticity Continuum damage mechanics

a b s t r a c t The field of computational constitutive modelling for engineering applications is an active research tread in academia. New advanced models and formulations are constantly proposed. However, when dealing with implementation aspects, often the main concern is to provide a minimal environment to show a certain model and its applications, with implementations made from scratch. Though advanced, usually such implementations lack of generality and are well-suited for a certain numerical method while not compatible with other ones, making it difficult to reuse the code in other contexts. The ObjectOriented Paradigm (OOP) to programming have been widely applied in the last years for the realization of numerous academic numerical simulation softwares, due to its fundamental properties of abstraction, inheritance and polymorphism that allow the creation of programs well-suited for an easy collaboration between developers with expertise in different fields of engineering and mechanics. As showed in this paper, the same properties can be effectively extended also to the constitutive aspects of a model. The application of the OOP to the constitutive modelling of a wide range of materials of engineering interest is investigated, aiming to the creation of a computational framework for constitutive models that is fully independent on the other components of a code and easy to expand. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction The last years have seen a wide spread of academic softwares based on the Object-Oriented Paradigm (OOP) devoted to the numerical simulation of different physical problems. The wellknown concepts of abstraction, inheritance and polymorphism of the OOP allow to achieve generality, reusability, extendibility and ease of maintenance of a code. These properties are of fundamental importance for academic softwares that, besides to be user-friendly, must also guarantee an easy collaboration of experts in different fields. The advantages of the OOP in the Finite Element Method (FEM) have been pointed out in a number of works (see, e.g., the works by Mackie [1,2] or the book by Nikishkov [3]). Among the numerous academic finite element codes based on the OOP, the first that have been proposed are the softwares SIC [4], CHARLY [5], CASTEM [6], OBJECT NAP [7], while more recent applications are represented by the softwares FEMOBJ [8], OOFEM [9] and KRATOS [10], for example. Besides finite element applications, the OOP has been successfully used to represent also different numerical methods, such as the Boundary Element Method (BEM) [11,12], the Generalized/Extended Finite Element Method (G/XFEM) [13,14], and ⇑ Corresponding author. E-mail address: [email protected] (L. Gori). http://dx.doi.org/10.1016/j.compstruc.2017.01.012 0045-7949/Ó 2017 Elsevier Ltd. All rights reserved.

Mesh-Free Methods (MFM) [15,16]. A recent comprehensive OOP code is the open-source software INSANE1 (INteractive Structural ANalysis Environment) [17], able to represent in a common framework different numerical methods like FEM, BEM [18], G/XFEM [19], and MFM [20], as well as different continuum descriptions, like classic and generalized ones [21,22]. A critical issue of most codes devoted to the analysis of physically non-linear problems is the effective implementation of constitutive models, as well as the treatment of non-linear systems solution. Constitutive models are usually represented inside a software in a matricial/vectorial form (the Voigt representation) using array objects, and that representation, in general, strongly depends both on the kind of analysis model of the problem (i.e., threedimensional, plane-stress, etc.) and on the numerical method. This dependence poses a serious limit to the generality and extendibility of a code since, in this case, an operation like the introduction of a new constitutive model requires the modification of parts of the code that are not strictly related to the constitutive modelling. A first attempt of generalization of the constitutive aspects in an object-oriented environment can be found in the work of Jeremic´ and Yang [23]. Taking advantage of the abstraction and inheritance mechanisms, the authors define a template class able to sythesize 1 The development code is available at the Git repository URL http://git.insane.dees. ufmg.br/insane/insane.git.

2

L. Gori et al. / Computers and Structures 187 (2017) 1–23

the main aspects of elasto-plastic materials. Peculiar constitutive models are then introduced in terms of inherited classes that define specific yield criteria, flow rules and hardening-softening rules. Such generalization is made possible by the use of a peculiar C++ library that allows to handle tensorial objects in a computational environment [24,25], easing the representation of solid mechanics tensorial equations. The concept of tensors representation in an object-oriented code can be found also in the work of Weller [26]. Again, the presence of tensor objects allows to represent continuum mechanics equations in a computational environment with a syntax that is as close as possible to their mathematical representation. Within this approach, the author proposes a C++ library suitable for the solution of fluid dynamic problems with the finite volume method. Regarding the solution of non-linear systems, Menétrey [27] showed that also the non-linear analysis procedures can be efficiently integrated into an object-oriented approach. These aspects have been further investigated by Dubois-Pèlerin [28] and Lages [29]. In the cited works the authors show that the modularization of linear and non-linear equations solving offered by the OOP allows an easy implementation of different equation solvers and solution control modes, and facilitates their combination in numerical applications. In this paper, a novel approach to computational constitutive modelling based on the OOP, that makes use of concepts previously discussed by Penna [30], is investigated. The initial approach proposed by Jeremic´ and Yang [23] for elasto-plastic materials is generalized in order to accomodate also elastic-degrading media, based on both single and multiple loading functions. The theoretical basis is represented by a unified formulation for constitutive models [31–37,30], i.e., a theoretical resource able to represent a large number of elasto-plastic and elastic-degrading models in a common framework, in terms of constitutive operators, loading functions and degradation/flow rules. From a theoretical point of view, the generality of the representation is guaranteed by the tensorial formalism of the involved constitutive equations, as illustrated by Rizzi and Carol [38]. Following concepts analogous to the ones discussed by Jeremic´ and Sture [24] and by Weller [26] for the handling of tensor objects in a computational environment, the constitutive models based on the unified formulation are introduced in the code directly with their tensor-based format rather than with the matricial/vectorial expressions of the Voigt representation, i.e., with a syntax that is as close as possible to their original mathematical representation. In this way the aforementioned unified formulation can be reproduced in the code taking full advantage of the object-oriented paradigm; the abstract general equations of the formulation are introduced in their tensorial format, and the peculiar constitutive models are derived within this scheme using the concept of inheritance. The proposed tensorbased approach to constitutive modelling results in a general framework for elasto-plastic and elastic-degrading media that, with respect to the aforementioned papers and to existent object-oriented codes, is able to guarantee a complete independence between the constitutive aspects and the other parts of the code, with specific emphasis on analysis models and numerical methods. In the first part of the paper, the basic aspects of the unified formulation for constitutive models are briefly recalled, and the derivation of known constitutive models within this formulation is showed. The second part focuses on the computational aspects of the proposed framework and on its implementation in the INSANE. The main aspects of the software are briefly presented, while the organization of the constitutive models framework, together with the relations that exist between the framework itself and the other parts of the code, is presented in details. Also the handling of tensor objects and the approach for non-linear system solving are discussed. Finally, numerical simulations are performed in order to emphasize the advantages of the proposed approach.

1.1. Notations Some standard notations used in the body of the paper are summarized here. With the symbols E; E and RN , the environment space (a three-dimensional Euclidean space), its associated vector space and the N-dimensional Euclidean space are indicated. The symbols ð
ei Þ and ð
ra Þ indicate, respectively, a basis of E and a basis of RN . Latin indexes are assumed to run from 1 to 3, while greek indexes
with assume values from 1 to N. Vectors are indicated as
x 2 E,
x ¼ xi
ei , while second-order and fourth-order tensors respectively
 E,
with x ¼ xij
ei 
ej , and by x
E
E
 E,
with ^2E by x 2 E ^ ¼ xijk‘
ei 
ej 
ek 
e‘ . The tensors of third-, fifth- and sixth-order x used in this paper are of mixed type, and are represented respec
 E,
~ 2 RN  E ~ ¼ xaij
r a 
ei 
tively by x with x ej , by N



 
x 2 R  E  E  E  E, with x ¼ xaijk‘ ra 
ei 
ej 
ek 
e‘ , and by
E
 RN  E
 E,
 ¼ xaijbk‘
ra 
 ¼ RN  E with x ei 
ej 
rb  x
ek 
e‘ The symbol  denotes both the standard dot product between vectors and the total contraction between tensors like, ^  y ¼ xijk‘ yk‘

¼ xi yi ; x  y
¼ xij yj
for example,
xy ei , x ei 
ej and the other possible combinations. The same symbol is used, with a slight abuse of notation, also for contractions with mixed order ~ ¼ xij yaij
r a , since there is no risk of confusion tensor, like x  y between the different indexes. With the symbol , the standard
¼ xi yj
ei 
ej or x  y ¼ xij yk‘
ei  tensorial product, as
xy
ej 
ek 
e‘ , is indicated. In case of mixed tensors, combinations ~y ~ ¼ xaij ybk‘
ra 
ei 
ej 
rb 
ek 
e‘ . are given, for example, by x In some applications the Voigt notation will be used to represent second-order and fourth-order tensors; once a certain coordinates system has been fixed, a generic second-order tensor with dimension three x can be represented by means of an array with nine components, indicated with the symbol fxg. In an analogus way, ^ can be represented a fourth-order tensor with dimension three x ^ . It should be noted that by means of a 9  9 matrix, indicated as ½x the provided dimensions refer to a general three-dimensional case; in different situations (e.g., plane-strain or plane-stress states), the size of arrays and matrices in Voigt representation is minor. The same symbols fg and ½ are also used in Section 3 to indicate, respectively, arrays and matrices in numerical equations. 2. Unified formulation for constitutive models As discussed in the introduction, the theoretical basis for the proposed computational framework is represented by a unified formulation for constitutive models that has been developed in the last years by a number of authors (see, e.g., [31–37,30]). The core of the formulation is constituted by a set of general equations that provide a common representation for different material behaviours; from such a general representation, specific constitutive models are defined in terms of peculiar constitutive operators, loading functions, and degradation/flow rules. In the following, a brief resume of the formulation presented in the cited papers is provided, with specific emphasis on its tensorial expression; the derivation within this formulation of some known constitutive models is also discussed. In a geometrically linear context, an elastic-degrading medium is characterized by the following total stress-strain relations

r ¼ E^ S  e;

1

e ¼ ðE^ S Þ  r

ð1Þ

^ S is the secant constitutive operator. These equations correwhere E spond to the assumption of an unloading-reloading process where the stiffness remains equal to the current secant one; in this case, a full unload leads to zero permanent strains. The time derivatives of the expressions in Fig. 1 result in the equations

3

L. Gori et al. / Computers and Structures 187 (2017) 1–23

r_ ¼ E^ S  e_ þ E^_ S  e ¼ r_ e þ r_ d ; 1

e_ ¼ ðE^ S Þ  r_ þ ðE^_ S Þ  r ¼ e_ e þ e_ d 1

ð2Þ

where the superscripts e and d indicate, respectively, the elastic and degrading parts of a quantity. The previous expressions can be resumed in the tangent relations 1

r_ ¼ E^ t  e_ ;

e_ ¼ ðE^ t Þ  r_

ð3Þ

^ t has the following general where the tangent constitutive operator E representation

^t ¼ E ^ S  z1  ðx ~y ~Þ E ðESijk‘

¼

1

 ðzab Þ

xbij

yak‘ Þ
ei 
ej 
ek 
e‘

N X

~ k_ b mðbÞ ¼
k_  m;

1 ^_ S Þ ¼ ðE

b¼1

N X

^ ðbÞ ¼
k_  M  k_ b M

ð5Þ

b¼1

each one represented by a fourth-order tensor. It can be easily   r holds. The different phases of ~ ¼M shown that the relation m the loading process are described in terms of N loading functions f a ðr; pÞ, that can be collected in a vector
f 2 RN , each one depending on the stress state r and on a generic set of internal variables p, that in a thermodynamical context are referred to as thermodynamical forces. As the inelastic strain rate, also the thermodynamical forces rate can be expressed in terms of a degradation rule as N X ðbÞ k_ b h

ð6Þ

b¼1

where the terms k_ b are inelastic multipliers, and the terms h are the directions of degradation of the set p. Within this approach, the operators appearing in Eq. (4) are represented by ðbÞ

^S  m ~ ¼ ðESijpq mbpq Þ
r b 
ei 
ej ~ :¼ E x

ð7Þ

^ S ¼ ðnav z ES Þ
r a 
ek 
e‘ ~ E ~ :¼ n y v zk‘ ^ S  mÞ ~  ðE ~ ¼ ðHab þ namn ES z :¼ H þ n

ð8Þ

mnrs

mbrs Þ
ra 
r b

ð9Þ

~ and H are the gradients of the loading functions, where the terms n defined by

~ :¼ n

@f a
r a 
ev 
ez ; @ rv z

H :¼ 

@f a ðbÞ @f h
r a 
r b ¼  a
r a 
r b @p @kb ð10Þ

The strain-based approach, on the contrary, relies on the additive decomposition for the stress rate r_ ¼ r_ e þ r_ d . In this case, the inelastic rate r_ d and the rate of the secant constitutive operator are defined in terms of the degradation rules

r_ d ¼

N X ~ ; k_ b mðbÞ ¼
k_  m b¼1

N X ðbÞ k_ b h

X ^ ðbÞ ¼
k_  M ^_ S ¼  k_ b M E N

ð11Þ

b¼1

where
k_ 2 RN is again a vector containing the N inelastic multipliers,  ~ is the third-order tensor containing the N directions of degradam

ð12Þ

b¼1

ð4Þ

~ is where
k_ 2 RN is the vector containing the N inelastic multipliers, m the third-order tensor containing the N directions of degradation mðbÞ of the strain degrading rate, each one represented as a second-order  is the fifth-order tensor containing the N directions of tensor, and M ^ ðbÞ of the inverse of the secant constitutive operator, degradation M

p_ ¼

ator, each one represented by a fourth-order tensor. It can be easily    e holds. As for the stress-based ~ ¼M shown that the relation m approach, a set of N loading functions f a ðe; p Þ describes the different phases of the loading process. In this case, each loading function depends on the strain state e and on a generic set of internal variables p . As the inelastic stress rate, also the thermodynamical forces rate p_  can be expressed in terms of a degradation rule as

p_  ¼

which specific expression depends on the chosen approach. The stress-based approach relies on the additive decomposition for the strain rate e_ ¼ e_ e þ e_ d . Within this approach, the inelastic rate e_ d and the rate of the inverse of the secant constitutive operator are defined in terms of the degradation rules

e_ d ¼

tion mðbÞ of the stress degrading rate, each one represented as a   is the fifth-order tensor containing second-order tensor, and M ^ ðbÞ of the secant constitutive operthe N directions of degradation M

ðbÞ where the terms k_ b are inelastic multipliers, and the terms h are the directions of degradation for the set p . Within this approach, the operators appearing in Eq. (4) are represented by

~  ¼ mbij
r b 
ei 
ej ~ :¼ m x  ~ ¼ nak‘
ra 
ek 
e‘ ~ :¼ n y

ð13Þ ð14Þ

z :¼ H ¼ Hab
ra 
r b

ð15Þ

~  and H are the gradients of the loading funcwhere the terms n tions, defined by

~  :¼ n

@f a
r a 
ek 
e‘ ; @ ek‘

H :¼ 

@f a ðbÞ @f
r a 
r b ¼  a
r a 
r b h @ p @kb ð16Þ

Within this general framework for elastic-degradation, damage models can be obtained once the secant material properties are assumed to depend on a reduced set of parameters, the damage ^ 0 Þ1 ; DÞ for the stress-based approach, and ^ S Þ1 ððE variables, as ðE

^ 0 ; D Þ for the strain-based approach. The symbols D and D , ^ S ðE E representing the sets of damage variables for the stress- and strain-based approaches, may indicate, in general, scalar, vectorial or tensorial damage variables [32,39]. For a damage model, the rate of the secant constitutive operator can be expressed as

^ S 1 1 ^_ S Þ ¼ @ðE Þ  D; _ ðE @D

^S ^_ S ¼ @ E  D _ E @D

ð17Þ

where with the symbol , a contraction operation compatible with the peculiar nature of the damage variables sets is indicated. The rates of the damage variables sets can be expressed in terms of the following degradation rules

_ ¼ D

N X k_ b MðbÞ ¼
k_  M; b¼1

_ ¼ D

N X k_ b MðbÞ ¼
k_  M

ð18Þ

b¼1

where
k_ 2 RN is the vector of inelastic multipliers while M and M contain, respectively, the N directions of degradation MðbÞ and MðbÞ of the sets of damage variables D and D . The directions of degradation of the damage variables sets and the ones of the secant constitutive operator are linked by the relations

^ S 1  ¼ @ðE Þ  M; M @D

^S   ¼ @ E  M M @D

ð19Þ

2.1. Some constitutive models It should be noted that the formulation described in the previous section represents a generalization of different constitutive models. Starting from the general expression for the tangent constitutive operator (Eq. (4)) different stress- or strain-based constitutive models can be obtained once the expressions of the secant

4

L. Gori et al. / Computers and Structures 187 (2017) 1–23

constitutive operator, of the set of loading functions, and of the directions of degradation are fixed. Some of them are briefly recalled here. Elasto-plastic models (see, e.g., the books of Lubliner [40] and Simo and Hughes [41]), for example, are usually derived within the stress-based approach, replacing the secant constitutive operator with the initial elastic one. In a mono-potential formulation, and considering the associativity in the strain-space, they are expressed by the tangent operator



   ^0  n  n  E ^0 E ^t ¼ E ^0    E ^0  n Hþn E

pffiffiffiffiffiffiffi 3J 2  ry ;

J2 ¼

1 dev dev r r ; 2

3 n ¼ pffiffiffiffiffiffiffi 2 3J 2

rdev ð21Þ

where J2 is the second invariant of the stress tensor, and ry the equivalent yield stress. Scalar damage models, on the contrary, are usually presented within the strain-based mono-potential formulation, and relies on a secant constitutive operator expressed as

^S ^_ S ðDÞ ¼ @ E D_ ¼ E ^ 0 D_ E @D

^ S ðD; E ^ 0 Þ ¼ ð1  DÞ E ^0; E

ð22Þ

where D is a scalar damage variable, assumed to vary from 0 ^0 (undamaged material) to 1 (completely damaged material), and E

^  ¼ E ^0; M

M ¼ 1;

^ 0  e ¼ r0 m ¼ E

ð23Þ

For a scalar damage model, a common choice for the loading function is represented by the following additive decomposition

f ðe; DÞ ¼ eeq ðeÞ  KðDÞ 6 0

ð24Þ

eeq

1  a þ aebðeeq K 0 Þ

eeq

8 pffiffiffiffiffiffiffiffi ee > > > qffiffiffiffiffiffiffiffiffi > > 0 > > > < q2w ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ 2w0 =E0 > > > > ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r > hP i > > 2 3 > ðkÞ : k¼1 ð< e >þ Þ



ð26Þ

ðMazars-LemaitreÞ ðSimo-JuÞ ð27Þ

ðLemaitre-ChabocheÞ ðMazarsÞ

^ 0  eÞ is the internal energy, E0 the initial Young’s where 2w0 ¼ e  ðE modulus, eðkÞ the k-th eigenvalue of the strain tensor, and < eðkÞ >þ ¼ ðeðkÞ þ jeðkÞ jÞ=2 its positive part. A generalization of the concept of scalar damage is represented, for example, by orthotropic damage models. In this case, the set of damage variables is represented in terms of an integrity tensor i, which expression in a system ðe1 ; e2 ; e3 Þ, alligned with the principal axes of the strain tensor e, is

0

i1

B i ¼ @0

is the initial constitutive operator. In this case, the operators described in the previous section are expressed as

_ k_ ¼ D;

K0 

where K 0 is a threshold value for the equivalent deformation, representing the onset of damage, and where a and b are parameters that define, respectively, the maximum allowed damage level and the damage evolution intensity. Different scalar damage models can be obtained choosing a peculiar equivalent deformation; the classic damage models of Mazars-Lemaitre [42], Simo-Ju [43], LemaitreChaboche [44] and Mazars [45,46] are defined by the equivalent deformations

ð20Þ

The classic Von Mises model (also known as J2-plasticity) is obtained with the loading function

f ðrÞ ¼

Dðeeq Þ ¼ 1 

0

0 i2 0

0

1

0

1  D1

C B 0 A ¼ @0 0 i3

0

0

1  D2

0

0

1  D3

1 C A

ð28Þ

where the ij are the integrity variables, assumed to vary from 1 (undamaged material) to 0 (completely damaged material), and the Dj are the damage variables, with opposite behaviour. A common representation for the secant consitutive operator in terms of the three integrity variables is ESijkl ¼ aijmn E0mnrs aklrs , with aijkl ¼ ðiik djl þ ijl dik Þ=2, resulting in

where eeq ðeÞ is a function depending only on the strain tensor, usually indicated as equivalent deformation, that represents the loading

condition of the continuum, while KðDÞ is an historical parameter that depends only on the damage variable and that is representative of the maximum level of deformation reached during the loading process. Such a loading function results in the tangent constitutive operator

where

  ^ t ¼ ð1  DÞ E ^ 0  1 r0  n ; E H  1 @KðDÞ @Dðeeq Þ H ¼ ¼ @D @ eeq

The degradation process is assumed to be ruled by three uncoupled loading functions, expressed as

n ¼

@ eeq ; @e ð25Þ

where DðKÞ is a prescribed evolution law for the damage like, for example, the following exponential damage law

A :¼

E0 ð1  m0 Þ ; ð1 þ m0 Þð1  2m0 Þ

f 1 ðe1 ; i1 Þ ¼ e1  Kði1 Þ; f 3 ðe3 ; i3 Þ ¼ e3  Kði3 Þ

B :¼

E0 m 0 ; ð1 þ m0 Þð1  2m0 Þ

C :¼

E0 2ð1 þ m0 Þ ð30Þ

f 2 ðe2 ; i2 Þ ¼ e2  Kði2 Þ; ð31Þ

where the ej are the eigenvalues of the strain tensor e, and K is, as in the case of a scalar damage model, an historical variable, depending

5

L. Gori et al. / Computers and Structures 187 (2017) 1–23

on the integrity values. Each Kðij Þ is representative of the maximum value reached by the correspondent eigenvalue ej during the loading process. The derivatives of the loading functions with respect to the strains are then expressed by

0

1 0 0

1

0

B C nð1Þ ¼ @ 0 0 0 A;

0 0 0

1

B C nð2Þ ¼ @ 0 1 0 A;

0 0 0

0

0 0 0

1

B C nð3Þ ¼ @ 0 0 0 A

0 0 0

0 0 1 ð32Þ

while, following the same considerations adopted for the scalar damage models, the operator H is given by

0 @i1

@ e1

B H ¼ B @0 0

0

11

0

@i2 @ e2

C 0 C A

0

@i3 @ e3

ð33Þ

where the ii ðei Þ are prescribed evolution laws for the integrity variables, that depend on the specific constitutive model. Recalling Eq. (18) and considering that the damage in each direction is assumed to be independent on the other directions, it follows that k_ b  i_b and that

0

ð1Þ

M

1

1 0 0 B C ¼ @ 0 0 0 A; 0 0 0

0

ð2Þ

M

1

0 0 0 B C ¼ @ 0 1 0 A; 0 0 0

0

ð3Þ

M

1

^ M

^S @E ¼ ; @ib

ðbÞ

m

^S @E ¼ e @ib

½K S  fUg ¼ fRg

ð36Þ

S

where ½K  is the global secant stiffness matrix of the system, fUg the nodal parameters vector, and fRg the vector of nodal dual parameters. In the presence of physical non-linearities, the solution of the problem consists in the search for a set of equilibrium configurations, parametrized in a quasi-static context by the pseudo-time. The procedure is based on an incremental-iterative method for which the pseudo-time is replaced by a finite discretization, by means of a set of increments (or steps); inside each step an iterative procedure is performed (predictor-corrector method). The incremental equilibrium equation at the iteration n of the step k is represented by k

½K t n1 fDUgkn ¼ Dkkn fPg þ fQgkn1

ð37Þ

0 0 0 B C ¼ @0 0 0A 0 0 1

where ½K t n1 is the tangent stiffness matrix at the iteration n  1 of the step k, depending in general on the current values of the field

ð34Þ

the iteration n of the step k; Dkkn is the increment of the load multiplier at the iteration n of the step k; fPg the vector of the nodal ref-

It can then easily shown that the directions of degradation of the secant operator and the ones of the inelastic stress rate are expressed as ðbÞ

parameters vector. It is worth to note that no reference has been made to a specific numerical method; indeed, each approximation function could be derived from the finite element method or from one of the mesh-free methods, for example. Replacing the aforementioned discretization in a weak form for static problems, the following matricial system follows

ð35Þ

3. Computational aspects The present section aims to illustrate the implementation aspects of the theoretical framework for constitutive models discussed in the previous section. As it will be shown, the use of an object-oriented approach results in a computational framework for constitutive models with an high level of generalization. The main guidelines that has been followed in the implementation were the independence of the constitutive model framework on the other components of the software like, for example, the analysis model (e.g., three-dimensional, plane-stress, plane-strain, etc.) and the numerical method (e.g., FEM, MFM, BEM, G/XFEM, etc.), and its modularity. The necessity of this independence should be viewed as follows: in case of any modification of the constitutive models framework, such as the introduction of new constitutive models, the presence of a strong link between the framework and the other components results in the need of modifications of parts of the software that are not strictly related to the constitutive aspects, in order to accomodate the new resources. Also the modularity of the framework aims to ease the task of implementation of new constitutive models; indeed, new models can take advantage of the general structure of the framework and can be defined by means of the inheritance mechanism. 3.1. Discrete models A discrete model is obtained once the state variable of the con
, is replaced at each tinuum problem, i.e., the displacement field u
ð
point
x by a finite approximation u xÞ ’ ½UfUg, where the matrix ½U contains the approximation functions, and fUg is the nodal

k

variables, fDUgkn is the vector of incremental nodal parameters at

erence loads, and fQ gkn1 the vector of the residual forces at the iteration n  1 of the step k, given by fQ gkn1 ¼ kkn1 fPg  fFgkn1 , where fFgkn1 is the vector of nodal forces equivalent to internal stresses. The stiffness matrix ½K t  and the vector fFg are assembled considering the contributions of each ‘‘element” (as it will be showed in the following, here the term element is used in a broad sense, with no restriction to FEM elements) composing the discrete model

Z ½K t el ¼

Xel

Z fFgel ¼

Xel

^ t ½B dX ½BT ½E

ð38Þ

½BT frg dV

ð39Þ

where, depending on the considered problem,

R Xel

may represent a

line, surface or volume integral over the ‘‘element”. The matrix ½B is the kinematical operator that allows to evaluate the strain vector feg at a certain point
x of the domain using the values of the nodal parameters that fall in the neighborhood of
x, while the stress vec^ S  as tor can be evaluated using the secant constitutive operator ½E ^ S feg (distinguished by the tangent one ½E ^ t  appearing in frg ¼ ½E Eq. (39)), or with a proper return-mapping algorithm. 3.2. INteractive Structural ANalysis Environment The INSANE system [17] is an open source software for computational mechanics. It is based on the Java language, and relies on the object-oriented paradigm. This, results in a robust software composed by a set of classes that interacts between them, and that are able to represent the different aspects of a numerical simulation with an high level of generalization. The high abstraction of the parts of the software eases the maintenance of the code, and its expansion. The software is composed by a set of interactive graphical applications allowing the operations of pre- and post-processing, and by a numerical core, responsible for the analysis of discrete models. As

6

L. Gori et al. / Computers and Structures 187 (2017) 1–23

Fig. 1. Numerical core.

showed in Fig. 1, with an UML diagram, the numerical core is composed by the interfaces Assembler, Model and Persistence, and by the abstract class Solution. The task of the Assembler (Fig. 2) is to mount the matricial system of the model, as the one of Eq. (36), returning the stiffness matrix and its partitions related to free and restrained degrees of freedom, and the vector of nodal dual parameters. Depending on the kind of approximation considered in the non-linear analysis, the stiffness matrix can be elastic (getC()), secant (getTotalC ()) and tangent (getIncrementalC()).2 In returning the vector fRg (Eq. (36)), the main task of the assembler is to mount the vector of nodal forces equivalent to internal stresses (getFp()). For problems modelled with the FEM, such interface is implemented by the class FemAssembler, that is extended by further classes in order to account for additional methods like MFM, BEM and GFEM (Fig. 2). The abstract class Solution provides a set of methods devoted to the solution of the matricial system expressed in Eq. (36). Different inherited classes allow the solution of linear and non-linear problems, both for static and dynamic analyses. The interface Model allows the representation of a generic discrete model in the numerical core of the software. It is composed by several lists of objects, each one representing a peculiar component of a discrete model, like nodes, elements, type of problem, type of analysis model and materials, for example. The interface Persistence is responsible for the processing of input and output data, which are persisted as XML (eXtensible Markup Language) files. Such class also guarantees the consistency of the data between the different components of the discrete model, applying the Observer-Observable design pattern. When an alteration in the state of an observed object (i.e., an object that extends the class java.util.observable) occours, the propagation mechanism of the alteration is triggered, and the observers (i.e., the objects implementing the interface java.util.observer) are notified to update themselves. In the INSANE system, as showed in Fig. 1, the role of observer is assumed by the interface Persistence, while the observed objects are the abstract classes Model and Solution. In the following sections, the structure of the proposed approach for constitutive modelling environment is widely discussed. Before present the framework for constitutive models, details on other parts of the code that are directly related to it are provided, with specific attention to the solution of static nonlinear problems and to the handling of tensorial quantities.

2 It is worth to note that in the described methods the symbol C is used to indicate the stiffness matrix, instead of K as in Eq. (36). This is justified by the fact that for a more general problem than the quasi-static one adopted here, Eq. (36) should account € þ ½B fUg _ þ ½C fUg ¼ fRg, where fUg _ and fUg € for additional terms, as ½A fUg represent, respectively, the first and second time derivatives of the nodal parameters. However, in the analytical expressions, the symbol K is maintained since it is traditionally adopted in the field of mechanics.

3.2.1. Non-linear problems solving The module Solution is devoted to the solution of both linear and non-linear systems of equations generated by each of the numerical methods embedded in the software. Following a standard methodology [47], non-linear systems in the form of Eq. (37) are analyzed decomposing additively the vector of incremental nodal parameters into two components as k

k

k

k

fDUgkn ¼ Dkkn fDU P gn þ fDU Q gn , such that ½K t n1 fDU P gn ¼ fPg and k

k

½K t n1 fDU Q gn ¼ fQ gkn1 . The different phases of the solution procedure adopted in the software are resumed by the pseudo-code of Fig. 3 [48]. As pointed out in the work of Lages et al. [29], the use of an object-oriented approach for the treatment of non-linear problems solution allows to obtain a certain modularity of the code, that makes possible to embed different solution methods in a same library, and also ease the expansion of the code. Despite the advantages discussed in the cited paper, it is observed that the proposed C++ library (NLS++: NonLinearSolver) seems to have a strong connection with the adopted numerical method (FEM, in that case) and with the constitutive aspects of the problem, hence preventing to take full advantage of the OOP in terms of independence of the library on the other components of the software. With respect to the library described by Lages et al. [29], the module Solution discussed in this section aims to generalize the concept of modularity of the implementation and to overcome the main drawbacks of the aforementioned library, especially in terms of independence on the other parts of the software; these aspects allow to use the module Solution to solve linear and non-linear problem that derives from the different numerical methods and constitutive models embedded in the software. This module is composed by the abstract class Solution and by the interfaces Step and IterativeStrategy. The class Solution (Fig. 4) is responsible for the representation inside the software of different solution procedures that can be adopted depending on the considered problem. Attention is here devoted to the class StaticEquilibriumPath. This class (Fig. 5) is responsible for the solution of static nonlinear problems. The method execute() implemented in such class triggers the solution algorithm represented in Fig. 3. Its attribute step, an instance of the interface Step, contains all the methods that are necessary to the execution of an incremental step of the non-linear analysis; the class StandardNewtonRaphson, detailed in Fig. 6, manages the solution with the Newton-Raphson method. The method execute() of this class is responsible for the activation of the iterative process inside each step. The attribute iterativeStrategy, an instance of the interface IterativeStrategy, specifies the control method adopted in the solution, and is responsible for the iterative predictorcorrector strategy defined in the algorithm of Fig. 3. Among the control methods implemented there are the load control, displacement control [47], different arc-length controls [49–53], and the generalized displacement control [48].

7

L. Gori et al. / Computers and Structures 187 (2017) 1–23

Fig. 2. Interface Assembler.

Fig. 3. Solution algorithm [48].

non-linear problems. Furthermore, it is emphasized that this module exhibits no direct connection with a peculiar numerical method or constitutive model, since the only contact with the other parts of the software is through the methods getIncrementalCuu() and getFp() of the interface Assembler. This aspect allows the use of this module for the solution of different problems independently on the peculiar numerical methods and constitutive models adopted in the analysis.

Fig. 4. Class Solution.

The implementation of the solution process for non-linear problems described in Fig. 3 is resumed by the sequence diagram of Fig. 7. The class Solution (i.e., the inherited class StaticEquilibriumPath in case of a static non-linear problem) initiates a loop over the steps of the incremental process while the method execute() triggers the iterative process for each step. The tangent stiffness matrix is provided3 by the Assembler, an instance of which is contained in the class StandardNewtonRaphson (Fig. 6). Once the tangent stiffness matrix has been obtained, the method getXPandXQ() allows to evaluate the incremental values k

k

fDU P gn and fDU Q gn of the state variable; the IterativeStrategy is then solicited to return the predicted value of the load factor increment at the first iteration (getPredictor()), or its correction for the other iterations (getCorrector()). The updating of both the load factor and the nodal parameters vector is performed by the method assignStepState(. . .), that receives as an input the evaluated load factor increment. Finally, the vector of residual forces is calculated using the vector of nodal forces equivalent to internal stresses, and the convergence is checked. As it can be observed in the class and sequence diagrams exposed in this section, the module Solution implemented in the software is highly modular, allowing an easy combination between different incremental methods and iterative strategies for the solution of various kind of static, dynamic, linear and

3.2.2. Tensors and arrays data structures As discussed in the introduction, the handling of matrices, vectors and tensors is an important aspect for the efficiency of a simulation code, and it can take advantage of the OOP [24]. The representation of tensor objects is also of fundamental importance for the effective implementation of constitutive models [23]. Interesting applications of tensors objects can be found in the work of Weller et al. [26] (which, however, doesn’t provide details about the implementation) and in the work of Jeremi and Sture [24]. In the latter, a C++ library for algebraic operations with matrices, vectors, and tensors is discussed, together with its practical application to solid mechanics problems. These three objects are represented within a common hierarchy structure; the classes devoted to matrices (matrix) and tensors (tensor) representation, inherit from the superclasses nDarray and nDarray_rep, while the class vector inherits from matrix. The Java library adopted in the present work differs from the aforementioned one, starting from the inheritance structure. In this case, matrices, vectors and tensors are represented, respectively, by the classes IMatrix, IVector, and Tensor, with no inheritance relations between them. The formers allow to represent, for example, the tangent stiffness matrix and the vector of nodal parameters appearing in Eq. (37), or the Voigt expressions of the constitutive operator and stress tensor appearing in Eq. (39) and (38), that in a plane-stress state are given by

0

rxx

rxy 3

This procedure, and the one for the evaluation of the vector of nodal forces equivalent to internal stresses, will be discussed in details in Section 3.2.3.

1

C ^ 0  feg!B frg ¼ ½E @ ryy A ¼

0 1 E0 B @ m0 1  m20 0

m0 0 1

0

0

1  m0

10

exx

1

CB C A@ eyy A

exy ð40Þ

8

L. Gori et al. / Computers and Structures 187 (2017) 1–23

Fig. 5. Class StaticEquilibriumPath.

Fig. 6. Class StandardNewtonRaphson.

The class Tensor instead, is used to represent tensor objects in the parts of the code devoted to constitutive modelling. As it will be shown, a conversion of an object between these classes is often required during the analysis; this issue will be discussed in a separate section (Section 3.2.5). The class IVector (Fig. 8(a)) possesses two attributes; an integer defining its size, and an array containing the vector components in double precision. In an analogous way, the class IMatrix (Fig. 8(b)) is characterized by three attributes: two integers defining the number of rows and columns, and a twodimensional array containing the matrix components in double

Fig. 7. Non-linear solution sequence diagram.

9

L. Gori et al. / Computers and Structures 187 (2017) 1–23

Fig. 8. Classes IVector and IMatrix.

one (Tensor(Tensor tensor)), an empy three-dimensional tensor with the specified order (Tensor(int order)), a zerothorder tensor with the specified value (Tensor(double t0)), and a tensor with specified order and dimensions, with the same dimension for all the indexes (Tensor(int order, int dim)) or with different dimensions for each index (Tensor(int order, int[] dim)). A tensor is internally stored with the attribute tensor, an instance of the class IMatrix; for example, a third-order tensor ~t with indexes dimension equal to three, is represented as the matrix Fig. 9. Class Tensor.

precision. For both the IVector and IMatrix classes, the numbering of the elements begins with zero. Some of the methods of these two classes are exposed in Fig. 8; there, are represented different constructors for their initialization, as well as different methods to get and set their components. Despite they are not explicitly exposed here, both the classes posses proper methods for the definition of standard algebraic operations like as sum and multiplication. Peculiar attention is devoted to the implementation of the class Tensor. As discussed in the cited papers, the definition of tensor objects in a code is of fundamental importance for the representation of constitutive equations that, for practical reasons of implementation simplicity and code clearness, should be expressed in a form that is as close as possible to their mathematical representation. The class Tensor (Fig. 9) possesses three attributes, an integer defining the order of the tensor, an array of integers that indicates the dimension of each index of the tensor, and an instance of the class IMatrix, that stores the elements of the tensor. Common tensors in solids mechanics, like the second-order stress ten^ have their indexes sor r or the fourth-order constitutive tensor E, running from 1 to 3; the implemented class however, allows the representation of generic tensors, with indexes dimension that may differ from 3. This characteristic is important,4 for example, for the representation of the third-order tensor ~ ¼ mbij
r b 
ei 
ej , appearing in Eq. (5), that contains the N direcm tions of degradation mðbÞ of the strain degrading rate in the stressbased representation (Section 2); while the indexes i and j run from 1 to 3, the dimension of the index b depends on the number of adopted directions of degradation. Such class is then structured to ideally accomodate tensors of the desired order and dimensions; however, at the current stage of implementation, only tensors of order from 0 to 6 are supported. In Fig. 9, the constructors of the class Tensor are also represented; they allow to initialize a tensor as a copy of an existing

4 Despite this aspect is not discussed in this paper, the possibility to use tensors with generic dimensions allows also to effectively handle continuum theories different from the Cauchy one as like, for example, micropolar continuum models [22].

In this way, the operations between tensors are internally performed using the operations already defined for the IMatrix class. The weak point of such operations is that they don’t match the mathematical syntax of the tensorial equations; hence, they must be filtered with properly defined method. First of all, the tensorial indexes numbering starts from 1 instead of 0, as in the Java arrays. As illustrated in Table 1 (for tensors with order from 0 to 3), proper getting and setting methods are defined to access the elements of a tensor, each one depending on the order of the tensor. These methods are able to map the indexes of the tensorial components to the ones of their matricial representation, in such a way as to make possible to use the IMatrix access methods. Taking advantage of the internal representation by means of an instance of the class IMatrix, the operations of sum and subtraction between two tensors of the same order and dimensions (Table 2), are simply reduced to the sum or subtraction of their representative matrices. Also the multiplication with a scalar quantity simply reduces to the scaling of the representative matrix with the informed scalar value (Table 2); the secant constitutive ^ S ¼ ð1  DÞ E ^ 0 appearing in Eq. (25), can then be operator E expressed in the code as secantOperator = elasticTensor. scale(1 - d).

Table 1 Access methods for tensor components. Get method

Set method

Order

getValue() getValue(int i) getValue(int i, int j)

setValue(double value) setValue(int i, double value) setValue(int i, int j, double value) setValue(int n, int i, int j, double value)

0 1 2

getValue(int n, int i, int j)

3

10

L. Gori et al. / Computers and Structures 187 (2017) 1–23

Table 2 Sum, subtraction and scaling. Method

Table 4 Tensorial product.

Description Return the sum of t1 and t2

add(Tensor t1, Tensor t2) add2(Tensor t2)

Sets the value of the current tensor to the sum of itself and t2 Return the sum of the current tensor and t2 Return the difference of t1 and t2

add3(Tensor t2) sub(Tensor t1, Tensor t2) sub2(Tensor t2)

Sets the value of the current tensor to the difference of itself and t2 Return the difference of the current tensor and t2 Returns the scalar multiplication of the factor s with the current tensor

sub3(Tensor t2) scale(double s)

Method

Mathematical operation

ijVECkl(Tensor b) mijVECnkl(Tensor b)

^ e‘ ei 
ej 
ek 
c ¼ a  b ¼ aij bk‘
~¼a b
 ~b ei 
ej 
en 
ek 
e‘ c¼a mij nk‘ em 


of development, only two methods of this kind have been implemented, since the formulation resumed in Section 2 just requires the product between second-order tensors and the product between third order tensors (Table 4). Using the methods exposed in this section, the general equation for the tangent constitutive operator (Eq. (4)), Etijk‘ ¼ ESijk‘  ðzab Þ1 xbij yak‘ , can be expressed in the code as secan tOperator.sub3(zInv.mnDOTmijnkl(x.mijVECnkl(y))).

Table 3 Tensors contraction. Method

Mathematical operation

iDOTi(Tensor b) ijDOTj(Tensor b) ijDOTijk(Tensor b) ijDOTijkl(Tensor b) mnDOTmijnkl(Tensor b)


¼a b
b c¼a i i
¼a b

c ¼ a  b ij j ei ~¼a b

c ¼ a  b ij ijk ek ^¼a b ek 
c ¼ab e‘ ij ijk‘
¼a b ^ c¼ab ej 
ek 
e‘ ei 
mn mijnk‘


The other fundamental operations that have to be considered by the class Tensor are the contraction between tensors, and the tensorial product. Regarding the contraction operation, a family of methods have been implemented, in such a way as to allow the application of the contraction operation to all the possible combinations of tensors of different order; some of the implemented methods, that perform the contraction of the current tensor with the informed one, are represented in Table 3. As showed in the code block of Fig. 10 with the operation
¼ a b
e , the implemented contraction doesn’t rely
c ¼ a  b ij j i directly on the algebraic operations of the IMatrix class, but rather on the access methods defined in the Tensor class (Table 1). It could be argued that the family of contraction methods implemented in the code (Table 3) should be made private, and accessed, for example, by a general public method reproducing the contraction operator , able to automatically select the proper private method depending on the order of the tensors involved in the operation. Despite this could lead to a cleaner code, the approach where the choice of the proper contraction method is leaved to the user is preferred, since it results in a more readable representation of the tensorial equations. In this way, the code exactly matches the index notation of the tensorial equations. Regarding the other fundamental operation, the tensorial product, the same approach have been followed. At the current state

Fig. 10. Tensor contraction ci ¼ aij bj .

3.2.3. Assembly of the tangent stiffness matrix and of the internal forces vector The procedure allowing to mount the tangent stiffness matrix and the vector of internal forces, used in Eq. (37), is resumed in the sequence diagram of Fig. 11. As showed in the diagram, the entire process is triggered by the class Solution; as already discussed for the diagram of Fig. 7, the class Solution initiates a loop over the steps of the incremental process, while its method execute() triggers the iterative process for each step. For each iteration n and each step k, the abstract class Assembler is solicited to k

provide the tangent stiffness matrix ½K t n1 (getIncrementalCuu ()). Such class, and the inherited classes that represent the different numerical methods, contains the main operations that allow to assemble the global stiffness matrix, making use of the contributions of the ‘‘elements” composing the discrete model (Eq. (39)). As discussed before, here the term element is used in a wide sense, since it may refer to the finite elements of a finite element model, to the quadrature cells of a mesh-free method, or to the boundary elements and internal cells of a boundary element model. These objects are represented by peculiar classes in the system, each one inherited from the same general abstract class (Element). Each ‘‘element” is then required to provide its tangent stiffness matrix; such object is evaluated by the numerical integration of Eq. (39), that can be expressed schematically as the following sum over the integration points

½K t el ¼

  X ^ t ½B J g W ip ½BT ½E ip ip ip

ip

ð41Þ

where W ip is a weighting factor, Jip the Jacobian representing the integration domain transformation, and g ip a certain geometrical property evaluated at the integration point. The abstract class ProblemDriver (called by the method getIncrementalC(Element e)) manages the integration procedure, that may vary depending on the peculiar problem analysed. For example, a physically non-linear problem requires the assembly of a material tangent stiffness, a problem with geometrical non-linearities need for a geometrical tangent stiffness, while a problem with nonconservative loads (the so-called follower loads) is characterized by a load stiffness; each one of these situations is represented by a proper class that extends the superclass ProblemDriver. A basic sketch of the integration procedure for a physically non-linear problem is represented is Fig. 12. In order to perform the numerical integration illustrated in Fig. 12, the ProblemDriver solicitates each integration point, represented by an instance of the abstract class Degeneration (Fig. 13), to return the value of the tangent constitutive operator, through the method mountCt().

11

L. Gori et al. / Computers and Structures 187 (2017) 1–23

Fig. 11. Tangent stiffness matrix and internal forces vector assembly sequence diagram.

expressed in Eq. (41). To do this, it is asked to mount the constitutive operator (initial, by the method mountC(), secant, by the method mountCs() or tangent, by the method mountCt()) evaluated at the integration point, which value is solicitated to an instance of the class ConstitutiveModel, illustrated in the following section. As showed in Fig. 11, an analogous procedure is repeated in order to mount the vector fFgkn of nodal forces equivalent to internal stresses, for each step k and iteration n. As for the tangent stiffness matrix, also this object is evaluated considering the contribution of the list of ‘‘elements” composing the discrete model, each one calculating numerically the integral of Eq. (38) as Fig. 12. Integration algorithm for physically non-linear problems.

fFgel ¼

  X W ip ½BT frg J ip g ip ip

The abstract class Degeneration (similar to the class GaussPoint originally defined by Mentrey and Zimmermann [27]) represents the degeneration of a certain solid model into a point. Each object of the class Degeneration is composed by a list of material points (instances of the class MaterialPoint, that represent, for example, a certain number of points discretizing a beam cross section) and by a Representation, that represents the integration point itself (with coordinates and weighting factor W ip ). Among the different kinds of degenerations, it is pointed out the presence of the class PrescribedDegeneration, which geometrical characteristics are explicitly specified (e.g., the thickness for a plane problem), and the class CrossSection, which characteristics depend on the set of material points discretizing the section of a beam. The object Degeneration is responsible for returning to the ProblemDriver its contribution to the numerical integration

ip

ð42Þ

It is worth to note that the operations resumed in the sequence diagrams of Figs. 7 and 11 don’t make use of the tensor objects described in Section 3.2.2, since all the involved objects are defined in terms of instances of the classes IVector and IMatrix. Tensor objects will be widely used in the next section, devoted to the proposed constitutive models framework. 3.2.4. Constitutive models framework The sequence diagram of Fig. 11 ends with two calls to the abstract class ConstitutiveModel with the methods mountCt () and mountDualInternalVariableVector() (abbreviated in Fig. 11 as mountDualIntVarVec()), that request such class to provide, respectively, the tangent constitutive operator and the internal stresses. The abstract class ConstitutiveModel encapsulates the different activities of a constitutive model; once all the necessary informations, like kind of analysis and material

12

L. Gori et al. / Computers and Structures 187 (2017) 1–23

Fig. 13. Abstract class Degeneration.

Fig. 14. Abstract class ConstitutiveModel.

properties, for example, are passed to the constitutive model, this is able to perform all the operations related to the evaluation of the constitutive operator and of the internal stresses. The general structure of the class is resumed in Fig. 14. The inherited class UnifiedConstitutiveModel is devoted to the representation of the constitutive models belonging to the unified formulation discussed in Section 2. More precisely, the inherited class UCMSingleLoadingFunction is devoted to the representation of the equations of the mono-potential approach (see, e.g., Eq. (20) for elasto-plastic models or Eq. (25) for scalar damage models), while the class UCMMultipleLoadingFunction contains the equations of the multi-potential formulation (see the general expression for the tangent operator of Eq. (4)). These classes contain the methods necessary to assemble the constitutive operator ^ t  (mountCt(), for the tangent one) and the vector of matrix ½E the internal stresses frg (mountDualInternalVariableVec tor()), at the considered integration point. As showed in Fig. 14, the class UnifiedConstitutiveModel contains an attribute of the class UnifiedConstitutiveModel Filter, the superclass of a set of classes called filters. Each filter is the representation of a peculiar constitutive model; the methods of the class UnifiedConstitutiveModel call directly a set of operations implemented in each filter, that allow to mount the ^ t and the vector of internal stresses r. While constitutive operator E

define the loading functions and their gradients (getLoadingFunctionPotential() and getHardeningSofteningPotential()), the directions of degradation (getInelasticPotential()), as well as the expressions of the secant constitutive operators (getSecantTensor()), as can be seen in Fig. 15. The procedures performed by the class UCMSingleLoadingFunction to mount the tangent constitutive operator are detailed in the code of Fig. 16; in order to illustrate the code, the tangent constitutive operator of Eq. (4) is considered here in its version for scalar damage models (Eq. (25)),  t 0   ^ ^ E ¼ ð1  DÞ E þ 1=H ðm  n Þ. The method mountCt() receives as inputs an instance of the class AnalysisModel (Section 3.2.5), an instance of the class Material (Section 3.2.6), and a map containing the values of the constitutive variables at the considered integration point. As already stated, the elements of the aforementioned equations are evaluated calling the methods of a proper instance of the class UnifiedConstitutiveModelFilter. The secant con^ 0 is obtained with the method stitutive operator ð1  DÞ E

getSecantTensor(. . .), the damage law derivative 1=H with ge tHardeningSofteningPotential(. . .), the loading function gradient n with getLoadingFunctionPotential(. . .), and finally the direction of degradation m with getInelasticPotential (. . .). At line 7, the tensorial product m  n is performed, using

the classes UCMSingleLoadingFunction and UCMMultipleLoadingFunction contains the main equations for the evaluation of the constitutive operator, each filter5 is responsible to

5 In this sense, the filters present some analogies with the Template3Dep class described by Jeremi and Yang [23], an object devoted to the representation of yield directions, flow directions and evolution laws. However, the approach of the cited authors is limited to elasto-plastic models, based on a single loading function. Furthermore, no details are provided on the relation of such class with the other parts of the code and on the eventual character of independence of the proposed template material models on, for example, different numerical methods and different analysis models.

Fig. 15. Class UnifiedConstitutiveModelFilter.

L. Gori et al. / Computers and Structures 187 (2017) 1–23

13

Fig. 16. Method mountCt() of UCMSingleLoadingFunction.

Fig. 17. Method mountDualInternalVariableVector() of UnifiedConstitutiveModel.

the method ijVECkl(Tensor b) of the class Tensor, while the whole tangent operator for the scalar damage model is assembled at line 19. At this point it is worth to note that the tangent constitutive operator (the variable tangentOperator) is evaluated internally as an instance of the class Tensor, while the method mountCt() returns to the class Degeneration the matricial representation of such tensor (i.e., its Voigt expression). This transformation is performed at line 20, with the method reduceToMatrix (Tensor tensor) of the class AnalysisModel (Section 3.2.5). The vector of internal stresses solicited by the class Degeneration is provided by the method mountDualInternalVariable Vector() implemented by the class UnifiedConstitutiveModel (Fig. 17). Such method takes as inputs a vector con-

taining the strain components, an instance of the class AnalysisModel, an instance of the class Material, and two maps containing the current and previous values of the constitutive variables at the considered integration point. The method, at line 3, solicitates the filter to update its constitutive variables (such as, for example, the strain tensor, the equivalent deformation, etc.); then, it evaluates the internal stresses, calling a proper returnmapping algorithm for elasto-plastic models, or using the secant constitutive operator for damage models. A partial representation of the filters classes organization is illustrated in Fig. 18. As for the class UnifiedConstitutiveModel, there are two families of filters, one devoted to the representation of mono-potential models (SLFConstitu-

14

L. Gori et al. / Computers and Structures 187 (2017) 1–23

Fig. 18. Class organization of the filters.

tiveModelFilter), like scalar damage models, and the other to the representation of multi-potential models (MLFConstitutiveModelFilter), like orthotropic damage models. The class SLFConstitutiveModelFilter is inherited by the class Elasto PlasticConstitutiveModelFilter, devoted to the representation of elasto-plastic models like the ones expressed by Eq. (20), and by the class IsotropicConstitutiveModelFilter, devoted instead to the representation of the scalar damage models defined by Eq. (25). Among the various mono-potential models contained in the framework, the ones recalled in Fig. 18 are the Von Mises elasto-plastic model [54], and the scalar damage models of Mazars-Lemaitre [42], Simo-Ju [43], Ju [55], Lemaitre-Chaboche [44], de Vree [56] and Mazars [45]. The classes MLFBiDissipati veConstitutiveModelFilter and OrthotropicConstituti veModelFilter inherit from the class MLFConstitutiveModelFilter; the former allows to represent bi-dissipative models like, for example, elasto-plasticity coupled with damage or scalar damage models with different damage variables for tension and compression (see, e.g., [57]), while the latter is devoted to the representation of orthotropic damage models and smeared crack models (see, e.g., [58]). In order to show the use of such filters, attention is focused on the classes IsotropicConstitutiveModelFilter and SLFICMLemaitreChaboche, where the latter is devoted to the representation of the scalar damage model of Lemaitre and Chaboche [44]. Among the methods collected in the class diagram of Fig. 15, the class IsotropicConstitutiveModelFilter implements the one devoted to the evaluation of the secant constitutive ^ S ¼ ð1  DÞ E ^ 0 (getSecantTensor(. . .), Fig. 19), that is operator E common to all the derived scalar damage models. A first use of the class Material can be observed; at line 2, such class is required to provide the material moduli that will be used by the proper analysis model to mount the elastic constitutive operator ^0. E The further methods of Fig. 15 are implemented by the class SLFICMLemaitreChaboche (Fig. 20). The method getInelas-

ticPotential(. . .) returns the direction of degradation ^ 0  e as in Eq. (23), getLoadingFunctionPotential m ¼ E (. . .) returns the gradient of the loading function n ¼ @ eeq =@ e as in Eq. (25), that for the Lemaitre-Chaboche model is expressed by n ¼ 1=ðeeq E0 Þ, and getHardeningSofteningPotential(. . .) returns the derivative of the damage law 1=H ¼ @Dðeeq Þ=@ eeq as in Eq. (25). As showed in Fig. 16, these quantities are used by the class UCMSingleLoadingFunction to mount the tangent constitutive operator. The class SLFICMLemaitreChaboche, as well as the other classes that inherit from IsotropicConstitutiveModelFil ter, contain further methods, that are used during the analysis. Among them there are the method getEquivalentStrain(. . .) (see Fig. 21), devoted to the evaluation of the equivalent strain qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi measure eeq (eeq ¼ 2w0 =E0 , for the Lemaitre-Chaboche model, as in Eq. (27)), and the method updateConstitutivesVariables (. . .) (Fig. 22), called in Fig. 17, that solicitates the filter to update its constitutive variables and to persist them. As showed in Fig. 22, in order to update the constitutive variables of the filter, the method unLoadLaw(. . .) (Fig. 23) is called. Such method is able to identify the different phases of the loading processes, making use of conditions on the loading function.

3.2.5. The AnalysisModel class The abstract class AnalysisModel is devoted to the representation inside the code of the different problems that can be analized, i.e., three-dimensional, plane-stress, plane-strain, etc. Its inherited classes characterize each problem in terms of number of degrees of freedom, kind and number of internal variables (i.e., strains) and dual internal variables (i.e., stresses). In this context, it is emphasized only its role in relation to the constitutive models framework. As already discussed in the previous sections, a difference exists in the representation of objects inside the proposed constitutive models framework and outside of it. The constitutive relations

Fig. 19. Method getSecantTensor(. . .) of IsotropicConstitutiveModelFilter.

15

L. Gori et al. / Computers and Structures 187 (2017) 1–23

Fig. 20. Methods of the class SLFICMLemaitreChaboche.

Fig. 21. Method getEquivalentStrain(. . .) of SLFICMLemaitreChaboche.

are expressed in the framework by means of tensor objects, with the advantages in terms of modularity and generalization discussed in the previous section. Outside the framework, as in conventional FEM codes, the matricial form of representation is recovered, for operations like numerical integration and assembly of stiffness matrices. The exchange of informations between the framework, that works with tensorial quantities, and the other parts of the code, that work with matricial objects, is guaranteed by the classes that inherit from AnalysisModel. Such classes indeed, are able to convert the tensorial representation of an operator into its Voigt expression, and vice versa. The method getActiveDualInter nalVariablesVector(. . .), for example, is able to convert the second-order stress tensor to its Voigt expression; for a planestress analysis model this results in the conversion

0

1

0

1

rxx rxx rxy 0 C B C r¼B @ ryx ryy 0 A ! frg ¼ @ ryy A rxy 0 0 0

ð43Þ

In an analogous way, the method reduceToMatrix(. . .), converts the fourth-order constitutive operator tensor into its Voigt expres-

sion; again, for a plane-stress problem, this results in a 3  3 matrix as

0 1 E B 0 ^0 ¼ ^ 0 ! ½E m0 E @ 1  m20 0

m0 0 1

0

0

1  m0

1 C A

ð44Þ

3.2.6. The Material class Since it has a strong connection with the constitutive models framework discussed in the previous sections, the abstract class Material is briefly presented here. The classes that inherit from the superclass Material are responsible for the representation of different materials inside the software; as already stated while commenting Fig. 19, one of the roles of such class is to provide the values of the material moduli, through the method getPs(), that can be used, for example, by an instance of the class AnalysisModel to mount a constitutive operator. Elasto-plastic materials and materials with damage must implement the interfaces Hardenable and Damageable, respectively. Focusing on the latter, such interface forces a class to implement the methods getDamage(. . .), that should return the damage

16

L. Gori et al. / Computers and Structures 187 (2017) 1–23

Fig. 22. Method updateConstitutivesVariables(. . .) of SLFICMLemaitreChaboche.

Fig. 23. Method unLoadLaw(. . .) of SLFICMLemaitreChaboche.

L. Gori et al. / Computers and Structures 187 (2017) 1–23

17

Fig. 24. Method getDamage(. . .) of DamageLawExponential.

level for a certain deformation, and getInelasticModulus(. . .), that should return instead the variation of the damage for a certain deformation, both of them called by the filters inside the method unLoadLaw(. . .) (Fig. 23). Such methods call analogous methods contained in the abstract class InelasticLaw, an instance of which may be contained in an class that inherits from Material. An exponential damage law like the one of Eq. (26), for example, is represented inside the code by the class DamageLawExponential, that inherits from InelasticLaw. In this case, the evaluation of the damage level with the method getDamage(. . .) is exposed in Fig. 24.

3.2.7. Comments on the independence of the constitutive models framework In the previous sections, a computational framework for constitutive models have been presented and discussed. Such framework is based on a theoretical unified formulation for constitutive models, briefly recalled in Section 2, that allows to represent different elasto-plastic and elastic-degrading media within a common set of equations. A peculiarity of the proposed framework is that it works directly with tensor objects, allowing to represent the equations of Section 2 in a way that is as close as possible to their mathematical format. As it will be commented in this section, such characteristics make the framework: highly modular and easy to expand, independent on the adopted numerical method, independent on the peculiar analysis model. Regarding the modularity of the framework, this naturally arises from the adopted theoretical formulation, combined with an object-oriented approach. As showed in the class organization of the filters (Fig. 18), for example, the unified formulation for constitutive models is well suited for an implementation in an OOP code. The common equations can be defined inside general classes, with the mechanism of abstraction, while specific constitutive models can be obtained by inheritance. Such scheme is easy to expand; if one has to introduce a new scalar damage model, for example, he has to implement just a new class that inherits from IsotropicConstitutiveModelFilter (Fig. 18), and that defines the methods getInelasticPotential(. . .), getLoadingFunctionPotential(. . .) and getHardeningSofteningPo tential(. . .) (Fig. 15) depending on the characteristics of the

new constitutive model, without requiring modifications in the superclasses. The independence of the framework on both the numerical method and on the analysis model are of fundamental importance; once a certain constitutive model has been implemented into the framework, it can be used directly with different numerical methods (such as FEM, BEM, G/XFEM and MFM), as well as with different analysis models (e.g., three-dimensional, plane-stress, planestrain, etc). Such independence is due to the tensorial representation adopted inside the framework and to the object-oriented approach followed in the other parts of the code. As showed in the sequence diagram of Fig. 11 for the evaluation of the tangent stiffness matrix of a model, the object-oriented approach allows to clearly separate the different tasks of the software. In the diagram, a specific numerical method is only defined at the assembler level, in terms of inherited classes (Fig. 2). Furthermore, the assembler is not directly linked to the constitutive models framework, since the flow of informations passes through different superclasses. The framework relates directly with the object Degeneration, a representation of an integration point, that is an entity that doesn’t depend on the adopted numerical method. The kind of analysis model influences the form of the Voigt representation ^ t  and of the tangent stiffness of the tangent constitutive operator ½E matrix of each ‘‘element” ½K t el . In this case the independence is due to the fact that the framework works with tensor objects, which representation is not influenced by a peculiar analysis model. As illustrated in Section 3.2.5, once the framework returns a constitutive operator in a tensor form, this is converted to its Voigt notation, suitable for the further operations, depending on the peculiar analysis model, by means of methods implemented by the inherited classes of AnalysisModel. The characteristics of modularity and independece on the numerical method are illustrated in the next section with the help of two different numerical simulations.

4. Numerical simulations In the previous sections, the main aspects of the proposed constitutive models framework have been discussed, and its advantages with respect to existent solutions have been also illustrated. In order to show the advantages of the proposed methodology, in the present section two different numerical simulations have been selected.

18

L. Gori et al. / Computers and Structures 187 (2017) 1–23

The first example consists in the numerical simulation of the three-point bending test of a plain concrete beam investigated by Petersson [59]. For the simulation, a finite element model is considered, in combination with different damage models. The aim of such example is to emphasize the benefits of the modularity of the framework, in the sense that a certain model, with defined mesh, loads and restraints, can be easily studied with different constitutive models with a minimum effort for the user. The purpose of the other example instead, is to emphasize the independence of the framework with respect to the adopted numerical method. With this intention, the bending test on an L-shaped plain concrete panel, investigated experimentally by Winkler et al. [60], is reproduced both with a finite element model and with an Element Free Galerkin model (EFG), using the same constitutive model (the same simulation have been performed, adopting the proposed framework, with a BEM model by Peixoto et al. [18], and with a G/XFEM model by Pinheiro et al. [61]). 4.1. Three-Point bending test The geometry of the beam is illustrated in Fig. 25(a). As in the paper of Petersson [59], the plain concrete is characterized by an elastic modulus E = 30,000 MPa, a Poisson’s modulus m ¼ 0:20, tensile and compressive uniaxial strengths ft = 3.33 MPa and fc = 33.30 MPa, a fracture energy Gf = 0.124 MPa, and a characteristic length of the material h = 40 mm. The finite element model (Fig. 25(b)) consists of 204 four-nodes quadrilateral elements in a plane-stress state, with a thickness of 50 mm. The loading process is driven by the generalized displacement control method [48], assuming a reference load P = 800 N, an initial loading factor increment of 0.02, and a tollerance for the convergence in displacement of 1  104 . In order to simulate the experimental test, five different constitutive models have been considered: four scalar damage models (from the superclass IsotropicConstitutiveModelFilter, Fig. 15) and a smeared-crack model (from the superclass Ortho

tropicConstitutiveModelFilter, Fig. 15). Regarding the scalar damage models, they are defined by the equivalent strain measures of Eq. (27), combined with an exponential damage load (Eq. (26)), which parameters are resumed in Table 5. The smeared crack model adopts the Carreira and Chu stress-strain law [62,63] for compression, and the Boone and Ingraffea stress-strain relation for tension [64]. The relation between the vertical displacement of the point where the load is applied and the load factor is presented, for the considered damage models, in Fig. 26. As can be observed, the different equilibrium paths are in good agreement with the experimental data obtained in [59], with the best results given by the simplified Mazars model both in terms of load factor peak and post-peak behaviour. The damage distribution is presented in Fig. 27 for the MazarsLemaitre and Mazars equivalent strain measures, for a state corresponding to a vertical displacement of the point of application of the load of 0.38 mm, that emphasizes that the Mazars-Lemaitre model, as well as the Lemaitre-Chaboche and Simo-Ju models, doesn’t distinguish between tensile and compressive behaviour, leading then to a damage initiation at both the sides of the beam model. Rather than the quality of the results, that obviously may be reproduced with other implementations, what is interesting to observe in the presented simulations is the easiness with which one can pass from a certain constitutive model to a different one. Once a numerical model has been defined, in terms of mesh, loads and restraints, different material behaviours may be simulated just by associating to the elements a different constitutive model, i.e., a different class from the tree illustrated in Fig. 15, without requiring further modifications of the numerical model. In Figs. 28, 29, the portion of the XML input file that refers to the definition of a certain element (element E87 in this case) for the three-point bending model is presented: the first one with the Mazars-Lemaitre model and the other with the Simo-Ju model. As it can be observed, the two inputs differ only for the child ConstitutiveModel

Fig. 25. Three-point bending test.

Table 5 Exponential damage law parameters.

a b K0

Lemaitre-Chaboche [44]

Simo-Ju [43]

Mazars-Lemaitre [42]

Mazars [46]

0.950 600

0.950 3

0.950 700

0.950 700

1:10  104

1:90  102

1:12  104

1:10  104

L. Gori et al. / Computers and Structures 187 (2017) 1–23

19

Fig. 26. Three-point bending – equilibrium paths.

Fig. 27. Damage distribution at dy = 0.38 mm.

Fig. 28. XML input file for the Mazars-Lemaitre model.

Fig. 29. XML input file for the Simo-Ju model.

(at line 5), that indicates the peculiar filter that has to be considered in the analysis; no further modifications of the input file are required. Besides the advantages for a developer in terms of ease of implementation of new constitutive models as discussed in Section 3, this example helps to illustrate also the advantages for a

user of the framework. With the proposed approach indeed, the numerical model preparation, that is usually the most time consuming task of an analysis, is optimized, in the sense that if one needs to investigate different material behaviours, he only has to made slight modifications in limited parts of the input files.

20

L. Gori et al. / Computers and Structures 187 (2017) 1–23

4.2. L-shaped panel The geometry of the simulated panel is illustrated in Fig. 30(a). As in the paper by Winkler et al. [60], the plain concrete is characterized by an elastic modulus E = 25,850 MPa, a Poisson’s modulus m ¼ 0:18, tensile and compressive uniaxial strengths ft = 2.70 MPa and fc = 31.00 MPa, a fracture energy Gf = 0.065 MPa, and a characteristic length of the material h = 28 mm. The panel is discretized with two different numerical models: a finite element model (Fig. 30(b)) consisting of 192 four-nodes quadrilateral elements, and a an Element Free Galerkin model (Fig. 30(c)) with the same nodal distribution as in the finite element model. Both the models are characterized by a plane-stress state, with a thickness of 100 mm. The loading process is driven by the generalized displacement control method [48], assuming a reference load q = 28 N/mm, loading factor increments of 0.02, and a tollerance for the convergence in displacement of 1  104 . Both the FEM and the EFG numerical models have been investigated using the same smeared crack constitutive model as adopted for the three-point bending test. The results of the analysis are presented, for both the numerical models, in Fig. 31, where the relation between the vertical displacement of the top left point and the load factor is presented (point A in Fig. 30(a)).

In the same figure, the experimental data obtained by Winkler et al. [60] are also reported. As it can be observed, both the simulations show results that are in good agreement with the experimental values, regarding the load factor peak and the post-peak behaviour. Both the models present an initial stiffness higher than the one observed in the experiment; however, this issue is common to all the simulations of such a test that can be found in the literature, and it can be observed also in the original paper of Winkler et al. [60]. Again, besides the quality of the results, it is interesting to discuss the advantages offered by the proposed methodology. This example indeed, shows directly the property of independence of the constitutive model with respect to the adopted numerical method. As widely discussed in Section 3, the peculiar theoretical formulation adopted by the framework, combined with an efficient object-oriented approach, allows to obtain a constitutive models library that relates to a numerical model only through the object Degeneration (Fig. 13), the representation of an integration point inside the code, a concept common to all the numerical methods implemented in the code. This makes possible to use a class of the tree of Fig. 15 independently with each one of the numerical methods that can be implemented in the code, with no need of specific adaptations. As can be observed in Figs. 32 and 33, both an element of a FEM model or a quadrature cell of an EFG model,

Fig. 30. L-shaped panel.

Fig. 31. L-shaped panel – equilibrium paths.

L. Gori et al. / Computers and Structures 187 (2017) 1–23

21

Fig. 32. XML input file for FEM model.

Fig. 33. XML input file for EFG model.

may share the same child ConstitutiveModel (at line 4), i.e., both the numerical model are able to use the same class describing, in this case, a smeared crack constitutive model. 5. Conclusions The present paper originated from the observation that, with a few exceptions commented in the introduction, current implementations of constitutive models present important issues. While new and more advanced constitutive models are constantly proposed, their computational aspects are often left aside, limiting them to the implementation of minimum working examples. This results in implementations that are often difficult to expand, in the sense that a different constitutive model usually need to be implemented from scratch, without the possibility to reuse existent code. Furthermore, such implementations are well-suited for a certain numerical method (e.g., the FEM), while their compatibility with other methods is not guaranteed. Another important consideration is that constitutive models are ofted defined in a code in their Voigt representation, that is strongly dependent on the peculiar analysis model that is being considered. In response to the aforementioned issues, and inspired in earlier works on the same theme, the present paper proposed a novel approach to computational constitutive modelling, strongly based on the object-oriented paradigm. As commented with more details in Section 3.2.7, the result is a computational framework for constitutive models that is: highly modular and easy to expand, independent on the adopted numerical method, independent on the peculiar analysis model. The first point has been made possible by the combination of the OOP with a unified formulation for constitutive models, widely investigated in the past, that allows the representation of a large family of models for engineering materials within a common set of equations. In this combination, the general equations that are shared by different models are defined in terms abstract objects, while specific models are obtained through the mechanism of inheritance. This should guarantee an easy expansion of the framework with different constitutive models. Also the second point is mainly due to the use of the OOP. In this case the independence

of the framework on the numerical method is achieved with an effective comunication between their respective codes; the flow of informations is transferred through the concept of integration point (Fig. 11), that is independent on both the framework and on the peculiar numerical method. In this way, the implementation of constitutive aspects should not concern on how the treatment of different numerical methods is structured inside the code. The first two points are also made possible by the use of tensor objects inside the framework. However, this aspect plays a more evident role in the third point, the independence on the analysis model. As it has been already discussed, current implementations usually define constitutive models directly in their Voigt notation, that is strongly dependent on a peculiar analysis model. This considerably limits the efficiency of the implementation, since the same constitutive model must be defined with different implementations, one for each analysis model; furthermore, the modification of a constitutive model involves the modification of a large part of code. The use of tensor objects instead, allows to express constitutive equations in a computational environment with a syntax that is as close as possible to their original mathematical representation. Beside the fact that such a representation eases the work of a developer, that can exactly reproduce the constitutive equations in the code as their are expressed mathematically, it also guarantees that the implemented constitutives models don’t depend on a peculiar analysis model, since a general tensorial expression of the equations is used in place of an analysis model-dependent Voigt expression. The main aspects of the implementation of the proposed constitutive models framework in a numerical simulation software have been investigated and presented with the aid of class diagrams and code blocks. Peculiar attention has been also devoted to the identification of the main procedures involved in the simulation process and to their representation in terms of sequence diagrams, aiming to point out the advantages offered by the proposed approach. Also the other parts of the code that directly relates with the framework have been discussed, with peculiar attention on the solution of non-linear problems, the handling of tensor objects and treatment of different analysis models. The latters have been shown to play a fundamental role for the achievement of the aforementioned results. In order to emphasize the results of the proposed approach to computational constitutive modelling, two different numerical simulations, performed with the software INSANE, have been

22

L. Gori et al. / Computers and Structures 187 (2017) 1–23

presented. The first one showed the capability of the framework to allow the use of a wide range of constitutive models for different simulations on a common discrete model, without requiring heavy user interventions, while the other one showed its compatibility with different numerical methods, i.e., the possibility to use a certain constitutive model belonging to the framework with any of the implemented numerical methods. As a final comment it can be stated that the proposed approach for computational constitutive models, as made clear along the paper, can virtually embedd unlimited constitutive models, as long as they can be included in the unified formulation that constitutes the theoretical basis of the framework; furthermore, such constitutive models can be independently applied to different analysis models and investigated with different numerical methods. Acknowledgements The first author (CNPq Scholarship – Brazil) and the other authors gratefully acknowledge the support of the brazilian research agencies FAPEMIG (in portuguese: Fundação de Amparo á Pesquisa do Estado de Minas Gerais – grant PMM-00669-15), CNPq (in portuguese: Conselho Nacional de Desenvolvimento Científico e Tecnológico – grant 308785/2014-2). The authors wish also to acknowledge the anonymous reviewers, which helpful comments led to a considerable improvement of the content of this paper. References [1] Mackie RI. Object-oriented methods and finite element analysis. United Kingdom: Saxe-Coburg Publications; 2002. [2] Mackie RI. Advantages of object-oriented finite-element analysis. Proc Inst Civil Eng – Eng Comput Mech 2009;162(1):23–9. http://dx.doi.org/10.1680/ eacm.2009.162.1.23. [3] Nikishkov GP. Programming finite elements in Java. United Kingdom: Springer-Verlag London; 2010. [4] Touzot G. S.i.c.1.1: Réflexion sur larchitecture des logiciels de modélisation. Tech. rep., Université de Technologie de Compiègne; 1986 [in French]. [5] De Saxcé G. Le projet charly: um logiciel de calcul par éléments finis et éléments frontières de seconde génération. In: Séminaire de génie logiciel, Division MSM, Université de Liège; 1987 [in French]. [6] Verpeaux P, Charras T, Millard A. Castem 2000: une approche moderne du calcul des structures. In: Fouet JM, Ladevèze P, Ohayon R, editors. Calcul des structures er intelligence artificielle, Pluralis; 1988 [in French]. [7] Forde BWR, Foschi RO, Stiemer SF. Object-oriented finite element analysis. Comput Struct 1990;34(3):355–74. http://dx.doi.org/10.1016/0045-7949(90) 90261-Y. [8] Zimmermann T, Dubois-Pèlerin Y, Bomme P. Object-oriented finite element programming: I. Governing principles. Comp Methods Appl Mech Eng 1992;98:291–303. http://dx.doi.org/10.1016/0045-7825(92)90180-R. [9] Patzák B, Bittnar Z. Design of object oriented finite element code. Adv Eng Softw 2001;32(10–11):759–67. http://dx.doi.org/10.1016/S0965-9978(01) 00027-8. [10] Dadvand P, Rossi R, Oñate E. An object-oriented environment for developing finite element codes for multi-disciplinary applications. Arch Comp Method E 2010;17(3):253–97. http://dx.doi.org/10.1007/s11831-010-9045-2. [11] Lage C. The application of object-oriented methods to boundary elements. Comp Methods Appl Mech Eng 1998;157(3):205–13. http://dx.doi.org/ 10.1016/S0045-7825(97)00235-1. [12] Wang W, Ji X, Wang Y. Object-oriented programming in boundary element methods using C++. Adv Eng Softw 1999;30(2):127–32. http://dx.doi.org/ 10.1016/S0965-9978(98)00050-7. [13] Bordas S, Nguyen PV, Dunant C, Guidoum A, Nguyen-Dang H. An extended finite element library. Int J Numer Methods Eng 2007;71(6):703–32. http://dx. doi.org/10.1002/nme.1966. [14] Chamrová R, Patzák B. Object-oriented programming and the extended finiteelement method. In: Proc. inst. civil eng. – eng. comput. mech., vol. 163-4; 2010. p. 271–8. doi:http://dx.doi.org/10.1680/eacm.2010.163.4.271. [15] Barbieri E, Meo M. A fast object-oriented matlab implementation of the reproducing kernel particle method. Comput Mech 2012;49(5):581–602. http://dx.doi.org/10.1007/s00466-011-0662-x. [16] Kanber B, Yavuz MM. Object-oriented programming in meshfree analysis of elastostatic problems. EAAS 2015;7(2):1–18. [17] INSANE Project. . [18] Peixoto R, Anacleto F, Ribeiro G, Pitangueira R, Penna S. A solution strategy for non-linear implicit BEM formulation using a unified constitutive modelling framework. Eng Anal Bound Elem 2016;64:295–310. http://dx.doi.org/ 10.1016/j.enganabound.2015.11.017.

[19] Alves PD, Barros FB, Pitangueira RLS. An object-oriented approach to the generalized finite element method. Adv Eng Softw 2013;59:1–18. http://dx. doi.org/10.1016/j.advengsoft.2013.02.001. [20] Silva RP. Análise não-linear de estruturas de concreto por meio do método Element Free Galerkin. Ph.D. thesis, UFMG Federal University of Minas Gerais, Belo Horizonte, Brazil; 2012 [in Portuguese]. [21] Gori L, Pitangueira RLS, Penna SS, Fuina JS. A generalized elasto-plastic micropolar constitutive model. Appl Mech Mater 2015;798:505–9. [22] Gori L, Penna SS, Pitangueira RLS. An enhanced tensorial formulation for elastic degradation in micropolar continua. App Math Model 2017;41:299–315. http://dx.doi.org/10.1016/j.apm.2016.08.025. [23] Jeremic´ B, Yang Z. Template elastic-plastic computations in geomechanics. Int J Numer Anal Met 2002;26(14):1407–27. http://dx.doi.org/10.1002/nag.251. [24] Jeremic´ B, Sture S. Tensor objects in finite element programming. Int J Numer Methods Eng 1998;41(1):113–26. http://dx.doi.org/10.1002/(SICI)1097-0207 (19980115)41:1<113::AID-NME277>3.0.CO;2-4. [25] Jeremic´ B, Runesson K, Sture S. Object-oriented approach to hyperelasticity. Eng Comp 1999;15(1):2–11. http://dx.doi.org/10.1007/s003660050002. [26] Weller HG, Tabor G, Jasak H, Fureby C. A tensorial approach to computational continuum mechanics using object-oriented techniques. Comput Phys 1998;12(6):620–31. http://dx.doi.org/10.1063/1.168744. [27] Menétrey P, Zimmermann T. Object-oriented non-linear finite element analysis: application to J2 plasticity. Comput Struct 1993;49(5):767–77. http://dx.doi.org/10.1016/0045-7949(93)90025-9. [28] Dubois-Pèlerin Y, Pegon P. Object-oriented programming in nonlinear finite element analysis. Comput Struct 1998;67(4):225–41. http://dx.doi.org/ 10.1016/S0045-7949(98)00054-6. [29] Lages E, Paulino G, Menezes I, Silva R. Nonlinear finite element analysis using an object-oriented philosophy – application to beam elements and to the Cosserat continuum. Eng Comp 1999;15(1):73–89. http://dx.doi.org/10.1007/ s003660050006. [30] Penna SS. Formulação multipotencial para modelos de degradação elástica: unificação teórica, proposta de novo modelo, implementação computational e modelagem de estruturas de concreto. Ph.D. thesis, UFMG – Federal University of Minas Gerais, Belo Horizonte, Brazil; 2011 [in Portuguese]. [31] Carol I, Rizzi E, Willam K. A unified theory of elastic degradation and damage based on a loading surface. Int J Solids Struct 1994;31(20):2835–65. http://dx. doi.org/10.1016/0020-7683(94)90072-8. [32] Rizzi E. Sulla localizzazione delle deformazioni in materiali e strutture. Ph.D. thesis, Dipartimento di Ingegneria Strutturale, Politecnico di Milano, Milano, Italy; 1995 [in Italian]. [33] Rizzi E, Carol I, Willam K. Localization analysis of elastic degradation with application to scalar damage. J Eng Mech-ASCE 1995;121(4):541–54. http:// dx.doi.org/10.1061/(ASCE)0733-9399(1995)121:4(541). [34] Carol I. Elastic degradation and damage: plasticity – like formulation, stiffness recovery and localization. In: Mecánica computacional, volume XVII. Number 4. Solid Mechanics (A), Tucumán Argentina; 1996. [35] Carol I, Willam K. Spurious energy dissipation/generation in stiffness recovery models for elastic degradation and damage. Int J Solids Struct 1996;33(20– 22):2939–57. http://dx.doi.org/10.1016/0020-7683(95)00254-5. [36] Carol I. New developments in elastic degradation and damage: anisotropic formulations and evolution laws in lag space, In: MECOM 99, Mendoza Argentina; 1999. [37] Hansen E, Willam K, Carol I. A two-surface anisotropic damage/plasticity model for plain concrete. In: Proceedings of fracture mechanics of concrete materials – Framcos-4 conferece, Paris, France; 2001. [38] Rizzi E, Carol I. A formulation of anisotropic elastic damage using compact tensor formalism. J Elast 2001;64(2):85–109. http://dx.doi.org/10.1023/ A:1015284701032. [39] Ju JW. Isotropic and anisotropic damage variables in continuum damage mechanics. J Eng Mech-ASCE 1990;116(12):2764–70. http://dx.doi.org/ 10.1061/(ASCE)0733-9399(1990)116:12(2764). [40] Lubliner J. Plasticity Theory. Dover Books on Engineering Series. Dover Pubblications; 2008. [41] Simo JC, Hughes TJR. Computational inelasticity. New York: Springer-Verlag; 1998. [42] Mazars J, Lemaitre J. Application of continuous damage mechanics to strain and fracture behavior of concrete. In: Shah SP, editor. Application of fracture mechanics to cementitious composites. NATO ASI series, vol. 94. Netherlands: Springer; 1984. p. 507–20. [43] Simo JC, Ju JW. Strain- and stress-based continuum damage models – I. Formulation. Int J Solids Struct 1987;23(7):821–40. http://dx.doi.org/10.1016/ 0020-7683(87)90083-7. [44] Lemaitre J, Chaboche JL. Mechanics of solid materials. Cambridge: Cambridge University Press; 1990. [45] Mazars J. Application de le mécanique de l’endommagement au comportement non lineaire et à la rupture du béton de structure. Ph.D. thesis, Université Pierre et Marie Curie – Laboratoire de Mécanique et Technologie, Paris, France; 1984 [in French]. [46] de Borst R, Gutiérrez MA. A unified framework for concrete damage and fracture models including size effects. Int J Fract 1999;95(1–4):261–77. http:// dx.doi.org/10.1023/A:1018664705895. [47] Batoz J-L, Dhatt G. Incremental displacement algorithms for nonlinear problems. Int J Numer Methods Eng 1979;14(8):1262–7. http://dx.doi.org/ 10.1002/nme.1620140811.

L. Gori et al. / Computers and Structures 187 (2017) 1–23 [48] Yang YB, Shieh MS. Solution method for nonlinear problems with multiple critical points. AIIA J 1990;28(12):2110–6. http://dx.doi.org/10.2514/3.10529. [49] Feng Y, Peric´ D, Owen D. A new criterion for determination of initial loading parameter in arc-length methods. Comput Struct 1996;58(3):479–85. http:// dx.doi.org/10.1016/0045-7949(95)00168-G. [50] Riks E. An incremental approach to the solution of snapping and buckling problems. Int J Solids Struct 1979;15(7):529–51. http://dx.doi.org/10.1016/ 0020-7683(79)90081-7. [51] Ramm E. Strategies for tracing the nonlinear response near limit points. In: Wunderlich W, Stein E, Bathe KJ, editors. Nonlinear finite element analysis in structural mechanics. Berlin, Heidelberg: Springer-Verlag; 1981. [52] Crisfield M. A fast incremental/iterative solution procedure that handles snapthrough. Comput Struct 1981;13(1):55–62. http://dx.doi.org/10.1016/00457949(81)90108-5. [53] Crisfield MA. An arc-length method including line searches and accelerations. Int J Numer Methods Eng 1983;19(9):1269–89. http://dx.doi.org/10.1002/ nme.1620190902. [54] Mises Rv. Mechanik der festen körper im plastisch-deformablen zustand. Nachr Ges Wiss Göttingen Mathematisch-Physikalische Klasse 1913;1913:582–92 [in German]. [55] Ju JW. On energy-based elastoplasticity damage theories: constitutive modeling and computational aspects. Int J Solids Struct 1989;25(7):803–33. http://dx.doi.org/10.1016/0020-7683(89)90015-2. [56] de Vree JHP, Brekelmans WAM, van Gils MAJ. Comparison of nonlocal approaches in continuum damage mechanics. Comput Struct 1995;55 (4):581–8. http://dx.doi.org/10.1016/0045-7949(94)00501-S.

23

[57] Comi C, Perego U. Fracture energy based bi-dissipative damage model for concrete. Int J Solids Struct 2001;38(36–37):6427–54. http://dx.doi.org/ 10.1016/S0020-7683(01)00066-X. [58] Rots JG. Computational modeling of concrete fracture. Ph.D. thesis, Delft University of Technology, Delft, The Nederlands; 1988. [59] Petersson PE. Crack growth and development of fracture zones in plain concrete and similar materials. 28 TVBM-1006, Division of Building Materials, Lund Institute of Technology, Lund, Sweden; 1981. [60] Winkler B, Hofstetter G, Lehar H. Application of a constitutive model for concrete to the analysis of a precast segmental tunnel lining. Int J Numer Anal Met 2004;28(7–8):797–819. [61] Pinheiro DCC, Barros FB, Pitangueira RLS, Penna SS. Métodos da partiç ao da unidade na análise n ao linear de estruturas. In: Proceedings of the XXXVII Ibero-Latin American congress on computational methods in engineering. Brasília, Brazil: ABMEC – Brazilian Association of Computational Methods in Engineering; 2016. [62] Carreira DJ, Chu KH. Stress-strain relationship for plain concrete in compression. J Am Concrete I 1985;82(6):797–804. [63] Carreira DJ, Chu KH. Stress-strain relationship for reinforced concrete in tension. J Am Concrete I 1986;83(1):21–8. [64] Boone T, Ingraffea AR. Simulation of the fracture process at rock interfaces. In: Proceedings of the fourth international conference in numerical methods in fracture mechanics; 1987.