1. Introduction

Several engineering problems can be described using partial differential equations relating field variables inside a particular domain. To obtain reasonable solutions it is used some numerical method since such problems have complex geometry and boundary conditions. In this context, it was developed the Finite Element Method (FEM), which is an efficient numerical tool to solve boundary value problems.

Although the FEM is a quite consolidated numerical technique in the study of several structural engineering problems, it presents some limitations related specially to the description of phenomena such as crack and damage propagation and large deformations. The nature of these phenomena leads to the modiﬁcation of the mesh in a very costly process.

Nowadays it is impossible to design innovative structures without FEM and the use of computational programs based on this method became easier due to the development of pre-processor and post-processor tools that provide interactive graphics resources.

Nonetheless, there are phenomena whose behavior cannot be satisfactorily described by conventional FEM and this fact has motivated the development of new strategies. Problems subjected to large deformations and to crack and damage propagation require modifications in discretization of the structure (remeshing) and methods such as the Generalized Finite Element Method (GFEM) have been developed to solve these issues.

The GFEM (Melenk and Babuška [^{1}]; Duarte *et al.* [^{2}]) can be considered as originated from the so-called meshless methods proposed in the 1990s. In spite of its theoretical bases be well established, there is an extensive area of research and of numerical experimentation to be investigated. According to Barros *et al.* [^{3}], GFEM is formulated in a way that the numerical simulation guarantees certain independence of the mesh of finite elements. The relative mesh independence can be observed by the possibility of introducing special functions on the numerical approximation, without modifying the mesh, and by the relative insensitivity to angular distortion of the elements.

The objective of this paper is to evaluate the use of linear and quadratic polynomial enrichment functions of the GFEM approach in nonlinear problems of concrete structures. The obtained results with GFEM and FEM strategies combined with an anisotropic constitutive model, based on the microplane theory, will be compared to each other. The validation of the results obtained with the GFEM approach is made through comparisons with experimental responses of concrete beams subjected to the three-point bending (Petersson [^{4}]) and of L-shaped panel by Winkler *et al.* [^{5}].

The outline of the paper is as follows. A brief explanation of GFEM is presented in Section 2. The Microplane Theory for anisotropic constitutive models, based in the formulation of Leukart [^{6}] and used in all simulations presented here, it is discussed in Section 3. In the Section 4, numerical examples are presented, and Section 5 is devoted to concluding remarks and discussion.

2. A summary explanation of GFEM

Initially proposed by Melenk [^{7}], the GFEM has the important characteristic of enriching the space of approximations originally constructed to the FEM with the introduction of special functions. Such functions are able to reproduce phenomena observed in specific regions of the analyzed domain.

Results presented by Strouboulis *et al.* [^{8}] and Duarte and Oden [^{2}] showed that the enrichment strategy with functions that reflect the local nature of the solution contributed to improve the quality of the approximation, with the introduction of a relatively small number of degrees of freedom.

In GFEM, the approximation functions are constructed as in the hp-Cloud Method (Duarte and Oden [^{9}]; [^{10}]), where cloud points are distributed arbitrarily and without links to each other, whereas in the GFEM a finite element is adopted for nodal positioning and structuring of the shape functions.

The use of the Partition of Unity (PU) functions (based on the method of Melenk and Babuška [^{1}]) over the mesh of finite elements and its enrichment carried out in a way analogous to the hp-Cloud Method make GFEM an unconventional approach of the FEM and allow it to be correlated with the Meshless Methods.

The strategy used in GFEM consists of enriching PU functions to define alternative shape functions. The standard functions of FEM (such as Lagrangian functions) facilitate the application of the GFEM and, differently from the meshless methods, directly verify the boundary conditions (Barros *et al.* [^{3}]).

In the GFEM the PU functions are constructed in regions called clouds ω_{j}. However, these clouds differ from those of the Meshless Methods and they are constituted by sets of finite elements sharing the nodal points x_{j}, according to Figure 1 (N_{j} (x) represents the PU functions).

For example, in ℝ^{1} the PU is obtained from linear functions associated with each cloud.

The enrichment functions are multiplied by the original PU, guaranteeing the improvement of the quality of the approximation. Aiming to clarify this strategy, it is considered a conventional mesh of finite element defined from a set of n nodal points
^{2}. It is defined a patch or cloud ω_{j} formed by all elements that share the nodal point x_{j}.

The set of interpolative Lagrangian functions associated with the node x_{j} defines the function N_{j} (x) whose support corresponds to the region ω_{j}, according to Figure 2(b).

