Modeling Particles Elements in Damaged Reinforced Concrete Structures

Abstract

In this paper, we introduce a finite element mesoscale modeling of damaged concrete structures, based on nodal positions. The mesoscale modeling consists of particle and fiber finite elements embedded in matrix finite elements. While the matrix elements represent the cement matrix, particle elements are used to simulate the coarse aggregates and fiber elements are used for reinforcement rebars. The embedded theory is used to immerse the reinforcement (both particle and fiber elements) without increasing the total number of degrees of freedom. This strategy does not require nodal coincidence, allowing randomly distribute the coarse aggregates. The materials nonlinear behavior is considered by a scalar damage model for the cement matrix and coarse aggregates, and one-dimensional elastoplastic model is used for the steel rebars. Four examples are presented, with good correlation between numerical and experimental results. It is shown that structures simulated with particulate elements could endure higher loads for the same displacement, although the maximum force is obtained in models without inclusion of particle elements.

Keywords:
Embedded particle; Positional FEM; Damage model; Elastoplasticity

Graphical Abstract

1 INTRODUCTION

Concrete is the most extensively used material in construction and building materials due the possibility of producing elements with different shapes and dimensions and the ability to combine low cost of production and execution. Besides, high mechanical strength is obtained when associated with steel rebars. As such, many efforts are directed to develop computational tools capable to evaluate concrete structures and its design.

The number of scientific papers that explore numerical methods associated to concrete is also very extensive. Most of the efforts grew after the 50s with the advent of computational mechanics. It is usual to employ continuous damage and fracture mechanics, associated with concepts from the theory of elasticity and plasticity. It began with the studies conducted by Hillerborg et al. (1976), Bažant (1976) and Swartz et al. (1978), which involved fracture mechanics applied to the complex behavior of concrete. These studies advanced with Bažant and Oh (1983) who proposed preferential cracking regions called cracking band or cracking process zone considering a constitutive model with softening.

Some contributions related to distributed cracking are found in De Borst and Nauta (1985) and Rots et al. (1985). Besides, some constitutive laws for concrete were proposed to evaluate its mechanical response with the distributed cracking models (Carreira and Chu, 1985, 1986; Boone et al., 1986). Concomitantly, models based on the continuous damage mechanics were proposed as an alternative to fracture mechanics. In these models, infinitesimal discontinuities due mechanical degradation are considered homogeneous at material point, and material degradation is represented by quantifying the defect densities (Lemaitre and Dufailly, 1987; Lemaitre, 1992).

The development of the continuous damage theory made possible to consider the overall reduction of stiffness when subjected to mechanical and thermal loads. It began with Kachanov (1958), who employed continuous damage variables. The reduction of mechanical properties due cracks was incorporated by Rabotnov (1969). However, only in the work of Lemaitre and Chaboche (1985) the formulation of continuous damage mechanics has been consistently described, based on the thermodynamics of irreversible processes. Nevertheless, Mazars (1984) proposed a simple yet efficient model based on equivalent strains. The mechanical degradation is represented by a scalar damage variable defined by the elongation strains of the strain tensor, uniformly reducing the stiffness in all directions.

The macroscale mechanical degradation is influenced by the internal structure of the concrete, observed only in mesoscale and microscale. Therefore, multiscale analyses stood out (Unger and Eckardt, 2011; Nguyen et al., 2011; Lloberas‐Valls et al., 2012; Etse et al., 2012; Rodrigues et al., 2020). In macroscale, in which most of classical models are formulated, the properties are considered homogeneous for different portions of the material. In mesoscale, employed in this work, the particle structure of concrete matrix has great relevance, in which the material heterogeneity leads to stress concentrations. Finally, at microscale, the hardened cement paste, composed of hydration products or hydrates, non-hydrated residual clinker and micropores, are the dominant parameters (Van Mier, 1997).

A 2D mesoscale model based on Representative Volume Element (RVE) is proposed by Borges et al. (2017) and Borges and Pituba (2017). The three-phase material is considered using the FEM. Round linear elastic inclusions are distributed into an elastoplastic matrix. The Interstitial Transition Zone (ITZ) is modeled by cohesive elements to simulate contact. A three-phase material is also employed by Chen et al. (2018) to evaluated the influence of both aggregate distribution and geometrical shapes using two-dimensional approach. The results indicated that the geometric parameters of coarse and fine aggregates do not significantly affect macroscopic tensile strength, but affect the macro compressive stress-strain curve. Zhou et al. (2019) investigated the influence of the aggregate shape on the mechanical properties of concrete. Their results showed that the aggregates size strongly influence cracking and stiffness. Forti et al. (2020) investigate failure behavior using a model that combines damage and plasticity on mortars for compressive softening, and cohesive fracture on interface elements for tensile softening and cracking. For an extensive review of geometry generation, particle distribution algorithm, constitutive relations and applications, the reader is referred to Thilakarathna et al. (2020). For mesoscale models and material behavior, the reader is referred to Caballero et al. (2006), Häfner et al. (2006), Wriggers and Moftah (2006), Kim and Al-Rub (2011), Du et al. (2014), Yılmaz and Molinari (2017), Zhou et al. (2017), Zhang et al. (2018), Li et al. (2019), Tang et al. (2019) and Meng et al. (2020).

The most recent discussion about multiscale modeling of concrete is seen in Rodrigues et al. (2020), who developed a 3D multiscale model, considering macro and mesoscale modeled independently yet coupled, based on Rodrigues et al. (2016) and Rodrigues et al. (2018). The mesoscopic scale is based on mesh fragmentation technique (MFT) (Sánchez et al., 2014). Interface elements (IEs) with high aspect ratio models the failure of the ITZ (Manzoli et al., 2012). The kinematic coupling is performed by a technique developed by Bitencourt Jr et al. (2015).

This notable advance in the study of concrete behavior is associated to the development of numerical methods. In mesoscale modeling of concrete structures different analysis methods stand out, such as Finite Element Modelling (FEM), Lattice Element Method (LEM), Rigid Body Spring Method (RBSM), Discrete Element Method (DEM) and Lattice Discrete Particle Method (LDPM). As an alternative to the FEM based on displacements, a FEM based on positions, called Positional FEM, has been proposed at the beginning of the 21st century, with the advantage of naturally consider the geometrical nonlinear effects (Bonet et al., 2000; Coda and Greco, 2004). Reinforced composite materials, such as concrete, are easily modelled by the Positional FEM employing embedded techniques. With this technique, the structure stiffness is composed by the contribution of the domain (concrete matrix), fiber (rebars) and particle elements (coarse aggregates), and does not require nodal coincidence between discretizations, maintaining the overall number of degrees of freedom (Vanalli et al., 2008; Radtke et al., 2010; Radtke et al., 2011; Sampaio et al., 2013; Paccola et al., 2015; Paccola and Coda, 2016; Ramos and Carrazedo, 2020).

In this work, we model the mechanical behavior of reinforced concrete structures considering the kinematic coupling of the cement matrix, the rebar and coarse aggregates using an embedded technique. The ITZ and voids distributed in the mortar are disregarded. Concrete physical nonlinearity is addressed by damage mechanics and the steel rebar nonlinearity is described by an elastoplastic model. The greatest advantage in our alternative mesoscale model is that the particles do not increase the number of degrees of freedom. Besides, large displacements can be naturally considered by the geometrical nonlinear formulation of the Positional FEM.

2 GEOMETRICALLY NONLINEAR FORMULATION OF COMPOSITE SOLIDS

In this section, the numerical model developed for the simulation of the mechanical behavior of reinforced concrete structures considering particle elements is exposed. The formulation is based on a positional version of the Finite Element Method, and the coupling is made by an embedded technique.

2.1 Geometric nonlinear equilibrium

The mechanical energy is composed by the strain energy ($U$), the potential energy of external forces ($P$) and the kinetic energy ($K$). The kinetic energy is disregarded for quasi-static applications, which is the case in this article. Therefore, considering the absence of dissipative potentials, the mechanical energy and its variation are written as:

$\Pi =U+Pand\delta \Pi =\delta U+\delta P$ (1)

Equilibrium is given by the Principle of Conservation of Mechanical Energy, or Principle of Stationary, defined as (Lanczos, 1986):

$\delta \Pi =\left(\frac{\partial U}{\partial Y_i}\right)\delta Y_i+\left(\frac{\partial P}{\partial Y_i}\right)\delta Y_i=0$ (2)

in which $Y_i$ is the position vector. Given the arbitrariness of $\delta Y_i$, the Eq. (2) can be rewritten as:

$g_i=\frac{\partial U}{\partial Y_i}+\frac{\partial P}{\partial Y_i}=F_i^{\mathrm{int}}-F_i^{ext}\cong 0_i$ (3)

Equation (3) is a system of nonlinear equations and represents the unbalanced forces between the internal and external forces. Since the solution depends on the current unknown position, Eq. (3) generates an unbalanced, residue vector $g_i$. Equilibrium or current configuration occurs when $g_i\left(Y_k\right)=0$. The current position is determined by resorting to Taylor series expansion for a trial position $Y_k^{trial}$, considering only the first order terms, resulting in:

$g_i\left(Y_k\right)=g_i\left(Y_k^{trial}\right)+{\frac{\partial g_i}{\partial Y_j}|}_{Y_k^{trial}}\Delta Y_j\cong 0_i$ (4)

Manipulation of Eq. (4) leads to:

$g_i\left(Y_k^{trial}\right)=-{\frac{\partial g_i}{\partial Y_j}|}_{Y_k^{trial}}\Delta Y_j=-{\frac{{\partial }^{2}U}{\partial Y_j\partial Y_i}|}_{Y_k^{trial}}\Delta Y_j=-H_{ji}\Delta Y_j=-H_{ij}\Delta Y_j$ (5)

where $\Delta Y_j$ is the correction of the trial positions and $H_{ij}$ refers to the Hessian or tangent stiffness matrix, which it is given by the second derivative of the strain energy in relation to the position vector, for conservative external forces.

Thus, the equilibrium configuration is determined through the Newton-Raphson iterative process, with $Y_k^{trial+1}=Y_k^{trial}+\Delta Y_k$, where the initial trial position corresponds to the undeformed configuration, that is, $Y_k^{trial}=X_k$. The convergence criterion is based on the norm of positions, in which the equilibrium condition is fully satisfied when $\sqrt{\Delta Y_k\Delta Y_k/X_l{X}_l}$ result less or equal to a previously stablished tolerance.

2.2 2D solid finite element kinematics

In this subsection, the continuum is discretized in elements with the introduction of the FEM approximations, seeking the internal forces and the Hessian matrix, required by the Newton-Raphson method. Two-dimensional solid finite elements are used for both the cement matrix and coarse aggregate. Considering a finite element, the mapping of its initial and current configuration is given by means of the shape functions, obtained from the product between the Lagrange polynomials, as follows:

$f_i^0=x_i={\varphi }_l\left({\xi }_1,{\xi }_2\right)X_{li}and{f}_i^1=y_i={\varphi }_l\left({\xi }_1,{\xi }_2\right)Y_{li}$ (6)

in which $f_i^0$ and $f_i^1$ maps of the initial and current configurations, respectively, to the dimensionless coordinates ${\xi }_1$ and ${\xi }_2$, ${\varphi }_l$ are the shape functions, and $X_{li}$ and $Y_{li}$ represent, respectively, the initial and current coordinates of node $l$ in direction $i$. In order to determine the quantities of interest, the Green-Lagrange nonlinear strain measure is defined, written as follows (Holzapfel, 2004):

$E_{ij}=\frac{1}{2}\left(A_{ki}{A}_{kj}-{\delta }_{ij}\right)$ (7)

in which $E$ is the Green-Lagrange strain tensor, $A$is the gradient of the change of configuration, also called deformation gradient, and $\delta$ is the Kronecker delta. The deformation gradient is composed by the mapping of the initial and current configuration as follows (Bonet et al., 2000; Coda and Greco, 2004):

$\nabla \overrightarrow{f}=\nabla \left({\overrightarrow{f}}^1\circ {\left({\overrightarrow{f}}^0\right)}^{-1}\right)=\nabla {\overrightarrow{f}}^1\cdot \nabla {\left({\overrightarrow{f}}^0\right)}^{-1}=A^1\cdot {\left(A^0\right)}^{-1}$ (8)

in which $A^0$ and $A^1$ are the gradients of the mappings from the dimensionless coordinates to the initial and current configurations, respectively, given by:

$A_{ij}^0=\frac{\partial f_i^0}{\partial {\xi }_j}={\varphi }_{l,j}{X}_{li}and{A}_{ij}^1=\frac{\partial f_i^1}{\partial {\xi }_j}={\varphi }_{l,j}{Y}_{li}$ (9)

The specific strain energy ($u_e$) is integrated into the initial volume of the solid ($V_0$) to find the strain energy ($U$), both written according to Eq. (10).

$u_e=E_{ij}{ℂ}_{ijkl}{E}_{kl}andU={\displaystyle {\int }_{V_0}{u}_{e}d{V}_0}$ (10)

in which $ℂ$ is the elastic constitutive tensor of the material, given by:

${ℂ}_{ijkl}=\mu \left({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk}\right)+\lambda {\delta }_{ij}{\delta }_{kl}$ (11)

where $\mu$ and $\lambda$ are Lamé constants which can be determined as a function of the longitudinal elastic modulus $\mathbb{E}$ and the Poisson's ratio $\nu$, according to the following expressions:

$\mu =\frac{\mathbb{E}}{2\left(1+\nu \right)}and\lambda =\frac{\nu \mathbb{E}}{\left(1+\nu \right)\left(1-2\nu \right)}$ (12)

Our formulation uses the second Piola-Kirchhoff stress tensor, energy conjugate of the Green-Lagrange strain, resulting in the constitutive law of Saint-Venant-Kirchhoff. Such constitutive relationship is similar to the generalized Hooke's law, linearly associating the respective stress and strain measurements, as expressed in Eq. (13):

$S_{ij}=\frac{\partial u_e}{\partial E_{ij}}={ℂ}_{ijkl}{E}_{kl}$ (13)

Particularly, the second Piola-Kirchhoff stress tensor is written for state of plane stress and plane strain by:

$S_{ij}=2\mu E_{ij}+\overline{\lambda }{E}_{kk}{\delta }_{ij}and{S}_{ij}=2\mu E_{ij}+\lambda E_{kk}{\delta }_{ij}$ (14)

in which $\overline{\lambda }=\nu \mathbb{E}/\left(1-{\nu }^2\right)$.

The internal forces developed in a finite element are defined as the strain energy derivative. Due to the complexity of the integral core, numerical integration is used. Integration is performed using Hammer quadrature. Thus, the internal forces of node $i$ for direction $j$ are written as:

$F_{ij}^{el}=\frac{\partial U^{el}}{\partial Y_{ij}}={\displaystyle {\int }_{V_0^{el}}\frac{\partial u_e\left(Y_{kl}\right)}{\partial Y_{ij}}d{V}_0^{el}}={\displaystyle \sum _{h=1}^n\frac{\partial u_e\left(Y_{kl}\left({\xi }_1^h,{\xi }_2^h\right)\right)}{\partial Y_{ij}}{J}_0\left({\xi }_1^h,{\xi }_2^h\right)w_h}$ (15)

in which $n$ denotes the number of integration points, with weights given by $w_h$, and $J_0$ is Jacobian of the initial configuration, calculated as $J_0=\det \left(A^0\right)$. The specific strain energy derivative can be further developed using Eq. (13), leading to:

$\frac{\partial u_e}{\partial Y_{ij}}=\frac{\partial u_e}{\partial E_{kl}}\frac{\partial E_{kl}}{\partial Y_{ij}}=S_{kl}\frac{\partial E_{kl}}{\partial Y_{ij}}$ (16)

Derivative of Green-Lagrange strain with respect to the positions is written according to Eq. (17). Notice that mixed notation is used here for readability.

$\frac{\partial E}{\partial Y_{ij}}=\frac{1}{2}\left({\left(A^0\right)}^{-t}\cdot \frac{\partial {\left(A^1\right)}^t}{\partial Y_{ij}}\cdot A^1\cdot {\left(A^0\right)}^{-1}+{\left(A^0\right)}^{-t}\cdot {\left(A^1\right)}^t\cdot \frac{\partial A^1}{\partial Y_{ij}}\cdot {\left(A^0\right)}^{-1}\right)$ (17)

in which the derivative of the gradient of the mapping of the current configuration with respect to the positions is given by:

$\frac{\partial A_{kl}^1}{\partial Y_{ij}}=\frac{\partial }{\partial Y_{ij}}\left({\varphi }_{m,l}{Y}_{mk}\right)={\varphi }_{m,l}\frac{\partial Y_{mk}}{\partial Y_{ij}}={\varphi }_{m,l}{\delta }_{im}{\delta }_{jk}={\varphi }_{i,l}{\delta }_{jk}$ (18)

The Hessian matrix is given by the second derivative of the strain energy, since external forces are considered conservative. Thus, performing the numerical integration in the initial volume of the solid, it is written as:

$H_{ijkl}^{el}=\frac{{\partial }^2{U}^{el}}{\partial Y_{ij}\partial Y_{kl}}={\displaystyle {\int }_{V_0^{el}}\frac{{\partial }^2{u}_e\left(Y_{op}\right)}{\partial Y_{ij}\partial Y_{kl}}d{V}_0^{el}}={\displaystyle \sum _{h=1}^n\frac{{\partial }^2{u}_e\left(Y_{op}\left({\xi }_1^h,{\xi }_2^h\right)\right)}{\partial Y_{ij}\partial Y_{kl}}{J}_0\left({\xi }_1^h,{\xi }_2^h\right)w_h}$ (19)

in which the second derivative of specific strain energy is calculated by deriving Eq. (16) with respect to positions, leading to:

$\frac{{\partial }^2{u}_e}{\partial Y_{ij}\partial Y_{kl}}=\frac{\partial E}{\partial Y_{ij}}:ℂ:\frac{\partial E}{\partial Y_{kl}}+S:\frac{{\partial }^{2}E}{\partial Y_{ij}\partial Y_{kl}}$ (20)

