Elastoplastic constitutive modeling for concrete: a theoretical and computational approach

Resumo This article presents a study of the plasticity model applicability to concrete in a theoretical framework that generalizes the formulation of constitutive models for physically nonlinear analysis of structures. In this sense, the theoretical framework for the computational implementation of the plasticity mathematical theory is described, detailing the models formulations capable to describe the inelastic behavior of concrete. The loading surfaces associated to Drucker Prager and Ottosen criterion are highlighted. Furthermore, the Cutting Plane return mapping algorithm, necessary to the integration of constitutive relations that govern the behavior of the material in the context of computational plasticity, is described. Finally, numerical simulations are presented, such as the direct tension loading and three-point bending tests. The results of these simulations are compared with those from the literature to verify the model stability and accuracy.


Introduction
A realistic solution for a structural problem involving concrete depends in large part on the choice of an appropriate constitutive model. The mechanical response of the concrete is complex and it seems unlikely that any phenomenological approach would be able to understand all the possible variations in the characteristics of the material. Even if a perfect model for the concrete could be built, it would be too complex to serve as a basis for the stress analysis of practical problems (CHEN; HAN [1]). However, intensive studies over the past decades have led to a better understanding of the constitutive behavior of quasi-brittle media. Research focused on modeling the mechanical behavior of concrete has led to formulations such as mathematical theory of plasticity, a necessary extension of the theory of elasticity which provides a more realistic approach about the behavior of the material. The theory of plasticity seeks to mathematically describe immediate and non-reversible deformations that occur in a solid body, i.e. the deformations that do not disappear completely when the causal forces are removed (CHEN; HAN [1]; LUBLINER [2]; SOUZA NETO et al. [3]). The description of the elastoplastic behavior of the concrete on the macroscopic point of view (designated as phenomenological behavior) for multiaxial stress states is based on the following assumptions: (1) The existence of an elastic field, that means, a region in which the material behaves as it is purely elastic without the appearance of permanent strains. The elastic domain is delimited by a flow function written in terms of yield stress. (2) The occurrence of inelastic deformation, that is, the deformations associated to stresses above the yield stress whose evolution can be described by a yielding rule. (3) The occurrence of the phenomenon of strain hardening of the material, in other words, the possibility of increasing flow stress, following the evolution of inelastic deformations. Rupture criteria for concrete are classified according to the num-ber of parameters that appear in the defining expressions. Simple models should be used, representing only those properties that are essential for the problem to be considered (CHEN; HAN [1]). To this purpose, this work aims to present a theoretical and computational framework necessary for the implementation of elastoplastic constitutive models, especially the models of Drucker Prager and Ottosen on INSANE computing system (Interactive Structural Analysis Environment), showing also the Cutting Plane return mapping algorithm, necessary for the integration of the constitutive relations governing the behavior of the material in the context of computational plasticity.

