Modeling Particles Elements in Damaged Reinforced Concrete Structures

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. Positional FEM, Damage model, Elastoplasticity.


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 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  and . 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 andAl-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. 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.

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.

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: Equilibrium is given by the Principle of Conservation of Mechanical Energy, or Principle of Stationary, defined as (Lanczos, 1986): in which i Y is the position vector. Given the arbitrariness of i Y δ , the Eq. (2) can be rewritten as: 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 i g . Equilibrium or current configuration occurs when ( ) 0 The current position is determined by resorting to Taylor series expansion for a trial position trial k Y , considering only the first order terms, resulting in: where j Y ∆ is the correction of the trial positions and ij H 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 where the initial trial position corresponds to the undeformed configuration, that is, trial The convergence criterion is based on the norm of positions, in which the equilibrium condition is fully satisfied when result less or equal to a previously stablished tolerance.

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: in which 0 i f and 1 i f maps of the initial and current configurations, respectively, to the dimensionless coordinates 1 ξ and 2 ξ , l φ are the shape functions, and li X and li Y 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): in which E is the Green-Lagrange strain tensor, A is the gradient of the change of configuration, also called deformation gradient, and δ 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): in which 0 A and 1 A are the gradients of the mappings from the dimensionless coordinates to the initial and current configurations, respectively, given by: The specific strain energy ( e u ) is integrated into the initial volume of the solid ( 0 V ) to find the strain energy (U ), both written according to Eq. (10).
in which  is the elastic constitutive tensor of the material, given by: where µ and λ are Lamé constants which can be determined as a function of the longitudinal elastic modulus  and the Poisson's ratio ν , according to the following expressions: ( ) ( )( ) and 2 1 1 1 2 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): Particularly, the second Piola-Kirchhoff stress tensor is written for state of plane stress and plane strain by: 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: in which n denotes the number of integration points, with weights given by h w , and 0 J is Jacobian of the initial configuration, calculated as The specific strain energy derivative can be further developed using Eq. (13), leading to: 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.
in which the derivative of the gradient of the mapping of the current configuration with respect to the positions is given by: Latin American Journal of Solids and Structures, 2021, 18(1), e345 6/24 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: in which the second derivative of specific strain energy is calculated by deriving Eq. (16) with respect to positions, leading to: 2 2 : : : The terms ij Y ∂ ∂ E ,  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).
is written as follows:

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: where 0 L and L correspond to the initial and current length, in which indices 1 and 2 refer to directions, and α and β 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: Since only linear truss element is used, solution is obtained analytically. Thus, the strain energy is written according to Eq. (24), with 0 0 0 Latin American Journal of Solids and Structures, 2021, 18(1), e345 7/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: Notice that E and S reduce to scalars in one-dimensional case. The Green-Lagrange strain derivative with respect to the nodal parameters is given by: The Hessian tensor is determined from strain energy second derivative with respect to the positions, given by: Using Eq. (25) in Eq. (27) and applying the product rule and chain rule, we have: 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 S , while the second derivative results in the material longitudinal elastic modulus  . In this sense, the Hessian is written as:

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: where p j X and p j Y are the initial and current positions of node p of the reinforcement in j direction, i φ are the shape functions applied to the dimensionless coordinates of node p , and ˆi j X and ˆi j Y 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 ( ).
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 ( trial j ξ ), and disregarding higher order terms, it is written: where j ∆ξ is the correction of dimensionless coordinates and ij L corresponds to a 2x2 matrix. The vector of dimensionless coordinates is then updated ( trial 1 trial k k k + ξ = ξ + ∆ξ ) and the process is repeated until a tolerance is achieved ( k k tol ∆ξ ∆ξ ≤ ). The resulting dimensionless coordinates must belong to the matrix dimensionless domain, that is, Thus, the internal forces of the composite are written as: 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: The Hessian matrix is also decomposed additively, including both reinforcement and matrix stiffness, as shown in the following expression: 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.

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 ( eq E ), as follows: where ( ) i E + 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: where variable lim E is the maximum elongation strain in the strain history. At the beginning of the process lim E is equivalent to the strain corresponding to the concrete tensile strength 0 d E . Since the formulation proposed here employs Green-Lagrange strain, it is converted as follows: in which 0 d ε represents the specific strain from a uniaxial tensile test.
The model also establishes two damage variables T D and C D due the concrete behavior related to tensile and compressive stresses. These variables describe the nonlinear stress-strain diagram for concrete, and are given by: where T A , C A , T B and C B are parameters of the Mazar's damage model.
For multiaxial stresses, damage is evaluated by linear combination of the variables, as follows: Coefficient T α and C α are calculated by decomposing the principal strains into two parts, according to Eqs. (43) and (44): From that, T α and C α are calculated as follows: where 0 1 T ≤ α ≤ ; 0 1 C ≤ α ≤ ; and 1 T C α + α = . The variable ( ) V E + represents the total state of elongation, given by the principal strains as follows: Stresses are found by the damaged elastic properties of the material, according to Eq. (47): where S  is the second Piola-Kirchhoff stress tensor,  is the fourth-order elastic constitutive tensor and E is the Green-Lagrange strain tensor.

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: where s S is stress acting on steel, sy S is the steel yield stress, K is the plastic hardening module and χ is the equivalent plastic strain measure.
In the one-dimensional case, sy S is determined from the Cauchy stress and the Piola-Kirchhoff stress, as follows: where s  is Young's modulus of steel, sy σ is the steel yield Cauchy stress and sy ε is the associated strain.
The increase in stress in the inelastic phase is given as follow: 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).

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.

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: m  = 20 GPa, m ν = 0.2, a  = 100 GPa, 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 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). 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. 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.    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).
where c  is the equivalent Young's modulus, m  and a  are the Young's modulus of the cementitious matrix and the coarse aggregates, m c  and a c are the volumetric rates of cement paste and aggregates, respectively. Numerical results of the obtained Young's modulus are presented in Table 2. 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.

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.  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 δ = -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. 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. We simulated 10 specimens with random distribution of the aggregates, seeking an equivalent elastic modulus of c  = 3000 kN/cm 2 and imposing a ratio of elastic modulus equal to m a   = 0.20. For the volumetric rate of aggregates of 30%, we obtained m  = 1990.21 kN/cm 2 and a  = 9951.05 kN/cm 2 . No changes to damage parameters C A , C B , T A and T B were made. Only 0 d ε was adjusted, to impose that degradation develops at same stress level in both materials.
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.  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. 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.

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 corotational 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 14 shows the diagrams obtained for compressive and tensile stresses. Two numerical analysis were simulated with 50 incremental steps of displacements δ = -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. In the first analysis, without particle elements, the numerical model properly converged until step 19 (imposed displacement of δ = -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 u F = 406.31 kN is achieved at step 19, and maximum force of max F = 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.  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 H δ = 5.972 cm, corresponding to step 28. A maximum force of max F = 449.97 kN was reached. Figure 17 shows the influence of the particle elements on the 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.

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.  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. 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. 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)  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.

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. 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.