The terms $\partial E/\partial Y_{ij}$, $ℂ$ and $S$ are defined respectively by Eqs. (17), (11) and (13). Only the second derivative of the Green-Lagrange strain related to positions remains, which is obtained by deriving Eq. (17). Thus, ${\partial }^{2}E/\partial Y_{ij}\partial Y_{kl}$ is written as follows:

$\frac{{\partial }^{2}E}{\partial Y_{ij}\partial Y_{kl}}=\frac{1}{2}\left({\left(A^0\right)}^{-t}\cdot \frac{\partial {\left(A^1\right)}^t}{\partial Y_{ij}}\cdot \frac{\partial A^1}{\partial Y_{kj}}\cdot {\left(A^0\right)}^{-1}+{\left(A^0\right)}^{-t}\cdot \frac{\partial {\left(A^1\right)}^t}{\partial Y_{kl}}\cdot \frac{\partial A^1}{\partial Y_{ij}}\cdot {\left(A^0\right)}^{-1}\right)$ (21)

2.3 Finite element kinematics of rebars

Steel rebars are modeled by fiber elements, represented by truss finite element. Therefore, mathematical developments regarding its positional formulation are briefly presented. The linear truss element is a particular case of an one-dimensional bar element with only axial stiffness, and two nodes with two degrees of freedom each. Thus, initial and current lengths can be determined from the nodal positions, according to the following expressions:

${\left({\overline{L}}_0\right)}^2={\left({\overline{X}}_{\alpha 1}-{\overline{X}}_{\beta 1}\right)}^2+{\left({\overline{X}}_{\alpha 2}-{\overline{X}}_{\beta 2}\right)}^{2}and{\left(\overline{L}\right)}^2={\left({\overline{Y}}_{\alpha 1}-{\overline{Y}}_{\beta 1}\right)}^2+{\left({\overline{Y}}_{\alpha 2}-{\overline{Y}}_{\beta 2}\right)}^2$ (22)

where ${\overline{L}}_0$ and $\overline{L}$ correspond to the initial and current length, in which indices 1 and 2 refer to directions, and $\alpha$ and $\beta$ denote the node of the element. A bar written over quantities is used to distinguish quantities related to truss elements. Therefore, the Green-Lagrange strain is calculated as the engineering nonlinear strain of a fiber immersed in the continuum, expressed as:

$\overline{E}=\frac{1}{2}\left(\frac{{\overline{L}}^2-{\overline{L}}_0^2}{{\overline{L}}_0^2}\right)$ (23)

Since only linear truss element is used, solution is obtained analytically. Thus, the strain energy is written according to Eq. (24), with ${\overline{V}}_0={\overline{A}}_0{\overline{L}}_0$.

${\overline{U}}^{el}={\displaystyle {\int }_{{\overline{V}}_0^{el}}{\overline{u}}_{e}d{\overline{V}}_0^{el}}={\overline{u}}_e{\overline{V}}_0^{el}={\overline{u}}_e{\overline{A}}_0{\overline{L}}_0$ (24)

As previously described, internal forces are defined as the strain energy derivative with respect to the positions, which for a truss element results in:

${\overline{F}}_{ij}^{el}=\frac{\partial {\overline{U}}^{el}}{\partial {\overline{Y}}_{ij}}={\overline{A}}_0{\overline{L}}_0\frac{\partial {\overline{u}}_e}{\partial {\overline{Y}}_{ij}}={\overline{A}}_0{\overline{L}}_0\frac{\partial {\overline{u}}_e}{\partial \overline{E}}\frac{\partial \overline{E}}{\partial {\overline{Y}}_{ij}}={\overline{A}}_0{\overline{L}}_0\overline{S}\frac{\partial \overline{E}}{\partial {\overline{Y}}_{ij}}$ (25)

where $\partial {\overline{u}}_e/\partial {\overline{Y}}_{ij}=\overline{S}$. Notice that $\overline{E}$ and $\overline{S}$ reduce to scalars in one-dimensional case. The Green-Lagrange strain derivative with respect to the nodal parameters is given by:

$\frac{\partial \overline{E}}{\partial {\overline{Y}}_{ij}}=\frac{{\left(-1\right)}^i}{{\overline{L}}_0^2}\left({\overline{Y}}_{2j}-{\overline{Y}}_{1j}\right)$ (26)

The Hessian tensor is determined from strain energy second derivative with respect to the positions, given by:

${\overline{H}}_{ijkl}^{el}=\frac{{\partial }^2{\overline{U}}^{el}}{\partial {\overline{Y}}_{ij}\partial {\overline{Y}}_{kl}}=\frac{\partial }{\partial {\overline{Y}}_{kl}}\left(\frac{\partial {\overline{U}}^{el}}{\partial {\overline{Y}}_{ij}}\right)$ (27)

Using Eq. (25) in Eq. (27) and applying the product rule and chain rule, we have:

${\overline{H}}_{ijkl}^{el}={\overline{A}}_0{\overline{L}}_0\left(\frac{{\partial }^2{\overline{u}}_e}{\partial {\overline{E}}^2}\frac{\partial \overline{E}}{\partial {\overline{Y}}_{kl}}\frac{\partial \overline{E}}{\partial {\overline{Y}}_{ij}}+\frac{\partial {\overline{u}}_e}{\partial \overline{E}}\frac{{\partial }^2\overline{E}}{\partial {\overline{Y}}_{ij}\partial {\overline{Y}}_{kl}}\right)$ (28)

For the Saint-Venant-Kirchhoff constitutive model, the first derivative of specific strain energy with respect to the Green-Lagrange strain corresponds to the Piola-Kirchhoff stress $\overline{S}$, while the second derivative results in the material longitudinal elastic modulus $\overline{\mathbb{E}}$. In this sense, the Hessian is written as:

${\overline{H}}_{ijkl}^{el}={\overline{A}}_0{\overline{L}}_0\left(\overline{\mathbb{E}}\frac{\partial \overline{E}}{\partial {\overline{Y}}_{kl}}\frac{\partial \overline{E}}{\partial {\overline{Y}}_{ij}}+\overline{S}\frac{{\partial }^2\overline{E}}{\partial {\overline{Y}}_{ij}\partial {\overline{Y}}_{kl}}\right)$ (29)

where $\partial \overline{E}/\partial {\overline{Y}}_{ij}$ is calculated according to Eq. (26) and ${\partial }^2\overline{E}/\partial {\overline{Y}}_{ij}\partial {\overline{Y}}_{kl}$ is given by:

$\frac{{\partial }^2\overline{E}}{\partial {\overline{Y}}_{ij}\partial {\overline{Y}}_{kl}}=\frac{\partial }{\partial {\overline{Y}}_{kl}}\left(\frac{\partial \overline{E}}{\partial {\overline{Y}}_{ij}}\right)=\frac{{\left(-1\right)}^i}{{\overline{L}}_0^2}\frac{\partial \left({\overline{Y}}_{2j}-{\overline{Y}}_{1j}\right)}{\partial {\overline{Y}}_{kl}}=\frac{{\left(-1\right)}^i{\left(-1\right)}^k}{{\overline{L}}_0^2}{\delta }_{jl}$ (30)

2.4 Kinematics coupling strategy

Fiber and particle elements are embedded within the matrix based on a kinematic strategy presented in Vanalli et al. (2008), Radtke et al. (2010), Radtke et al. (2011), Sampaio et al. (2013), Paccola et al. (2015), Paccola and Coda (2016) and Ramos and Carrazedo (2020). It consists of writing the nodal positions of the reinforcement elements according to the nodal positions of the matrix elements. Thus, the strain energy of the reinforcement contributes in the strain energy of the matrix elements based on its shape functions. It requires that the dimensionless coordinates of the reinforcement nodes in which matrix element is contained are known. But first, notice that the initial and current positions are determined as follows:

${\overline{X}}_j^p={\varphi }_i\left({\xi }_1^p,{\xi }_2^p\right){\widehat{X}}_{ij}and{\overline{Y}}_j^p={\varphi }_i\left({\xi }_1^p,{\xi }_2^p\right){\widehat{Y}}_{ij}$ (31)

where ${\overline{X}}_j^p$ and ${\overline{Y}}_j^p$ are the initial and current positions of node $p$ of the reinforcement in $j$ direction, ${\varphi }_i$ are the shape functions applied to the dimensionless coordinates of node $p$, and ${\widehat{X}}_{ij}$ and ${\widehat{Y}}_{ij}$ are the initial and current positions of the matrix element nodes $i$ where the reinforcement node is contained, in $j$ direction. Notice that, in this section, the matrix quantities are written using circumflex ($\widehat{\odot }$).

The first expression of Eq. (31) leads to a nonlinear system of equations, required to find the dimensionless coordinates. The solution is given by the Newton-Raphson method. Thus, expanding in a Taylor series over a tentative dimensionless coordinate (${\xi }_j^{trial}$), and disregarding higher order terms, it is written:

${\overline{X}}_i={\varphi }_l\left({\xi }_1^{trial},{\xi }_2^{trial}\right){\widehat{X}}_{li}+\left({\frac{\partial {\varphi }_l\left({\xi }_1,{\xi }_2\right){\widehat{X}}_{li}}{\partial {\xi }_j}|}_{\left({\xi }_1^{trial},{\xi }_2^{trial}\right)}\right)\Delta {\xi }_{j}or{\overline{X}}_i={\overline{X}}_i^{trial}+L_{ij}\Delta {\xi }_j$ (32)

