Evaluation of a proposed model for concrete at mesoscopic scale

Abstract

This work deals with numerical modeling of mechanical behavior in quasi-brittle materials, such as concrete. For this propose, a two-dimensional meso-scale model based on RVE existence is presented. The material is considered as a three-phase material consisting of interface zone (ITZ), matrix and inclusions - each constituent modeled by an independent constitutive model. The Representative Volume Element (RVE) consists of inclusions idealized as circular shapes symmetrically and non-symmetrically placed into the specimen. The interface zone is modeled by means of cohesive contact finite elements. The inclusion is modeled as linear elastic and matrix region is considered as elastoplastic material. Our main goal here is to show a computational homogenization-based approach as an alternative to complex macroscopic constitutive models for the mechanical behavior of the brittle materials using a finite element procedure within a purely kinematical multi-scale framework. Besides, the fundamental importance of the representing dissipative phenomena in the interface zone to model the complex microstructural responses of materials like concrete is focused in this work. A set of numerical examples, involving the microcracking processes, is provided in order to illustrate the performance of the proposed modeling.

Keywords:
homogenization; quasi-brittle materials; cohesive contact finite element; concrete; plasticity

Resumo

Este trabalho trata da modelagem numérica do comportamento mecânico em materiais quase-frágeis, tal como o concreto. Para este fim, um modelo 2D de escala mesoscópica baseado no conceito de Elemento de Volume Representativo (EVR) é apresentado. O material é considerado como composto por três fases consistindo de zona de interface, matriz e inclusões, onde cada constituinte é modelado independentemente. O EVR consiste de inclusões idealizadas como de forma circular dispostas de maneira simétrica e não simétrica. A zona de interface é modelada por meio de elementos finitos coesivos de contato. A inclusão é modelada como sendo um material elástico linear, já a matriz é considerada como material elastoplástico. Nosso principal objetivo é mostrar que uma formulação baseada na homogeneização computacional é uma alternativa aos modelos constitutivos macroscópicos complexos para o comportamento mecânico de matérias frágeis usando um procedimento baseado no Método dos Elementos Finitos no âmbito de uma teoria multi-escala. Além disso, o trabalho foca na fundamental importância em representar os fenômenos dissipativos na Zona de Transição para obter uma resposta microestrutural de um material complexo como o concreto. Uma série de exemplos envolvendo processos de microfissuração é apresentada de modo a ilustrar o desempenho da modelagem proposta.

Palavras-chave:
homogeneização; materiais quase-frágeis; elemento finito coesivo; concreto; plasticidade

1. Introduction

Recently, different composite materials have been used in Structural Engineering, as example the concrete. These materials are made by two or more different materials in order to create new materials or to improve the mechanical behavior of well known materials.

In general, as the concrete is a composite material, it presents a very complex mechanical behavior, which is very difficult to be modeled (see Pituba and Fernandes [¹], Brancherie and Ibrahimbegovic [²], Zhu et al. [³] and others). Initially, the constitutive theories called as phenomenological could represent satisfactorily the mechanical behavior of such materials, as example the Continuous Damage Mechanics that can provide sophisticated constitutive models to simulate the mechanical behavior of heterogeneous materials, as the concrete, presenting satisfactory results (see Pituba et al. [⁴], Pereira Jr at el. [⁵], Pituba [⁶]). But in order to better represent the mechanical behavior of such complex materials, these kind of constitutive models require a sophisticated formulation as well as a great number of parameters, sometimes difficult to identify. Then, in the last decades, theories that analyze the structure in different scales have been proposed, where the constitutive response of the material is obtained by analyzing the structure in different lengths and/or time. Therefore, the mechanical behavior of the material at micro-scale is connected to the structural behavior at macro-scale, leading to a more accurate behavior of the continuum, Péric et al. [⁷] and Miehe and Koch [⁸]. In this context, some works had been developed to model the mechanical behavior of the concrete (see Gitman [⁹], Wriggers and Moftah [¹⁰] and López et al. [¹¹])

In summary, the elastic macroscopic response is modeled by a multi-scale constitutive theory based on homogenization techniques, where the stress and strain vectors for a particular point of the macro-continuum are defined as the volumetric average of their respective microscopic fields at the RVE (Representative Volume Element) related to that point of the macro- continuum. In this model, the micro-structure elastic response depends on the choice of kinematical restrictions imposed to the cinematically admissible displacements field of the RVE.

One of the main advantages of multi-scale modeling is that the physical phenomena can be better evaluated, as the properties of each material of the microstructure can be considered as well as the imperfections and voids and then take these informations to a bigger scale of observation. If the analysis is performed at macro level, the concrete behaves as a continuum material, but its microstructure is composed by several materials, presenting different mechanical behaviors.

The phenomena treated by conventional theories, in fact is a macroscopic reflection of what happens at microstructure level. Thus, when analyzing heterogeneous materials, specially the concrete, more efficient constitutive models can be obtained if its microstructure is observed and a multi-scale modeling considered, where adopting simple constitutive models at the microstructure, complex phenomena can be reproduced at the macrostructure, Pituba and Souza Neto [¹²].

In this work only the mechanical behavior of the concrete at microstructure level is considered in order to validate qualitatively the proposed model. The presented formulation is proposed in the context of the multi-scale analysis, where the RVE must be defined as well as homogenization techniques, according to the multi-scale formulation developed in Fernandes et al. [¹³] e Fernandes et al. [¹⁴] in which the Finite Element Method is used in the RVE modeling. In the proposed model the Fracture Mechanics as well as the Plasticity Theory have been considered to model the dissipative phenomena in the interface zone as well as inside the matrix taking into account the geometry and properties of the materials defined at the concrete microstructure. The proposed model is an alternative to the complex phenomenological constitutive models used to represent the behavior of heterogeneous quasi-brittle materials. In future work the proposed model will be coupled to a macro-continuum formulation in order to perform a full coupled multi-scale analysis.

In what follows a brief description of the proposed model for the concrete is presented in section 2, where the homogenization techniques, the constitutive models based on the Contact and Fracture Mechanics, the Mohr-Coulomb model adopted to represent the matrix behavior, as well as a contact-cohesive finite element used to model the interface zone have been discussed. In section 3 numerical examples are analyzed to show the potentialities and limitations of the proposed model. Finally, in section 4, final considerations have been discussed.

