A meshless approach to non-local damage modelling of concrete

Engineering Analysis with Boundary Elements 79 (2017) 62–74

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

A meshless approach to non-local damage modelling of concrete a,b

Behzad V. Farahani P.M.G.P. Moreirab a b

a,b,⁎

, J. Belinha

a

MARK

a

, F.M. Andrade Pires , António J.M. Ferreira ,

Faculty of Engineering, University of Porto, Dr. Roberto Frias Str., S/N, 4200-465 Porto, Portugal Institute of Science and Innovation in Mechanical and Industrial Engineering, Dr. Roberto Frias Str., 400, S/N, 4200-465 Porto, Portugal

A R T I C L E I N F O

A BS T RAC T

Keywords: Non-linear, Constitutive damage model Meshless method Concrete RPIM

A non-linear continuum damage mechanics model for concrete constructions is analysed using a radial point interpolation meshless method (RPIM). The fundamental mathematical relations and the material model are fully characterized. The 2D plane stress RPIM formulation is extended to a rate-independent standard (local) damage model considering both tension and compression static states. Additionally, in this work, the local damage formulation is modified considering a non-local constitutive damage criterion with regard to a Helmholtz free energy potential. Here, the internal variational fields, such as local and non-local damage variables, are determined by a return-mapping damage algorithm. Due to the non-linear nature of the phenomenon, a displacement controlled Newton-Raphson iterative approach is adopted to attain the non-linear damage solution. In the end, the performance of the proposed non-local damage model is evaluated using an experimental test of a notched-three-point bending beam available in the literature. The obtained solution shows that the meshless methods are capable to effectively analyse concrete structures assuming a non-linear non-local continuum damage model.

1. Introduction Nowadays, the industrial structural design relies mainly on the finite element method (FEM) [1] analysis to obtain efficiently accurate solutions. This status quo can be explained with the vast number of commercial FEM software packages available in the market and with the robustness and reliability of the FEM. However, this fact cannot hide some of the FEM drawbacks, such as the solution dependency on the element mesh discretization (due to the mesh-based interpolation) or the reduced continuity of its shape functions. Numerous challenging fields in computational mechanics involve the study of the numerical non-linear local damage solution of concrete materials using FEM formulations, see e.g. [2–9]. Distinct aspects of the continuum damage mechanics field have been investigated by researchers to attain the experimental solution [9,10]. The standard –local– damage mechanics theory is applicable to analyse several problems as introduced by Kachanov [11]. Different types of material experiencing brittle or ductile behaviours could be studied through continuum damage mechanics theory [see i.e. Krajcinovic et al. [12,13]; Resende et al. [14]; Voyiadjis et al. [2]]. The foregoing authors have precisely focused on the local damage mechanics in concrete structures. Afterwards, mathematical relations for rate-independent damage mechanics associated with the local ⁎

damage were presented by Crisfield [15], Cervera et al. [3–5,16]. The continuum rate-independent damage formulation considered in this study respects Crisfield’s and Cervera’s hypothesis. Basically, the degradation of the constitutive model due to the presence of tensile and compressive enforced displacement states include various principal stress terms in standard damage models [8]. As a fact, the standard local constitutive models are inappropriate whenever strong strain softening is encountered. Thus, the governing differential relationships might lose ellipticity. Numerically, this situation appears itself by spurious mesh sensitivity of finite element computations. As the mesh is refined, the strain starts to localize into a narrow band whose width depends on the element size and tends to zero. Hence, the corresponding relation of load-displacement always experiences snapback for a sufficiently fine mesh, and the total energy dissipated by fracture converges to zero [17]. In FEM studies, the most trustful approach to tackle the aforementioned disturbance, is to regulate the post-peak slope of the stressstrain curve as a function of the element size. Consequently, the energy dissipated in a band of cracking elements will be independent from the bandwidth. It is possible to find in the literature more refined strategies guaranteeing the corresponding objectivity (so-called localization limiters) including higher-order gradient models, see e.g. [18–22] and also

Corresponding author at: Faculty of Engineering, University of Porto, Dr. Roberto Frias Str., S/N, 4200-465 Porto, Portugal. E-mail address: [email protected] (J. Belinha).

http://dx.doi.org/10.1016/j.enganabound.2017.04.002 Received 9 August 2016; Received in revised form 5 April 2017; Accepted 10 April 2017 0955-7997/ © 2017 Elsevier Ltd. All rights reserved.

Engineering Analysis with Boundary Elements 79 (2017) 62–74

B.V. Farahani et al.

2. Meshless method

for viscoplastic regularization [23]. The concept of non-local averaging is sufficient for the localization limiters and it can be applied to any kind of constitutive model. The idea of non-local continuum models was proposed by Eringen [24] and later on it was developed for the strain-softening materials by Bazant [25]. Afterwards, Pijaudier-Cabot proposed an improvement leading to the non-local damage theory [26]. Its early extension was established into various approaches of non-local models for damage and fracture mechanics by Jirásek [17]. Furthermore, non-local plasticity combined with the finite element method was presented by Stromberg et al. [27]. Recently, there are other works available in the literature regarding the application of non-local damage formulations, such as: on nanocomposites and composite laminated panels [28,29]; concrete materials [30]; complex microstructures [31], dynamic ductile fracture problems combined with an Extended FEM (XFEM) viscoplasticity model [32] and also in manufacturing simulations [33]. It is possible to find some experimental tests reporting the behaviour of concrete materials, such as the softening response of concrete under monotonic uniaxial tension test [34], the behaviour of concrete under compressive enforced displacement [35], the response of concrete under biaxial stress states [36] and the three-point-bending tests on single-edge notched beams [37]. Various demanding isotropic non-linear damage models for concrete structures were analysed with the FEM formulations, such as linear elastic models, see e.g. [9,38,39], rate-dependent models [3,5], viscous-damage models [16] and, more recently, elasto-plastic damage models for crack propagation using XFEM formulations [40]. However, the state-of-the-art lacks a reliable work on meshless methods combined with a non-linear damage constitutive model, particularly whenever the non-local model is required. Therefore, this work aims to fulfil a gap in the RPIM state-of-theart, presenting an extensive and complete numerical study of the RPIM regarding the analysis of continuum damage mechanics theory. In the literature, it is possible to find relevant research works combining meshless methods with damage mechanics [41,42]. Nevertheless, those works use particle methods, such as the smoothed particle hydrodynamics (SPH) and reproducing kernel particle methods (RKPM). These formulations, in opposition to RPIM, lack the delta Kronecker property. The efficient aspects of the non-local and standard damage models are fully addressed. Since the authors fully developed their original code to analyse the non-linear phenomenon here proposed, all the relevant attributes of the RPIM can be studied with detail and validated for the proposed damage constitutive model, such as: the size of the influence domain; the integration scheme; the internal fields and optimized damage variables; and the global efficiency. The present formulation is suited for static or quasi-static structural applications, since inertia is not accounted in this continuum local damage model. This work is organized as follows; first, in Section 2, the Radial Point Interpolation Method (RPIM) is presented, including the integration scheme, the nodal connectivity, radial point interpolators. Additionally, the 2D elastic solid mechanics formulation for the plane stress state combined with the RPIM, assuming the Galerkin weak form, is shown. In Section 3, first, the rate-independent standard damage theory is formulated based on a Helmholtz free energy function. Furthermore, the foregoing damage formulization is extended to the non-local fashion. Besides, the non-linear return-mapping algorithm of the numerical implementation is also introduced in Section 3. In Section 4, notched-three-point-bending concrete beams are considered and solved with the proposed computational approach. Convergence studies are performed aiming to accomplish the optimum non-local damage characteristics. Then, the results obtained are compared to the experimental solution. The manuscript ends with the final remarks and conclusion, in Section 5.