where $\Delta {\xi }_j$ is the correction of dimensionless coordinates and $L_{ij}$ corresponds to a 2x2 matrix. The vector of dimensionless coordinates is then updated (${\xi }_k^{trial+1}={\xi }_k^{trial}+\Delta {\xi }_k$) and the process is repeated until a tolerance is achieved ($\sqrt{\Delta {\xi }_k\Delta {\xi }_k}\le tol$). The resulting dimensionless coordinates must belong to the matrix dimensionless domain, that is, $0\le {\xi }_1\le 1$, $0\le {\xi }_2\le 1$ and $0\le 1-{\xi }_1-{\xi }_2\le 1$. This will only happen if the reinforcement node is contained in the current matrix element.

Equation (33) establishes the strain energy of the composite, where $\widehat{U}$ and $\overline{U}$are the contribution of the matrix and reinforcement strain energies, respectively.

$U=\widehat{U}\left({\widehat{Y}}_{ij}\right)+\overline{U}\left({\overline{Y}}_{mn}\left({\widehat{Y}}_{ij}\right)\right)$ (33)

Thus, the internal forces of the composite are written as:

$F_{ij}^{\mathrm{int}}=\frac{\partial \widehat{U}}{\partial {\widehat{Y}}_{ij}}+\frac{\partial \overline{U}}{\partial {\widehat{Y}}_{ij}}=\frac{\partial \widehat{U}}{\partial {\widehat{Y}}_{ij}}+\frac{\partial \overline{U}}{\partial {\overline{Y}}_{mn}}\frac{\partial {\overline{Y}}_{mn}}{\partial {\widehat{Y}}_{ij}}={\widehat{F}}_{ij}+{\overline{F}}_{mn}\frac{\partial {\overline{Y}}_{mn}}{\partial {\widehat{Y}}_{ij}}$ (34)

Since the internal forces of both matrix and reinforcement elements are already known, calculated by Eqs. (15) and (25), only the coupling term must be defined, given by:

$\frac{\partial {\overline{Y}}_{mn}}{\partial {\widehat{Y}}_{ij}}=\frac{\partial }{\partial {\widehat{Y}}_{ij}}\left({\varphi }_l\left({\xi }_1^m,{\xi }_2^m\right){\widehat{Y}}_{\ln }\right)={\varphi }_l\left({\xi }_1^m,{\xi }_2^m\right)\frac{\partial {\widehat{Y}}_{\ln }}{\partial {\widehat{Y}}_{ij}}={\varphi }_l\left({\xi }_1^m,{\xi }_2^m\right){\delta }_{il}{\delta }_{jn}$ (35)

The Hessian matrix is also decomposed additively, including both reinforcement and matrix stiffness, as shown in the following expression:

$H_{ijkl}=\frac{{\partial }^2\widehat{U}}{\partial {\widehat{Y}}_{ij}\partial {\widehat{Y}}_{kl}}+\frac{{\partial }^2\overline{U}}{\partial {\widehat{Y}}_{ij}\partial {\widehat{Y}}_{kl}}$ (36)

The first part of Eq. (36) is calculated according to Eq. (19). The second part is expanded in terms of the specific strain energy, as follows:

$\frac{{\partial }^2{\overline{U}}^{el}}{\partial {\widehat{Y}}_{ij}\partial {\widehat{Y}}_{kl}}={\displaystyle {\int }_{{\overline{V}}_0^{el}}\frac{{\partial }^2{\overline{u}}_e}{\partial {\widehat{Y}}_{ij}\partial {\widehat{Y}}_{kl}}d{\overline{V}}_0^{el}}={\displaystyle {\int }_{V_0^{el}}\left(\begin{array}{l}\frac{{\partial }^2{\overline{u}}_e}{\partial {\overline{Y}}_{mn}\partial {\overline{Y}}_{mn}}\frac{\partial {\overline{Y}}_{mn}}{\partial {\widehat{Y}}_{ij}}\frac{\partial {\overline{Y}}_{mn}}{\partial {\widehat{Y}}_{kl}}+\frac{{\partial }^2{\overline{u}}_e}{\partial {\overline{Y}}_{mn}\partial {\overline{Y}}_{op}}\frac{\partial {\overline{Y}}_{mn}}{\partial {\widehat{Y}}_{ij}}\frac{\partial {\overline{Y}}_{op}}{\partial {\widehat{Y}}_{kl}}+\\ \frac{{\partial }^2{\overline{u}}_e}{\partial {\overline{Y}}_{op}\partial {\overline{Y}}_{mn}}\frac{\partial {\overline{Y}}_{op}}{\partial {\widehat{Y}}_{ij}}\frac{\partial {\overline{Y}}_{mn}}{\partial {\widehat{Y}}_{kl}}+\frac{{\partial }^2{\overline{u}}_e}{\partial {\overline{Y}}_{op}\partial {\overline{Y}}_{op}}\frac{\partial {\overline{Y}}_{op}}{\partial {\widehat{Y}}_{ij}}\frac{\partial {\overline{Y}}_{op}}{\partial {\widehat{Y}}_{kl}}\end{array}\right)d{V}_0^{el}}$ (37)

With the proper contributions of reinforcement and matrix elements to internal forces and Hessian matrix of the two-dimensional solid, the geometric nonlinear solution strategy based on the Newton-Raphson iterative incremental method can be applied to seek the equilibrium configuration of the reinforced composite.

3 MECHANICAL DEGRADATION

3.1 Continuous damage mechanics

The nonlinear mechanical behavior of concrete is addressed with theory of continuous damage mechanics, based on the thermodynamics of irreversible processes. In this study, the damage model proposed by Mazars (1984) is adopted, since it requires a low number of parameters, which can be easily obtained from experimental or theoretical uniaxial stress-strain diagrams. The elastic properties are penalized to represent the mechanical degradation due crack growth. Thus, this formulation is based on the following hypotheses: (i) residual strains are not considered; (ii) isotropic damage, represented by a scalar variable; and (iii) damage only occurs due elongation strains.

Mazar’s damage model proposes to measure the elongation state at a given point in the continuum by a scalar variable, represented by the equivalent strain ($E_{eq}$), as follows:

$E_{eq}=\sqrt{{\left(E_1\right)}_{+}^2+{\left(E_2\right)}_{+}^2+{\left(E_3\right)}_{+}^2}$ (38)

where ${\left(E_i\right)}_{+}$ corresponds to the positive components of the principal values of the strain tensor.

The damage increases when a reference strain is exceeded. Thus, the damage criterion can be defined as:

$f\left(E_{eq},D\right)=E_{eq}-E_{\lim }\left(D\right)\le 0$ (39)

where variable $E_{\lim }$ is the maximum elongation strain in the strain history. At the beginning of the process $E_{\lim }$ is equivalent to the strain corresponding to the concrete tensile strength$E_{d0}$. Since the formulation proposed here employs Green-Lagrange strain, it is converted as follows:

$E_{d0}\left({\epsilon }_{d0}\right)=0.5{\epsilon }_{d0}{}^2+{\epsilon }_{d0}$ (40)

in which ${\epsilon }_{d0}$ represents the specific strain from a uniaxial tensile test.

The model also establishes two damage variables $D_T$ and $D_C$ due the concrete behavior related to tensile and compressive stresses. These variables describe the nonlinear stress-strain diagram for concrete, and are given by:

$D_T=1-\frac{E_{d0}\left(1-A_T\right)}{E_{eq}}-\frac{A_T}{e^{B_T\left(E_{eq}-E_{d0}\right)}}and{D}_C=1-\frac{E_{d0}\left(1-A_C\right)}{E_{eq}}-\frac{A_C}{e^{B_C\left(E_{eq}-E_{d0}\right)}}$ (41)

where $A_T$, $A_C$, $B_T$ and $B_C$ are parameters of the Mazar’s damage model.

For multiaxial stresses, damage is evaluated by linear combination of the variables, as follows:

$D={\alpha }_T{D}_T+{\alpha }_C{D}_C$ (42)

Coefficient ${\alpha }_T$ and ${\alpha }_C$are calculated by decomposing the principal strains into two parts, according to Eqs. (43) and (44):

${\left(E_i^T\right)}_{+}=E_{i}if{E}_i>0and{\left(E_i^T\right)}_{+}=0if{E}_i\le 0$ (43)

${\left(E_i^C\right)}_{-}=-E_{i}if{E}_i<0and{\left(E_i^C\right)}_{-}=0if{E}_i\ge 0$ (44)

From that, ${\alpha }_T$ and ${\alpha }_C$ are calculated as follows:

${\alpha }_T={\displaystyle \sum _i{\left(E_i^T\right)}_{+}/{\left(E_V\right)}_{+}}and{\alpha }_C={\displaystyle \sum _i{\left(E_i^C\right)}_{-}/{\left(E_V\right)}_{+}}$ (45)

where $0\le {\alpha }_T\le 1$; $0\le {\alpha }_C\le 1$; and ${\alpha }_T+{\alpha }_C=1$. The variable ${\left(E_V\right)}_{+}$ represents the total state of elongation, given by the principal strains as follows:

${\left(E_V\right)}_{+}={\displaystyle \sum _i{\left(E_i^T\right)}_{+}}+{\displaystyle \sum _i{\left(E_i^C\right)}_{-}}$ (46)

Stresses are found by the damaged elastic properties of the material, according to Eq. (47):

${\tilde{S}}_{ij}=\left(1-D\right){ℂ}_{ijkl}{E}_{kl}$ (47)

where $\tilde{S}$ is the second Piola-Kirchhoff stress tensor, $ℂ$ is the fourth-order elastic constitutive tensor and $E$ is the Green-Lagrange strain tensor.

3.2 Rebar elastoplasticity

Ductile materials undergo permanent strain in response to applied forces, non-reversible when unloaded. Plasticity is mathematically defined by three basic hypotheses: (i) definition of the yield strength, the transition from elastic to plastic behavior; (ii) definition of a flow rule, that describes the evolution of plastic strain; and (iii) definition of a hardening law, which establishes the stress-strain relationship in the plastic stage.

Additive decomposition of the strains can be considered if only moderate strains are allowed. In this paper, the nonlinear behavior of steel is given by an elastoplastic model with positive linear isotropic hardening, in which yielding function is given by:

$f\left(S_s,\chi \right)=\left|S_s\right|-\left(S_{sy}+K\chi \right)\le 0$ (48)

where $S_s$ is stress acting on steel, $S_{sy}$ is the steel yield stress, $K$ is the plastic hardening module and $\chi$ is the equivalent plastic strain measure.

In the one-dimensional case, $S_{sy}$ is determined from the Cauchy stress and the Piola-Kirchhoff stress, as follows:

$S_{sy}=\frac{{\sigma }_{sy}}{\left({\epsilon }_{sy}+1\right)}or{S}_{sy}=\frac{{\sigma }_{sy}{\mathbb{E}}_s}{\left({\sigma }_{sy}+{\mathbb{E}}_s\right)}$ (49)

where ${\mathbb{E}}_s$ is Young's modulus of steel, ${\sigma }_{sy}$ is the steel yield Cauchy stress and ${\epsilon }_{sy}$ is the associated strain.

The increase in stress in the inelastic phase is given as follow:

$dS=\mathbb{E}{s}_{}\left[dE-{\mathbb{E}}_{s}dE/\left({\mathbb{E}}_s+K\right)\right]=\left[{\mathbb{E}}_{s}K/\left({\mathbb{E}}_s+K\right)\right]dE$ (50)

where ${\mathbb{E}}_{s}K/\left({\mathbb{E}}_s+K\right)$ is the tangent elastoplastic module corresponding to hardening phase.

More details related to the theory of plasticity in the scope of computational methods can be consulted in the works of Simo and Hughes (1998), Chen and Han (2007) and Souza Neto et al. (2011).

4 EXAMPLES

Four examples are provided. The first shows the adopted kinematic coupling adopted for matrix and particle elements, applied to evaluation of the Young’s modulus of concrete. The next three examples were select to validate the proposed formulation. In the second example, a beam with and without inclusion of particle elements to simulate coarse aggregate is performed. The third example is a column with geometric imperfection and eccentric load, outlining physical and geometrical nonlinearities. In the fourth example, a reinforced concrete frame built by a monolithic set of two beams and two columns is simulated with coarse aggregates. The examples aim to explore the applicability of the proposed formulation in the mechanical modeling of usual reinforced concrete structures.

4.1 Evaluation of the Young’s modulus of concrete

This example shows a simple methodology to numerically estimate the Young’s modulus of concrete, based on only two phases: cement paste and coarse aggregates. The physical parameters for cement and aggregates are: ${\mathbb{E}}_m$ = 20 GPa, ${\nu }_m$ = 0.2, ${\mathbb{E}}_a$ = 100 GPa, ${\nu }_a$ = 0.0. Square specimens of sides 100 cm are simulated by a finite element mesh composed of 200 triangular elements with cubic approximation and 961 nodes, discretizing the matrix elements. Figure 1 shows where the specimens are restrained, as well its geometry and discretization. A displacement of ${\delta }_c$ = 0.1 cm was imposed in the upper side. Since only small displacements are imposed, it is admitted physical linear behavior. Thus, no correction is required for the Young’s modulus of concrete due to the constitutive law. The standard test method for the determination of the compressive strength of cylindrical concrete specimens is presented by ASTM C39/C39M-20 (2020).

Figure 1
Geometry and finite element mesh adopted for the specimen.

Three volume fractions of aggregates were adopted (20%, 25% and 30%). Each coarse aggregate was modeled by 8 particle triangular elements with 9 nodes, randomly distributed in the cement matrix. This was made possible due embedded technique, which does not require nodal coincidence between matrix and particle meshes. Geometric details regarding coarse aggregates elements are shown in Figure 2.

Figure 2
Details of the aggregates included in the matrix.

We also adopted three discretization for each volumetric rate of aggregates to observe the influence of the aggregates dimensions. The adopted aggregate discretizations are listed in Table 1 and shown in Figures 3, 4 and 5.

Table 1
Particulate specimen characteristics.

Volumetric rate of Aggregate (%)	Number of aggregates	Aggregate volume (cm3)	Larger particle size (cm)	Number of particle elements	Number of nodes
20	20	100	14.142	160	180
20	200	10	4.472	1600	1800
20	2000	1	1.414	16000	18000
25	25	100	14.142	200	225
25	250	10	4.472	2000	2250
25	2500	1	1.414	20000	22500
30	30	100	14.142	240	270
30	300	10	4.472	2400	2700
30	3000	1	1.414	24000	27000

Figure 3
Aggregate distribution for volumetric rate of aggregate of 20%: (a) 160 particles, (b) 1600 particles and (c) 16000 particles.

Figure 4
Aggregate distribution for volumetric rate of aggregate of 25%: (a) 160 particles, (b) 1600 particles and (c) 16000 particles.

Figure 5
Aggregate distribution for volumetric rate of aggregate of 30%: (a) 160 particles, (b) 1600 particles and (c) 16000 particles.

Results are compared with the homogenization models proposed by Voigt, Reuss, Hansen, Hirsch and Counto, whose formulation are transcribed respectively by Eqs. (51) to (55) (Mehta and Monteiro, 2014; Li et al. 2019).

${\mathbb{E}}_c={\mathbb{E}}_m{c}_m+{\mathbb{E}}_a{c}_a$ (51)

$\frac{1}{{\mathbb{E}}_c}=\frac{c_m}{{\mathbb{E}}_m}+\frac{c_a}{{\mathbb{E}}_a}$ (52)

${\mathbb{E}}_c=\left[\frac{{\mathbb{E}}_m{c}_m+\left(1+c_a\right){\mathbb{E}}_a}{\left(1+c_a\right){\mathbb{E}}_m+{\mathbb{E}}_a{c}_m}\right]{\mathbb{E}}_m$ (53)

$\frac{1}{{\mathbb{E}}_c}=\frac{1}{2}\left[\left(\frac{c_a}{{\mathbb{E}}_a}+\frac{c_m}{{\mathbb{E}}_m}\right)+\left(\frac{1}{c_a{\mathbb{E}}_a+c_m{\mathbb{E}}_m}\right)\right]$ (54)

$\frac{1}{{\mathbb{E}}_c}=\left[\frac{1-\sqrt{c_a}}{{\mathbb{E}}_m}\right]+\frac{1}{\left[\frac{1-\sqrt{c_a}}{\sqrt{c_a}}\right]{\mathbb{E}}_m+{\mathbb{E}}_a}$ (55)

where ${\mathbb{E}}_c$ is the equivalent Young’s modulus, ${\mathbb{E}}_m$ and ${\mathbb{E}}_a$ are the Young’s modulus of the cementitious matrix and the coarse aggregates, $c_m$ and $c_a$ are the volumetric rates of cement paste and aggregates, respectively. Numerical results of the obtained Young’s modulus are presented in Table 2.

Table 2
Comparison of the Young’s modules.

Volumetric rate of Aggregate (%)	Aggregate volume (cm3)	Equivalent modulus of elasticity (GPa)
		Numerical	Voigt	Reuss	Hansen	Hirsch	Counto
20	100	26.829	36.000	23.810	26.154	28.662	28.045
	10	31.134
	1	35.432
25	100	29.175	40.000	25.000	28.000	30.769	30.000
	10	35.172
	1	39.583
30	100	32.081	44.000	26.316	30.000	32.934	32.055
	10	39.613
	1	43.670

Results range from Reuss to Voight homogenization models, and, as the volume of the aggregate decreases, increasing the number of particle elements, the numerical results ranges from median to Voight. That result is expected, since the relationship between aggregate dimensions and dimensions of the structural element significantly influences the stiffness of the structure.

4.2 Reinforced concrete beam with concentrated load

In this example, we model a double support beam with concentrated load applied in the center of the upper surface, proposed and measured by Mazars (1984). We aim to validate the numerical implementation of Mazar’s damage model associated with the composite solid model based on FEM. To do so, we evaluate two cases: a beam composed by a concrete matrix with embedded steel rebar (matrix and fiber elements), and a beam composed by a cement matrix, embedded coarse aggregates and embedded steel rebar (matrix, particle and fiber elements). The geometric properties of the experimental test done by Mazars (1984) are shown in Figure 6.

