A direct FEM approach to model mesoscale concrete and connect non-matching meshes in multiscale analysis

Uma abordagem direta pelo MEF para modelar concreto em mesoescala e conectar malhas não conformes em análises multiescala

Abstracts

Abstract

The mechanical degradation of concrete structures is a phenomenon dependent on the material heterogeneity observed at mesoscale. As the mechanical degradation is a localized phenomenon, structural members and structures may be simulated using the concurrent multiscale analysis technique. Thus, only the most critical regions are modeled in mesoscale, reducing the computational cost compared to the simulation of the entire structure at this scale. This work presents two contributions in concurrent multiscale analysis. The first contribution introduces an alternative representation of the mesoscale interfacial transition zone (ITZ) of the concrete together with a strategy that allows modeling particles (coarse aggregates) without degrees of freedom. The resulting ITZ representation allows the simulation of more realistic discrete cracks in concrete modeling. The second contribution uses particle-like elements without degrees of freedom as coupling elements to model non-matching meshes between different media. The proposed coupling technique does not add degrees of freedom and does not use penalty or Lagrange Multipliers methods. Experimental and numerical results are used in order to validate the proposed multiscale formulation regarding concrete specimen simulations.

Keywords:
mesoscale concrete; interfacial transition zone; coupling element; multiscale; non-matching meshes

Resumo

A degradação mecânica de estruturas de concreto é um fenômeno dependente da heterogeneidade do material observada em mesoescala. Como a degradação mecânica é um fenômeno localizado, elementos estruturais e estruturas podem ser simulados usando a técnica de análise multiescala concorrente. Assim, apenas as regiões mais críticas são modeladas em mesoescala, reduzindo o custo computacional em comparação com a simulação de toda a estrutura nessa escala. Este trabalho apresenta duas contribuições para a estratégia de análise multiescala concorrente. A primeira contribuição apresenta uma representação alternativa da zona de transição interfacial da mesoescala (ITZ) do concreto juntamente com uma estratégia que permite modelar partículas (agregados graúdos) sem graus de liberdade. A representação da ITZ resultante permite a modelagem de fissuras discretas mais realistas na simulação de concreto. A segunda contribuição usa elementos semelhantes a agregados sem graus de liberdade como elementos de acoplamento para modelar malhas não conformes entre diferentes meios. A técnica de acoplamento proposta não adiciona graus de liberdade e não utiliza métodos de penalização ou multiplicadores de Lagrange. Resultados experimentais e numéricos são usados para validar a formulação multiescala proposta aplicada em simulações de amostras de concreto.

Palavras-chave:
concreto em mesoescala; zona de transição interfacial; elemento de acoplamento; multiescala; malhas não conformes

1 INTRODUCTION

It is well known the fact that composites mix two or more materials in order to obtain a new material with physical and economic properties better than those presented separately by each constituent material. Particle composites are characterized by the random presence of particles with different sizes and shapes immersed inside a continuous matrix. Some particle composites are glass reinforced with mica flakes, brittle polymers reinforced with rubberlike particles, and concrete [¹].

Concrete is a particulate composite widely used around the world and improving knowledge of its mechanical behavior is very important. In the analysis of structures, it is common to admit that concrete is a homogeneous material, but this hypothesis may be not sufficient when considering its nonlinear mechanical behavior. This happens because the nonlinear behavior is strongly influenced by the heterogeneity of the concrete. When the objective is to perform an accurate analysis of the mechanical behavior of structures at the macroscale, the representation of concrete at the mesoscale proves to be quite adequate [²], [³].

At the mesoscale, concrete is analyzed as a composite of particles (called coarse aggregates), mortar and the ITZ between these phases. The consideration of ITZ is important because it plays an important role in the degradation process due to its lower strength than the other components [⁴], [⁵]. The representation of concrete heterogeneity makes the mechanical model gain physical meaning [⁶] improving the mechanical behavior understanding and reducing the dependence of macroscale phenomenological constitutive models’ calibration.

Numerical mesoscale modeling uses simple constitutive models for each phase of the composite, allowing the representation of its complex responses [²]. Among the behaviors that can be captured are the location of cracks in the least resistant paths (ITZ) [⁷], [⁸], the size effects on structural strength [⁹], [¹⁰], and the random nature of the structural response [¹¹], [¹²].

Mesoscale concrete simulation techniques are classified according to the domain representation. Some modeling strategies represent the phases of the concrete in a discrete way, mainly using the Lattice Model, [¹³ ]–[¹⁵] and the particle model [¹⁶ ]–[¹⁸]. Other models represent the phases of concrete using continuous strategies and the finite element method (FEM) predominates. In the present study, we are interested in the FEM strategies that use interface elements. The interface elements represent the possible paths for cracks within the matrix and also the ITZ. Strategies that use zero thickness interface elements combined with cohesive models [²], [⁶], [¹⁹] and high aspect ratio interface elements combined with a damage model [⁷], [²⁰]–[²²] are commonly adopted. The interface elements are advantageous as they do not require the tracking step of the crack path, which is necessary in enrichment techniques [²³], [²⁴] or the XFEM approach [²⁵].

In order to consider the heterogeneity in mesoscale models, aggregates can be numerically generated [²], [¹¹], [¹³] or placed through imaging techniques that collect information on real structures [²⁶]–[²⁸]. When a statistical analysis of heterogeneity is desired (e.g., Monte Carlo simulations [⁹], [¹¹], [²⁹], [³⁰]) numerous samples must be provided by changing the position or shape of the coarse aggregates. A limitation of the above-mentioned mesoscale concrete FE-based simulation techniques is that matrix and particles are simultaneously created in a unique mesh. Thus, for each new numerical case (to analyze the same specimen) a completely new mesh must be created with new boundary conditions numbering.

Paccola and Coda [³¹] developed the direct FEM approach for particulate reinforced solids that allows representing matrix and particles by different meshes and with non-coincident nodes. In this strategy, particles are embedded inside the continuous mesh of the matrix material and do not add degrees of freedom to the problem. The strategy allows the matrix mesh and boundary conditions to remain unchanged, regardless of the position and shape of particles. Thus, it is possible to reuse the matrix mesh and the boundary conditions when changing the particles mesh for new particulate specimen analyses. Paccola and Coda [³¹] showed the efficiency of the direct FEM approach for elastic problems and Ramos et al. [³²] extended this representation to model reinforced concrete structures using a distributed damage model, however without representing the ITZ. The formulation presented by Paccola and Coda [³¹] is developed in the context of the FEM positional approach [³³], [³⁴]. Thus, the formulation is geometrically exact and allows more precise analyzes than the models formulated with the hypothesis of small displacements. It is worth mentioning that this particulate embedding strategy [³¹] has been extensively tested for fiber reinforcements in references [³⁵ ]–[³⁸].

The computational and memory costs in mesoscale numerical simulations are high, limiting analyses to small samples when simple workstations are employed. In this sense, the use of concurrent multiscale techniques [²], [³⁹] - that divide the structure into subdomains employing different discretization scales - is justified. In concurrent multiscale techniques, only subdomains in the most critical regions are modeled with an explicit representation of material’s heterogeneity. The other subdomains are modeled considering equivalent homogeneous material. The representation of concrete through concurrent multiscale analysis is successfully used by [²], [²¹], [³⁹], [⁴⁰], and hierarchical modeling is adopted by [⁴¹], for example. The two major challenges of concurrent multiscale technique are coupling the scales correctly and determining the domain location that must be discretized on the mesoscale [²].

In concurrent multiscale analysis, the connection between the meshes of different scales can be made through transition elements. In this case, the dimension of the finite elements near the interface between scales decreases in size until they reach the dimension of the elements of the smallest scale with coincident nodes. However, the quality of the problem response in this region can be impacted by the presence of distorted elements [⁴²]. Coupling techniques that allow direct connection between non-matching meshes avoid the presence of distortions and reduce the number of elements, since this transition region no longer exists. Some strategies use penalty methods [⁴³], [⁴⁴] to couple the subdomains. The use of penalty does not add degrees of freedom but may result in ill-conditioned systems. Other strategies use Lagrange multipliers to enforce continuity as, for example, methods of Mortar [⁴⁵ ]–[⁴⁷] and of Arlequin [⁴⁸], [⁴⁹]. A disadvantage of Lagrange multipliers is the increase of degrees of freedom. A mixed approach (Nitsche’s method [⁵⁰]), in which the Lagrange Multipliers are replaced by their physical representation and an extra penalty-like term, is applied to handle non-matching meshes in isogeometric analysis [⁵¹], [⁵²]. Another strategy adopts coupling non-conventional elements created from elements of each subdomain in the region to be coupled [⁵³], [⁵⁴]. In this case, the adjacent elements are modified in order to convert the non-matching interface into a matching one and, as mentioned by Bitencourt et al. [⁴⁴], some additional computational cost to build these elements appears.