In this work, an advanced discretization meshless technique, the Radial Point Interpolation Method (RPIM), is exerted [43,44]. The RPIM is an interpolator meshless method and forces the nodal connectivity using the influence domain concept. The background integration mesh (required to integrate numerically the integro-differential equations) is constructed using a background regular lattice filled with integration points respecting the Gauss Legendre quadrature principle. In meshless methods the nodal cloud discretizing the problem domain does not form a mesh, because no previous information regarding the relation between each node is required to construct the interpolation functions [41].

2.1. Numerical integration and nodal connectivity Similar to FEM, meshless methods are classified in discrete numerical approaches. In meshless methods the domain is discretized with a nodal distribution X = {x1 , x2 , …, xN } ∈ Ω Λ xi ∈ 2 , being N the total number of nodes discretising the problem domain. Then, a background integration mesh, is constructed, Q = {x1 , x2 , …, xQ } ∈ Ω Λ xI ∈ 2 , being Q the total number of integration points discretising the problem domain. Within the RPIM, the integration mesh is completely independent on the nodal mesh and it is required to numerically integrate the integro-differential relations governing the studied phenomenon, the weak-form of Galerkin. As it is well-known, in the FEM formulation, the background integration mesh is constructed using the geometry of the elements and it can be accurately defined because the order of the polynomial shape functions, obtained from the nodal arrangement of the assumed finite element, is known. Hence, it is feasible to determine the number of integration points existing in any integration cell in advance accurately [1,44,45]. However, since in meshless methods the degree of shape function is undisclosed, this predefinition is out of question. Nevertheless, Wang et al. [46] suggested an integration scheme acceptable for RPIM formulations. First, the problem domain is divided in a regular lattice, forming a quadrilateral grid, with no voids or overlaps. Then, each quadrilateral is filled with integration points respecting the Gauss-Legendre quadrature principle [43]. The literature shows that each integration cell should contain 3 × 3 integration points for a nodal arrangement as the one presented in Fig. 1. Regarding the nodal connectivity, the RPIM uses the influencedomain concept [43]. This concept implies that each integration point xI has to search in the problem domain for the n nearest nodes. Therefore, each integration point xI will possess its own influencedomain. This concept will permit to construct the shape functions of each integration point xI . This issue is described with detail in the literature, see e.g. [43,44]. The performance and accuracy of the final solution is strongly influenced by the size of the influence domains. Thus, it is important to allocate approximately the same number of nodes on each influence domain [43]. In the literature, it is possible to find some meshless

Fig. 1. - An integration cell with 3×3 integration points in the discrete model [45].

63

Engineering Analysis with Boundary Elements 79 (2017) 62–74

B.V. Farahani et al.

distance is directionless, ri (xj , yj ) = rj (xi , yi), i.e. Rij = Rji , matrix R is symmetric. A unique solution is obtained if the inverse of the radial moment matrix R exists, {a b}T = G−1 { u z }T . The solvability of this system is usually guaranteed by the requirements rank ( p )=m ≤ n [53]. In this study, the influence-domain will always possess enough nodes to largely satisfy the previously mentioned condition. It is possible to obtain the interpolation with

works proposing an optimal value for the number of nodes inside the influence-domain, the majority suggest that for 2D analyses n should be between 9 and 16 nodes, [43–48]. 2.2. Radial Point Interpolators The RPIM shape functions are obtained using the Radial Point Interpolators (RPI), which combine radial basis functions with polynomial basis functions. In this work, only simplified two-dimensional domains Ω ⊂ 2 are studied, therefore it considers an interpolation function uh (x) defined in an influence-domain ΩI ⊂ Ω of an interest point xI ∈2 . The nodal set of the domain is defined in the twodimensional space by X = {x1 , x2 , …, xN } ∧ xi ∈2 . Thus, the influence domain of an interest point xI is defined by {x1 , x2 , …, xn } ⊂ X , where n is the number of nodes in the influence-domain of xI . The RPI constructs the interpolation function uh (x) capable to pass through all nodes within the influence-domain, meaning that since the nodal function value is assumed to be ui at the node xi , ui=u (xi ), consequentlyuh (xi )=u (xi ). Using a radial basis function r (x) and polynomial basis function p (x), the interpolation functionuh (x) can be defined at the interest point xI ∈Q (not necessarily coincident with any xi ∈X ), n

uh (xI ) =

⎧ u⎫ ⎧ u⎫ uh (xI ) = {r (xI )T ; p (xI )T } G−1 ⎨ z ⎬ = {Φ (xI )T ; Ψ (xI )T } ⎨ z ⎬ ⎩ ⎭ ⎩ ⎭

where the interpolation function vector Φ(xI ) is identified byΦ(xI ) = { φ1 (xI ) φ2 (xI ) … φn (xI )} and the residual vector Ψ(xI ), with no relevant physical meaning, is expressed as:Ψ(xI ) = { ψ1 (xI ) ψ2 (xI ) … ψm (xI )}. The RPI test functions Φ(xI ) depend uniquely on the distribution of scattered nodes [43]. Previous works [43,46,51] show that RPI test functions possess the Kronecker delta property. Since the acquired RPI test functions have a local compact support it is possible to assemble a well-conditioned and a banded stiffness matrix. If a polynomial basis is included, the RPI test functions have reproducing properties and possess the partition of unity property [43].

m

∑ ri (xI ) ai + ∑ pj (xI ) bj = r (xI )T a+p (xI )T b = u (xI ) i =1

j =1

2.3. Equation System (1) The discrete equation system is obtained using the Galerkin weakform, which can be expressed as,