2. Proposed model for the concrete

The proposed formulation represents the mechanical behavior of a particular point of the macro-continuum, which can be an integration point of a finite element defined at the macro-structure, as example. By solving the macrostructure problem, this point is subjected to a strain vector which is imposed to its microstructure and then the stresses and constitutive tensor related to that point of the macro-continuum can be computed after solving the microstructure problem. For that, the material microstructure is defined as RVE, whose dimensions are not important, but the distribution and proportionality of the materials composing the microstructure affect its behavior. Adopting concepts of volumetric average and of energy equivalence between the macro and micro-continuum, different values for the homogenized stresses and constitutive tensor can be obtained according to the multi-scale model adopted, which depends on the boundary conditions adopted for the RVE. Note that in this work, the material microstructure is analyzed in the context of multi-scale analysis, where different RVEs subjected to a strain vector have been analyzed, but a full coupled multi-scale analysis of a structure is not presented.

Therefore, to simulate the concrete mechanical behavior, a RVE is used to represent the meso-scale, whose discretization by the Finite Element Method is shown in Figure 1. The aggregates are considered approximately circular, where elastic triangular finite elements are defined, while the matrix can present elastoplastic behavior governed by the Mohr-Coulomb model. Besides, cohesive-contact elements are used to model the interface zone in order to simulate the opening and/or closure of fractures that occur mainly in this region leading to dissipative phenomena during the fracture process of the concrete microstructure. Therefore, in what follows the proposed model is used to simulate the mechanical behavior of the conventional concretes.

Figure 1
Multi-scale analysis scheme

2.1 Proposed formulation for meso-scale

The RVE is considered as a continuum medium, so that the stress concept is valid at micro-scale.. The macroscopic quantities for strain ε(x,t) and stress σ(x,t) at a point x of the macro-continuum are defined as the volumetric average of their respective field ε_μ = ε_μ (y,t) or σ_μ = σ_μ (y,t) over the RVE, considering all points y of the RVE related to the point x. Thus, for an arbitrary instant t the following expressions are defined:

$\mathit{\epsilon }\left(x,t\right)=\frac{1}{V_{\mathit{\mu }}}{\int }_{{\mathrm{\Omega }}_{\mathit{\mu }}}^{}{\mathit{\epsilon }}_{\mathit{\mu }}\left(y,t\right)dV$ (1)

$\mathit{\sigma }\left(x,t\right)=\frac{1}{V_{\mathit{\mu }}}{\int }_{{\mathrm{\Omega }}_{\mathit{\mu }}}^{}{\mathit{\sigma }}_{\mathit{\mu }}\left(y,t\right)dV$ (2)

Equations (1) and (2) represent the macroscopic or homogenized values for strain and stress, as a microscopic filed have been transformed into a macroscopic quantity by means of a homogenization technique. Besides, the microscopic stress can be written in terms of the microscopic strain, as follows:

${\mathit{\sigma }}_{\mathit{\mu }}\left(y,t\right)=f_y\left({\mathit{\epsilon }}_{\mathit{\mu }}\left(y,t\right)\right)$ (3)

where f_y is the constitutive functional, defined in this work by the Mohr-Coulomb model. Moreover, the microscopic strain ε_μ can be written in terms of the microscopic displacement filed u_μ of the RVE, as follows:

${\mathit{\epsilon }}_{\mathit{\mu }}\left(y,t\right)={\nabla }^S{u}_{\mathit{\mu }}(y,t)$ (4)

where ${\nabla }^S$ is the symmetric gradient operator of the displacement field u.

Without loss of generality, the microscopic displacement filed u _μ can be defined as the sum of three parts:

$u_{\mu }\left(y,t\right)=u\left(x,t\right)+\overline{u_{\mu }}\left(y,t\right){+\tilde{u}}_{\mu }\left(y,t\right)$ (5)

being the first one constant representing a rigid body motion coincident to the macroscopic displacement u (x, t) related to the point x, the second one is obtained from the macroscopic strain ε as follows:

$\overline{u_{\mathit{\mu }}}\left(y,t\right)=\mathit{\epsilon }\left(x,t\right)y$ (6)

which varies linearly with the coordinate y, and a displacement fluctuation field _{^{ũμ (y, t)}} . Thus, Equation (5) can be written as:

$u_{\mathit{\mu }}\left(y,t\right)=\mathit{\epsilon }\left(x,t\right)y{+\tilde{u}}_{\mathit{\mu }}\left(y,t\right)$ (7)

In Equation (7) the part εy varies linearly with y resulting from the multiplication of the macroscopic strain ε of the EVR, which is constant, by the coordinates of the point y. In the case of having uniform microscopic displacement ε_μ, the displacement fluctuation ũ_μ is null. In the RVE the following relations for the microscopic strain ε_μ and the microscopic strain fluctuation ${\mathit{\epsilon }}_{\mathit{\mu }}$ have to be satisfied:

${\mathit{\epsilon }}_{\mathit{\mu }}={\nabla }^S{u}_{\mathit{\mu }}(y,t)$ (8)

${\tilde{\mathit{\epsilon }}}_{\mu }={\nabla }^S{\tilde{u}}_{\mathit{\mu }}(y,t)$ (9)

Considering Equations (7) to (9) the microscopic strain can also be written as:

${\mathit{\epsilon }}_{\mathit{\mu }}\left(y,t\right)=\mathit{\epsilon }\left(x,t\right){+\tilde{\mathit{\epsilon }}}_{\mathit{\mu }}\left(y,t\right)$ (10)

After some manipulations (Fernandes et al. [¹³]), Equation (10) can be written in terms of velocity, where a microscopic strain velocity is cinematically admissible if:

${\dot{\mathit{\epsilon }}}_{\mathit{\mu }}\left(y,t\right)={\nabla }^s{\dot{u}}_{\mathit{\mu }}=\dot{\mathit{\epsilon }}\left(x,t\right)+\dot{{\mathit{\epsilon }}_{\mathit{\mu }}}\left(y,t\right)\forall {\dot{u}}_{\mathit{\mu }}\in v_{\mathit{\mu }}$ (11)

where ν_µ is the space of cinematically admissible displacements of the RVE. More details can be found in Fernandes et al. [¹³].