Figure 6
Geometric properties of the beam.

Only half domain was modeled due symmetry, see Figure 7. The following elastic properties were adopted for the first simulation: ${\mathbb{E}}_c$ = 3000 kN/cm² and ${\nu }_c$ = 0.2 for the concrete, and ${\mathbb{E}}_s$ = 19600 kN/cm² and ${\nu }_s$ = 0.0 for the steel reinforcement. For Mazar’s damage model, the following parameters were used: ${\epsilon }_{d0}$ = 1.15x10^-5, $A_C$ = 1.40, $B_C$ = 1850, $A_T$ = 0.80 and $B_T$ = 20000.

Figure 7
Beam discretization.

3520 triangular elements with cubic approximation were used to model the matrix, and 240 truss linear elements were used to model the reinforcement. A vertical displacement of $\delta$ = -1.1 mm was applied, divided into 100 steps.

For the second simulation, a random distribution of particles elements to represent the coarse aggregate was introduced. The volumetric rate of aggregate is 30% of the cement matrix volume, leading to 528 coarse aggregates distributed randomly, each with 8 triangular elements with quadratic approximation. The aggregate maximum dimension is 1.414 cm, and no overlapping is allowed. Thus, the particulate phase was discretized by 4224 particle elements. The mesh used in this simulation is shown in Figure 8.

Figure 8
Distribution of particles in the domain.

Before evaluate the beam, we first performed numerical compression tests on 10x10 cm specimens to stablish equivalent Young’s modulus for the cement matrix and for the coarse aggregates, following the procedures described in the previous example. The cement matrix was discretized by 200 triangular elements with cubic approximation, and the coarse aggregate was discretized by 240 triangular elements. Figure 9 shows the adopted finite element mesh and the problem geometry.

Figure 9
Finite element mesh adopted for the specimen.

We simulated 10 specimens with random distribution of the aggregates, seeking an equivalent elastic modulus of ${\mathbb{E}}_c$ = 3000 kN/cm² and imposing a ratio of elastic modulus equal to ${\mathbb{E}}_m/{\mathbb{E}}_a$ = 0.20. For the volumetric rate of aggregates of 30%, we obtained ${\mathbb{E}}_m$ = 1990.21 kN/cm² and ${\mathbb{E}}_a$ = 9951.05 kN/cm². No changes to damage parameters $A_C$, $B_C$, $A_T$ and $B_T$ were made. Only ${\epsilon }_{d0}$ was adjusted, to impose that degradation develops at same stress level in both materials. Thus, ${\epsilon }_{d0}$= 1.733x10-4 was adopted for the cementitious matrix and ${\epsilon }_{d0}$ = 3.467x10^-5 for the coarse aggregate.

Numerical results are compared to the experimental measurements of Mazars (1984) and the numerical results of Santos (2015), shown in Figure 10. Two results are shown - with and without particle elements introduced in the matrix. There is a slight reduction on the beginning stage of the cracking process, and higher loads for the same displacement in the nonlinear stage in the simulation with particle elements. This is probably provided by the local stiffening around the particle elements, which spread damage to adjacent regions, see Figure 11.

Figure 10
Comparison of results: force vs displacements.

Figure 11
Distribution of the damage variable.

The damage variable was calculated at the points of integration and extrapolated to the element nodes using the least squares. In Figure 11, the evolution of damage that indicate the presence of cracks can be visualized accordingly. Notice that in the case of the homogenized beam, the propagation of damage follows main cracks. In the case of the particulate beam, damage deviate particles and create new crack surfaces. Figure 12 presents a close up in the lower half of central portion of the beam, and the evolution of the damage process. One may notice the steel rebar highlighted in red.

Figure 12
Damage evolution.

Thus, it is clear that the kinematic coupling of matrix and particle elements provided good results for nonlinear mechanical analysis of composite materials without homogenization.

4.3 Column with eccentric load and geometric imperfection

In this example, we numerically evaluate a column with eccentric load, considering an initial geometric imperfection. The experimental evaluation was made by Espion (1993), which is an important benchmark, often used for validation of nonlinear numerical models. This laboratory assessment was also modeled by Parente Jr et al. (2014) using a co-rotational formulation and concrete models provided by Eurocode 2 (2004). Theoretical models for concrete are also presented by ACI 318 (2014). Another numerical result is provided by Matias et al. (2017) to validate the Mazar’s damage model. Geometric properties, boundary conditions, adopted finite element mesh and the position of the nodes for the application and reading of forces and displacements are shown in Figure 13.

Figure 13
Geometric properties and discretization adopted for the column.

The column was discretized by 960 triangular finite elements with cubic approximation. The reinforcements were modeled with 900 linear truss elements. The Young’s modulus of concrete and steel are ${\mathbb{E}}_c$ = 3360 kN/cm² and ${\mathbb{E}}_s$= 21000 kN/cm2, respectively. The Poisson’s ratio of concrete and steel are assumed ${\nu }_c$ = 0.2 and ${\nu }_s$ = 0.0, respectively. Perfect elastoplastic behavior was adopted for steel rebars, with yield stress of ${\sigma }_{sy}$ = 46.5 kN/cm². The parameters of Mazar’s damage model were adopted following Parente Jr et al. (2014), adjusted to Eurocode 2 (2004), resulting in the following parameters: ${\epsilon }_{d0}$ = 8.65x10^-5, $A_C$ = 1.20, $B_C$ = 1500, $A_T$ = 0.50 and $B_T$ = 9000. Figure 14 shows the diagrams obtained for compressive and tensile stresses.

Figure 14
Uniaxial behavior of concrete: (a) compressive stresses and (b) tensile stresses.

As in previous example, we numerically simulated 10 specimens with aggregates randomly distributed in cement matrix, seeking an equivalent elastic modulus of ${\mathbb{E}}_c$ = 3360 kN/cm² and imposing a ratio of elastic modulus equal to ${\mathbb{E}}_m/{\mathbb{E}}_a$ = 0.20. For the volumetric rate of aggregates of 30%, we obtained the following values: ${\mathbb{E}}_m$ = 2240 kN/cm² and ${\mathbb{E}}_a$ = 11200 kN/cm².

Two numerical analysis were simulated with 50 incremental steps of displacements $\delta$ = -0.012 cm, observing the reaction force in the same point. The first analysis employed only concrete for matrix and trusses elements for rebars. The second analysis employed a cement matrix, coarse aggregates as particles, and rebars. Figure 15 shows the results of both analyses as well results provided by references.

Figure 15
Equilibrium trajectory of the column.

In the first analysis, without particle elements, the numerical model properly converged until step 19 (imposed displacement of $\delta$ = -0.228 cm). After that, there was a sudden reduction in the stiffness due to damage localization at the base of the column. This behavior is shown in Figure 16, in which the deformed configuration was enlarged by 10 times. Structural collapse is observed in step 20 at the bottom of the column. An ultimate force of $F_u$ = 406.31 kN is achieved at step 19, and maximum force of $F_{\max }$ = 454.00 kN is found at step 15. Table 3 provides a comparison of the maximum force achieved by the models proposed here and the references.

Figure 16
Deformed configuration and damaged column state.

Table 3
Analysis of the maximum force supported by the column.

References	Maximum force (kN)	Relative difference (%)	Relative difference (%)
		(without particles)	(with particles)
Without particles	454.00	-	0.90
With particles	449.97	0.89	-
Espion (1993)	450.19	0.84	0.05
Parente Jr et al. (2014)	460.56	1.44	2.35
Matias et al. (2017)	465.01	2.43	3.34

In the second analysis, the equilibrium path was smooth, and no numerical instabilities pointed out the collapse mechanism. Particle elements distributed the damage in the matrix, and the analysis converged a while longer. It was possible to impose a displacement of ${\delta }_H$ = 5.972 cm, corresponding to step 28. A maximum force of $F_{\max }$ = 449.97 kN was reached. Figure 17 shows the influence of the particle elements on the damage evolution at the base of the column.

Figure 17
Damage evolution at the base of the column.

Results obtained by the proposed model and references are quite close, validating the formulation and its accuracy. Besides, calibrating the damage parameters by normative instruments has proven to be viable in the absence of experimental data. In this example, both physical and geometric nonlinearities take place simultaneously. Therefore, the model is proven to reliable and consistent.

4.4 Reinforced concrete frame

In this fourth example, a reinforced concrete frame composed by two beams and two columns is modeled. The structure was experimentally tested by Vecchio and Emara (1992), with total height of 440 cm and span of 350 cm. Geometric characteristics and the adopted finite element mesh are shown in Figure 18. Figure 19 shows the cross sections for both beams and columns.

Figure 18
Geometric properties and discretizations of the frame.

Figure 19
Cross sections of the frame.