In this work, the direct FEM embedding approach for particulate reinforced solids developed by Paccola and Coda [³¹] is extended to represent the coarse aggregate in a mesoscale model of concrete and to connect non-matching meshes in concurrent multiscale analysis. The developed mesoscale concrete model includes the mortar (matrix mesh), coarse aggregate and ITZ phase. The interesting property of the direct FEM embedding approach to be explored is its capacity of a selective switching among immersed meshes (coarse aggregate). The present work considers also the exact nonlinear geometric behavior in the formulation. Rabczuk and Belytschko [⁵⁵] observed that it is important to consider the geometric nonlinearity in modeling the degradation of reinforced concrete frame corners.

The main contributions of the present study can be summarized as: (i) development of a strategy for creating the ITZ around coarse aggregates without degrees of freedom embedded in the mortar matrix, (ii) non-matching meshes (macroscale and mesoscale meshes) are connected by the direct FEM embedding [³¹], (iii) the degrees of freedom of the fine mesh used to describe mortar condenses the coarse aggregate immersed discretization eliminating its degrees of freedom and (iv) the natural invariance of the mesh numbering that allows analyzing different random distributions of aggregates, with different volumetric fractions, without the need to change the matrix mesh.

It is interesting to mention that a damage model [⁵⁶] is applied here to all interface elements (including ITZ elements), so that they become responsible for the degradation of the material. The advantage of using high aspect ratio interface elements is that they are integrated as regular elements and allow representing discrete crack formation. Moreover, the couple of macroscale and mesoscale meshes in concurrent multiscale analysis by the direct FEM embedding does not add degrees of freedom and does not result in ill-conditioning, as the same mechanical properties of the material are used.

The paper is organized as follows: Section 2 presents the positional geometric nonlinear approach, including the direct FEM embedding strategy. Section 3 shows the mesoscale modeling including the coarse aggregate definition, the interface finite elements (including ITZ elements), the employed damage model and the definition of matrix-aggregate interface regions. In Section 4, some examples are presented in order to validate and show the range of applications of the proposed model. Finally, the main conclusions are drawn in Section 5.

2 GEOMETRICALLY NONLINEAR FORMULATION OF COMPOSITE SOLIDS

This section presents the formulation of the FEM approach used in this work and the particle inclusion strategy via kinematics compatibility (embedding strategy). This strategy is necessary to represent the coarse aggregates without adding degrees of freedom and to create the coupling elements for non-matching meshes.

2.1 Geometric nonlinear equilibrium

The implemented computational code uses the positional approach of FEM for which the unknowns are nodal coordinates instead of displacements. The positional FEM is a naturally non-linear geometric version of the FEM that firstly appear in Bonet et al. [³³] and was generalized for various applications in [³¹], [³⁴], [³⁷], [³⁸], for example. The elastic version of the positional FEM is based on the Variational Principle of Mechanics [⁵⁷], using positions as the main variable. This work is restricted to static problems, thus the mechanical energy ($\Pi$), Equation 1, includes the potential energy of conservative external forces ($P$) and the strain energy ($U$):

$\Pi =U+P$ (1)

A general mechanical energy variation using positions $\overrightarrow{Y}$ as variables is expressed in Equation 2.

$\frac{\partial \Pi }{\partial \overrightarrow{Y}}\cdot \delta \overrightarrow{Y}=\frac{\partial U}{\partial \overrightarrow{Y}}\cdot \delta \overrightarrow{Y}+\frac{\partial P}{\partial \overrightarrow{Y}}\cdot \delta \overrightarrow{Y}=0$ (2)

Considering that $\delta \overrightarrow{Y}$ is arbitrary and recalling the force/position energy conjugate, the nonlinear equilibrium equations, as described by [³¹], [³⁷], [³⁸], are written in Equation 3.

$g_i=\frac{\partial U}{\partial Y_i}+\frac{\partial P}{\partial Y_i}=F_i^{int}\left(\overrightarrow{Y}\right)-F_i^{ext}=0_i$ (3)

in which $F_i^{int}$ is the internal force (a function of the nodal position vector $\overrightarrow{Y}$) and $F_i^{ext}$ is the external force. Vector $g_i$ is understood as an unbalanced mechanical residuum that is not null if, for an iterative solver, the assumed position $\overrightarrow{Y}$ is a trial position ${\overrightarrow{Y}}^{trial}$. The adopted iterative solver is the Newton-Raphson method that is written from a truncated Taylor expansion of $g_i$, i.e.:

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

From which one finds the position correction

$\Delta Y_j=-{\left({\left.\frac{\partial g_i}{\partial Y_j}\right|}_{{\overrightarrow{Y}}^{trial}}\right)}^{-1}{g}_i\left({\overrightarrow{Y}}^{trial}\right)$ (5)

in which $\partial g_i/\partial Y_j$ is called Hessian matrix. The position is updated as ${\overrightarrow{Y}}^{trial}={\overrightarrow{Y}}^{trial}+\Delta \overrightarrow{Y}$ until $‖\Delta \overrightarrow{Y}‖/‖\overrightarrow{X}‖$ becomes smaller than a tolerance, in which $\overrightarrow{X}$ is the initial position vector.

To solve Equation 5 it is necessary to calculate the internal force and the Hessian matrix for a given specific strain energy potential. As this study is worried with concrete, only small strains take place, thus the Saint-Venant-Kirchhoff potential is adopted to represent the elastic behavior and its potential is given in Equation 6.

$u\left(E\right)=\frac{1}{2}{E}_{ij}{\mathfrak{C}}_{ijkl}{E}_{kl}$ (6)

In which ${\mathfrak{C}}_{ijkl}$ is the elastic constitutive tensor (similar to the generalized Hooke’s law) and $E_{ij}$ is the Green-Lagrange strain. The energy conjugate stress of the Green-Lagrange strain is the second Piola-Kirchhoff stress, given by:

$S_{ij}=\frac{\partial u}{\partial E_{ij}}={\mathfrak{C}}_{ijkl}{E}_{kl}$ (7)

The constitutive tensor is also related with the strain energy in Equation 8.

${\mathfrak{C}}_{ijkl}=\frac{\partial u}{\partial E_{ij}\partial E_{kl}}$ (8)

For plane stress, the object of this study, the specific strain energy is written as:

$u\left(E\right)=\frac{G}{1-\nu }\left\{E_{11}^2+E_{22}^2+2\nu E_{11}{E}_{22}+(1-\nu )(E_{12}^2+E_{21}^2)\right\}$ (9)

in which $G$ is the shear elastic modulus and $\nu$ is the Poisson ratio.

From the specific strain energy in Equation 9, one writes the internal force and Hessian matrix as follows:

$F_i^{int}=\frac{\partial U}{\partial Y_i}={\int }_{V_0}\frac{\partial u}{\partial Y_i}d{V}_0={\int }_{V_0}{S}_{kl}\frac{\partial E_{kl}}{\partial Y_i}d{V}_0$ (10)

$H_{ij}=\frac{{\partial }^{2}U}{\partial Y_i\partial Y_j}={\int }_{V_0}\frac{\partial E_{oz}}{\partial Y_i}{\mathfrak{C}}_{ozkl}\frac{\partial E_{kl}}{\partial Y_j}+S_{mn}\frac{{\partial }^2{E}_{mn}}{\partial Y_i\partial Y_j}d{V}_0$ (11)

More details regarding the positional FEM basic equations can be seen in references [³¹], [⁵⁸].

2.2 Aggregates embedding via kinematic compatibility

In this subsection, based on studies dedicated to linear elasticity [³¹], we present the necessary equations and the adopted strategy to immerse aggregates meshes inside the matrix meshes and to couple non-matching meshes. To simplify the description of the kinematics compatibility steps, the finite elements of the coarse aggregates or coupling elements (non-matching meshes) are simply called reinforcement finite elements. The mesh that represents the matrix in the case of mesoscale analysis and the meshes that represent the two scales in the case of coupling non-matching meshes are simply called the matrix.

The approach consists of writing the nodal positions of the reinforcement elements using the shape functions of the matrix elements in which they are immersed. To do this, it is necessary to determine the dimensionless coordinates $({\xi }_1^P,{\xi }_2^P)$ of the reinforcement’s nodes $P$ inside the matrix elements at the initial configuration using the initial position mapping, as:

$X_{Pi}={\varphi }_{\mathcal{l}}({\xi }_1^P,{\xi }_2^P)X_{\mathcal{l}i}$ (12)

where ${\varphi }_{\mathcal{l}}$ are shape functions associated with node $\mathcal{l}$, $X_{\mathcal{l}i}$ are the initial coordinates ($i$ direction) of the matrix finite element for which a reinforcement node $P$ belongs and $X_{Pi}$ is the coordinates of the reinforcement node.

Equation 12 is solved for $({\xi }_1^P,{\xi }_2^P)$ iteratively using the Newton-Raphson method. The reinforcement node belongs to the matrix element if results are in the intervals $0\le {\xi }_1^P\le 1$, $0\le {\xi }_2^P\le 1$ and $0\le 1-{\xi }_1^P-{\xi }_2^P\le 1$. It is noteworthy that the method produces accurate responses and presents fast convergence. More details can be found in [³¹], [³²]. Once the dimensionless coordinate $({\xi }_1^P,{\xi }_2^P)$ is achieved the current position of reinforcement nodes $Y_{Pi}$ are written as a function of the current position of matrix nodes, Equation 13.