A set of enrichment functions, I_{j}, named local approximation functions, is composed by q_{j} linearly independent functions defined to each node x_{j} with support on the cloud ω_{j}:

At the end of the process, the shape functions ϕ_{ji} (x) of GFEM, shown on Figure 3(b), associated with the node x_{j} are built through the enrichment of the PU functions by the components of the set I_{j}.

Thus, according to equation (3), ϕ_{ji} (x) (Figure 3(b)) can be obtained by the product between the basic functions that form the PU (Figure 2b) and the enrichment functions (Figure 3a).

The functions in equation (2) can be polynomial or not depending on the problem analyzed. The use of the functions of FEM as the PU simplifies the implementation and avoids, according to Barros *et al.* [^{3}], problems related to the numerical integration and to the imposition of the boundary conditions. Thus, a generic approximation

where u_{j} and b_{ji} are nodal parameters associated with standard (N_{j}) and GFEM (N_{j} (x) × L_{ji} (x)) shape functions, respectively.

Furthermore, aiming to minimize round-off errors, Duarte *et al.* [^{2}] suggest that a transformation should be performed over the L_{ji} (x) functions, if they are of polynomial type. In such case, the coordinate x is replaced as follows:

in which h_{j} is the diameter of the smallest circle that circumscribed the cloud ω_{j}.

It is obtained the product function that presents the characteristics of the local approximation function while inherits the compact support of the PU.

Among some advantages of the GFEM when compared to the standard FEM, it is possible to mention the enrichment of the approximations to treat specific problems with functions specially defined to this purpose. Thus, the enrichment can be performed in regions of interest in which the behavior is more pronounced without elevating computationally the analysis.

3. Anisotropic constitutive models based on the microplane theory

The adaptation of Microplane Theory to concrete structures occurred from the association between solid structure of the heterogeneous material (cementitious matrix with aggregates of different particle sizes) and the existence of multiple plans of discontinuities, positioned at the interfaces of its grains.

This association is quite pertinent because of the occurrence and propagation of microcracks in different directions that lead to an inelastic response of the material. Such propagation generally occurs at the existing interface between the cementitious matrix and the aggregates, as shown in Figure 4.

The formulations of the Microplane Models generally follow three main steps: the projection of the strains in the microplanes, the definition of the constitutive laws and the homogenization process, as can be seen in Figure 5.

The model proposed by Leukart [^{6}], and used in all simulations of this paper adopts: (step 1) a decomposition of the macroscopic strain tensor into its volumetric (ε^{V}) and deviatoric (ε^{D}) components (V-D split); (step 2) that the damage process is the main dissipation mechanism which describes the degradation on the material and that the degradation is evaluated through a single equivalent strain combined with a single damage law; (step 3) that the free energy on the microplanes exists and its integral over all microplanes is equal to the macroscopic free energy of Helmholtz.

Wolenski [11] generalizes the computational implementation of Leukart [^{6}] proposition, in order to allow any microplane equivalent strain measure and any damage law. Such an improvement has been implemented in the context of the Unified Computational Environment, proposed by Penna [^{12}] and Gori *et al.* [^{13}], on the **INSANE**. This system has been expanded by Alves *et al.* [^{14}] with the enclosing the standard version of GFEM formulation with minimum impact in the code structure. Based on this expansion an object oriented design of GFEM to physically nonlinear analysis has been extended by Monteiro *et al.* [^{15}], being used in all simulations of this paper. Further details about the **INSANE** system can be found in **INSANE** Project [^{16}].

The numerical simulations presented in this paper use one of the options of the unified environment for microplane models of the **INSANE** system. Specifically, the simulations uses volumetric-deviatoric strain split proposed by Leukart [^{6}] and the equivalent strain defined by de Vree *et al.* [^{17}], according to:

where ε^{V} is volumetric part of the strain tensor, ε_p is the p component of the deviatoric strain tensor, k_{1} and k_{2} are material parameters that relate to tensile and compression strength of concrete and η_{Vree} is equivalent strain measure adopted on the different damage laws such as the exponential, polynomial, and linear laws defined by equations (7), (8) and (9), respectively:

where d^{mic} is the damage measure, κ is the current equivalent strain (equation 6), κ_{0} and κ_{u} are material parameters that specifies a limit for κ referring, respectively, to the beginning and end of the damage process, while E_{0} is the Young’s modulus.

The parameter α is the maximum material degradation, β is the parameter governing the shape of the post-peak branch, f_{e} is the equivalent stress related to the material strength limit. They are dimensionless numerical parameters of the constitutive model.