As already mentioned, the microscale is represented by the RVE, being the FEM the numerical method used to solve the RVE equilibrium problem. The RVE characteristics as its dimensions, elastic properties and constitutive models adopted for the materials are the same for all RVEs defined at macroscale. The RVE solution is obtained, i.e., the computation of displacements, internal forces, stress and constitutive tensors, for all finite elements, are obtained when the convergence of the equilibrium problem is achieved according to the adopted tolerance. But in order to solve the RVE equilibrium problem boundary conditions in terms of displacement fluctuations must be imposed to the RVE Then, the numerical response can varies according to the boundary condition adopted. Thus, in what follows, the following topics will be discussed: RVE equilibrium problem, Hill-Mandel Principle, Stress homogenization, Boundary conditions imposed to the RVE, Homogenized Constitutive Tangent Tensor.

To simplify the model the inclusions domain ${\mathrm{\Omega }}_{\mathit{\mu }}^i$ and the matrix domain ${\mathrm{\Omega }}_{\mathit{\mu }}^m$ will be considered together as the solid domain ${\mathrm{\Omega }}_{\mathrm{\mu }}^{\mathrm{s}}$. Neglecting the inertia forces and that the RVE is subjected to the body force b = b (y, t) and to surface force field te = te (y, t) acting along the boundary, the Principle of Virtual Displacements establishes that the RVE is in equilibrium if, and only if, the stress field σμ over Ωμ satisfies the classic variational equation of the elasticity:

$\underset{{\text{Ω}}_{\mu }^S}{\overset{}{{\displaystyle \int }}}{\sigma }_{\mu }\left(y,t\right):{\nabla }^{S}dV-\underset{{\text{Ω}}_{\mu }^S}{\overset{}{{\displaystyle \int }}}b\left(y,t\right).dV+\underset{{\text{Ω}}_{\mu }^v}{\overset{}{{\displaystyle \int }}}{\sigma }_{\mu }\left(y,t\right):{\nabla }^{S}dV-\underset{{\text{Ω}}_{\mu }^v}{\overset{}{{\displaystyle \int }}}b\left(y,t\right).dV-\underset{\partial {\text{Ω}}_{\mu }^{}}{\overset{}{{\displaystyle \int }}}{t}^e\left(y,t\right).dA=0\forall \in v_{\mu }$ (12)

The works of Hill and Mandel (Giusti et al. [¹⁵]) have established the Macro Homogeneity Principle which defines that the macroscopic stress power in a arbitrary point of the macrocontinuum must be equal to the volumetric average of the microscopic stress power over the RVE related to that point for any movement cinematically admissible of the RVE. (Giusti et al. [¹⁵]). Considering Equation (11), assuming ${\dot{\tilde{u}}}_{\mathit{\mu }}=$ and considering that the voids are in equilibrium, after some manipulations (see details in Fernandes et al. [¹³] and Fernandes et al. [¹⁴]), we can conclude that the Hill-Mandel principle is valid if, and only if, the following integrals are nulls:

$\underset{{\partial \mathrm{\Omega }}_{\mu }^{}}{\overset{}{\int }}{t}^e\left(y,t\right).{\dot{\tilde{u}}}_{\mathit{\mu }}dA=0\forall {\dot{\tilde{u}}}_{\mathit{\mu }}\in {\mathit{v}}_{\mathit{\mu }}$ (13)

$\underset{{\partial \mathrm{\Omega }}_{\mu }^S}{\overset{}{\int }}b\left(y,t\right).{\dot{\tilde{u}}}_{\mathit{\mu }}dV=0\forall {\dot{\tilde{u}}}_{\mathit{\mu }}\in {\mathit{v}}_{\mathit{\mu }}$ (14)

Considering Equations (8) and (10) and writing σ_μ σ_μ as σ_μ = f_y (ε_μ), where ^_fy is the constitutive functional, the following Equation in terms of displacement fluctuation can be obtained to represent the equilibrium problem of the solid part of the RVE:

$\underset{{\mathrm{\Omega }}_{\mathit{\mu }}^S}{\overset{}{\int }}{f}_y(\mathit{\epsilon }\left(x,t\right)+{\nabla }^s\widetilde{u_{\mathit{\mu }}}\left(y,t\right)):{\nabla }^{S}dV=0\forall \in {\mathit{v}}_{\mathit{\mu }}$ (15)

Finally the formulation is completed by the appropriated choice of the space ${\mathit{v}}_{\mathit{\mu }}$, with the choice of the kinematical restrictions to be imposed to the RVE. Thus, the microscopic equilibrium problem consists of, given the macroscopic strain tensor ε, finding the field ${{\tilde{u}}_{\mathit{\mu }}\in v}_{\mathit{\mu }}$ such that for each instant t, the Equation (15) is satisfied. As η is an arbitrary field η, after the RVE domain discretization into finite elements, whose domain is referred as ${\mathrm{\Omega }}_{\mu }^h$, the following incremental microscopic equilibrium equation must hold for a load increment in time $\Delta t_n=t_{n+1}-t_n$ and a domain discretization h, finding the displacement fluctuation ${\tilde{u}}_{\mathit{\mu }(n+1)}={\tilde{u}}_{\mathit{\mu }\left(n\right)}+{∆\tilde{u}}_{\mathit{\mu }\left(n\right)}$:

$G_h^{n+1}={\int }_{{\mathrm{\Omega }}_{\mathit{\mu }}^h}^{}{B}^T{f}_y({\mathit{\epsilon }}_{n+1}+B{\tilde{u}}_{\mathit{\mu }\left(n+1\right)})dV=0$ (16)

where B is the global matrix relating strain and displacement, ${\mathrm{\Omega }}_{\mu }^h$ is the RVE discretized domain. If the load increment is non-linear, Equation (16) is solved by applying the Newton-Raphson Method which consists of finding the fluctuation correction $\mathit{}\mathit{\delta }{\tilde{\mathit{u}}}_{\mathit{\mu }}^{i+1}$ for iteration i+1, such that:

$F^i+K^i\mathit{\delta }{\tilde{u}}_{\mathit{\mu }}^{i+1}=0$ (17)