$Y_{Pi}={\varphi }_{\mathcal{l}}({\xi }_1^P,{\xi }_2^P)Y_{\mathcal{l}i}$ (13)

In order to find the reinforcement’s internal force and the Hessian matrix contributions regarding the matrix mesh one applies the chain rule. Thus, it is necessary to write the derivative of the current positions of the reinforcement nodes regarding the current position of matrix nodes as:

$\frac{\partial Y_{Pi}}{\partial Y_{\beta \alpha }}=\frac{\partial Y_{\mathcal{l}i}}{\partial Y_{\beta \alpha }}{\varphi }_{\mathcal{l}}({\xi }_1^P,{\xi }_2^P)={\delta }_{\beta \mathcal{l}}{\delta }_{\alpha i}{\varphi }_{\mathcal{l}}({\xi }_1^P,{\xi }_2^P)={\delta }_{\alpha i}{\varphi }_{\beta }({\xi }_1^P,{\xi }_2^P)$ (14)

in which ${\delta }_{\alpha i}$ is the Kronecker delta.

The total strain energy $U$ of the problem is the sum of the energies of the reinforcement $U_R$ and the matrix $U_{mat}$, that is:

$U=U_R+U_{mat}$ (15)

Using Equation 14, the internal force $F_{\beta \alpha }^{int}$ at node $\beta$ in direction $\alpha$ of the matrix mesh is given by:

$F_{\beta \alpha }^{int}=\frac{\partial (U_{mat}+U_R)}{\partial Y_{\beta \alpha }}=\frac{\partial U_{mat}}{\partial Y_{\beta \alpha }}+\frac{\partial U_R}{\partial Y_{\left(P\right)i}}\frac{\partial Y_{\left(P\right)i}}{\partial Y_{\beta \alpha }}=F_{\beta \alpha }^{mat}+{\varphi }_{\beta }({\xi }_1^P,{\xi }_2^P)F_{\left(P\right)\alpha }^R$ (16)

where the parenthesis in index $\left(P\right)$ denotes that there is no sum regarding it. Thus, the total internal force is the sum of two parts, the natural internal force of the matrix elements ($F_{\beta \alpha }^{mat}$) and the transferred (to matrix nodes) reinforcement’s internal force $F_{P\alpha }^R$.

From Equation 15 the Hessian matrix is achieved, Equation 17, by the second derivative of the total strain energy regarding the matrix’s nodes, i.e.:

$H=\frac{{\partial }^2(U_{mat}+U_R)}{\partial Y_{\beta \alpha }\partial Y_{\xi \gamma }}=H_{\beta \alpha \xi \gamma }^{mat}+H_{\beta \alpha \xi \gamma }^R$ (17)

for which, using twice the chain rule, results:

$H_{\beta \alpha \xi \gamma }^R=\frac{{\partial }^2{U}_R}{\partial Y_{\left(P\right)j}\partial Y_{\left(P\right)i}}\frac{\partial Y_{\left(P\right)j}}{\partial Y_{\xi \gamma }}\frac{\partial Y_{\left(P\right)i}}{\partial Y_{\beta \alpha }}+\frac{{\partial }^2{U}_R}{\partial Y_{\left(\bar{P}\right)j}\partial Y_{\left(P\right)i}}\frac{\partial Y_{\left(\bar{P}\right)j}}{\partial Y_{\xi \gamma }}\frac{\partial Y_{\left(P\right)i}}{\partial Y_{\beta \alpha }}+\frac{{\partial }^2{U}_R}{\partial Y_{\left(P\right)j}\partial Y_{\left(\bar{P}\right)i}}\frac{\partial Y_{\left(P\right)j}}{\partial Y_{\xi \gamma }}\frac{\partial Y_{\left(\bar{P}\right)i}}{\partial Y_{\beta \alpha }}+\frac{{\partial }^2{U}_R}{\partial Y_{\left(\bar{P}\right)j}\partial Y_{\left(\bar{P}\right)i}}\frac{\partial Y_{\left(\bar{P}\right)j}}{\partial Y_{\xi \gamma }}\frac{\partial Y_{\left(\bar{P}\right)i}}{\partial Y_{\beta \alpha }}$ (18)

where $P$ and $\bar{P}$ are reinforcement’s nodes (no summation), $\beta$ and $\xi$ are matrix’s nodes of one or two different elements and $\alpha$ and $\gamma$ represent directions.

It is important to mention that for overlapping aggregates, as in [³¹], the modulus of elasticity of the aggregates is taken as the difference between the real modulus of elasticity and the modulus of elasticity of the matrix. When using the formulation for non-matching meshes, the coupling elements have the real elastic modulus.

2.3 Coupling non-matching meshes

In this subsection, the proposed non-matching meshes coupling strategy is described.

Consider the uncoupled non-matching meshes of subdomains ${\Omega }_1$ (coarse mesh) and ${\Omega }_2$ (refined mesh) that share the boundary interface ${\Gamma }_{12}$, see Figure 1. Elements $e_1$ and $e_2$ and nodes 1, 2, 3 and 4 of Figure 1 are part of the subdomain ${\Omega }_1$ while all other nodes and elements belongs to subdomain${\Omega }_2$. The nodes that have the same coordinates, but belongs to different subdomains are independent, i.e., uncoupled when originally generated.

Figure 1
Subdomains of non-matching meshes.

The strategy consists in writing the positions of element nodes of ${\Omega }_2$ that are on ${\Gamma }_{12}$ as a function of positions and shape functions of element nodes of ${\Omega }_1$ where the referred nodes are immersed. In this way, the scales become strongly coupled. Considering the subdomains of Figure 1, the elements $e_3$ through $e_6$ (finer mesh) that have nodes belonging to the interface ${\Gamma }_{12}$, i.e., nodes 5 through 7, are converted into coupling elements. The idea of creating coupling elements starting from the elements of one of the subdomains has been proposed by [⁴⁴] in a different way of the one presented here, because we choose the smaller scale as reference and does not use penalty techniques [⁴⁴], [⁵⁹].

In Figure 2 the coupling elements are detached in order to promote a direct relation of the non-matching mesh coupling and the aggregate embedding technique [³¹] described in section 2.2. The coupling elements are now identified by $l_1$, $l_2$, $l_3$ and $l_4$ and can be understood as an aggregate whose nodes are written in terms of the nodes of the subdomains ${\Omega }_1$ and ${\Omega }_2$. However, in the subdomain ${\Omega }_2$ there are nodes that coincide, that is, the nodes $b$, $d$ and $f$ are simple nodes 8, 9 and 10 of the finer mesh. So, unlike the embedded aggregates presented in Paccola and Coda [³¹], here only part of the nodes of the coupling element has its degree of freedom written as a function of the degrees of freedom of other elements. In the case considered in Figure 2, only the nodes $a$, $c$ and $e$ fall inside elements $e_1$ and $e_2$, thus their internal force and Hessian matrix are written as a function of the coarse mesh nodes positions by Equations 16 and 18, respectively, i.e., there is no penalty constants. Furthermore, it is important to highlight that the coupling elements are identical to elements $e_3$ through $e_6$ of Figure 1, which are eliminated. In addition to the strong coupling between subdomains, the coupling elements fulfill the structural function of the original elements.

Figure 2
Creating non-matching mesh coupling elements.

When meshes with very different scales are coupled some spurious stress concentrations appear [²], but in the proposed strategy this effect is less important as the coupling elements have the same elastic properties of the homogeneous macro scale and are not allowed to degenerate, as they are far from the damaged region. As a final remark, it is interesting to note that the generation of coupling elements takes place entirely in the pre-processing stage.

It is important to mention that the same strategy used for overlapping aggregates, is applied to the case of overlapping meshes [³¹].

3 MESOSCALE MODELING

In this section, the main points necessary to carry out the mesoscale modeling are presented. They are: (i) the strategy of creation and positioning the coarse aggregates (particles); (ii) the geometric generation and adopted constitutive model of the interface elements; and (iii) the proposed strategy for the identification and setting of physical parameters for ITZ representation. In order to simplify the modeling description, three different nomenclatures will be defined to label mesh combinations:

• Overlapped mesh (OM): In this model, the aggregate’s mesh is immersed in the matrix mesh, without adding degrees of freedom, there are both interface elements (mortar) and ITZ. This model is the complete mesoscale strategy proposed here.

• Overlapped mesh without ITZ (OM - no ITZ): In this model, the aggregate’s mesh is immersed in the mesh of the matrix, without adding degrees of freedom. The interface elements are generated allowing cracking propagation, but there are no ITZ properties.

• Single mesh (SM): In this model, the aggregate’s mesh is represented in the same mesh generation as the matrix, i.e., it uses a conforming mesh generator before creating the interface elements that allows crack propagation. This model is equivalent to the mesoscale concrete model presented by [⁷]. The only difference is that the model also considers the geometric nonlinearity of the problem.

3.1 Modeling of coarse aggregate