The finite element mesh was composed of 832 triangular elements with cubic approximation. Reinforcements were modeled by 800 linear truss elements for longitudinal reinforcements and 774 linear truss elements for transverse reinforcements. The following material properties were used: Young’s modulus of ${\mathbb{E}}_c$ = 3040 kN/cm² and ${\mathbb{E}}_s$ = 19250 kN/cm² for concrete and steel, respectively, Poisson's ratio of ${\nu }_c$ = 0.20 and ${\nu }_s$ = 0.0 for concrete and steel. The damage parameters for concrete were obtained in Nogueira (2010): ${\epsilon }_{d0}$ = 6.5x10^-5, $A_C$ = 0.9717, $B_C$ = 1204, $A_T$ = 0.9205 and $B_T$ = 10390. The constitutive behavior of the reinforcement is represented by an elastoplastic model with positive isotropic hardening, in which ${\sigma }_{sy}$ = 41.8 kN/cm² represents the yield stress of the material, $K$= 1925 kN/cm2 the tangent modulus and ${\sigma }_{su}$ = 59.8 kN/cm² the tension rupture of steel.

A second analysis employing particle elements to simulate the coarse aggregate was modeled with 144000 triangular finite elements with linear approximation. Geometry properties, dimensions and volumetric rate are similar to the previous examples. The distribution of particles in the frame is shown in Table 4.

Table 4
Particulate specimen characteristics.

Structural component	Total volume (cm3)	Volume of aggregates (cm3)	Number of aggregates	Number of finite elements	Number of nodes
Left column	528000	158400	5280	42240	47520
Right column	528000	158400	5280	42240	47520
Upper beam	372000	111600	3720	29760	33480
Lower beam	372000	111600	3720	29760	33480

As in previous example, to attain an equivalent Young’s modulus for concrete, imposing a ratio of elastic modulus equal to ${\mathbb{E}}_m/{\mathbb{E}}_a$ = 0.20, and a volumetric rate of aggregates of 30%, we obtained the following elastic modulus for the mortar matrix and coarse aggregates: ${\mathbb{E}}_m$= 1438.05 kN/cm2 and ${\mathbb{E}}_a$ = 7190.23 kN/cm².

There are two load steps. In the first step, a vertical force of 700 kN were applied gradually over 10 loading substeps at the top of the left column. In the second step, substeps of horizontal displacements increments of 0.05 cm were gradually imposed. The analysis run up to 3.30 cm (66 steps) for the model without inclusions and 3.10 cm (64 steps) for the model with particle elements. After that, damage severely affect the bearing capacity of the structure, and the model was not able to converge anymore.

Numerical results were compared with the experimental values of Vecchio and Emara (1992) and the numerical results of La Borderie et al. (1992) and Nogueira (2010). Nogueira (2010) presented two results - one using Bernoulli kinematics (B) and other considering Timoshenko model with the contribution of transverse reinforcement and dowel action (TSD). The horizontal force vs horizontal displacement at the top of the left column are shown in Figure 20.

Figure 20
Equilibrium trajectory of the frame.

The results obtained with the proposed model are in agreement with other numerical models, especially with the model of La Borderie et al. (1992). We also provide a comparative analysis of the ultimate force, as seen in Table 4. The proposed with and without particle elements resulted in 329.30 kN and 333.43 kN, respectively, near the experimental result obtained by Vecchio and Emara (1992).

Table 4
Analysis of the maximum force supported by the frame.

References	Maximum force (kN)	Relative difference (%)	Relative difference (%)
		(without particle)	(with particle)
Without particles	333.43	-	1.25
With particles	329.30	1.24	-
B - Nogueira (2010)	367.60	10.25	11.63
TSD - Nogueira (2010)	375.70	12.68	14.09
La Borderie et al. (1992)	306.90	7.96	6.80
Vecchio and Emara (1992)	332.30	0.34	0.91

Figure 21 shows the mapping of the damage variable with particle elements for steps 20 and 66, associated to the horizontal displacements of 1.00 cm and 3.30 cm, respectively. Notice that at step 66, the frame is about to collapse. Besides, both models (with and without particles) provides quite close results.

Figure 21
Damage to the frame.

5 CONCLUSIONS

In this paper the simulation of structures considering the inclusion of particle elements with the embedded technique has been successfully tested for physical and geometric nonlinear behaviors. The mechanical responses obtained with the proposed method were validated with several results of other research groups, showing accurate description and excellent results without over increasing computational cost, as the number of degrees of freedom is maintained with the embedeed particle elements. Based on the analysis of the results, the following conclusions are presented:

1. The embedded technique proved to be a reliable tool in the simulation of the mesoscale concrete constituent phases. The generation of independent meshes provides versatility to the simulation, enabling random inclusions distributed in the domain. With this, it is possible to represent coarse aggregates along the cementitious matrix.

2. The developed model can be used to estimate an equivalent Young’s modulus of the concrete, considering its phases. However, the relationship between the aggregate dimensions and the specimen is still a concern.

3. In the second example, higher loads are registered for the same displacement in the nonlinear stage of the simulation with particle elements. This behavior can be justified by the local stiffening provided by the coarse aggregate, which redistribute damage in the finite elements around it.

4. Decrease in maximum force was noticed in the last two examples when considering inclusions. In the third example, the maximum force decreased 0.89%, from 454.00 to 449.47 kN. Meanwhile, in the fourth example, the maximum force decreased about 1.24%, from 333.43 to 329.30 kN. The results are quite close to the experimental reference. When inclusions are considered, a difference of 0.05% and 0.91% is found in the third and fourth examples for the maximum force, respectively.

5. The damage model is capable of representing the preferred cracking zones due to the location of the scalar damage variable. When particles are included in the simulation, new paths of damage surrounding the aggregate are developed, indicating that the model can be improved if ITZ is included.

Therefore, it was possible to simulate the concrete at mesoscale, considering coarse aggregates and rebars in the mortar matrix. Results shown the applicability, reliability and new possibilities of the proposed model in the simulation of reinforced concrete structures subject to nonlinearities. In the next step, ITZ will be included, as well as the geometry, size fractions (grading) and distribution of coarse aggregate will be discussed.

Acknowledgements

The research supported by the Brazilian National Council for Scientific and Technological Development (CNPq 133981/2018-5, 310564/2018-2 and 428762/2018-2) is gratefully acknowledged. This study was also financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001.

References