where F is the force vector and K the tangent rigidity matrix of the RVE. After computing the correction $\mathit{\delta }{\tilde{u}}_{\mathit{\mu }}^{i+1}$ defined in Equation (17), the next step is to obtain the displacement fluctuation field to be considered at iteration i+1 given by: ${{\tilde{u}}_{\mathit{\mu }}^{i+1}=\tilde{u}}_{\mathit{\mu }}^i+\mathit{\delta }$.

The homogenized stress is computed from Equation (2), considering that the RVE is composed by voids and a solid part (matrix and agregates) ${\mathrm{\Omega }}_{\mu }={\mathrm{\Omega }}_{\mu }^S\cup {\mathrm{\Omega }}_{\mu }^v$, resulting into:

$\mathit{\sigma }=\mathit{\sigma }\left(x,t\right)=\frac{1}{{\mathrm{V}}_{\mathit{\mu }}}\underset{{\mathrm{\Omega }}_{\mathit{\mu }}^S}{\overset{}{\int }}{\mathit{\sigma }}_{\mathit{\mu }}\left(y,t\right)dV+\frac{1}{V_{\mathit{\mu }}}\underset{{\mathrm{\Omega }}_{\mathit{\mu }}^v}{\overset{}{\int }}{\mathit{\sigma }}_{\mathit{\mu }}\left(y,t\right)dV$ (18)

The RVE equilibrium problem is completed with the choice of the kinematical restrictions to be imposed to the RVE, leading to different classes of multi-scale models and consequently to different numerical results (Peric et al. [⁷]). In this work only periodic displacement fluctuations is considered. For that, each RVE side ${\mathrm{\Gamma }}_i^{+}$ whose normal direction is $n_i^{+}$, must correspond to an equal side ${\mathrm{\Gamma }}_i^{-}$ with normal direction $n_i^{-}$, being $n_i^{+}={-n}_i^{-}$. Similarly, for each point ${,y}^{+}$ defined on ${\mathrm{\Gamma }}_i^{+}$ must exist a point ${,y}^{-}$ on the side ${\mathrm{\Gamma }}_i^{-}$.To have periodic displacement fluctuation on the boundary, for every pair of points ${(y}^{+}{,y}^{-})$ the following relation must be verified:

${\tilde{u}}_{\mathit{\mu }}\left(y^{+},t\right)={\tilde{u}}_{\mathit{\mu }}\left(y^{-},t\right)\forall \left\{y^{+}{,y}^{-}\right\}\in \partial {\mathrm{\Omega }}_{\mathit{\mu }}$ (19)

2.2 Cohesive fracture model and plasticity

Pituba et al. [¹²] have proposed an extension of a cohesive fracture law presented in Cirak et al. [¹⁶] in order to deal with damage process leading to the complete failure of microstructures in ductile media. In general way, this model has been developed to represent the cracking process where traction is still possible to be transmitted between fracture lips. The proposed model describes the finite-deformation irreversible cohesive law. The cohesive free energy is given by:

$\mathrm{\Phi }=\mathrm{\Phi }\left({\delta }_n,{\delta }_s,\mathit{q}\right)$ (20)

Where, δn is the normal opening displacement due to mode I; δs is the sliding opening displacement due to mode II and q is the internal variable that describes the inelastic processes related to decohesion.

It is possible to assume that the deformation due to sliding opening process is a scalar value independent of the direction of sliding on the cohesive surface, thus δs = |δs|, therefore the behavior has an isotropic characteristic and the cohesive law is written introducing an effective opening displacement expressed by:

$\delta =\sqrt{{\beta }^2{\delta }_S^2+{\delta }_n^2}$ (21)

The parameter β assumes different values (from 0 to 1) to the sliding and normal opening displacements given a weight ratio between the sliding and normal directions. On the other hand, the φ free energy potential depends of δ, and the cohesive law is expressed as:

$\mathit{t}=\frac{t}{\delta }({\beta }^2{\mathit{\delta }}_{\mathit{S}}+{\delta }_n\mathit{n})$ (22)

Where, n is the unit normal to the cohesive surface; δs is the sliding opening vector located on the cohesive surface, t is the cohesive traction on the crack; t is a scalar effective traction.

On the other hand, the released cohesive energy in the microstructure of the material proposed in this work (Equation (20)) is given by:

$\mathrm{\Phi }=e{\sigma }_c{\delta }_c\left[1-e^{⌊-\left(1+\frac{\delta }{{\delta }_c}\right)⌋}\right]$ (23)

Where the law for the scalar effective traction for the loading cases is obtained from Equation (23) as:

$t=\frac{\partial \varphi }{\partial \delta }={\sigma }_c{e}^{-\delta /{\delta }_c}if\delta ={\delta }_{max}and\dot{\delta }\ge 0$ (24)

For the scalar effective traction for the unloading cases is proposed a law considering an elastic behavior, i. e., without residual effective opening displacement as follows:

$t=\frac{t_{max}}{{\delta }_{max}}\delta if\delta <{\delta }_{max}or\dot{\delta }<0$ (25)

Where e is the e-number, σ_c is the maximum tension cohesive normal traction and δc is a characteristic opening displacement that indicates a critical opening. Thus, β, σ_c and δ_c are parameters of the cohesive model. Besides, $\dot{\delta }$ is the opening displacement rate. The effective stress versus effective opening displacement cohesive laws for loading and unloading cases are presented in Figure 2a and Figure 2b, respectively.

Figure 2
Cohesive law written in terms of effective stress versus effective opening displacement: a) Loading case; b) Unloading case

Accordingly to Ortiz and Pandolfi [¹⁷], there is a relation between the cohesive law and the critical energy released rate (_^GC ) for crack propagation in the microstructure in our case. Assuming the direction 1 as the direction on the fracture surface and towards to the its propagation, _^GC can be written as:

$G_c=\underset{0}{\overset{R}{\int }}t.{\delta }_{,1}d{x}_1$ (26)

Where R is the cohesive zone length. The Equation (26) can also be defined as:

$G_c=\underset{0}{\overset{\infty }{\int }}t.{\delta }_{,1}d{x}_1={\varphi }_{\infty }$ (27)

For the cohesive law presented in this work, using Equation (24), the critcial energy released rate is given by:

$G_c=e{\sigma }_c{\delta }_c$ (28)