Notice that ai is the non-constant coefficient of ri (xI ) and bj is the non-constant coefficient for pj (xI ). The integers n and m are the number of nodes inside the influence-domain of the interest point xI . The vectors are defined as aT = {a1, a2, …, an}, bT = {b1, b 2, …, bm}, r (x)T = {r1 (x), r2 (x), …, rn (x)} and p (x)T = { p1 (x), p2 (x), …, pm (x)}, in which xi = (xi , yi). This work uses the Multiquadrics Radial Basis Function (MQ-RBF) [43,47,48,49,50], which can be defined by p ri (xI )=s (diI ) = (diI2 + c 2 ) , where diI is a distance between the interest point xI = (xI , yI ) and the node xi = (xi , yi), where diI = (xi − xI )2 + ( yi − yI )2 . The parameters c and p are the MQ-RBF shape parameters, which are fixed parameters determined in previous works [45,46]. The variation of these parameters can affect the performance of the MQ-RBFs. In the work of Wang and Liu [46,47], it was shown that the optimal values are c = 1.42 and p = 1.03, which are the values used in this work. The original RPI formulation requires a complete polynomial basis function, following the Pascal’s triangle. For instances, for the 2D space the quadratic polynomial basis is defined as p (xi )T = {1,xi , yi , xi2 , xi yi , yi2,…}. However, it was shown in previous RPI research works [43,51,52] that using a simple constant basis increases the RPI formulation efficiency. Thus, in this work only the constant basis is considered p (xi ) = {1}, for which the number of monomial terms is defined by m = 1. The coefficients ai and bj in Eq. (1) are determined by enforcing the interpolation to pass through all n nodes within the influencedomain [43]. The interpolation at the k th node is defined by, n

uh (xk , yk ) =

i =1

∫Ω δεT σ d Ω = ∫Ω δuT b d Ω + ∫Γ δuT t d Γ + δuT q

j =1

⎧ ux (xI ) ⎫ ⎬= uh (xI ) = ⎨ ⎩ u y (xI )⎭

n

⎡ φi (xI )

∑⎢ i =1

⎣ 0

0 ⎤ ⎧ u (xi )⎫ ⎬= ⎥⎨ φi (xI )⎦ ⎩ v (xi ) ⎭

n

∑ Hi (xI ) u (xi ) i =1

(6)

Besides, as a result of Eq. (6), it is possible to developed the strain field equation:

⎡ ∂ 0⎤ ⎢ ∂x ⎥ n ∂ ⎥ ⎢ ε (xI ) = ⎢ 0 ∂y ⎥ uh (xI ) = Luh (xI ) = L ∑ Hi (xI ) u (xi ) i =1 ⎢∂ ∂⎥ ⎢⎣ ∂y ∂x ⎥⎦

(2)

The inclusion of the following polynomial term is an extra-requirement n that guarantees unique approximation, ∑i =1 pj (xi , yi) ai=0, j = 1,2, …, m [43,52]. Thus, the computation of the shape functions can be written in a matrix form as

⎧ a ⎫ ⎧ u⎫ ⎡ R P ⎤ ⎧ a⎫ ⎧ u ⎫ ⎢ T ⎥ ⎨ ⎬ = ⎨ ⎬ ⟺G ⎨ b ⎬ = ⎨ z ⎬ ⎣ P Z ⎦ ⎩ b⎭ ⎩ a ⎭ ⎩ ⎭ ⎩ ⎭

(5)

Vector b is the body force vector, t is the external traction force vector applied on a natural surface boundary Γ and q represents the external force vector applied locally. The infinitesimal domain dΩ is defined asd Ω=dx. dy. dz . In the RPIM, the weak form has local support, which means that the discrete system of equations is developed firstly for every influence-domain. After that, the local system of equations is assembled to form the global system of equations, which is solved afterwards. The RPIM trial function is given by Eq. (4), thus for each degree of freedom, it is possible to n n write:uxh (xI ) = ∑i =1 φi (xI ) ux (xi )and u yh (xI ) = ∑i =1 φi (xI ) u y (xi ). In the abovementioned equations, φi (xI ) represents the RPIM interpolation function while ux (xi ) and u y (xi ) are the nodal parameters of the i th node belonging to the nodal set defined in the influence domain of the interest point xI . Subsequently, it is reasonable to generate a more general equation:

m

∑ ri (xk , yk ) ai + ∑ pj (xk , yk ) bj = uk , k =1, 2, …, n

(4)

n

ε (xI ) =

∑ Bi (xI ) u (xi ) = ∑ i =1

(3)

Matrix G is the complete moment matrix, matrix Z is a null matrix defined as Z ij =0, ∀{{i , j}∈ {i , j}≤m} and the null vector a can be represented by a i=0, ∀{i∈i ≤ m}, being  the set of natural numbers. The vector for function values is defined as ui=u (xi ),∀{i∈i ≤ n}. The p radial moment matrix R is represented as Rij = ri (xj ) = (dij2 + c 2 ) and the polynomial moment matrix P is defined Pki = pk (xi ). Since the

n

i =1

⎡ ∂φi (xI ) ⎢ ∂x ⎢ ⎢ 0 ⎢ ∂φi (xI ) ⎢⎣ ∂y

(7)

0 ⎤⎥

∂φi (xI ) ⎥ ⎧ ux (xi ) ⎫ ⎬ ⎨ ∂y ⎥ u (x ) ⎩ y i ⎭ ∂φi (xI ) ⎥ ⎥ ∂x ⎦

(8)

The matrixB is the deformation matrix,B = LH . As mentioned before, the stress field is a function of the strain vector. Thus, the developed relation of the stress vector for an interest point (xI ) could be written, using the matrix notation, as σ (xI )=Cε (xI ). The material constitutive matrix C is defined as, 64

Engineering Analysis with Boundary Elements 79 (2017) 62–74

B.V. Farahani et al.

⎡1 ν ⎢ν 1 E C= ⎢ (1+ν )(1−ν ) ⎢ 0 0 ⎣

0 ⎤ 0 ⎥ ⎥ 1−ν ⎥ 2 ⎦

tion of all the elastic material properties such as Young’s modulus, Poisson’s ratio, shear modulus and other relevant properties. Considering the damage parameters, it is possible to infer that the corresponding model could predict mixed mode failure. Consequently, the equivalent effective tensile and compressive norms are obtained after splitting the stress as follows:

(9)

The Young’s modulus and Poisson’s ratio are defined by E and ν , respectively, and e the thickness of the plate. To evaluate the stiffness matrix, first it is necessary to present the general integration of the weak formulation for any interest point(xI ): δuT

⎧b ⎫

⎧t ⎫

τ+ =

⎧q ⎫

− σoct

Being A the domain area and ΓC the curve in which t is applied. The linear system of equations based on Eq. (10) is represented as,δuT (Ku−fb −ft −fq )=0 , leading to: u = K−1 ( fb+ft +fq ). Since the RPI test functions possess the delta Kronecker property, the essential boundary conditions are directly imposed in the global stiffness matrix. 3. Rate-independent damage formalism In this section, the continuum damage models used in this work are presented. First, the continuum local damage model is described and then, the numerical extension to the non-local model is shown. Due to its simplicity, in this study, isotropic materials with isotropic damage are considered. Indeed, the response of the material is independent on time, as well as the influence of the temperature and humidity are here neglected. The material behaviour is assumed as brittle-elastic.

g+(τ +, r +) = τ +−r +≤0

(17)

g−(τ −, r −) = τ −−r −≤0

(18)

Eq. (18) is Drucker-Prager cone for compressive state, r + and r − are current damage thresholds controlling the expansion of the damaged surface. In the beginning of the damage routine, r0+ and r0−are initialized as material properties associated to the uniaxial maximum strengths. The kinematic of the tensile and compressive internal variables are defined based on the consistency conditions for loading through KuhnTucker equations as follows [54]:

∂G+(r +) = Ġ +(r +)≥0 ∂r +

(19)

∂G−(r −) r −̇ = τ −̇ ≥0, ḋ − = r −̇ = Ġ −(r −)≥0 ∂r −

(20)

r +̇ = τ +̇ ≥0, ḋ + = r +̇

3.1. Standard local damage model The theory of the continuum damage mechanics relies on the definition of the effective stress concept connected to the equivalent strain. Meaning that the strain value related to the damage state, when the stress σ applied, is equivalent to the strain obtained from the undamaged state under the effective stressσ . Considering this principle, the effective stress tensor is expressed as follows:

G+

G−

where and present the proper monotonic increasing function, which is captured from experimental tests [3]. Therefore,

r + = max(r0+, max (τ +)), d + = G+(r +) r− =

(11)

max(r0−,

max (τ −)), d − = G−(r −)

σi pi ⊗pi

i =1

r0+ exp (A+ (1 − r + / r0+)) if r +≥r0+ r+

(23)

− dlocal = G−(r −)=1 −

r0− r−

(24)

∑

σi pi ⊗pi

i =1

and in compression and A+ in The local damage coefficients tension will be determined in each numerical example. Regarding A+ , there exists a parameter, so-called “characteristic length, lch ” , playing a substantial role in its determination. In the literature, it is possible to find semi-empirical equations to determine the characteristic length [55–57]. Precisely, Faria et al. [3,6] [62] reported that the characteristic length could be defined as follows:

(12)

(13)

Notice that ⊗ operates as the dyadic product and the sign of (+) and (-) indicates the tensile and compressive entities, respectively. The Macaulay brackets . indicates that the expression is equal to the enclosed value when positive and zero if it is negative. Besides σ can be determined as σ =σ − σ . Assuming σ = σ ++σ −, it is necessary to present the Cauchy second order stress tensor σ as the final relation characterizing the full stress tensor for the isotropic model. Thus,

σ =(1 − d +) σ ++(1 − d −) σ −

(1−A− ) − A−exp (B−(1 − r − / r0−)) if r −≥r0−

A−

3

σ − = ⟩σ ⟨ =

(22)

+ dlocal = G+(r +) = 1 −

3

∑

(21)

Then, local damage variables are computed within the following relations [3]:

where : operates as the tensor multiplication. The tensors σ and ε are the second order tensors of the effective stress and strain, respectively. Moreover, the linear elastic constitutive tensor,C , is calculated based on Eq. (9). Indeed, it is necessary to split the effective stress tensor into tensile and compressive components. Therefore, the following relations are applicable [3–5,16]:

= σ =

(16)

− τoct

and are the octahedral normal and shear stresses Tensors that should be obtained from σ − tensor, Eq. (13). Furthermore, K is a material property depending on the ratio between the biaxial and uniaxial compressive strengths for concrete materials, obtained as in the literature [3]: K = 2 (β −1)(2β −1). The constant β is typically assumed as β = 1.16 for concrete materials [3]. The damage criteria for tensile and compressive states are expressed as [54]:

(10)

σ+

(15)

− − 3 (Kσoct + τoct )

τ− =

∫A eBT CB d A u = δuT ∫A eH T ⎨⎩ bxy ⎬⎭ dA+δuT ∫ΓC eH T ⎨⎩txy ⎬⎭ d ΓC +δuT ⎨⎩ qyx ⎬⎭

σ =C:ε

σ + : C−1 : σ +

lch =

B−

A+ H +(2+A+ )

(25)

H + is

dependent on the material properties obtained Notice that within the following equation:

(14)

H+ =

In Eq. (14), the damage parameters in tension and compression are represented as d + andd −. As mentioned before, both values are scalar and they should be limited as 0 ≤ d +, d −≤1. Notice that thermomechanical consideration related to the non-negativeness of the dissipation + − expresses the following condition, {d ̇ , d ̇ }≥0 [3]: Both damage parameters have the potential to control the degrada-

( f +) 2 (r0+)2 = 0 + + 2G f 2EG f

(26)

So that G +f represents the tensile fracture energy per unit area and f 0+ is the uniaxial tensile strength available for the material in the literature. The corresponding conditions must be satisfied in the 1 problem analysis: A+ ≥0 and lch≤ H + . The importance of lch is to control the maximum size of the finite 65

Engineering Analysis with Boundary Elements 79 (2017) 62–74

B.V. Farahani et al.

elements (in the FEM) or influence domain (in meshless methods). It means that if the size of elements in FEM or discretization in meshless method increases, the softening branch of the response undergoes an almost vertical slope and perhaps the fracture procedure is more brittle.

Table 1 Weight function for non-local damage process.

3.2. Non-local damage model

J =1

wd+/− +/− wtotal

w=1 w=−

w=2

(dIJ )2 +1 μ2

(dIJ )3 (d )2 −3 IJ2 +1 μ3 μ

+ − − + σ = (1 − dnon − local ) σ +(1 − d non − local ) σ

σ+

(29)

σ−

and are the effective stress tensor in tension and Tensors compression obtained from Eqs. (12) and (13), respectively. To calculate the total damage parameter, the equivalent von-Mises effective stress σ∼ of the new Cauchy stress tensor, derived from Eq. (29), is used. σ∼ σ∼ = [1 − (dtotal ) I ] σ∼ → (dtotal ) I = 1 − ∼ (30) σ where d presents the total damage value calculated for each interest point xI ∈Ω , σ∼ is the corrected equivalent von-Mises effective stress which is based on Eq. (11) [15]. Since the current damage model is based on the uncoupled elasticity theory proposed by Faria et al. [58], the Helmholtz free energy potential is applied, see e.g. [58,59]. 3.3. Numerical scheme of the constitutive law The numerical approach here proposed relies on the strain-driven formalism of the corresponding constitutive model. Box 1 represents the return-mapping algorithm for rate-independent elastic damage to evaluate Cauchy stress tensor with regard to a given strain tensor for any pseudo-time stepping scheme using RPIM meshless formulations. Indeed, this algorithm trusts on the mathematical relations for standard and non-local damage models. The return-mapping algorithm evolves the non-linear Newton-Raphson method to acquire the non-linear damage solution. In the same way, the principal effective stress tensor is adopted instead of the general form. Besides, the computational implementation of the non-local damage is demonstrated in Box 2.

(27)

Box 1. Computational algorithm to evaluate the damage parameter for rate-independent continuum damage models

The number of neighbours is represented by n and wJ presents their weights using the spatial position with respect to xI . Adding all the n +/− = ∑J =1 wJ , the non-local damage variable neighbours’ weight as wtotal associated to the interest point xI is evaluated for tensile and compressive states as; +/− dnon − local =

0th 2nd

After the damage localization stage, it is possible to calculate the Cauchy second order stress tensor based on the following relation – developed fashion of Eq. (14) [15]:

n +/− +/− ) J ×wtotal ∑ (dlocal

Weight function

3rd

The previously presented standard formulation of damage theory was local. Now, in this subsection, the non-local damage model is described with detail. Since, the standard local constitutive models are inappropriate for materials exhibiting strong strain softening, the localization process is more suitable to tackle the mentioned numerical difficulty. A simple computational solution is to adjust post peak slope of stress-strain curve as a function of the element size [17]. For a specific interest point xI ∈Ω , first, the local tensile and compressive damage must be calculated within the rate-independent continuum damage mechanics formulations as explained previously, Eqs. (23) and (24). Afterwards, a circle with a radius, μ , centred at the interest point, is specified. This circle includes the allowable integration points (neighbours) participating in damage localization process, see Fig. 2. The radius of this circle is obtained by: μ = s. h which is dependent on the average distance between nodes, h . Consider the nodes discretised in specific divisions along x and y directions, h is computed as L D h = divisions along x = divisions along y , in which L and D are the specimen dimensions along x and y directions. Additionally, s is a normalization parameter having an influence on μ . Subsequently, to evaluate the distance between the interest point xI = (xI , yI ) and its neighboursxJ = (xJ , yJ ), the Euclidean norm is applied:dIJ = (xI − xJ )2 + ( yI − yJ )2 <μ, as shown in Fig. 2. Within the non-local damage models, a weight function must be used to smooth the local damage surrounding the intended interest point xI . Consequently, three weight functions were considered in this work, Table 1. First, the corresponding Jth neighbours (the ones surrounding the interest point xI and satisfying dIJ <μ) are selected. Then, these neighbour nodes will contribute to the localization procedure. Thus, the final non-local weighted damage for the Ith interest point is computed with the following multiplication:

(wd+/−) I =

Order

Step n = 0 , Initialization. (i) Set rn+ = r0+ , rn− = r0− dn+=0 , dn−=0 Step n + 1 (ii) Evaluate εn +1 Then, compute σn + 1 = C : εn + 1 Principal stress and directions (iii) Split σn + 1 into σ n++ 1 and σ n−+ 1 − − Obtain σoct , τoct + (iv) Calculate τn+1 and − τn+1 (v) Check the tensile and compressive damage If τn++1 > rn+ ,set

(28)

rn++1 = max {rn+, τn++1}

Fig. 2. - A schematic view of the interest and neighbour integration points in three point bending beam.

66

Stress norms Scalar damage variables Trial strain, Eq. (8) Effective stress of damage routine, Eq. (11)

Effective tensile and compressive stress, Eqs. (12) and (13) Tensile and compressive equivalent stress, Eqs. (15) and (16)

Update tensile damage threshold, Eq. (21)

Engineering Analysis with Boundary Elements 79 (2017) 62–74

B.V. Farahani et al.

Fig. 3. A notched-three-point bending beam (a) geometry (b) a regular mesh and (c) an irregular mesh with 2441 nodes.

Fig. 4. - RPIM nodal distributions for three-point-bending beam corresponds to (a) 295 nodes, (b) 505 nodes, (c) 771 nodes and (d) 2441 nodes.

If τn−+1 > rn− , set = max {rn−, τn−+1} (vi) Update of damage variables + (dlocal )n +1 = G+(rn++1)

rn−+1

− (dlocal )n +1 = G−(rn−+1)

(v) For the Ith interest point; (vi) Loop over all Jth integration points (neighbours) in the vicinity of xI ; - Check the distance dIJ < μ ; If yes, it is the acceptable neighbour point. Use the weight function (w) based on Table 1;

Update compressive damage threshold, Eq. (22)

Local Tensile damage parameter, Eq. (23) Local Compressive damage parameter, Eq. (24)

n

+/− Update wtotal = ∑J =1 wJ ;Update wd+ and wd−.

End of the local damage test. (vii) Go to Box 2, obtain the non-local tensile and compressive damage parameters; (viii) Compute the final Cauchy stress tensor (Non-local State); + + − − σn +1 = (1 − (dnon − local )n +1) σ n+1+(1 − (d non − local )n +1) σ n+1

(ix)

(27′)

If no, the neighbour point is outside the circle, go for the next integration point. (vii) Calculation of the non-local damage for xI ; +/− (dnon − local ) I =

(29′)

Calculation of total damage parameter;

wd+/− +/− wtotal

(28′)

(viii) Updating the vector of non-local tensile and compressive damage for all interest points; End of test.

σ∼ dtotal = 1 − ∼ (30′) σ (x) Update the total damage vector for all the integration points; End of test.

4. Verification of the non-local damage model The main purpose of the present manuscript is to present in detail the extension of the RPIM to a continuum damage model for concrete materials. Hence, no large deformation or instability phenomena (such as snap-through and snap-back effects) are considered here. Thus, the present damage formulation deals exclusively with the material nonlinearity. Therefore, a well-known benchmark example (the three-point bending beam) was selected from the literature, in order to validate the presented non-linear formulation. In this section, a notched-three-point bending beam associated with the concrete material is investigated using the proposed rate-independent damage model combined with the RPIM formulation. The results are compared to the experimental result derived from literature [37].

Box 2. The numerical algorithm to obtain non-local damage parameters.

Step n = 0, Initialization. − + = 0; (i) Set wtotal = 0 and wtotal + − (ii) Set wd = 0 and wd = 0 ; (iii) Defining s∈[0. 5, 2. 1]; (iv) Radius of the circle to include the neighbours μ=s. h ; Step n+1, loop over all integration points.

67

Engineering Analysis with Boundary Elements 79 (2017) 62–74

B.V. Farahani et al.

Table 2 μ values based on different s obtained for distinct RPIM meshes.

s

μ (m ) Mesh 1

Mesh 2

Mesh 3

0.5 0.7 0.8 1.1 1.6 2.1

0.0082 0.0115 0.0131 0.0181 0.0263 0.0345

0.0062 0.0086 0.0099 0.0135 0.0197 0.0259

0.0049 0.0069 0.0079 0.0108 0.0158 0.0207

h (m)

h1=0.0164

h2=0.0123

h3=0.0099

Fig. 6. - Response of load P- deflection on point A with regard to 2nd-order weight function corresponds to (a) mesh 1, (b) mesh 2 and (c) mesh 3.

4.1. Model's parameters definition Fig. 3 illustrates the three-point bending beam with the essential boundary conditions. Taking advantage of the symmetric geometry, only half of the beam is considered. As it is visible in Fig. 3(b), the nodal distribution is optimized to capture the potential stress concentration on the beam mid-span. Regarding the domain discretization, the beam is divided into two sections: [0, 3L/8] and [3L/8, L/2], each one showing a distinct discretization level, Fig. 3(b). In this example, distinct nodal distributions are considered in order to capture the influence of the nodal discretization on the obtained solution. The beam presents a total length of L = 0.788(m ) and a width of D = 0.102(m ). The thickness of the beam is e = 0.102(m ) and the initial

Fig. 5. - Response of load P- deflection on point A with regard to 0th-order weight function corresponds to (a) mesh 1, (b) mesh 2 and (c) mesh 3.

68

Engineering Analysis with Boundary Elements 79 (2017) 62–74

B.V. Farahani et al.

Fig. 8. -Variation of damage vs. effective strain for 0- order weight function; (a) mesh 1, (b) mesh 2 and (c) mesh 3. Fig. 7. - Response of load P- deflection on point A with regard to 3rd-order weight function corresponds to (a) mesh 1, (b) mesh 2 and (c) mesh 3.

enforcement, imposed compressively on point A in the Oy direction, is identified as δ . The present analysis initiates with a convergence study on distinct RPIM nodal distribution patterns to examine the capability of the meshless method formulation extended for non-local damage model. At the end, it will be possible to obtain the optimal damage characteristics and then apply them in further benchmarks.

crack length in the mid-span of the beam is equal to a 0=0.051(m ). The material properties and damage characteristics of concrete are settled according to the work represented by Malvar et al. [37]: E = 21.7(GPa ), ν = 0.2 , f 0+ =2.4(MPa ), f 0− =29(MPa ) and G +f =30(N / m ). Based on the study conducted by Faria et al. [3] the corresponding compressive parameters defining the compressive damage behaviour of concrete used in Eq. (24) are: A− =1.00 andB−=0.89 . Moreover, the tensile damage character is identified as A+ =0.001. Consequently, another tensile property so-called H + is obtained based on Eq. (26) as: H + = 4.244(m−1). Following that, the characteristic length, associated with the element size and affecting the degradation of the damage phenomenon, is calculated based on Eq. (25) as: lch=1.1296 × 10−4 (m ). Concerning the algorithm configuration, the total displacement

4.2. Convergence study The present test analysis in compression is performed for several nodal discretizations, shown in Fig. 4, considering δ = 0.25(mm ). All numerical results obtained with the non-local damage model here presented are compared with the experimental solution reported by Malvar et al. [37] for a load P versus the deflection on point A (as seen in Fig. 3). Moreover, the non-local damage variation, in terms of 69

Engineering Analysis with Boundary Elements 79 (2017) 62–74

B.V. Farahani et al.

Fig. 9. Variation of damage vs. effective strain for the 2nd- order weight function; (a) mesh 1, (b) mesh 2 and (c) mesh 3.

Fig. 10. Variation of damage vs. effective strain for the for 3rd-order weight function; (a) mesh 1, (b) mesh 2 and (c) mesh 3.

effective strain, is compared to the one derived from the local model. Thus, the non-local damage formulation is employed to analyse the first three meshes for distinct s values ranging from 0.5 to 2.1 [45,60]. Moreover, three distinct weight functions presented on Table 1 are considered. Moreover, μvariation in terms of s for any of three nodal meshes was evaluated, see Table 2. Considering 0th-order weight function (see Table 1), the graph correlating the load P with the correspondent deflection on point A, for different s values, is presented in Fig. 5. Subsequently, assuming 2ndorder and 3rd-order weight functions, the new studies were performed and the response of P-A for various s values are presented in Fig. 6 and Fig. 7, respectively. Considering a particular integration point, xc , located at the crack tip, it is rational to acquire its total damage variation in terms of effective strain using the proposed non-local model. Once more, the analysis was performed for different weight functions. Hence, total damage is plotted at the particular integration point near the crack tip. Naturally, the exact position of xc is distinct for three RPIM meshes

shown in Fig. 4. Nevertheless, the variable magnitudes are comparable. Considering different s values, the obtained outcome is plotted in Fig. 8, Fig. 9 and Fig. 10 for 0th , 2nd and 3rd-order weight functions, respectively. The load-displacement diagrams demonstrated in Fig. 5, Fig. 6 and Fig. 7 visualise the influence of s variable, leading to distinct μ magnitudes. It is observable from the corresponding figures that large s values produce numerical solutions far from the experimental result. Additionally, it is confirmed that higher order weight functions permit to obtain RPIM solutions close to the experimental result. Considering all the above-mentioned description, the normalization parameter can be optimized ass=0.8 leading to produce the results consistently close to the experiment. Regarding the behaviour of weigh function order, it is recognisable from foregoing diagrams that the results are getting far away from the experimental solution when the 0th-order weight function was assumed. On the other hand, the RPIM solution are closer to the

70

Engineering Analysis with Boundary Elements 79 (2017) 62–74

B.V. Farahani et al.

Fig. 13. Response of P-deflection on point A obtained from irregular mesh compared to the regular mesh and experimental solution.

Table 3 Displacement enforcement stages.

Displacement enforcement (mm)

First stage

Second stage

Third stage

1.2841E−01

1.7526E−01

3.1581E−01

4.3. High density discretization analysis The main goal of the current study is to observe the efficiency and performance of the numerical algorithm when the optimum values previously reported are considered. So, material properties and damage characteristics were the same used for the convergence studies. It must be remarked that the total displacement enforcement was applied here as δ = 0.32(mm ). The average distance between nodes was adopted as h = 0.0049(m ) and μ was evaluated as:μ = 0.0039(m ). So, the load displacement response (load P - deflection on point A) was captured for this analysis and compared with the experimental result [37], as illustrated in Fig. 11(a). Furthermore, the variation of damage versus effective strain in the case of non-local and local damage at xc is presented in Fig. 11(b). Referring to Fig. 11, it is observable that the RPIM solution agrees very well with the experimental one. Therefore, this agreement confirms that the chosen methodology was efficient and practical so that the success of the proposed damage model and the optimization procedure was fulfilled. Overall, the load-displacement curves obtained with the proposed numerical approach for four nodal distributions considered in this study, are represented in Fig. 12.

Fig. 11. The final analysis results correspond to (a) response of P-deflection on point A and (b) damage variation versus effective strain.

Fig. 12. Comparison of P-deflection on point A response for different RPIM nodal distributions with the experimental curve.

4.4. Study of irregular mesh

experimental solution when the 2nd and 3rd orders are considered. Thus, the 0th- weight function order is discarded. In the same way, it is perceptible that the 2nd-order weight function increases the peak load point leading to delay the damage phenomenon. Notice that the damage initiates after reaching the ultimate elastic value where the concavity of the curves changed. Therefore, it is rational to abandon the 2nd-order weight function. Considering all aforementioned explanations, only the 3rd-order weight function will be considered in further studies. Consequently, the current analysis ends up with the last RPIM discretization as shown in Fig. 4(d), assuming optimal non-local damage parameters: s=0.8 and 3rd-order weight function.

In this subsection an irregular mesh is considered. Once more, taking advantage of the symmetry of the beam, only half of the problem was analysed. Thus, the irregular nodal distribution discretizing half of the problem domain is presented in Fig. 3(c), in which the nodes were distributed on the plate randomly. Seven irregular random meshes, as the one shown in Fig. 3(c), were constructed and analysed. All the material properties considered in the previous analyses were once more assumed. Additionally, all irregular meshes were solved assuming: s = 0.8 and 3rd-order weight function. The average solution of the seven irregular meshes is plotted in Fig. 13. It is observable that the obtained average curve has an acceptable agreement with the similar regular mesh with 2441 nodes and particularly the experimental curve. 71

Engineering Analysis with Boundary Elements 79 (2017) 62–74

B.V. Farahani et al.

Fig. 15. Equivalent von-Mises stress profile for three-point-bending beam at (a) first stage, (b) second stage and (c) third stage of enforced displacement.

Fig. 14. - Total damage contour for three-point-bending beam at (a) first stage, (b) second stage and (c) third stage of enforced displacement.

5. Conclusions

4.5. Graphical representation of damage variables

In this work, the Radial Point Interpolation Method (RPIM) was combined with a rate-independent continuum damage mechanics model. The standard damage mechanics fashion was extended to a non-local formulation, suited for non-linear analyses of concrete structures. The essential internal and damage variables were obtained within a return-mapping damage algorithm where the displacement was enforced for each pseudo-time stepping scheme. The predictive ability of the non-local model was compared with the local solution, which is inappropriate whenever strong strain softening is encountered. The chosen methodology presents several advantages: (a) using the principal effective stress instead of general stress tensor permits to get rid of the repetitive computations; (b) it can be implemented for both local and non-local constitute damage models producing results very close to the experimental solution; (c) it has the potential to study the micro-scale concrete structures explicitly; (d) it is capable to produce accurate and stable results with irregular mesh due to its insensitivity to the nodal arrangements. Furthermore, the present work focused on the damage characteristics influencing the response of load in terms of deflection and

Considering results previously acquired within the regular mesh presented in Fig. 4(d), in this subsection the damage profile and total equivalent von-Mises stress distribution are discussed. Since the proposed non-linear algorithm is incremental (displacement increments are imposed in the mid-span of the beam), it is rational to show the progress of the mentioned variable fields for specific incremental displacement levels. Three vertical displacement levels were selected, see Table 3. It was verified that the compressive damage remained null during all displacement enforcement stages, consequently, the total damage has the same value as the tensile one. The obtained contours, related to the damage profiles based on the particular displacement enforcement stages are presented in Fig. 14. Similarly, the total equivalent vonMises mapped stress is demonstrated in Fig. 15. These profiles confirm the prediction of damaged region near the crack tip. Furthermore, the crack tip zone presents larger damage values, as expected. Besides, it is perceptible an increasing trend of total effective stress value when the enforcement grows, particularly in the cracked region. 72

Engineering Analysis with Boundary Elements 79 (2017) 62–74

B.V. Farahani et al.

1986;53(4):791. [20] Vardoulakis I, Aifantis EC. A gradient flow theory of plasticity for granular materials. Acta Mech . 1991;87(3–4):197–217. [21] DE BORST R. Simulation of strain localization: a reappraisal of the cosserat continuum. Eng Comput . 1991;8(4):317–32. [22] Pamin JK. Gradient-dependent plasticity in numerical simulation of localization phenomena [12-Dec]. TU Delft, Delft University of Technology; 1994. [23] Needleman A. Material rate dependence and mesh sensitivity in localization problems. Comput Methods Appl Mech Eng . 1988;67(1):69–85. [24] Eringen AC. A unified theory of thermomechanical materials. Int J Eng Sci . 1966;4(2):179–202. [25] Bazant ZP. Imbricate Continuum and its Variational Derivation. J Eng Mech 1984;110(12):1693–712. [26] Pijaudier-Cabot G, Bazant ZP. Nonlocal damage theory. J Eng Mech 1987;vol. ASCE 113:1512–33. [27] Strömberg L, Ristinmaa M. FE-formulation of a nonlocal plasticity theory. Comput Methods Appl Mech Eng . 1996;136(1–2):127–44. [28] Silani M, Talebi H, Hamouda AM, Rabczuk T. Nonlocal damage modelling in clay/ epoxy nanocomposites using a multiscale approach. J Comput Sci . 2015. [29] Patel BP, Gupta AK. An investigation on nonlocal continuum damage models for composite laminated panels. Compos Part B Eng . 2014;60:485–94. [30] He W, Wu Y-F, Xu Y, Fu T-T. A thermodynamically consistent nonlocal damage model for concrete materials with unilateral effects. Comput Methods Appl Mech Eng . 2015;297:371–91. [31] Boeff M, Gutknecht F, Engels PS, Ma A, Hartmaier A. Formulation of nonlocal damage models based on spectral methods for application to complex microstructures. Eng Fract Mech . 2015;147:373–87. [32] Broumand P, Khoei AR. X-FEM modeling of dynamic ductile fracture problems with a nonlocal damage-viscoplasticity model. Finite Elem Anal Des . 2015;99:49–67. [33] Abiri O, Lindgren L-E. Non-local damage models in manufacturing simulations. Eur J Mech - A/Solids . 2015;49:548–60. [34] Gopalaratnam VS, Shah SP. Softening response of plain concrete in direct tension. J Proc 1985;82(3):310–23. [35] Karsan ID, Jirsan JO. Behavior of concrete Under compressive loadings. J Struct Devision 1969;95(ST12):2543–63. [36] Kupfer H, Hilsdorf HK, Rusch H. Behavior of concrete Under biaxial stresses. J Am Concr Inst 1969;66(8):656–66. [37] Malvar J, Warren G. Fracture energy for three-point bend tests on single-edge notched beams [no. March]. Nav Civ Eng Lab 1988:1–28. [38] Khan AR, Al-Gadhib AH, Baluch MH. Elasto-damage Model for High Strength Concrete Subjected to Multiaxial Loading. Int J Damage Mech 2007;16(3):361–98. [39] Labadi Y, Hannachi NE. Numerical Simulation of Brittle Damage in Concrete. Strength Mater 2005:268–81. [40] Roth S-N, Léger P, Soulaïmani A. A combined XFEM–damage mechanics approach for concrete crack propagation. Comput Methods Appl Mech Eng . 2015;283:923–55. [41] Li S, Liu WK (Wing K. Meshfree particle methods. New York: Springer-Verlag; 2004. [42] Danielson KT, Hao S, Liu WK, Uras RA, Li S. Parallel computation of meshless methods for explicit dynamic analysis. Int J Numer Methods Eng . 2000;47(7):1323–41. [43] Belinha J. Meshless Methods in Biomechanics: Bone Tissue Remodelling Analysis, 16. Netherlands: Springer; 2014. [44] Vasheghani Farahani B, Manuel Berardo J, Drgas R, César de Sá JMA, Ferreira AJM, Belinha J. The axisymmetric analysis of circular plates using the radial point interpolation method. Int J Comput Methods Eng Sci Mech 2015;16(6):336–53. [45] Farahani BV. The radial point interpolation meshless method extended to axisymmetric plates and non-linear continuum damage mechanics. University of Porto; 2015. [46] Wang JG, Liu GR. A point interpolation meshless method based on radial basis functions. Int J Numer Methods Eng 2002;54(11):1623–48. [47] Wang JG, Liu GR. On the optimal shape parameters of radial basis functions used for 2-D meshless methods. Comput Methods Appl Mech Eng 2002;191(23– 24):2611–30. [48] Nguyen VP, Rabczuk T, Bordas S, Duflot M. Meshless methods: a review and computer implementation aspects. Math Comput Simul 2008;79(3):763–813. [49] Belinha J, Dinis LMJS, Natal Jorge RM. Analysis of thick plates by the natural radial element method. Int J Mech Sci 2013;76:33–48. [50] Belinha J, Dinis LMJS, Jorge RMN. Composite laminated plate analysis using the natural radial element method. Compos Struct 2013;103:50–67. [51] Dinis LMJS, Jorge RMN, Belinha J. Analysis of 3D solids using the natural neighbour radial point interpolation method. Comput Methods Appl Mech Eng 2007;196(13–16):2009–28. [52] Dinis LMJS, Natal Jorge RM, Belinha J. Analysis of plates and laminates using the natural neighbour radial point interpolation method. Eng Anal Bound Elem 2008;32(3):267–79. [53] Wang JG, Liu GR, Wu YG. A point interpolation method for simulating dissipation process of consolidation. Comput Methods Appl Mech Eng 2001;190(45):5907–22. [54] Simo JC, Ju JW. Strain- and stress-based continuum damage models - I. Formulation. Int J Solids Struct 1987;23(7):821–40.

particularly the damage variable against the effective strain. These corresponding characters were studied for the three-point-bending beam test in different nodal discretizations, and then some non-local damage parameters have been optimized, such as s and the weight function order. Consequently, these optimum values were used to analyse the same beam model considering a denser nodal distribution. The obtained numerical results were very close to the experimental solution. Nevertheless, the most interesting analysis is related with the irregular mesh study, in which the problem domain was discretized in several random irregular meshes. The numeric outcome demonstrated a suitable agreement with the similar regular case and the experimental result. Thus, this study proves that the developed meshless damage model is insensitive to the nodal distribution. Acknowledgments The authors truly acknowledge the funding provided by Ministério da Educação e Ciência– Fundação para a Ciência e a Tecnologia (Portugal), under grant PD/BD/114095/2015, SFRH/BPD/75072/ 2010 and SFRH/BPD/111020/2015, and by project funding UID/ EMS/50022/2013. Dr. Moreira acknowledges POPH - QREN-Tipologia 4.2 Promotion of scientific employment funded by the ESF and MCTES. Authors gratefully acknowledge the funding of Project NORTE-010145-FEDER-000022 - SciTech - Science and Technology for Competitive and Sustainable Industries, cofinanced by Programa Operacional Regional do Norte (NORTE2020), through Fundo Europeu de Desenvolvimento Regional (FEDER). References [1] Zienkiewicz OC, Taylor RL. The finite element method, 4th ed. London: McGrawHill Book Company Europe; 1994. [2] Voyiadjis G, Taqieddin Z. Elastic plastic and damage model for concrete materials: Part I-Theoretical formulation. Int J Struct Chang Solids 2009;1(1):31–59. [3] Cervera M, Oliver J, Manzoli O. A Rate-dependent Isotrpic Damage Model for Seismic Analysis of Concrete Dams. Earthq Eng Struct Dyn 1996;25:987–1010. [4] Oliver J, Cervera M, Oller S, Lubliner J. Isotropic Damage Models and Smeared Crack Analysis of Concrete, In: Proceedings of the 2nd Conference on Computer Aided Analysis and Design of Concrete Structures, pp. 945–957; 1990. [5] Cervera M, Oliver J, Faria R. Seismic evaluation of concrete dams via continuum damage models. Earthq Eng Struct Dyn 1995;24(9):1225–45. [6] He W, Wu Y-F, Liew KM, Wu Z. A 2D total strain based constitutive model for predicting the behaviors of concrete structures. Int J Eng Sci . 2006;44(18– 19):1280–303. [7] Yu RC, Ruiz G, Chaves EWV. A comparative study between discrete and continuum models to simulate concrete fracture. Eng Fract Mech 2008;75(1):117–27. [8] Lee J, Fenves GL. A return-mapping algorithm for plastic-damage models: 3-D and plane stress formulation | Jeeho Lee; Gregory L. Fenves | digital library booksc. Int J Numer Methods Eng 2001;50(2):487–506. [9] Tao X, Phillips DV. A simplified isotropic damage model for concrete under bi-axial stress states. Cem Concr Compos 2005;27(6):716–26. [10] Wu JY, Li J, Faria R. An energy release rate-based plastic-damage model for concrete. Int J Solids Struct 2006;43(3–4):583–612. [11] Kachanov LM. Introduction to Continuum Damage Mechanics. Dordrecht, the Netherlands: Martinus Nijhoff Publishers; 1986. [12] Krajcinovic D, Fonseka G. the Continuum Damage Theory of Britle Materials, Part 1: general Theroy. ASME J Appl Mech 1981;48:809–15. [13] Krajcinovic D, Fonseka G. the continuum damage theory for brittle material. ASME J Appl Mech 1983;50:355–60. [14] Resende L, Martin J. A progressive damage continuum model for granular materials. Comput Methods Appl Mech Eng 1984;42:1–18. [15] Crisfield MA. Non-linear Finite Element Analysis of Solids and Structures, 2. Wiley; 1996. [16] Faria R, Oliver J, Cervera M. A strain-based plastic viscous-damage model for massive concrete structures. Int J Solids Struct 1998;35(14):1533–58. [17] Jirásek M. Nonlocal models for damage and fracture: Comparison of approaches. Int J Solids Struct . 1998;35(31–32):4133–45. [18] Aifantis EC. On the Microstructural Origin of Certain Inelastic Models. J Eng Mater Technol . 1984;106(4):326. [19] Schreyer HL, Chen Z. One-Dimensional Softening With Localization. J Appl Mech .

73

Engineering Analysis with Boundary Elements 79 (2017) 62–74

B.V. Farahani et al.

[58] Faria R, Oliver J. A rate dependent plastic-damage constitutive model for large scale computations in concrete structures. CIMNE Monogr 1993(17). [59] Shao JF, Jia Y, Kondo D, Chiarelli AS. A coupled elastoplastic damage model for semi-brittle materials and extension to unsaturated conditions. Mech Mater . 2006;38(3):218–32. [60] Farahani BV, Belinha J, Andrade Pires FM, Ferreira AJM, Moreira PMGP. Extending a radial point interpolation meshless method to non-local constitutive damage models. Theor Appl Fract Mech 2016;no. 85:84–98.

[55] Ghaffarzadeh H, Mansouri A. Numerical modeling of wave propagation using RBFbased meshless method [in]. Future Trends Struct, Civil, Environ Mech Eng 2013. [56] Wu CT, Hu W, Guo Y. the Recent update of LS-DYNA meshfree and advanced FEM for manufacturing application. Livermore Softw Technol Corp LS-DYNA Forum, Bamb 2014. [57] Ghaffarzadeh H. Assessment of effective parameters on Rbf-based meshfree method for wave propagation modeling [in]. Latest Trends Eng Mech Struct Eng Geol 2014.

74