These formulations are detailed in Wolenski [^{11}] and they were adopted on the numerical simulations presented in Wolenski *et al.* [^{18}].

4. Numerical simulations

In this section some numerical simulations are presented aiming to illustrate and to validate discussions about the use of the anisotropic constitutive model, which is able to represent the behavior of concrete structures together with GFEM approach.

These simulations also allow illustrating the use of the GFEM according to the characteristics of the analyzed problem, providing different ways of investigating the problem and performing the relevant numerical simulations, with the choice of the enrichment functions, the nodes to be enriched, the combination of different functions in the same problem, among other possibilities.

For all simulations, the approximation functions of enrichment, with monomials expressed in coordinates x and y, are defined by:

The solution procedure to nonlinear problems described here is the Newton-Raphson Algorithm and other characteristics of each simulation are detailed in the respective item. The geometry and boundary conditions are illustrated throughout each problem, since the aiming is to use distinct types of elements combined with different enrichment strategies. The results obtained in the simulations are compared with the experimental ones presented in literature.

*4.1 Three-point bending*

Petersson [^{4}] experimentally studied concrete beams subjected to three-point bending. The experimental results obtained by the author were used by Monteiro at al. [^{15}] and extended here to compare with the numerical results of the simulations performed with GFEM and FEM approaches.

The material parameters obtained experimentally by Petersson [^{4}] and adopted to the numerical simulations are shown in Table 1.

Young's modulus | Poisson ratio | Uniaxial yield stress | Fracture energy |

E_{0} = 30000 MPa |
ν_{c} = 0,20 |
σ_{t} ? 3,0 MPa |
G_{f} ? 0,130 N/mm |