Obviously, the _^GC for conventional modeling is developed with phenomenological constitutive models applied on the homogenized macrocontinuum. In the present work, the concept of fracture energy is closely related to that which occurs in the ITZ of the microstructure at mesoscale level. When there are nucleation and growth processes of the microcracks, it may generate a localization process leading to a nucleation of fracture on the macrocontinuum. This process leads to a relation between fracture energies on macrocontinuum (homogenized material) and on mesoscale of the material, as described in Oliver et al. [¹⁸].

Before the nucleation process, it is assumed the existence of stiffness between the lips of the future fracture situated between triangule finite elements. This stiffness is simulated by another parameter of the proposed model called penalty factor (λ_p ). In a practical view, high values for this parameter are adopted in order to obtain a accurate approach. This procedure ensures that the future fracture be kept closed until the separation criterion is reached and, at the same time, guarantees the physical admissibility of the entire process. The penalty factor is, therefore, a stiffness imposed to the closuring of the crack.

In general way, that strategy intends to create stiffness between the nodes of the embedded cohesive contact finite elements in the matrix zone in order not to allow penetration of the surfaces of the fracture. On the other hand, in tension regimes, this penalty factor effectively replaces the initial rigid part of the cohesive law for a linear response given by Equation (26). In order to detect the cohesive contact phenomenon, the concept of the opening displacement gap between the Gauss points of the cohesive contact finite element is adopted.

$t={\lambda }_p\delta if{\lambda }_p\delta \le {\sigma }_c$ (29)

The finite elements used in this work are composed of two surfaces which are coincident in the undeformed configuration of the RVE. The cohesive contact finite element is defined as an element with four nodes and its geometry is compatible with the two triangule finite elements used to model the matrix and aggregate zones. The formulation of the cohesive contact finite element is presented in Pituba and Souza Neto [¹²] and Pituba et al. [¹⁹].

On the other hand, to deal with plastic strains presented on the macrostructure of the concrete when the material is subjected mainly to compression stress, the well-know Mohr-Coulomb model is used to represent the mechanical behavior of the cement matrix. Therefore, this is another dissipative process to be modeled in the microstructure of the material together the possible microcracking nucleation that occurs mainly in the ITZ. In case of predominant tension regime, the proposed modeling evidences the microcracking process in the ITZ as the most important dissipative phenomenon, mainly in the initial loading stages. This phenomenon is simulated by the cohesive contact finite elements placed in the ITZ. Obviously, the microcracking process in the matrix zone that occurs in the softening regime of the macrostructure of the material is also important leading to the need of the insertion of cohesive contact finite elements in the matrix zone. However, this embedded finite elements can generates numerical instabilities, mainly in the peak stress regime. On the other hand, in predominant compression regimes, the yielding process of the matrix is understood as principal phenomenon in conjunction with the microcracking process in the ITZ.

The modeling of the microcracking process in the ITZ is fundamental to estimate the collapse of the concrete microstructure, mainly in predominant tension regimes. This assertion is based on results obtained by Pituba and Souza Neto [¹²], even considering elastic behavior for the matrix and aggregates, the fracture process in the ITZ modeled by cohesive contact finite elements together with the geometry of the aggregates have allowed the simulation of homogenized plastic macroscopic strains in unloading and reversal loadings situations, leading to the capture of the unilateral effect of the concrete.

The Mohr-Coulomb model as well as the contact and cohesive fracture model has been implemented on the computational code developed to analyze RVEs of heterogeneous materials submitted to a state of macroscopic strain. The macroscopic strain is divided in increments and imposed to the RVE and the stresses are corrected by the constitutive models. In case of Mohr-Coulomb model, with the stress values it is possible to update the yielding surface of the constitutive model. More details can be obtained in Souza Neto et al. [²⁰].

On the other hand, the equations proposed are valid to the three-dimensional cases. For plane stress cases, Souza Neto et al. [²⁰] suggests a strategy to approach to the stress plane case using constitutive equations of the plane stress state restraining the stresses in Gauss points of the finite element. As the stresses σ₁₃ and σ₂₃ are nulls, it is added σ₃₃ in the computational implementation, where its value has to be zero, considering the direction 3 normal to the middle surface of the finite element. However, σ₃₃ is related to the strain ε₃₃ that depends of others stresses on the plane of the finite element. Thus, for a given initial elastic strain (ε_^e)33 , an algorithm for the axisymmetric case of stresses is applied in order to obtain the value for σ₃₃ . If the convergence of the process is achieved (|σ₃₃ |≤є_tol), the updating of the stress is valid, otherwise it is necessary to update the strain (ε_^e)33 . The correction procedure is performed using the Newton-Raphson algorithm and a self-consistent linearization of the problem leading to a convergence with quadratic rate for the numerical solution. This is a very interesting aspect because in multi-scale analyses the computational cost has a high value.

3. Results and discussions

In order to evaluate the application of the proposed formulation in quasi-brittle materials as concrete, some numerical simulations based on mechanical behavior of the microstructure of the concrete are performed. Initially, RVEs with dimensions l x l and thickness l/10 containing inclusions placed into a matrix are generated. Obviously, the inclusions (aggregates) and matrix have different mechanical properties evidencing the heterogeneous characteristic of the medium submitted to plane stress states. For the matrix zone, an elastopplastic behavior is assumed following the Mohr-Coulomb criterion with the parameter values given by, Assad et al. [²¹]: Young´s modulus E is 20 GPa and Poisson ration ν is 0.2, friction angle and dilatation angle are φ = 5º and Ψ = 10%, respectively. The aggregates are considered elastic media with E = 35 GPa and ν = 0.26, Mehta and Monteiro [²²].

For the ITZ, in the situations where the fracture process has been evaluated, cohesive contact finite elements have been used. The parameters for the cohesive law are given by: λ_{^{p = 200000 N/mm3}} , β = 0.7, σc=0.09 MPa and δc = 0.02mm, Ortiz e Pandolfi [¹⁷], Oliver et al. [¹⁸], Pituba e Souza Neto [¹²].

3.1 Influence of the fracture process in the ITZ