The coarse aggregates are artificially generated with random positions, dimensions and shapes. A specific computational code was developed for the generation and positioning of aggregates using the take-and-place method [⁷], [²¹], [²²], [⁶⁰ ]–[⁶³] allowing generating aggregates in the form of regular polygons with different numbers of sides. It is possible to use any granulometric distribution, whether experimental or theoretical, to define the amount of coarse aggregates in each size range of the studied cases. For complex packing strategies readers are invited to consult references [²⁶], [²⁷], [⁶⁴], [⁶⁵], for instance.

After generating all the geometries of the coarse aggregates, the finite element mesh of each aggregate is created and the mechanical properties are assigned. For the above defined OM and OM - no ITZ models it is noteworthy that, as aggregates overlap the mortar matrix, the value to be assigned to their Young modulus must be the difference between their real value and the mortar’s one, as done in [³¹].

3.2 Modeling of interface solid finite element

The interface elements are represented by finite elements of high aspect ratio and small thickness. Manzoli et al. [⁶⁶] showed that when the smallest height of these elements approaches zero, at the limit, their mechanical behavior is equivalent to that observed in the Continuous Strong Discontinuity Approach [⁶⁷], [⁶⁸], allowing the crack propagation modeling by using an isotropic damage model.

The interface elements are generated using an adapted version of the mesh fragmentation technique [⁵⁶]. This fragmentation technique has 3 steps depicted in Figure 3. Step 1 generates a regular mesh for the entire sample, Step 2 reduces the dimensions of all elements leaving a gap between them, which in Step 3 are filled by inserting pairs of interface elements. In reference [⁵⁶] triangular elements sizes are reduced keeping the element centroid fixed for all elements, thus all sides of elements are reduced even for boundary edges and corners. This imposes a small reduction in size of the analyzed solid. Here, for edges we adopt the segment center and for corners the element sides intersection as fixed references, respectively, eliminating the sizing reduction, see Figure 3.

Figure 3
Steps for creating high aspect ratio interface elements.

It is opportune to mention that some interface elements are transformed in ITZ elements by changing their physical properties, but the appropriate localization of ITZ inside the domain is necessary and described in item 3.4.

3.3 Damage model

The degradation occurs only in the interface elements (including ITZ elements) and considers a tensile damage model (orthogonal direction of the interface) [⁷], [⁵⁶]. Thus, other elements inside the mesoscale (usually called bulk elements) are elastic with mechanical properties of mortar or aggregates. In order to write the adopted damage model, it is necessary to calculate the Cauchy stress from the second Piola-Kirchhoff stress, see Equation 7, by the well-known formulae:

$\sigma =\frac{1}{J}\mathit{A}\cdot \mathit{S}\cdot {\mathit{A}}^t$ (19)

where $J$ is the determinant of the gradient of the change of configuration function.

The criteria for the existence of damage $\varphi$ is given by:

$\varphi ={\sigma }_{nn}-q\left(r\right)\le 0$ (20)

where ${\sigma }_{nn}$ is the Cauchy stress normal to the base of the interface element, $r$ is a stress-like internal variable and the function $q\left(r\right)$ defines the softening law. Dividing all terms of $\varphi$ by $(1-d)$, in which $d$ is called damage variable, the damage criterion for the effective stress is given by:

$\bar{\varphi }={\bar{\sigma }}_{nn}-r\le 0$ (21)

in which $r$ is defined by $r=q\left(r\right)/(1-d)$ and controls the size of the elastic domain in the space of the effective stress. Isolating $d$, the damage variable is given by:

$d=1-\frac{q\left(r\right)}{r}$ (22)

In the damage evolution, for the loading-unloading conditions, one must respect the relations given by:

$\bar{\varphi }\le 0$, $\dot{r}\ge 0$ and $\dot{r}\bar{\varphi }=0$ (23)

To ensure the well-known consistency condition one writes:

$\dot{r}\dot{\bar{\varphi }}=0$ if $\varphi =0$ (24)

Under these conditions, for a pseudo-time $t$ associated with the loading process, the variable $r$ is always given by the highest value between the current and its initial value $r_0$ given by the tensile strength $f_t$ of the material. So, it can be written as:

$r=\mathit{max}\left({\bar{\sigma }}_{nn}\left(s\right),r_0\right)\left.\right|s\in \left[0,t\right]$ (25)

To represent the softening, the same exponential law already used successfully by Rodrigues et al. [⁷] for mesoscale concrete is chosen and given by:

$q\left(r\right)=f_t\mathit{exp}\left(\frac{f_t^2}{G_{f}E}h\left(1-r/f_t\right)\right)$ (26)

where $h$ is the smallest height of the interface element and $G_f$ is the fracture energy for mode I.

The damage implementation uses the implicit-explicit integration scheme (IMPL-EX) developed by Oliver et al. [⁶⁹]. In this strategy the value of the damage variable used in each step depends on the linear extrapolation of the variable $r$ obtained in the previous step. Table 1 presents the algorithm implemented in this study.

Table 1
IMPL-EX scheme for the tension damage model at step n + 1.