• ACI 318. (2014). Building Code Requirements for Structural Concrete (ACI 318M-14) and Commentary (ACI 318RM-14). Farmington Hills.
• ASTM C39/C39M-20. (2020). Standard Test Method for Compressive Strength of Cylindrical Concrete Specimens. ASTM International, West Conshohocken.
• Bažant, Z. P. (1976). Instability, ductility, and size effect in strain-softening concrete. ASCE Journal of the Engineering Mechanics Division 102:331-344.
• Bažant, Z. P., Oh, B. H. (1983). Crack band theory for fracture of concrete. Matériaux et construction 16:155-177.
• Bitencourt Jr, L. A., Manzoli, O. L., Prazeres, P. G., Rodrigues, E. A., Bittencourt, T. N. (2015). A coupling technique for non-matching finite element meshes. Computer Methods in Applied Mechanics and Engineering 290:19-44.
• Bonet, J., Wood, R. D., Mahaney, J., Heywood, P., (2000). Finite element analysis of air supported membrane structures. Computer Methods in Applied Mechanics and Engineering 190:579-595.
• Boone, T. J., Wawrzynek, P. A., Ingraffea, A. R. (1986). Simulation of the fracture process in rock with application to hydrofracturing. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts 23:255-265.
• Borges, D. C., Quaresma, W. G., Fernandes, G. R., Pituba, J. J. C. (2017). Evaluation of a proposed model for concrete at mesoscopic scale. IBRACON Structures and Materials Journal 10(5):1087-1099.
• Borges, D. C., Pituba, J. J. C. (2017). Analysis of quasi-brittle materials at mesoscopic level using homogenization model. Advances in Concrete Construction 5(3):221-240.
• Caballero, A., López, C. M., Carol, I. (2006). 3D meso-structural analysis of concrete specimens under uniaxial tension. Computer Methods in Applied Mechanics and Engineering 195(52):7182-7195.
• Carreira, D. J., Chu, K. H. (1985). Stress-strain relationship for plain concrete in compression. Journal Proceedings 82:797-804.
• Carreira, D. J., Chu, K. H. (1986). Stress-strain relationship for reinforced concrete in tension. Journal Proceedings 83:21-28.
• Chen, W. F., Han, D. J. (2007). Plasticity for structural engineers, J. Ross Publishing.
• Chen, H., Xu, B., Mo, Y. L., Zhou, T. (2018). Behavior of meso-scale heterogeneous concrete under uniaxial tensile and compressive loadings. Construction and Building Materials 178:418-431.
• Coda, H. B., Greco, M. (2004). A simple FEM formulation for large deflection 2D frame analysis based on position description. Computer Methods in Applied Mechanics and Engineering 193:3541-3557.
• De Borst, R., Nauta, P. (1985). Non-orthogonal cracks in a smeared finite element model. Engineering Computations 2:35-46.
• Du, X., Jin, L., Ma, G. (2014). Numerical simulation of dynamic tensile-failure of concrete at meso-scale. International Journal of Impact Engineering 66:5-17.
• Espion, B. (1993). Benchmark examples for creep and shrinkage analysis computer programs: creep and shrinkage of concrete. In: Rilem Proceedings. Chapman & Hall.
• Etse, G., Caggiano, A., Vrech, S. (2012). Multiscale failure analysis of fiber reinforced concrete based on a discrete crack model. International Journal of Fracture 178:131-146.
• Eurocode 2. (2004). 1-1. 2004. BS EN 1992: 2004 - Eurocode 2: Design of concrete structures, Part 1-1: General rules and rules for buildings. British Standards Institution, London.
• Forti, T., Batistela, G., Forti, N., Vianna, N. (2020). 3D Mesoscale Finite Element Modelling of Concrete under Uniaxial Loadings. Materials 13(20):4585.
• Häfner, S., Eckardt, S., Luther, T., Könke, C. (2006). Mesoscale modeling of concrete: Geometry and numerics. Computers & structures 84(7):450-461.
• Hillerborg, A., Modéer, M., Petersson, P. E. (1976). Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cement and Concrete Research 6:773-781.
• Holzapfel, G. A. (2004). Nonlinear Solid Mechanics: A Continuum Approach for Engineering Science, Chichester, John Wiley & Sons Ltd.
• Kachanov, L. M., (1958). On the time to failure under creep conditions, izv. AN SSSR 8:1-8.
• Kim, S. M., Al-Rub, R. K. A. (2011). Meso-scale computational modeling of the plastic-damage response of cementitious composites. Cement and Concrete Research 41(3):339-358.
• La Borderie, C., Mazars, J., Pijaudier-Cabot, G. (1992). Response of plain and reinforced concrete structures under cyclic loadings. ACI 134:147-172.
• Lanczos, C. (1986). The Variational Principles of Mechanics, New York, Dover.
• Lemaitre, J. (1992). A Course on Damage Mechanics, Berlim, Springer.
• Lemaitre, J., Chaboche, J. C. (1985). Mécanique des matériaux solids. Paris, Dunod-Bordas.
• Lemaitre, J., Dufailly, J. (1987). Damage measurements. Engineering Fracture Mechanics 28:643-661.
• Li, K. Q., Li, D. Q., Li, P. T., Liu, Y. (2019). Meso-mechanical investigations on the overall elastic properties of multi-phase construction materials using finite element method. Construction and Building Materials 228:116727.
• Lloberas‐Valls, O., Rixen, D. J., Simone, A., Sluys, L. J. (2012). Multiscale domain decomposition analysis of quasi‐brittle heterogeneous materials. International Journal for Numerical Methods in Engineering 89:1337-1366.
• Manzoli, O. L., Gamino, A. L., Rodrigues, E. A., Claro, G. K. S. (2012). Modeling of interfaces in two-dimensional problems using solid finite elements with high aspect ratio. Computers & Structures 94:70-82.
• Matias, B. S., Parente Jr, E., Araújo, T. D. P. (2017). Use of Mazars constitutive model in static nonlinear analysis of reinforced concrete flat frames. In: XXXV Ibero-Latin American Congress on Computational Methods in Engineering, Florianopolis.
• Mazars, J. (1984). Application de la mécanique de l’endommagement au comportement non linéaire et à la rupture du béton de structure. Ph.D. Thesis, Université Paris 6, France.
• Mehta, P. K., Monteiro, P. J. M. (2014). Concrete: Microstructure, Properties, and Materials (in Portuguese), São Paulo, IBRACON.
• Meng, Q. X., Lv, D., Liu, Y. (2020). Mesoscale computational modeling of concrete-like particle-reinforced composites with non-convex aggregates. Computers & Structures 240:106349.
• Nguyen, V. P., Lloberas-Valls, O., Stroeven, M., Sluys, L. J. (2011). Homogenization-based multiscale crack modelling: from micro-diffusive damage to macro-cracks. Computer Methods in Applied Mechanics and Engineering 200:1220-1236.
• Nogueira, C. G. (2010). Development of mechanical, reliability and optimization models for application in reinforced concrete structures. Ph.D. Thesis (in Portuguese), University of São Paulo, Brazil.
• Paccola, R. R., Coda, H. B. (2016). A direct FEM approach for particulate reinforced elastic solids. Composite Structures 141:282-291.
• Paccola, R. R., Piedade Neto, D., Coda, H. B. (2015). Geometrical non-linear analysis of fiber reinforced elastic solids considering debounding. Composite Structures 133:343-357.
• Parente Jr, E., Nogueira, G. V., Meireles Neto, M., Moreira, L. S. (2014). Material and geometric nonlinear analysis of reinforced concrete frames. IBRACON Structures and Materials Journal 7:879-904.
• Rabotnov, Y. N. (1969). Creep problems in structural members. North-Holland.
• Radtke, F. K. F., Simone, A., Sluys, L. J. (2010). A computational model for failure analysis of fibre reinforced concrete with discrete treatment of fibres. Engineering fracture mechanics 77:597-620.
• Radtke, F. K. F., Simone, A., Sluys, L. J. (2011). A partition of unity finite element method for simulating non‐linear debonding and matrix failure in thin fibre composites. International Journal for Numerical Methods in Engineering 86:453-476.
• Ramos, É. S., Carrazedo, R. (2020). Cross-section modeling of the non-uniform corrosion due to chloride ingress using the positional finite element method. Journal of the Brazilian Society of Mechanical Sciences and Engineering 42(10):1-18.
• Rodrigues, E. A., Manzoli, O. L., Bitencourt Jr, L. A. (2020). 3D concurrent multiscale model for crack propagation in concrete. Computer Methods in Applied Mechanics and Engineering 361:112813.
• Rodrigues, E. A., Manzoli, O. L., Bitencourt Jr, L. A., Bittencourt, T. N. (2016). 2D mesoscale model for concrete based on the use of interface element with a high aspect ratio. International Journal of Solids and Structures 94:112-124.
• Rodrigues, E. A., Manzoli, O. L., Bitencourt Jr, L. A., Bittencourt, T. N., Sánchez, M. (2018). An adaptive concurrent multiscale model for concrete based on coupling finite elements. Computer Methods in Applied Mechanics and Engineering 328:26-46.
• Rots, J. G., Nauta, P., Kuster, G. M. A., Blaauwendraad, J. (1985). Smeared crack approach and fracture localization in concrete. HERON 30:1-48.
• Sampaio, M. S. M., Paccola, R. R., Coda, H. B. (2013). Fully adherent fiber-matrix FEM formulation for geometrically nonlinear 2D solid analysis. Finite Elements in Analysis and Design 66:12-25.
• Sánchez, M., Manzoli, O. L., Guimarães, L. J. (2014). Modeling 3-D desiccation soil crack networks using a mesh fragmentation technique. Computers and Geotechnics 62:27-39.
• Santos, S. B. S. (2015). An isotropic scalar damage model application in the analysis of reinforced concrete structures, Master Thesis (in Portuguese), University of Brasília, Brazil.
• Simo, J. C., Hughes, T. J. R. (1998). Computational Inelasticity, New York, Springer‐Verlag.
• Souza Neto, E. A., Peric, D., Owen, D. R. (2011). Computational methods for plasticity: theory and applications. John Wiley & Sons.
• Swartz, S. E., Jones, G. L., Hu, K. K. (1978). Compliance monitoring of crack growth in concrete. Journal of the Engineering Mechanics Division 104:789-800.
• Tang, L., Zhou, W., Liu, X., Ma, G., Chen, M. (2019). Three-dimensional mesoscopic simulation of the dynamic tensile fracture of concrete. Engineering Fracture Mechanics 211:269-281.
• Thilakarathna, P. S. M., Baduge, K. K., Mendis, P., Vimonsatit, V., Lee, H. (2020). Mesoscale modelling of concrete - a review of geometry generation, placing algorithms, constitutive relations and applications. Engineering Fracture Mechanics 231:106974.
• Unger, J. F., Eckardt, S. (2011). Multiscale modeling of concrete. Archives of Computational Methods in Engineering 18:341-393.
• Vanalli, L., Paccola, R. R., Coda, H. B. (2008). A simple way to introduce fibers into FEM models. Communications in Numerical Methods in Engineering 24:585-603.
• Van Mier, J. G. M. (1997). Fracture processes of concrete: assessment of material parameters for fracture models, CRC Press.
• Vecchio, F. J., Emara, M. B. (1992). Shear deformations in reinforced concrete frames. ACI Structural Journal 89:46-56.
• Wriggers, P., Moftah, S. O. (2006). Mesoscale models for concrete: Homogenisation and damage behaviour. Finite Elements in Analysis and Design 42(7):623-636.
• Yılmaz, O., Molinari, J. F. (2017). A mesoscale fracture model for concrete. Cement and Concrete Research 97:84-94.
• Zhang, J., Wang, Z., Yang, H., Wang, Z., Shu, X. (2018). 3D meso-scale modeling of reinforcement concrete with high volume fraction of randomly distributed aggregates. Construction and Building Materials 164;350-361.
• Zhou, R., Song, Z., Lu, Y. (2017). 3D mesoscale finite element modelling of concrete. Computers & Structures 192:96-113.
• Zhou, Y., Jin, H., Wang, B. (2019). Modeling and mechanical influence of meso-scale concrete considering actual aggregate shapes. Construction and Building Materials 228:116785.

Edited by

Editor:

Marco L. Bittencourt.

Publication Dates

• Publication in this collection
  10 Mar 2021
• Date of issue
  2021

History

• Received
  08 Sept 2020
• Reviewed
  03 Nov 2020
• Accepted
  24 Dec 2020
• Published
  15 Jan 2021