In this section, RVEs containing 4 and 8 inclusions with circular shape placed in the matrix zone are analyzed. Figure 3 presents the RVE with 4 inclusions that represent 12% of volume fraction. For the RVE discretization, 798 triangule finite elements and 64 cohesive contact finite elements are used when included in the analyses. Besides, Figure 4 presents the RVE with 8 inclusions that represent 14% of volume fraction. For the RVE discretization, 1184 triangule finite elements and 128 cohesive contact finite elements are used when included in the analyses. In other analysis, perfect bonding between aggregates and matrix has been considered.

Figure 3
Representative volume element with 4 inclusions

Figure 4
Representative volume element with 8 inclusions

The RVEs have been submitted to a total macroscopic strain ε_{^{x = 0,0001}} e ε_{^{y= -0,00001}} divided in 20 increments. The distortional strain has been considered null. The homogenized stresses obtained in the analyses for x-direction versus the macroscopic strain in the same direction are plotted in Figures 5 and 6. Note that the numerical responses considering perfect bonding or fracture process in the ITZ are the same in the initial steps of loading. Nevertheless, when the dissipative processes take place those numerical responses diverge. The consideration of fracture process in the ITZ evidences an important contribution in the non-linear behavior of the material at microscale level. Therefore, the cohesive fracture model decreases the stiffness of the RVE beyond the yielding limit. The numerical responses for both RVEs are quite similar due to the very close values for the volumetric fractions of inclusions.

Figure 5
Tension normal homogenized stress in the x-direction versus elongation macroscopic strain imposed in the x-direction of the RVE with 4 inclusions

Figure 6
Tension normal homogenized stress in the x-direction versus elongation macroscopic strain imposed in the x-direction of the RVE with 8 inclusions

In order to visualize the stress distribution in the x-direction inside the RVE, Figures 7 and 8 illustrate the impact on the numerical response when considering the fracture process in the ITZ. Figure 7 shows the stress distribution for the RVE containing 4 inclusions considering perfect bonding between inclusions and matrix whereas Figure 8 represents the RVE considering the fracture process in the ITZ. Note in Figure 7 that the stresses are perfectly transmitted submitting high levels of stresses at the inclusions. Therefore, as the aggregates present stiffness values higher than the cimenticious matrix, they increase the homogenized stiffness of the material, see Figure 5.

Figure 7
Normal stress distribution in the x-direction inside the RVE with 4 inclusions without considering the fracturing process in the ITZ

Figure 8
Normal stress distribution in the x-direction inside the RVE with 4 inclusions considering fracturing process in the IT

On the other hand, the stress values in the aggregates presented in Figure 8 are smaller when compared to the Figure 7 due to the dissipative processes presented in the cohesive contact finite elements. Thus, the impact of the aggregates in the homogenized stiffness of the material is smaller, as shown in Figure 5.

3.2 Influence of the boundary conditions

On the Multi-scale formulation, the homogenized response is influenced by the boundary conditions as presented in section 2.1. For the numerical analysis presented here, a new RVE model proposed by Nguyen et al. [²³] is used. This RVE presents a randomly distribution of aggregates with various dimensions, as shown in Figure 9. The use of this RVE is due to no existence of symmetry in any axes. This property will be important in the analyses of the section 3.3. The volume fraction of inclusions in this RVE is 35.6%. For the RVE discretization, 520 triangule finite elements and 95 cohesive contact finite elements are used when included in the analyses.

Figure 9
Representative volume elements proposed by Nguyen et al. [²³]

Figure 10 shows the homogenized responses of stress in x-direction when using linear and periodic boundary conditions on the RVE with 4 inclusions presented in Figure 4. Figure 11 presents the homogenized responses for the RVE shown in Figure 9. Both analyses consider the fracture process in the ITZ.

Figure 10
Tension normal homogenized stress in the x-direction versus elongation macroscopic strain imposed in the x-direction of RVE with 4 inclusions considering different boundary conditions

Figure 11
Tension normal homogenized stress in the x-direction versus elongation macroscopic strain imposed in the x-direction of the RVE presented in Nguyen et al. [²³] considering different boundary conditions

The use of linear and periodic boundary conditions did not influence the responses presented in Figure 10 due to the symmetric disposition of aggregates in the RVE. But, when the RVE does not present symmetry, as shown in Figure 9, the homogenized response is influenced by the fracture process in the ITZ. This influence has not significance at the microscopic level. Nevertheless, this influence has to be very important when performing multi-scale analysis.

3.3 Anisotropic behaviour

In order to check how the proposed modeling is capable to obtain complex responses on the macrostructure, as the anisotropic behavior, as instance, a set of analyses has been performed using the RVEs described in Figures 3 and 4. Initially, the macro strains ε_{^{x = 0.0001}} and ε_{^{y= - 0.00001}} have been imposed. Soon after, ε_{^{y = 0.0001}} e ε_{^{x = - 0.00001}} have been applied. The numerical results are presented in Figures 12 and 13 expressed by homogenized stress versus macro strains. In Figure 12 is presented the numerical results for RVE with 4 inclusions and Figure 13 presents the numerical results for RVE with 8 inclusions.

Figure 12
Tension normal homogenized stress in the x-direction versus elongation macroscopic strain imposed in the RVE with 4 inclusions

Figure 12 shows that the mechanical behavior in different directions has been the same. This happens because the RVE has symmetry related to the axes x and y. But in Figure 13, the mechanical behavior in different directions presents significative differences when the fracture nucleation in the ITZ takes place. In this last case, the non-symmetric distribution of the inclusions in the RVE is responsible for this phenomenon. Therefore, the geometry of the inclusions and the consideration of the fracture process in the ITZ lead to simulate a complex macroscopic phenomenon using simple constitutive models at mesoscale of the material.

Figure 13
Tension normal homogenized stress in the x-direction versus elongation macroscopic strain imposed in the RVE with 8 inclusions

3.4 Homogenized responses in tension and compression regimes

In this section is evaluated the fracture process in the ITZ for RVEs in predominant tension and compression regimes. The RVE presented in Figure 9 is used. The first analysis consists in the application of the macro strain simulating a predominant tension regime given by: ε_{^{x = 0.0001}} and ε_{^{y = - 0.00001}} , where γ_xy is approximately null. After that, a predominant compression regime has been considered with the same magnitude before. Therefore, the macro strains are: ε_{^{x = - 0.0001}} and ε_{^{y = 0.00001}} . The Figures 14 and 15 present results of the homogenized stress in x-direction versus imposed macro strain in x-direction.