INPUT: $E_{\left(n+1\right)},\mathrm{}{\overline{\sigma }}_{n{n}_{\left(n\right)}},\mathrm{}{r}_{\left(n\right)},\mathrm{}\Delta r_{\left(n\right)}$
(i) Compute the effective second Piola-Kirchhoff stress tensor
${\overline{S}}_{\left(n+1\right)}=\mathfrak{C}:E_{\left(n+1\right)}$
(ii) Compute the Cauchy effective stress tensor
${\overline{\sigma }}_{\left(n+1\right)}=\frac{1}{J}A\cdot {\overline{S}}_{\left(n+1\right)}\cdot A^t$
(iii) Compute the effective stress normal to the base of the element
${\overline{\sigma }}_{n{n}_{\left(n+1\right)}}=\overrightarrow{n^t}\cdot {\overline{\sigma }}_{\left(n+1\right)}\cdot \overrightarrow{n}$
(iv) Compute loading-unloading conditions
If ${\overline{\sigma }}_{n{n}_{\left(n+1\right)}}\le r_{\left(n\right)}\mathrm{}\mathrm{}⟶\mathrm{}\mathrm{}{r}_{\left(n+1\right)}=r_{\left(n\right)}$
If ${\overline{\sigma }}_{n{n}_{\left(n+1\right)}}>r_{\left(n\right)}\mathrm{}\mathrm{}⟶\mathrm{}\mathrm{}{r}_{\left(n+1\right)}={\overline{\sigma }}_{n{n}_{\left(n+1\right)}}$
(v) Compute the strain-like internal variable increment
$\Delta r_{\left(n+1\right)}=r_{\left(n+1\right)}-r_{\left(n\right)}$
(vi) Compute the explicit linear extrapolation of the strain-like internal variable
${\tilde{r}}_{\left(n+1\right)}=r_{\left(n\right)}+\frac{\mathrm{\Delta }{r}_{\left(n\right)}}{\mathrm{\Delta }{t}_{\left(n\right)}}\mathrm{\Delta }{t}_{\left(n+1\right)}$,
where $\mathrm{\Delta }{\mathrm{}t}_{\left(n+1\right)}=t_{\left(n+1\right)}-t_{\left(n\right)}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{\Delta }{t}_{\left(n\right)}=t_{\left(n\right)}-t_{\left(n-1\right)}$
(vii) Update the damage variable
${\tilde{d}}_{\left(n+1\right)}=1-\frac{q\left({\tilde{r}}_{\left(n+1\right)}\right)}{{\tilde{r}}_{\left(n+1\right)}}$
(viii) Compute the approximate nominal Cauchy effective stress tensor
${\tilde{S}}_{\left(n+1\right)}=\left\{\begin{array}{c}{\left(1-{\tilde{d}}_{\left(n+1\right)}\right){\overline{S}}_{\left(n+1\right)}if\overline{\sigma }}_{n{n}_{\left(n\right)}}>0\\ {\overline{S}}_{\left(n+1\right)}if{\overline{\sigma }}_{n{n}_{\left(n\right)}}\le 0\end{array}\right.$
(ix) Compute the approximate constitutive tensor
${\tilde{\mathfrak{C}}}_{\left(n+1\right)}=\left\{\begin{array}{c}\left(1-{\tilde{d}}_{\left(n+1\right)}\right)\mathfrak{C}\mathfrak{}\mathfrak{}\mathfrak{}\mathfrak{}if{\overline{\sigma }}_{n{n}_{\left(n\right)}}>0\\ \mathfrak{C}\mathfrak{}\mathfrak{}\mathfrak{}\mathfrak{}\mathfrak{}\mathfrak{}\mathfrak{}\mathfrak{}\mathfrak{}\mathfrak{}\mathfrak{}\mathfrak{}if{\overline{\sigma }}_{n{n}_{\left(n\right)}}\le 0\end{array}\right.$
OUTPUT: ${\tilde{\mathfrak{C}}}_{\left(n+1\right)},{\mathrm{}\tilde{\sigma }}_{\left(n+1\right)},{\mathrm{}\overline{\sigma }}_{n{n}_{\left(n+1\right)}},{\mathrm{}r}_{\left(n+1\right)},\mathrm{}\Delta r_{\left(n+1\right)}$

Once the damage variable has been determined, the solution of the linear system of equations that allows finding the equilibrium of the structure is started. The values of the second Piola-Kirchhoff stress tensor and the elastic constitutive tensor penalized by the damage variable are applied in Equations 10 and 11. Once this is done, the linear system of equations in Equation 4 can be solved.

3.4 Determination of matrix-aggregate interfaces (ITZ)

This subsection presents the proposed strategy to represent the ITZ around the coarse aggregates in the mesoscale modeling of concrete. The strategy is applied in cases where the above defined overlapping mesh (OM) is adopted. The proposed strategy transforms interface elements near the edges of aggregates into ITZ elements by changing their mechanical properties [¹³].

To describe the proposed strategy, we consider Figure 4. Let ${\Omega }_1$ be the matrix geometry and ${\Omega }_2$ be the aggregate geometry (immersed in ${\Omega }_1$). Around ${\Omega }_2$, pseudo finite element pairs are created as shown (in green) in Figure 4a. The fragmented mesh generated for ${\Omega }_1$ is illustrated in Figure 4b. The matrix interface elements that have the centroid inside those pseudo finite elements receive ITZ mechanical properties turning into ITZ elements, highlighted in Figure 4c. To identify if the centroid belongs to a pseudo element the same algorithm used in section 2.2 (Equation 12) is employed.

Figure 4
Scheme for defining the interface elements that form the ITZ: (a) Matrix, aggregate and pseudo finite element geometries. (b) Subdomains mesh. (c) Interface elements with ITZ property.

The pseudo finite elements are created from the aggregate faces using its normal vector and the adopted ITZ region thickness (H). This height is an important variable for the accuracy of the proposed strategy. A study regarding this thickness dimension is carried out in the next section.

4 NUMERICAL EXAMPLES

In this section, 4 examples of numerical modeling are presented. The first example studies the influence of matrix’s mesh refinement, ITZ presence and its thickness. In the second example, using the results obtained in the previous one, different samples are studied to identify the effect of aggregates’ volume fractions. The third example validates the proposed non-matching meshes coupling strategy to link different scales. The last example studies a concrete beam modeled by the proposed concurrent multiscale and compares results with experimental data.

4.1 Mesoscale uniaxial tension tests – sensitivity analysis

This example tests the geometries and boundary conditions shown in Figure 5. Displacement control $\delta$ is adopted as depicted in the Figure 5 and the sample thickness is 50 mm. Regular pentagon-shaped aggregates are immersed in the central region of samples. The size of the aggregates is defined as a function of the diameter (d) of the circle that inscribes the pentagons. We adopted $d=60mm$ and $d=24mm$ for cases with 1 and 5 aggregates, respectively. The regions for which displacements are prescribed have the same elastic properties of the central region, but only the last one presents interface elements, i.e., allows degradation.

Figure 5
Geometry and boundary conditions: sample with 1 and with 5 aggregates.

The same set of analyzes is repeated for samples with 1 and 5 aggregates. This is done using the models described in section 3, i.e., Overlapped mesh (OM), Overlapped mesh without ITZ (OM - no ITZ) and Single mesh (SM). All simulations of this example were performed by dividing the prescribed displacement into 400 steps. For larger amount of steps, solutions are stable.

For the SM model, the matrix elements have Young modulus and Poisson ratio of $20GPa$ and $0.2$, respectively. The interface elements have the same Young modulus but null Poisson ratio. The aggregates have a Young modulus of $40GPa$ and a Poisson ratio of $0.2$. The damage properties of interface elements are $f_t=2MPa$ and $G_f=0.04N/mm$. ITZ elements have the same elastic properties of interface elements, but damage properties of $f_t=1.0MPa$ and $G_f=0.02N/mm$.

For the OM and OM - no ITZ models the adopted aggregates’ Young modulus is $20GPa$ as, due to material’s superposition, it is necessary to discount the matrix’s Young modulus. Finally, for the OM - no ITZ model, the adopted interface elements’ damage properties are $f_t=1.5MPa$ and $G_f=0.03N/mm$.

All meshes are defined as a function of the aggregate diameter ($d$), thus the side of a finite element $l_{FE}$ (matrix or aggregate) has an approximate dimension of a fraction of $d$. A mesh with $l_{FE}\cong d/20$ for the SM strategy is used as reference for all analyzed cases.

In Figure 6 we use the model OM – no ITZ to verify the influence of the mesh density on the global behavior of the analyzed problem. We keep a $d/20$ mesh to represent the aggregate and meshes $d/12$ and $d/20$ to model the matrix. As expected, we conclude that the OM – no ITZ model does not have good accuracy along the post peak branch, but it indicates that a $d/20$ mesh for the matrix captures well the peak value. This mesh is used to verify the behavior of the OM model when varying the height (H) of the pseudo-elements used to define the ITZ region.

Figure 6
Force-Displacement curves: Comparison of SM model and OM - no ITZ model responses for sample with (a) 1 and (b) 5 aggregates.

The achieved results for the OM model using different H values are shown in Figure 7. We remember that, for the adopted matrix mesh $l_{FE}\cong d/20$. As one may observe, only the value $H=0.5l_{FE}$ does not accurately represent the force peak value and for values larger then $l_{FE}$ results do not present significant variation. As expected, the post-peak behavior (softening) shows less agreement with the reference value. This is because the ITZ in the SM model is a straight line, while in the OM model it has a zig-zag geometry. It should be noted that both results are numerical and that the real representation of the ITZ geometry depends on experimental observations that are beyond the scope of this study. This comment becomes evident when comparing the last two cases in Figure 7.

Figure 7
Force-Displacement curves: Study of the height of pseudo finite elements for samples with (a) 1 and (b) 5 aggregates.

As in the OM – no ITZ model there is no reduction in fracture energy or tension strength, the nucleation of fractures occurs in regions different from those expected experimentally [¹³], [⁷⁰] see Figures 8 and 9. Despite the similar aspects of fracture paths (Figure 8) for the OM model with $H=0.5l_{FE}$ and $H=1.0l_{FE}$, the number of fracture paths allowed by the first case is smaller, justifying the stiff behavior of the model with $H=0.5l_{FE}$ (Figure 7) and also the difference between the fracture paths of Figure 8. As the identification of the ITZ elements is done by their centroid position in relation to the pseudo-element, the value $H=0.5l_{FE}$ is under the lower limit that guarantees a good representation of the continuity of the fracture paths, see the big difference of behaviors in Figure 9. Due to the small difference between behaviors of the OM formulation (Figure 7) for H values above $1.0l_{FE}$, it is understood that when $H\ge 1.0l_{FE}$ the ITZ provides a complete crack path. Results also confirm that nucleation occurs in regions very close to the aggregate interface.

Figure 8
Sample with 1 aggregate: cracks for different used models.

Figure 9
Sample with 5 aggregates: cracks for different used models.

Keeping $H=1.0l_{FE}$, in Figure 10 we show the influence of the matrix mesh on the behavior of the OM formulation response. The following values are adopted for the $l_{FE}$ in the matrix mesh: $d/30$, $d/20$; $d/12$ and $d/6$. A constant mesh with $d/12$ is adopted to model aggregates. For the coarsest mesh ($d/6$) we observed that the ITZ elements are very large and the start of nucleation is anticipated when the number of aggregates is 5. In the case that has only 1 aggregate, this behavior is masked because of the simple stress distribution. When finer discretizations are considered, the peak load and the post-peak behavior are a little more rigid than the response of the SM model, confirming the results and discussions regarding Figure 7. There is also a tendency of smaller sensitivity of results between $d/20$ and $d/30$ meshes, since convergence is not expected in a nonlinear problem with softening.

Figure 10
Force-Displacement curves: Study of matrix mesh refinement for sample with (a) 1 and (b) 5 aggregates.

The aggregate immersion strategy in the OM model does not significantly influence the nonlinear behavior of the analyzed problems, as the aggregate discretization is not used to determine the size of the interface elements or the ITZ region. This can be verified by the small variability of the results presented in Figure 11, in which we keep the matrix mesh in $d/20$ and vary the aggregate’s mesh using $d/20$, $d/12$, $d/6$ and $d/3$.

Figure 11
Force-Displacement curves: Study of aggregate mesh refinement for sample with (a) 1 and (b) 5 aggregates.

4.2 Mesoscale uniaxial tension tests for different volumes of aggregates

In this example, for OM and SM models, using the same specimens of the previous example (same boundary conditions), we analyze the effect of the aggregates’ volume fractions. Three cases are studied with the following aggregate volume fractions: 30%, 35% and 40%. These volume fractions are chosen because according to Ma et al. [⁶⁴] the volume of coarse aggregates in conventional concrete is in the range between 30% and 40% of the total volume.

The shapes (polygons with 5, 6, 7 and 8 sides) and positions of the aggregates are chosen randomly, while their sizes are determined using the Fuller's theoretical granulometric curve [⁶¹] with diameters between 5 mm and 10 mm.

Meshes with approximate $l_{FE}$ dimension $d/12$ are used for the matrix elements of the OM model and for all elements of the SM model, for which $d=7.5mm$ (average aggregates’ dimension). Based on the previous example results, achieved for the OM model, a discretization corresponding to $d/3$ is adopted for the aggregates and the height of the pseudo finite element to define the ITZ is adopted as $H=1.0l_{FE}$. All simulations of this example were performed by dividing the prescribed displacement into 800 steps. For larger amount of steps, solutions are stable.

Observing the force - displacement results of Figure 12, we conclude that the alternative formulation proposed here (OM) has similar results to the SM formulation for both peak values and post-peak behavior. The stiffness differences are very small, leading to the conclusion that the models are mechanically equivalent in situations closer to real composites. In addition, as the matrix’s node numbering does not change for different volume fractions, the OM formulation has a simpler mesh generation strategy and boundary conditions application.

Figure 12
Force-Displacement curves: samples analyzed using SM model and OM model with different aggregate volumes.

As can be seen in Figure 13, the fracture paths have many similarities, that is, the concentration of stresses in the SM and OM models (that lead to the nucleation and propagation of fractures) are very close, revealing the equivalence in quality of the formulations.

Figure 13
Cracks in samples analyzed using SM model and OM model with different aggregate volumes.

In what follows, a study is carried out to identify the influence of random aggregate distribution on response curves and crack pattern. Three different distributions of aggregates were generated for the volume fraction of 35% and analyzed by the OM and SM models. The mesh refinements were the same described previously and in the case of the OM model the matrix mesh does not change. The results presented in Figures 12 and 13 with volume fraction of 35% are also included in this new set of analyses and called SM1 and OM1. The new samples were named SM2, SM3 SM4, OM2, OM3 and OM4.

As expected, Figure 14 confirms that the position of the aggregates influences the response curves for both models. Figure 15 shows, as expected, that the aggregates position influences the starting point and paths of cracks. Furthermore, it was evident that the cracks obtained for the same distribution of aggregates in the OM model and in the SM model were close in the four studied cases.

Figure 14
Force-Displacement curves: samples analyzed using OM (a) and SM (b) models with 35% of aggregate volume fraction and different distributions.

Figure 15
Samples with different aggregate distributions for OM and SM models: Cracks patterns with volume fraction of 35%.

It is important to mention that in this example the processing time was measured. The OM model takes 39% more time than the SM model. However, the intention of the developing the OM model is the possibility to analyze multiple random aggregate distributions, with different volume fractions, which is very difficult using the SM model, doing to the pre-processing step.

4.3 Coupling non-matching meshes for different boundary conditions

The third example validates the proposed non-matching meshes coupling strategy to link different scales. It is considered an elastic specimen (no degradation) subjected to traction, compression, and shear stresses. The problem is modeled using the SM formulation (unique mesh) and the OM formulation adapted to consider different meshes’ scales in adjacent regions.

Figure 16 shows the geometry and boundary conditions of the problem, whose thickness is 50 mm. The prescribed displacements are $\delta =0.20mm$ and $\delta =-0.20mm$ vertically and $\delta =0.20mm$ horizontally for tension, compression and shear-like cases, respectively.

Figure 16
Geometry and boundary conditions: (a) tension, (b) compression and (c) shear tests.

The SM discretization has 2500 elements and 1326 nodes, the non-matching (OM) discretization has two regions, one with the same mesh density as the SM case and the other with a coarser mesh. The coupling strategy, see section 2, results in 1298 elements and 714 nodes. In all cases, the Young modulus is 20000 MPa and the Poisson ratio is 0.20.

In Figure 17, the normal stress and vertical displacement fields of the tensioned case are shown. As one can see, vertical displacements are in very good agreement and the relative difference in stress field has a maximum of 1.5%.

Figure 17
Tension test.

In Figure 18, the vertical displacement and normal stress fields for the compressed case are presented. As expected, the relative difference between the SM and OM representations are the same as the tension case, i.e., 1.5%.

Figure 18
Compression test.

In Figure 19, the horizontal displacement and normal stress fields in the vertical direction for the sample subjected to shear-like displacements are shown. It is important to note that in the contact region between non-conforming meshes, nodal results are not averaged, considering that the nodes are not coincident. In this sense, the stresses in the well-discretized region of the OM model (mesoscale) do not present significant differences in relation to the stresses of the SM model, resulting in the possibility of performing very precise nonlinear analyses. In the region of coarse mesh (OM) it is not expected an accurate stress distribution since physical non-linear analyses are not to be performed there. Finally, the displacement distribution is practically coincident for the two models.

Figure 19
Shear test.

From the previous discussion, using the proposed non-matching approach is a good alternative to differentiate mesoscale and macro scale representations, reducing the computational cost and keeping the accuracy at the region in which degradation may occur.

4.4 Multiscale three-point bending test

The last example studies a concrete beam modeled by the proposed concurrent multiscale model and compares results with experimental data. From the above results and the existence of experimental data, the SM strategy is not used to describe aggregates.

A concrete beam subjected to a three-point bending test is simulated using the proposed OM model to represent mesoscale and results are compared with experimental data given by Grégoire et al. [⁷¹]. The coupling between macroscale and mesoscale is performed using (i) a mesh refinement in transitions’ region, resulting in a conforming mesh and (ii) the proposed non-matching meshes coupling strategy of section 2. In both analyses, the mesoscale is modeled using aggregates that do not add degrees of freedom and with ITZ generated using the OM strategy proposed here.

The beam geometry and boundary conditions are shown in Figure 20. The beam has a notch on the lower edge in the central region with a length equal to one fifth of the beam height and a width of 2 mm. The notch induces stress concentration and the crack propagation, thus only this region needs to be represented in mesoscale. The thickness of the beam is 50 mm and displacement control is applied at the midpoint of the upper edge, where the reaction force F is numerically measured.

Figure 20
Geometry and boundary conditions.

For both analyzed cases, the coarse aggregates are generated in the form of regular polygons with 5, 6, 7 and 8 sides and 40% of volume fraction. They are represented in the size range limited by a minimum diameter of 5 mm and a maximum diameter of 10 mm. The size of the coarse aggregates follows the granulometric distribution given by the average experimental curve of [⁷¹].

The concrete’s Young modulus (macro scale) is the same presented by [⁷¹]. The properties of mortar, coarse aggregates and interfaces were adopted with the same values estimated by [⁷], [²¹]. The values adopted for the mechanical properties are presented in Table 2. We recall that the numerical value of the aggregate Young modulus is 19.8 GPa, due to the superposition of mortar and aggregate meshes.

Table 2
Mechanical properties of each material.

Material	Young modulus (GPa)	$\nu$	$G_f$ (N/mm)	$F_t$ (MPa)
Concrete	37	0.2	-	-
Coarse aggregate	50	0.2	-	-
Mortar	30.2	0.2	-	-
Coupling element	30.2	0.2	-	-
ITZ	30.2	0	0.05	2.6
Int. matrix-matrix	30.2	0	0.1	5.2

The mean aggregate diameter ($d=7.5mm$) is used to define the mesoscale mesh refinement. The mortar mesh is formed by elements with sides of the order of $d/20$. The height of the pseudo-elements used to define the ITZ is $H=1.25l_{FE}$. The aggregate mesh is formed by elements with sides of the order of $d/3$. The macroscale mesh is formed by elements with sides of the order of 4 mm. All simulations of this example were performed by dividing the prescribed displacement into 800 steps. For larger amount of steps, solutions are stable.

The curves shown in Figure 21 relate the reaction force F with the crack mouth opening displacement (CMOD). The achieved numerical responses are very close to the average experimental curve of [⁷⁰]. The main difference between results (peak) is due to the fact that we adopted the properties chosen by [⁷] and [²¹] which do not take into account the randomness of the positioning and orientation of the aggregates, whose detailed assessment is beyond the objective of the present study.

Figure 21
Comparison of the numerical results with those obtained experimentally by [⁷¹]: reaction force-CMOD curves.

Comparing the results of cases (i) and (ii), in Figure 21, one observes that they are practically identical, confirming that the meshing coupling strategy works properly for concurrent multiscale analysis. One can also observe in Figure22 that the crack patterns of both cases are very similar, thus the stress field provided by the non-matching strategy is considered accurate.

Figure 22
Numerical crack patterns obtained: (a) beam without and (b) beam with use of coupling element.

From the previous observations, we conclude that the use of the OM strategy (ITZ and interface elements) in conjunction with the strategy of coupling non-matching meshes can be used in general problems of multiscale concrete modeling.

5 CONCLUSIONS

In this work, a computational framework to numerically model concurrent multiscale analysis of particulate composite is implemented. Two main contributions are detached: The first introduces an alternative representation of the mesoscale interfacial transition zone (ITZ) of the concrete, improving a previously developed strategy that allows modeling aggregates without degrees of freedom. The main advantage of the proposed ITZ representation, together with the model of representation of aggregates embedded in the matrix, is the possibility of representing discrete cracks that pass on the face of the aggregates, a characteristic observed experimentally in conventional concrete. Another advantage of mesoscale modeling that comes with this improvement is the possibility of performing parametric studies of coarse aggregates while preserving the matrix mesh and boundary conditions, which simplifies the preprocessing stage. The second contribution consists in using immersed particle-like elements without degrees of freedom as coupling elements of non-matching meshes. The proposed coupling elements do not add degrees of freedom to the problem and does not use penalization technics.

From the performed analyzes, including comparisons with other numerical techniques and with experimental results, we conclude that the proposed computational framework can be used in general multiscale analyzes of concrete structures. Future developments should consider reliability analysis involving the randomness of parameters: geometry, relative quantity, and physical properties of the constituent media.

From the good results using the proposed two-dimensional simulations (plane stress state) it is also proposed as future developments its extension to the three-dimensional case.

ACKNOWLEDGEMENTS

The authors acknowledge the financial support of Sao Paulo Research Foundation (FAPESP), Brazil - Grant #2020/05393 – 4 and the Coordination for the Improvement of Higher Education Personnel (CAPES) finance code 001.

• Financial support: São Paulo Research Foundation (FAPESP), Brazil - Grant #2020/05393 – 4 and the Coordination for the Improvement of Higher Education Personnel (CAPES) finance code 001.

• Data Availability:

  none.
• How to cite: W. H. Vieira, H. B. Coda, and R. R. Paccola, “A direct FEM approach to model mesoscale concrete and connect non-matching meshes in multiscale analysis,” Rev. IBRACON Estrut. Mater., vol. 17, no. 1, e17108, 2024, https://doi.org/10.1590/S1983-41952024000100008

REFERENCES

• ¹
  I. M. Daniel and O. Lshai, Engineering Mechanics of Composite Materials, vol. 2. London, UK: Oxford Univ. Press, 2006.
• ²
  J. F. Unger and S. Eckardt, “Multiscale modeling of concrete,” Arch. Comput. Methods Eng., vol. 18, pp. 341–393, 2011.
• ³
  P. Chen, J. Liu, X. Cui, and S. Si, “Mesoscale analysis of concrete under axial compression,” Constr. Build. Mater., vol. 337, pp. 127580, 2022.
• ⁴
  B. D. Barnes, S. Diamond, and W. L. Dolch, “The contact zone between portland cement paste and glass “aggregate” surfaces,” Cement Concr. Res., vol. 8, pp. 233–243, 1978.
• ⁵
  R. Zimbelmann, “A contribution to the problem of cement-aggregate bond,” Cement Concr. Res., vol. 15, pp. 801–808, 1985.
• ⁶
  C. M. López, I. Carol, and A. Aguado, “Meso-structural study of concrete fracture using interface elements. I: numerical model and tensile behavior,” Mater. Struct., vol. 41, pp. 583–599, 2008.
• ⁷
  E. A. Rodrigues, O. L. Manzoli, L. A. G. Bitencourt Jr., and T. N. Bittencourt, 2D mesoscale model for concrete based on the use of interface element with a high aspect ratio,” Int. J. Solids Struct., vol. 94-95, pp. 112–124, 2016.
• ⁸
  Y. Xia, W. Wu, Y. Yang, and X. Fu, “Mesoscopic study of concrete with random aggregate model using phase field method,” Constr. Build. Mater., vol. 310, pp. 125199, 2021.
• ⁹
  X. Wang, Z. Yang, and A. P. Jivkov, “Monte Carlo simulations of mesoscale fracture of concrete with random aggregates and pores: a size effect study,” Constr. Build. Mater., vol. 80, pp. 262–272, 2015.
• ¹⁰
  H. Zhou, H. Zhou, X. Wang, W. Cao, T. Song, and Q. Qiao, “Static size effect of recycled coarse aggregate concrete: experimental study, meso-scale simulation, and theoretical analysis,” Structures, vol. 34, pp. 2996–3012, 2021.
• ¹¹
  D. Tal and J. Fish, “Stochastic multiscale modeling and simulation framework for concrete,” Cement Concr. Compos., vol. 90, pp. 61–81, 2018.
• ¹²
  S. Naderi, W. Tu, and M. Zhang, “Meso-scale modelling of compressive fracture in concrete with irregularly shaped aggregates,” Cement Concr. Res., vol. 140, pp. 106317, 2021.
• ¹³
  E. Schlangen and J. G. M. Van Mier, “Simple lattice model for numerical simulation of fracture of concrete materials and structures,” Mater. Struct., vol. 25, pp. 534–542, 1992.
• ¹⁴
  H. K. Man and J. G. Van Mier, “Damage distribution and size effect in numerical concrete from lattice analyses,” Cement Concr. Compos., vol. 33, pp. 867–880, 2011.
• ¹⁵
  B. B. Aydin, K. Tuncay, and B. Binici, “Simulation of reinforced concrete member response using lattice model,” J. Struct. Eng., vol. 145, pp. 1–17, 2019.
• ¹⁶
  Z. P. Bažant, M. R. Tabbara, M. T. Kazemi, and G. Pijaudier-Cabot, “Random particle model for fracture of aggregate or fiber composites,” J. Eng. Mech., vol. 116, pp. 1686–1705, 1990.
• ¹⁷
  M. Yip, Z. Li, B.-S. Liao, and J. Bolander, “Irregular lattice models of fracture of multiphase particulate,” Int. J. Fract., vol. 140, pp. 113–124, 2006.
• ¹⁸
  G. Cusatis, A. Mencarelli, D. Pelessone, and J. Baylot, “Lattice Discrete Particle Model (LDPM) for failure behavior of concrete. II: calibration and validation,” Cement Concr. Compos., vol. 33, pp. 891–905, 2011.
• ¹⁹
  R. Zhou and Y. Lu, “A mesoscale interface approach to modelling fractures in concrete for material investigation,” Constr. Build. Mater., vol. 165, pp. 608–620, 2018.
• ²⁰
  E. A. Rodrigues, O. L. Manzoli, L. A. Bitencourt, T. N. Bittencourt, and M. Sánchez, “An adaptive concurrent multiscale model for concrete based on coupling finite elements,” Comput. Methods Appl. Mech. Eng., vol. 328, pp. 26–46, 2018.
• ²¹
  E. A. Rodrigues, O. L. Manzoli, and L. A. Bitencourt Jr., “3D concurrent multiscale model for crack propagation in concrete,” Comput. Methods Appl. Mech. Eng., vol. 361, pp. 1–33, 2020.
• ²²
  E. A. Rodrigues, M. Gimenes, L. A. Bitencourt, and O. L. Manzoli, “A concurrent multiscale approach for modeling recycled aggregate concrete,” Constr. Build. Mater., vol. 267, pp. 121040, 2021.
• ²³
  J. Oliver, M. Caicedo, E. Roubin, A. E. Huespe, and J. A. Hernández, “Continuum approach to computational multiscale modeling of propagating fracture,” Comput. Methods Appl. Mech. Eng., vol. 294, pp. 384–427, 2015.
• ²⁴
  E. Roubin, A. Vallade, N. Benkemoun, and J. B. Colliat, “Multi-scale failure of heterogeneous materials: a double kinematics enhancement for embedded finite element method,” Int. J. Solids Struct., vol. 52, pp. 180–196, 2015.
• ²⁵
  X. Du, L. Jin, and G. Ma, “Numerical modeling tensile failure behavior of concrete at mesoscale using extended finite element method,” Int. J. Damage Mech., vol. 23, pp. 872–898, 2014.
• ²⁶
  Y. Huang, Z. Yang, W. Ren, G. Liu, and C. Zhang, “3D meso-scale fracture modelling and validation of concrete based on in-situ X-ray computed tomography images using damage plasticity model,” Int. J. Solids Struct., vol. 67-68, pp. 340–352, 2015.
• ²⁷
  Ł. Skar and J. Tejchman, “Experimental investigations of fracture process in concrete by means of X-ray micro-computed tomography,” Strain, vol. 52, pp. 26–45, 2016.
• ²⁸
  Z. Yang et al., “In-situ X-ray computed tomography characterization of 3D fracture evolution and image-based numerical homogenization of concrete,” Cement Concr. Compos., vol. 75, pp. 74–83, 2017.
• ²⁹
  Y. Huang, D. Yan, Z. Yang, and G. Liu, “2D and 3D homogenization and fracture analysis of concrete based on in-situ X-ray computed tomography images and Monte Carlo simulations,” Eng. Fract. Mech., vol. 163, pp. 37–54, 2016.
• ³⁰
  H. Zhang, Y. J. Huang, M. Lin, and Z. J. Yang, “Effects of fibre orientation on tensile properties of ultra high performance fibre reinforced concrete based on meso-scale Monte Carlo simulations,” Compos. Struct., vol. 287, pp. 115331, 2022.
• ³¹
  R. R. Paccola and H. B. Coda, “A direct FEM approach for particulate reinforced elastic solids,” Compos. Struct., vol. 141, pp. 282–291, 2016.
• ³²
  É. S. Ramos, R. Carrazedo, and R. R. Paccola, “Modeling particles elements in damaged reinforced concrete structures,” Lat. Am. J. Solids Struct., vol. 18, pp. 1–24, 2021.
• ³³
  J. Bonet, R. D. Wood, J. Mahaney, and P. Heywood, “Finite element analysis of air supported membrane structures,” Comput. Methods Appl. Mech. Eng., vol. 190, pp. 579–595, 2000.
• ³⁴
  H. Coda and M. Greco, “A simple FEM formulation for large deflection 2D frame analysis based on position description,” Comput. Methods Appl. Mech. Eng., vol. 193, pp. 3541–3557, 2004.
• ³⁵
  L. Vanalli, R. R. Paccola, and H. B. Coda, “A simple way to introduce fibers into FEM models,” Commun. Numer. Methods Eng., vol. 24, pp. 585–603, 2008.
• ³⁶
  F. K. F. Radtke, A. Simone, and L. J. Sluys, “A computational model for failure analysis of fibre reinforced concrete with discrete treatment of fibres,” Eng. Fract. Mech., vol. 77, pp. 597–620, 2010.
• ³⁷
  M. S. Sampaio, R. R. Paccola, and H. B. Coda, “Fully adherent fiber-matrix FEM formulation for geometrically nonlinear 2D solid analysis,” Finite Elem. Anal. Des., vol. 66, pp. 12–25, 2013.
• ³⁸
  R. R. Paccola, D. N. Piedade, and H. B. Coda, “Geometrical non-linear analysis of fiber reinforced elastic solids considering debounding,” Compos. Struct., vol. 133, pp. 343–357, 2015.
• ³⁹
  O. Lloberas-Valls, D. J. Rixen, A. Simone, and L. J. Sluys, “Multiscale domain decomposition analysis of quasi-brittle heterogeneous materials,” Int. J. Numer. Methods Eng., vol. 89, pp. 1337–1366, 2012.
• ⁴⁰
  O. Lloberas-Valls, D. Rixen, A. Simone, and L. Sluys, “On micro-to-macro connections in domain decomposition multiscale methods,” Comput. Methods Appl. Mech. Eng., vol. 225-228, pp. 177–196, 2012.
• ⁴¹
  H. Su, J. Hu, and H. Li, “Multi-scale performance simulation and effect analysis for hydraulic concrete submitted to leaching and frost,” Eng. Comput., vol. 34, pp. 821–842, 2018.
• ⁴²
  A. Sellitto, R. Borrelli, F. Caputo, A. Riccio, and F. Scaramuzzino, “Methodological Approaches for kinematic coupling of non- matching finite element meshes,” Procedia Eng., vol. 10, pp. 421–426, 2011.
• ⁴³
  A. Pantano and R. C. Averill, “A penalty-based finite element interface technology,” Comput. Struc., vol. 80, pp. 1725–1748, 2002.
• ⁴⁴
  L. A. Bitencourt, O. L. Manzoli, P. G. Prazeres, E. A. Rodrigues, and T. N. Bittencourt, “A coupling technique for non-matching finite element meshes,” Comput. Methods Appl. Mech. Eng., vol. 290, pp. 19–44, 2015.
• ⁴⁵
  B. I. Wohlmuth, “A mortar finite element method using dual spaces for the Lagrange multiplier,” SIAM J. Numer. Anal., vol. 38, pp. 989–1012, 2001.
• ⁴⁶
  B. P. Lamichhane and B. I. Wohlmuth, “Mortar finite elements for interface problems,” Comput., vol. 72, pp. 333–348, 2004.
• ⁴⁷
  H. Fang, D. Zhang, Q. Fang, L. Cao, and M. Wen, “An efficient patch-to-patch method for coupling independent finite element subdomains with intersecting interfaces,” Comput. Methods Appl. Mech. Eng., vol. 388, pp. 114209, 2022.
• ⁴⁸
  H. B. Dhia, “Multiscale mechanical problems: the Arlequin method,” Comp. Rendus Académ. Sci. Ser. IIB Mech.-Phys.-Astron., vol. 326, pp. 899–904, 1998.
• ⁴⁹
  H. B. Dhia and G. Rateau, “The Arlequin method as a flexible engineering design tool,” Int. J. Numer. Methods Eng., vol. 62, pp. 1442–1462, 2005.
• ⁵⁰
  J. Nitsche, “Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind,” Abh. Math. Semin. Univ. Hambg., vol. 36, pp. 9–15, 1971.
• ⁵¹
  A. Apostolatos, R. Schmidt, R. Wüchner, and K.-U. Bletzinger, “A Nitsche-type formulation and comparison of the most common domain decomposition methods in isogeometric analysis,” Int. J. Numer. Methods Eng., vol. 97, pp. 473–504, 2014.
• ⁵²
  V. P. Nguyen, P. Kerfriden, M. Brino, S. P. Bordas, and E. Bonisoli, “Nitsche’s method for two and three dimensional NURBS patch coupling,” Comput. Mech., vol. 53, pp. 1163–1182, 2014.
• ⁵³
  H. Kim, “Interface element method: treatment of non-matching nodes at the ends of interfaces between partitioned domains,” Comput. Methods Appl. Mech. Eng., vol. 192, pp. 1841–1858, 2003.
• ⁵⁴
  J. Zhang and C. Song, “A polytree based coupling method for non-matching meshes in 3D,” Comput. Methods Appl. Mech. Eng., vol. 349, pp. 743–773, 2019.
• ⁵⁵
  T. Rabczuk and T. Belytschko, “Application of particle methods to static fracture of reinforced,” Int. J. Fract., vol. 137, pp. 19–49, 2006.
• ⁵⁶
  O. L. Manzoli, M. A. Maedo, L. A. Bitencourt, and E. A. Rodrigues, “On the use of finite elements with a high aspect ratio for modeling cracks in quasi-brittle materials,” Eng. Fract. Mech., vol. 153, pp. 151–170, 2016.
• ⁵⁷
  C. Lanczos, The Variational Principles of Mechanics New York, NY, USA: Dover, 1986.
• ⁵⁸
  H. B. Coda, Introdução ao Método dos Elementos Finitos Posicionais: Sólidos e Estruturas - Não Linearidade Geométrica e Dinâmica São Carlos, Brazil: EESC-USP, 2018.
• ⁵⁹
  O. L. Manzoli, M. Tosati, E. A. Rodrigues, and L. A. G. Bitencourt Jr., “Computational modeling of 2D frictional contact problems based on the use of coupling finite elements and combined contact/friction damage constitutive model,” Finite Elem. Anal. Des., vol. 199, pp. 103658, 2022.
• ⁶⁰
  Z. M. Wang, A. K. Kwan, and H. C. Chan, “Mesoscopic study of concrete I: generation of random aggregate structure and finite element mesh,” Comput. Struc., vol. 70, pp. 533–544, 1999.
• ⁶¹
  P. Wriggers and S. Moftah, “Mesoscale models for concrete: homogenisation and damage behaviour,” Finite Elem. Anal. Des., vol. 42, pp. 623–636, 2006.
• ⁶²
  Z. Qian, E. J. Garboczi, G. Ye, and E. Schlangen, “Anm: a geometrical model for the composite structure of mortar and concrete using real-shape particles,” Mater. Struct., vol. 49, pp. 149–158, 2016.
• ⁶³
  S. C. Seetharam et al., “A mesoscale framework for analysis of corrosion induced damage of concrete,” Constr. Build. Mater., vol. 216, pp. 347–361, 2019.
• ⁶⁴
  H. Ma, W. Xu, and Y. Li, “Random aggregate model for mesoscopic structures and mechanical analysis of fully-graded concrete,” Comput. Struc., vol. 177, pp. 103–113, 2016.
• ⁶⁵
  Z. Zhang, X. Song, Y. Liu, D. Wu, and C. Song, “Three-dimensional mesoscale modelling of concrete composites by using random walking algorithm,” Compos. Sci. Technol., vol. 149, pp. 235–245, 2017.
• ⁶⁶
  O. Manzoli, A. Gamino, E. Rodrigues, and G. Claro, “Modeling of interfaces in two-dimensional problems using solid finite elements with high aspect ratio,” Comput. Struc., vol. 94-95, pp. 70–82, 2012.
• ⁶⁷
  J. Oliver, M. Cervera, and O. Manzoli, “Strong discontinuities and continuum plasticity models: the strong discontinuity approach,” Int. J. Plast., vol. 15, pp. 319–351, 1999.
• ⁶⁸
  J. Oliver, “On the discrete constitutive models induced by strong discontinuity kinematics and continuum constitutive equations,” Int. J. Solids Struct., vol. 37, pp. 7207–7229, 2000.
• ⁶⁹
  J. Oliver, A. E. Huespe, S. Blanco, and D. L. Linero, “Stability and robustness issues in numerical modeling of material failure with the strong discontinuity approach,” Comput. Methods Appl. Mech. Eng., vol. 195, pp. 7093–7114, 2006.
• ⁷⁰
  J. Mazars, “A description of micro - and macroscale damage of concrete structures,” Eng. Fract. Mech., vol. 25, pp. 729–737, 1986.
• ⁷¹
  D. Grégoire, L. B. Rojas-Solano, and G. Pijaudier-Cabot, “Failure and size effect for notched and unnotched concrete beams,” Int. J. Numer. Anal. Methods Geomech., vol. 37, pp. 1434–1452, 2013.

Edited by

Editors: Osvaldo Manzoli, Guilherme Aris Parsekian.

Data availability

none.

Publication Dates

• Publication in this collection
  19 June 2023
• Date of issue
  2024

History

• Received
  05 Dec 2022
• Accepted
  11 Apr 2023