Figure 6 depicts the geometric data of the beam and the mesh of 115 four-node quadrilaterals elements (with 4×4 Gauss Quadrature). The numerical simulations are performed with the conventional FEM (260 degrees of freedom and with no enrichment - P_{0}) and with GFEM applying linear (P_{1}) and quadratic (P_{2} enrichment functions highlighted at the nodes of the figure (384 degrees of freedom).

Table 2 presents the numerical parameters to the different damage laws and they are grouped according to the equivalent strain defined by de Vree *et al.* [^{17}].

Exponential damage law | ||

α = 0,960 | β = 30000,00 | κ_{0} = 0,0002 |

Polynomial damage law | ||

f_{e} = 5,95 |
E_{0} = 30000,00 |
κ_{0} = 0,000385 |

Linear damage law | ||

κ_{u} = 0,00460 |
κ_{0} = 0,000190 |

In possession of such parameters, the nonlinear analyses have been performed under plane stress conditions and with the adoption of the generalized displacement control method (Yang and Shieh [^{19}]), with initial load factor equal to 0,020, tolerance to convergence of 1 × 10^{-4}(×100%) = 0,010% in relation to the norm of the incremental displacements vector, a reference load of P = 800N and a secant approximation to the constitutive tensor.

Figures 7, 8 and 9 show the equilibrium paths, describing the vertical displacement of the node 10, together with the experimental results obtained by Petersson [^{4}], using different damage laws, according to equations (7), (8) and (9). The equilibrium path GFEM-P1+P2 refers to the trajectory obtained using GFEM approach and the equilibrium path denominated FEM-P0 refers to the trajectory obtained by the FEM approach.

From Figures 7, 8 and 9 it is possible to observe that the GFEM approach provided stability to the equilibrium paths and agreement with those experimentally obtained by Petersson [^{4}]. The enrichment strategy improved the solution around the region where the nonlinear phenomenon happens and it does not require modification in the neighboring elements of that region, as it would be the case to FEM if selective *h* or *p* refinements are applied.

Additionally, it is recognized that the standard FEM analysis presented a worse description of the equilibrium path due to the poorer discretization adopted. A ﬁner mesh around the notch could provide a better solution. In such case, special attention would be required to avoid problems related to the transition from the bigger to the smaller elements. Another strategy could be using higher order elements on that same region, but the inclusion of irregular nodes lead to a constrained approximation.

Figure 10 shows the number of iterations to achieve equilibrium for each load step. This figure allows inferring that the number of iterations GFEM does not significantly vary, being smaller in some situations and greater in others.

*4.2 L-shaped panel*

The numerical simulations of a L-shaped panel are presented to discuss the influence of mesh refinement and enrichment functions on the numerical response performed with GFEM and FEM approaches, as well as the use of the Microplane Constitutive Model by Leukart [^{6}].

Such simulations are compared with the result presented by Winkler *et al.* [^{5}] that performed experimental tests on concrete panel, according to the geometry and boundary conditions presented in Figure 11(a).

The numerical simulations are performed with triangular finite elements (T3) in the following way: conventional FEM with 410 elements (Figure 11(c)) (476 degrees of freedom and with no enrichment - P_{0}) and GFEM with 278 elements applying enrichments P_{1} in 142 nodes and P_{2} in 23 nodes from the center corner of the panel, as illustrated in Figure 11(b) (1082 degrees of freedom).

It was adopted enrichment functions P_{2} for an area with a higher occurrence of damage, in order to achieve the same effect of mesh refinement via FEM in this region. A conventional FEM analysis (Figure 11(c)) was applied considering all nodes without enrichment (P_{0} ).

The material parameters experimentally obtained by Winkler *et al.* [^{5}] and adopted to the numerical simulations are shown in Table 3.

Young's modulus | Poisson ratio | Uniaxial yield stress | Fracture energy |

E_{0} = 25850 MPa |
ν_{c} = 0,18 |
σ_{t} ? 2,70MPa |
G_{f} ? 0,065 N/mm |

From Table 3 it was possible to obtain the numerical parameters used in the simulations, according to Table 4.

In possession of such parameters, the nonlinear analyses have been performed under plane stress conditions and it is adopted generalized displacement control method, with initial load factor equal to 0,03, tolerance to convergence of 1 × 10^{-4} (×100%) = 0,010% in relation to the norm of the incremental displacements vector, reference load of q = 32 N/mm and a secant approximation to the constitutive tensor.

Figure 12 shows the numerical results of the displacement control point shown in Figure 11(a), together with the experimental results obtained by Winkler *et al.* [^{5}].

Both results show good concordance with the experimental ones obtained by Winkler *et al.* [^{5}]. However, it is noted better stability along of the inelastic regime obtained with GFEM approach and the equilibrium path was closer to the experimental result than the equilibrium path obtained with standard FEM. Only the insertion of enrichment functions into the coarser mesh, without the need for refinement, was able to improve the result in relation to the FEM.

To illustrate the evolution of panel degradation throughout the analysis, the Figures 13 and 14 shows the damage for both meshes 1 (GFEM) and 2 (FEM), respectively.

These figures show that the propagation of damage represents the expected behavior, since the degradation begins in a concentrated manner in the center of the panel and propagates horizontally along of the panel. In this sense, the mesh with a greater refinement presented a behavior slightly closer to that obtained by Winkler *et al.* [^{5}], for the experimental and numerical cases, as illustrated in Figure 15.

The substitution of finite element mesh refinement by GFEM approach proved to be efficient for the problem in question, indicating the flexibility and feasibility of this feature to the physically nonlinear analysis.

5. Final remarks

In this paper a summary of GFEM formulation and of the Anisotropic Constitutive Models, based in Microplane Theory, were made. It was analyzed a concrete beam subject a three-point bending and a L-Shaped Panel, and a computational framework for nonlinear analysis by GFEM approach combined with the Microplane Constitutive Model was used.

The obtained results were compared with the experimental ones provided by Petersson [^{4}] and Winkler *et al.* [^{5}] and they showed good agreement in all simulations performed to different damage laws. It was possible to verify the versatility of the GFEM approach in **INSANE** system because of the application of different polynomial enrichment functions combined with different damage laws of the Microplane Model by Leukart [^{6}].

In both simulations (Sections 4.1 and 4.2) the enrichment strategy under the GFEM approach provided stability to the equilibrium paths. The polynomial enrichment strategy provided by GFEM allows improving the quality of the approximation in a very simple and straightforward way, without overloading the analysis with a very refined mesh.

All the simulations demonstrated that GFEM are able to reproduce the results of FEM and even improve the solutions just applying enrichment functions in some nodes of the meshes, allowing varied analysis in which the reﬁnement of the mesh (that can lead to numerically induced strain localization) is not necessary to achieve better solutions.

New investigations can be performed aiming to verify the numerical stability of GFEM to the nonlinear analysis, mainly when its application requires enrichments of higher order or a large number of finite elements. Additionally, it is possible to apply other constitutive models and also to apply the GFEM global-local implementation to nonlinear analysis, which is under development.

These aforementioned research themes are among the future works of our research group, aiming to empower the INSANE computational platform to solve wider range of the solid mechanics problems. Finally, a data set along within put files of this work can be found at https://doi.org/10.6084/m9.figshare.5505898 [^{20}], in order to reproduce the results here presented.