Figure 14
Tension normal homogenized stress in the x-direction versus elongation macroscopic strain imposed in the x-direction of the EVR of Nguyen on the tension predominant regime

Figure 15
Compression normal homogenized stress in the x-direction versus compression macroscopic strain imposed in the x-direction of the EVR of Nguyen on the compression predominant regime

Figure 14 shows the important contribution of the fracture process in predominant tension regime leading to a decreasing of the homogenized stiffness and strength of the material. Figure 15 shows that the fracture process is not so important in predominant compression regimes. This assertion is based on the mechanical behavior of the cohesive contact finite elements placed in the ITZ. The stress transmission lost between the surfaces of the fracture is evident when the value of the cohesive traction increases, but this is not happen in many cohesive contact finite elements in RVE submitted to predominant compression regime. Therefore, many cohesive contact finite elements are submitted to compression loading conditions, where the contact law is activated. For the visualization proposes, Figure 15 presents positive signals for compressive homogenized stress and macro strains.

3.5 Alternative modeling of the ITZ

Some authors consider the ITZ as a low resistance zone surrounding the aggregates with mechanical properties less stiffened than the mechanical properties of the matrix. For instance, Ramesh [²⁴] and Yang [²⁵] consider a average value for the Young modulus of the ITZ as a percentage value of the Young modulus for the matrix. In this section, a RVE containing one inclusion placed on the central region is analyzed, see Figure 15. The volume fraction of inclusions is 19.6%. For the RVE considering the ITZ modeled by cohesive contact finite elements, 612 triangule finite elements and 36 cohesive contact finite elements have been used in the analysis. For the RVE considering a band for the ITZ, 612 triangule finite elements have been used.

The homogenized results for both proposed modeling have been compared. For the thickness of the band zone, a value of 40 μm has been adopted following Mehta and Monteiro [²²]. The authors suggest values between 20 and 50 μm.

The same values used in the last sections for the parameters of the constitutive models are used here. The macroscopic strain imposed to the RVE is given by: ε_{^{x = 0.0001}} and ε_{^{y = - 0.00001}} , where γ_xy is approximately null. For the green region of the Figure 16 a Young modulus with 50% of the value given to matrix has been used. Figure 17 presents the homogenized stress in x-direction versus macroscopic strain imposed in the x-direction.

Figure 16
ITZ modeled by a triangular finite element band

Figure 17
Tension normal homogenized stress in the x-direction versus elongation macroscopic strain imposed in the x-direction of RVE with 1 inclusion

Note that model containing a band for the ITZ initially presents a lower homogenized stiffness when compared to the model with ITZ modeled by cohesive contact finite elements. However, when the loading increases, after yielding and fracture processes have begun, the EVR with cohesive elements presents a lower homogenized stiffness. This is due to the consideration of cohesive contact finite elements which plays a important role in collapse regimes.

4. Conclusions

In this work a model to simulate the concrete mechanic behavior, in the context of multiscale modeling, has been proposed, being the numerical analyses restricted to the material mesostructure. Using simple constitutive models and defining accordingly the geometry of the different phases of the RVE, some important macroscopic phenomena have been represented.

By considering cohesive-contat finite elements at the concrete transition zone, the loss of rigidity in the homogenized response could be evidenced without presenting total loss of the stress propagation, as expected. As this kind of element does not influence the homogenized response while the element opening is not activated, this model proved to be viable, if compared to others strategies to simulate the effect of the concrete transition zone, as the definition of less rigid finite elements around the agregates.

On the other hand, the consideration of the plasticity model in the multi-scale analysis proved to be satisfactory to overcome the problem of having rigid responses for predominant compression regime, as discussed in Pituba and Souza Neto [¹²]. But considering only the plasticity model for the mortar without defining the cohesive-contact finite element is not satisfactory, because as the aggregates are assumed to have elastic behavior, they present bigger rigidity modifying the homogenized response. Therefore, in this work is shown the importance of considering the dissipative phenomena at Transition Zone for better represent the concrete mechanical behavior. On the other hand, a limitation of the proposed model is not being able to reproduce the softening behavior for predominant tension regimes when occurs the microcracking process inside the mortar. Besides, the definition of cohesive-contact elements in the matrix can lead to instabilities of the numerical response for predominant tension regimes. To overcome these difficulties, in a future work these cohesive-contact elements could be replaced for the high aspect ratio elements developed by Rodrigues et al. [²⁶]. But for concrete structures, in service regimes, the proposed model has shown to be a viable tool to perform multi-scale analysis of structures. The authors believe that this strategy will be able to represent the softening behavior for tension regimes.

Although the proposed model presents some limitations, in the numerical examples analyzed in this work, it has represented very well the concrete mechanical behavior, as it has captured complex phenomena by adopting simple constitutive models, what encourage us to proceed with this research. The proposed model will be considered for identification of quantitative responses for the concrete, as well as for full coupled multi-scale analyses of concrete structures, based on the works developed in Fernandes et al. [¹³] and Fernandes et al. [¹⁴]. Moreover, it is important the development of a formulation considering the local damage phenomena at microstructure, what can lead to a fracture at macrostructure, as discussed in Sanchéz et al. [²⁷].

5. Acknowledgments

The authors wish to thank CAPES Foundation (Ministry of Education of Brazil) and CNPq (National Council for Scientific and Technological Development).

6. References

• ^[1]
  PITUBA, J. J. C.; FERNANDES, G. R.. An anisotropic damage for the concrete. Journal of Engineering Mechanics - ASCE, v. 137, 2011; p. 610-624.
• ^[2]
  BRANCHERIE, D.; IBRAHIMBEGOVIC, A. Novel anisotrtopic continuum-discrete damage model capable of representing localized failure of massive structures. Part I: theoretical formulation and numerical implementation. Engineering Computations, v. 26, 2009; p. 100-127.
• ^[3]
  ZHU, Q.; KONDO, D.; SHAO, J.; PENSEE, V. Micromechanical modelling of anisotropic damage in brittle rocks and application. International Journal of Rock Mechanics and Mining, v. 45, 2008; p. 467-477.