Plasticity for concrete
The classical theory of plasticity was originally developed for the study of metals and some of its proposed fundamentals are not safe for other engineering materials such as concrete. However, they still have some similarities, particularly in pre-failure regime. For example, the concrete exhibits a nonlinear stress-strain behavior during loading and has a substantial irreversible deformation under unloading regimen. Especially under compressive loads with confining pressure, the concrete may show some ductile behavior. Thus, concrete irreversible strains are induced by microcracks and can be treated by the theory of plasticity (CHEN; HAN [1]). A variety of constitutive models were then proposed in order to mathematically reproduce the material stress-strain relationships for different load conditions. The majority adopts a phenomenological approach, i.e., models that describe the material from the macroscopic point of view, neglecting the microscopic mechanisms, and considering the material media as continuous and homogeneous. The approach of plasticity falls into this category. A constitutive model suitable for concrete structures requires a complete description of the material behavior, as that one showed in Figure [  From a macroscopic point of view, the classical plasticity can simulate the behavior of the concrete particularly in the pre-peak regime, such as nonlinearity of the stress-strain curve and the irreversible strains after loading. Many papers have been presented by researchers seeking to adapt the classical theory of plasticity in order to get a better representation for concrete (PARK; KIM [5]; GRASSL; JIRASEK [6]). In the concrete plasticity modeling, it is important to observe some characteristics such as sensitivity to hydrostatic pressure, not associative flow rule, compatible inelastic law, and tensile strength limit. Some of these models are mathematically highly complex, making them undesirable for many engineering applications, especially the analysis of simple structural elements. The failure surface definition according to phenomenological models is performed by using a yielding parameter. The rupture surface of the concrete can be generally expressed by: where, I 1 is the first invariant of the stress tensor and J 2 and J 3 are respectively the second and the third principal invariant of the deviatoric stress tensor. These invariants are generally represented by the following expressions: wherein, σ 1 , σ 2 and σ 3 are the values of the principal stresses and e s ij are the components of the deviatoric stress tensor given by: The explicit form of the failure function for concrete is defined by experimental data, since concrete resistance tests are well documented in the literature (CHEN; HAN [1]). Such functions shall have the following characteristics: (1) Represent a smooth convex surface, with the exception of its apex; (2) The meridians of its surface are parabolic, and open towards the negative hydrostatic axis; (3) Failure curve is roughly triangular to tensile stresses and low compressive stresses, becoming more circular as the compressive stress increases.

The Drucker-Prager Model
The criteria proposed by D. C. Drucker and W. Prager in 1952 requires that the plastic flow occurs when the second invariant of the deviatoric stress tensor, J 2 , and the hydrostatic pressure reaches a critical combination. The function that models the criteria of Drucker-Prager is given by: where η is a constant related to the material and the function k. In case the material undergoes an isotropic hardening, the function k relates to the uniaxial stress-strain curve (CHEN; HAN [1]) and can be defined by: where σ is a function of the hardening internal variable α. The yielding surface in the principal stress space is represented by a circular cone whose axis is the axis of symmetry of the hydrostatic pressure. The Drucker-Prager yield surface is shown in Figure [2].

The Ottosen Model
The rupture surface of four parameters (α, β, and c 1 and c 2 ) for concrete, proposed by OTTOSEN [7] involves the strain invariants I 1 , J 2 , and the loading angle, θ. Its smoothness, convexity and its curved meridians that have a gradual transition from an almost triangular shape to an almost circular in the deviatoric plan, as the hydrostatic pressure increases, make this criterion suitable for failure simulation of concrete structures ( Figure [3]). The mathematical representation of the criteria is given by: Where σ c is a function of the hardening internal variable, α and β are material parameters. The parameter λ depends on two other parameters (c 1 and c 2 ) and is given by: Still in the equation (9), the angle involved in this criteria is given by: There are several propositions for the determination of the four

Elastoplastic constitutive modeling for concrete: a theoretical and computational approach
There are several propositions for the determination of the four parameters (α, β, c 1 e c 2 ) of the Ottosen model. The INSANE system has implemented the calibration proposals of OTTOSEN [9], CEB-FIP Model Code [10] and DAHL [11]. According to OTTOSEN [9], the four parameters can be determined based on the following tests: (1) f c -uniaxial compressive strength (θ = 60°); (2) f t -uniaxial tensile strength (θ = 0°); (3) f bc ≅ 1,16 f c -biaxial compressive strength (θ = 0°); (4) -the triaxial stress state on the compressive meridian (θ=60°). The values obtained for the parameters from these tests depend on the average tensile and compression ratio k = f tm / f cm . Table [1] shows some of the most commonly used values from this calibration. The parameters for intermediate values of k can be obtained by interpolation. Another way to obtain the model parameters is through the expressions recommended by the CEB-FIP Model Code [10], which also make use of the relation k = f tm / f cm .
This calibration allows obtaining the parameters to any values of k automatically, provided that the compressive strength values do not exceed prescribed values. DAHL's [11] proposal is based on the observation that the CEB's recommendations are in agreement with the experimental results only for low strength concrete, thus suggesting a way to obtain the coefficients using only the mean resistance to compression of the concrete (f cm ). NETO [12] describes the hardening as a process that is physically connected to increased dislocation density (geometric defect in atomic arrangement). For many real materials, the yielding stress limit of the material is dependent on a measure of accumulated plastic strain. In uniaxial model, after reaching the yielding, the stress-strain curve continues to grow (in hardening) or decreasing (in the case of softening) causing a variation in yielding stress during plastic flow. In models of two and three dimensions, the hardening is characterized by changes in the set of internal variables α during plastic yielding. These changes can generally affect the size, shape and orientation of the yielding surface, defined by Φ (σ, α) = 0. Depending on the type of material, the stress-strain curves can have different forms, being convenient to idealize some of these behaviors. Figure [4] illustrates three models commonly used to describe materials that have an elastoplastic behavior (MALAVOLTA [13]). Case (a) corresponds to the perfect plasticity model with a material having an elastic portion with modulus of elasticity E, and after yielding, the material remains with a constant level of stress as the strain increases. Model (b) is called bilinear elastoplastic, where the first slope corresponds to the elastic portion and, after reaching the yielding, begins a new line with an inclination H associated with the material hardening, corresponding to the plastic domain. Case (c) is the nonlinear elastoplastic model, which once reached the yielding, the hardening shall be described by a nonlinear law. In this paper has been used a linear and potential hardening laws. The linear hardening law is a generally used law, and it is defined by:

Hardening laws
The potential law, proposed by BOUCHARD et al [14] and VAZ JR; ROJAS [15], is characterized by:  where a e b are material parameters and k is the internal variable related to the hardening or softening.

Cutting plane return mapping algorithm
The elastoplastic analysis requires the integration of the constitutive law, so that the elastic and plastic portions of the total strain increment are obtained, which is an iterative process, due to the fact that the elastoplastic module is a function of plastic deformation. For the incremental nature of the numerical models of plasticity, a return mapping algorithm, able to obtain the update of the stresses needed to balance internal forces in the nonlinear analysis, must be used. The development of efficient schemes for the integration of constitutive relations in a numeric context is still subject of recent research in the world, mainly because of its importance in engineering problems involving plastic deformation. There are a variety of methods of integration with different levels of complexity (TAQIEDDIN [16]). Conceptually, the idea of the algorithm is quite simple and consists of an explicit process that includes the first elastic equations from the stress in the previous step, to obtain the trial stress on the current step. One of the great advantages of Cutting Plane is the fact that there is no need to compute gradients of yielding function and hardening law, as this task can be extremely cumbersome for complex plasticity models. The general case of this scheme involves the following steps: (1) Assume the existence of a plastic load, i.e., , such that λ > 0. Define the residual plastic flow R n+1 and the yield condition.
(3) Update the state variables and the consistent parameters The algorithm convergence to the final value of the state variables is obtained in a quadratic rate. These quadratic convergence rates are achieved here in spite of the simplicity provided by the method, which end up making the cutting plane algorithm very attractive for large-scale calculations in more elaborate models, mainly in the explicit codes that do not require the solution of a global system of equilibrium equations.

Formulation of constitutive models
Constitutive models typically have a proper notation and, although in many cases they keep similarities, the lack of unity of formulations prevents a generic and objective computational implementation. The constitutive models framework proposed by PENNA [17] and GORI et al. [18] presents an expansion of the theoretical framework proposed by CAROL et al. [19], being able to contemplate various constitutive models -elastoplastic or elastic degradation; isotropic, orthotropic, or anisotropic -formulated with one or more loading functions. Next, is presented the design for the elastoplastic constitutive models following the theoretical proposed framework. The plots necessary to the description of each model are explained indicating the correlation between the original form and this work objectives.

Drucker-Prager Model
The mathematical representation of the Drucker-Prager criterion is given by the function: For the case of isotropic hardening, the function σ (α), used for the Elastoplastic constitutive modeling for concrete: a theoretical and computational approach determination of parameter k (given by [7] equation) is given by: (19) in which σ y is the initial yield stress, α is the accumulated plastic flow, and H is the strain hardening modulus. The Drucker-Prager model is non-associative (F ≠ Q). The Drucker-Prager non-associative law is obtained by adopting for the plastic potential function, a similar function to the yielding function, wherein the angle of friction, ϕ, is replaced by dilatancy angle, ψ. (20) where ψ < ϕ and is an additional constant of material. Therefore, F and Q gradients represented by tensors m and n are given by: The inelastic modulus, obtained from a law associated with hardening or softening phenomenon, is given by: The application of Drucker-Prager model should take into account the existing singularity on the yielding surface, its apex. Therefore, it must be used an alternative solution strategy for the implementation of the constitutive relations integration algorithm. Various methods have been proposed in the context of yielding surfaces with singularities like corners and vertices, such as SIMO; HUGHES [20] and SOUZA NETO et al. [3]. When the return occurs at the vertex, the yielding function (equation [18]) and plastic flow potential (equation [20]) should be changed to (SZABÓ; KOSSA [21]):

Ottosen model
The Ottosen criteria is given by the following yielding function where in σ c = (σ y + Hα). The derivative of F for an isotropic material may be obtained by chain rule to: in which the invariant stress derivatives are: where, δ ij is the Krönecher delta, s ij are the tensor components of deviatoric stress and t ij is the quadratic tensor of the deviatoric stress. The derivatives from the yielding function in relation to invariants are:  the model adopted is associative (n ij = m ij ). The inelastic module associated with hardening or softening phenomenon, is given by:

Diametral compression test
The diametral compression test is commonly used to determine the tensile strength of concrete and consists of applying diametrically opposed loads on a cylindrical specimen in order to produce an indirect traction in its central region. In this sense, a plasticity criteria can be used to determine the failure stress. Therefore, the Drucker-Prager's model with the internal approach of the Mohr-Coulomb surface was adopted. The model geometry and loading conditions are specified in Figure [5]. The material parameters are shown in Table [2] and were based on the study made by CECÍLIO [22]. The finite elements model ( Figure [5-b]) consists of 404 quadrilateral four-nodes elements in plane strain state, with a thickness of 300 mm and 2x2 points for integration. For the nonlinear analysis, direct displacement control method was adopted with increment of 0,00002179 mm, controlling the horizontal displacement of the highlighted node in Figure [5] with tolerance for convergence of 5 × 10 -3 and reference load of P = 60 kN. The result of stresses in the specimen's center, while the material is in the elastic region may be obtained analytically depending on the applied load, P, on the diameter D and on the length L, by the equation: The simulation was not able to describe the inelastic behavior of the specimen. However, the simulation could represent the material behavior in the elastic range. Figure [6] shows the instant when the plastic strain appears in the specimen (step 111 of 250 steps), represented by the accumulated plastic flow variable, illustrating the elastic limit of the material. The stress distribution is shown in Figure [7], where the condition of the indirect tensile test is clearly observed. The tension stress limit in the elastic range, according to equation [35] on the presented conditions, is:

Direct tension of a concrete plate (Dogbone-shaped panel)
This example shows a finite element model to simulate the experimental behavior obtained from a direct tensile test on a concrete flat specimen reinforced with fibers. Due to the symmetry, only a quarter of the plate was discretized. The boundary conditions and finite element mesh adopted are presented in Figure [10]. The plate was discretized with 12 quadrilateral eight-nodes elements in plane stress state, with 3 x 3 points integration.
The material is considered a cement matrix composite reinforced with fibers with 2% of vinyl polyvinyl acetate (PVA), according to PEREIRA et al. [23]. The material parameters of the experiments are given in Table [3]. For the numerical simulation, the Ottosen's model with linear hardening law was adopted, using the generalized displacement control method, with initial load factor of 1,0 and tolerance for

Figure 8
Distribution of normalized stress σ x along the y axis

Figure 9
Distribution of normalized stress σ y along the y axis

Figure 10
Dogbone-shaped panel -Finite element mesh, geometry and image of supports used during the test performed by PEREIRA et al [23]. Adapted from PEREIRA et al. [8] convergence of 1× 10 -4 . The Stress x Strain curve is shown in Figure [ 11] in comparison with the experimental results obtained by PEREIRA et al. [23].
The numerical model has showed capable to simulate the behavior of the material and was able to represent the experimental results with good accuracy. Figure [12 -a] shows the deformed mesh and it's observed that the orthogonal lines along the longitudinal axis of the plate remained parallel after deformation. The pattern of evolution of displacements, after reaching the yield stress obtained with the Ottosen's model, implemented in this work, can be seen in Figure [12 The results are in excellent agreement with the experimental and numerical results presented by PEREIRA et al. [8] (Figure [12-c]). In order to attest the convergence behavior, 4 mesh were adopted, with 3, 12, 24 and 128 quadrilateral eight-nodes elements. The meshes are schematically represented in Figure [ 13]. Figure [14] shows the load factor-displacement curves. The results present the convergence of the solution and indicates no approximation errors related to the discretization and the mesh refinement. However, more tests should be performed in order to attest the mesh sensitivity and the general behavior of the model under highly refined meshes.

Reinforced concrete beam
In   Elastoplastic constitutive modeling for concrete: a theoretical and computational approach The Ottosen's criteria with potential hardening law for concrete was adopted. The materials parameters are given in Table [4]. For steel, the elastoplastic criterion of von Mises was adopted, assuming perfect plasticity and perfect bond between steel and concrete. The materials parameters are given in Table [5].
The Figure [16] shows the discrete model, consisting of 132 quadrilateral eight-nodes elements to represent the concrete, and 22 onedimensional three-nodes elements to represent the reinforcement. In the analysis, the direct displacement control method was adopted, incrementing of 0,001 mm the horizontal displacement of the supported right node, tolerance for convergence of 1 × 10 -4 and reference load P = 1,0N. The model was analyzed considering plane stress conditions.

Figure 15
Reinforced concrete beam Load-Displacement curves corresponding to the vertical displacement of the point where the load is applied are shown in Figure [17]. The numerical results are compared to the experimental values presented by MAZARS and PIJAUDIER-CABOT [24]. The graphic shows good agreement between the experimental results and the results obtained with Ottosen's model, even that the model presents a higher initial stiffness and a lower yielding load when compared with the experimental values.

Conclusion
In this paper an elastoplastic constitutive model applied to concrete, emphasizing the criteria of Drucker-Prager and Ottosen, has been presented. In addition to the constitutive models, also has been presented equations for the implementation of the Cutting-plane return mapping algorithm, required for the integration of the constitutive relations governing the behavior of the material. The constitutive models have been implemented in the computational system INSANE (INteractive Structural ANalysis Environment), according to the theoretical and computational environment for constitutive models developed by PENNA [17] and GORI et al. [18]. Classical models of associated and non-associated plasticity were easily incorporated in the theoretical framework, varying only the return algorithm according to the needs of each model.
Numerical simulations presented in order to illustrate, validate and emphasize the individual characteristics of each model. From the results presented, the following considerations can be made: i. Models showed appropriate behavior, and the responses were consistent; ii. Models showed appropriate behavior using the Cutting Plane return mapping algorithm; iii. By analyzing all numerical simulations presented, the constitutive models representing the elastoplastic behavior of concrete showed a good correlation among numerical, experimental and analytical results.

Figure 16
Finite element mesh

Figure 17
Numerical results obtained with Ottosen's model compared with experimental results of reinforced concrete beam performed by MAZARS and PIJAUDIER-CABOT [24]