• ^[4]
  PITUBA, J. J. C.; DELALIBERA, R. G.; RODRIGUES, F. S. Numerical and statistical analysis about displacements in reinforced concrete beams using damage mechanics. Computers and Concrete, an International Journal, v. 10, n. 3, 2012; p. 307-330.
• ^[5]
  PEREIRA JR, W. M.; ARAÚJO, D. L.; PITUBA, J. J. C. Numerical analysis of steel-fiber-reinforced concrete beams using damage mechanics. Ibracon Structures and Materials Journal, v. 9, n. 2, 2016; pp. 153-191.
• ^[6]
  PITUBA, J. J. C. A damage model formulation: unilateral effect and RC structures analysis. Computers and Concrete, an International Journal, v. 15, n. 5, 2015; p. 709-733.
• ^[7]
  PERIC, D.; DE SOUZA NETO, E. A.; FEIJÓO, R. A.; PARTOVI, M.; CARNEIRO MOLINA, A. J. On micro-to-macro transitions for multiscale analysis of heterogeneous materials: unified variational basis and finite element Implementation. International Journal for Numerical Methods in Engineering, v. 87, 2011; p. 149-170.
• ^[8]
  MIEHE, C.; KOCH, A. Computational micro-to-macro transitions of discretized microstructures undergoing small strains. Archive of Applied Mechanics, v. 72, 2002; p. 300-317.
• ^[9]
  GITMAN, I. M. Representative volumes and multi-scale modeling of quasi-brittle materials. PhD Thesis, Technische Universiteit Delft, 2006.
• ^[10]
  WRIGGERS, P.; MOFTAH, S. O. Mesoscale models for concrete: homogenization and damage behavior. Finite Elements in Analysis and Design, v. 42, 2006; p. 623-636.
• ^[11]
  LÓPEZ, C. M.; CAROL, I.; AGUADO, A. Meso-structural study of concrete fracture using interface elements. I: numerical model and tensile behavior. Material and Structures, v. 41, 2008; p. 583-599.
• ^[12]
  PITUBA, J. J. C.; SOUZA NETO, E. A. Modeling of unilateral effect in brittle materials by a mesoscopic scale approach. Computers and Concrete, An International Journal, v.15, 2015; p. 1-25.
• ^[13]
  FERNANDES, G. R.; PITUBA J. J. C.; DE SOUZA NETO, E. A. Multi-scale modelling for bending analysis of heteregeneous plates by coupling BEM and FEM. Engineering Analysis with Boundary Elements, v. 51, 2015; p. 1-13.
• ^[14]
  FERNANDES, G. R.; PITUBA J. J. C.; DE SOUZA NETO, E. A. FEM/BEM formulation for multi-scale analysis of stretched plates. Engineering Analysis with Boundary Elements, v. 54, 2015; p. 47-59.
• ^[15]
  GIUSTI, S. M.; BLANCO, P. J.; DE SOUZA NETO, E. A.; FEIJÓO, R. A. An assessment of the Gurson yield criterion by a computational multi-scale approach. Engineering Computations, v. 26, 2009; p. 281-301.
• ^[16]
  CIRAK, F.; ORTIZ, M.; PANDOLFI, A.. A cohesive approach to thin-shell fracture and fragmentation. Computer Methods in Applied Mechanics Engineering, v. 194, 2005; p. 2604-2618.
• ^[17]
  ORTIZ, M.; PANDOLFI, A. Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis. International Journal for Numerical Methods in Engineering, v. 44, 1999; p. 1267-1282.
• ^[18]
  OLIVER, J.; CAICEDO, M.; ROUBIN, E.; HERNADÉZ, J. A.; HUESPE, A. Multi-scale (FE²) analysis of materials failure in cement/aggregate-type composite structure. In: EURO-C, 2004, St. Anton am Alberg, XX. Anais. Computational Modelling of Concrete Structure, Londres: CRC PRESS, 2014, v.1, p. 39-49.
• ^[19]
  PITUBA, J. J. C.; FERNANDES, G. R.; SOUZA NETO, E. A. Modelling of cohesive fracture and plasticity processes in composite microstructures. Journal of Engineering Mechanics, v. 142, 2016, pp. 04016069-1 - 04016069-15.
• ^[20]
  SOUZA NETO, E. A.; PERIÉ, D.; OWEN, D. R. J. Computational Methods for Plasticity: Theory and Applications. Swensea: Wiley, 2008.
• ^[21]
  ASSAAD, J. J.; HARB, J.; MAALUF, Y. Measurement of yield stress of cement pastes using the direct shear test. Journal of Non-Newtonian Fluid Mechanics, v. 214, 2014; p. 18-27.
• ^[22]
  MEHTA, P. K.; MONTEIRO, P. J. M. Concreto: Estrutura, Propriedades e Materiais. 1ª. ed. São Paulo: PINI, 2008.
• ^[23]
  NGUYEN, V.P.; LLOBERAS VALLS, O.; STROEVEN, M.; SLUYS, L.J. On the existence of representative volumes for softening quasi-brittle materials - a failure zone averaging scheme”, Computer Methods in Applied Mechanics Engineering, v. 199, 2010; p. 3026-3036.
• ^[24]
  RAMESH, G.; SOTELINO, E.D.; CHEN, W.F. Effect of transition zone on elastic moduli of concrete materials. Cement and Concrete Research, v. 26, 1996; p. 611-622.
• ^[25]
  YANG, C.C. Effect of transition zone on elastic moduli of mortar. Cement and Concrete Research, v. 28, 1998; p. 727-736.
• ^[26]
  RODRIGUES, E. A.; MANZOLI, O. L.; BITENCOURT JR, L. A. G.; BITTENCOURT, T. N. 2D mesoscale model for concrete based on the use of interface element with a high aspect ratio, International Journal of Solids and Structures, v. 94-95, 2016; p. 112-124.
• ^[27]
  SÁNCHEZ, P. J., BLANCO, P. J., HUESPE, A. E. AND FEIJÓO, R. A. Failure-Oriented Multi-scale Variational Formulation: Micro-structures with nucleation and evolution of softening bands, Computer Methods in Applied Mechanics and Engineering, v. 257, 2013; p. 221-247.

Publication Dates

• Publication in this collection
  Sept 2017

History

• Received
  06 Mar 2016
• Accepted
  15 Jan 2017