Elastoplastic constitutive modeling for concrete: a theoretical and computational approach

Modelagem constitutiva elastoplástica para o concreto: uma abordagem teórica e computacional

Abstract

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.

Keywords:
computational plasticity; constitutive models; hardening laws; return mapping

Resumo

Este artigo aborda modelos de plasticidade aplicados ao concreto em um framework teórico que visa generalizar a formulação de modelos constitutivos para a análise fisicamente não lineare de estruturas. Nesse sentido, é descrito o arcabouço teórico para a implementação computacional da teoria matemática da plasticidade, detalhando as formulações de modelos capazes de descrever o comportamento inelástico do concreto. As superfícies de carga associadas aos critérios Drucker Prager e Ottosen são destacadas. Além disso, é descrito o algoritmo de mapeamento de retorno “Cutting Plane”, necessário à integração de relações constitutivas que governam o comportamento do material no contexto da plasticidade computacional. Finalmente, são apresentadas simulações numéricas, como os testes de compressão diametral, tração direta e flexão de três pontos. Os resultados dessas simulações são comparados com os da literatura para verificar a estabilidade e a precisão do modelo.

Palavras-chave:
plasticidade computacional; modelos constitutivos; leis de encruamento; mapeamento de retorno

1. 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 [¹]).

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 [¹]; LUBLINER [²]; SOUZA NETO et al. [³]).

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 number 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 [¹]). 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.

2. 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 [¹]).

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 [1], wherein the pre-failure (hardening) and post-failure (softening) behaviours are displayed.

Figure 1
Uniaxial behavior of concrete (CHEN; HAN [⁴])

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 [⁵]; GRASSL; JIRASEK [⁶]).

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:

$f\left(I_1,J_2,J_3\right)=0$ (1)

where, I₁ is the first invariant of the stress tensor and J₂ and J₃ are respectively the second and the third principal invariant of the deviatoric stress tensor. These invariants are generally represented by the following expressions:

$I_1={\sigma }_1+{\sigma }_2+{\sigma }_3$ (2)

$J_2=\frac{1}{2}{s}_{ij}{s}_{ji}$ (3)

$J_3=\frac{1}{3}{s}_{ij}{s}_{jk}{s}_{ki}$ (4)

wherein, σ₁, σ₂ and σ₃ are the values of the principal stresses and e s_ij are the components of the deviatoric stress tensor given by:

$s=\mathit{\sigma }-\frac{1}{3}\left(tr\left(\mathit{\sigma }\right)\right)I$ (5)

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 [¹]). 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.

2.1 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₂, and the hydrostatic pressure reaches a critical combination. The function that models the criteria of Drucker-Prager is given by:

$\mathrm{\Phi }\left(\mathit{\sigma },c\right)=\sqrt{J_2}+\eta p\left(\mathit{\sigma }\right)-k$ (6)

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 [¹]) and can be defined by:

$k=\left(\eta +\frac{1}{\sqrt{3}}\right)\sigma \left(\alpha \right)$ (7)

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

Figure 2
The Drucker-Prager yield surface in principal stress space. (SOUZA NETO et al. [³])

2.2 The Ottosen Model

The rupture surface of four parameters (α, β, and c₁ and c₂) for concrete, proposed by OTTOSEN [⁷] involves the strain invariants I₁, J₂ , 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]).

Figure 3
Representation of the Ottosen yield criterion in the deviatory planes using the Haigh-Westergaard stress space, for (a) low and (b) high hydrostatic stress. (PEREIRA et al.[⁸])

The mathematical representation of the criteria is given by:

$\mathrm{\Phi }\left(\mathit{\sigma }\right)=[\alpha J_2+{\sigma }_c{\left(\lambda \sqrt{J_2}+\beta I_1)\right]}^{\frac{1}{2}}-{\sigma }_c$ (8)

Where σ_c is a function of the hardening internal variable, α and β are material parameters. The parameter λ depends on two other parameters (c₁ and c₂) and is given by:

$\mathrm{\lambda }=\left\{\begin{array}{c}{c}_1\cos \left[{\cos }^{-1}(c_2\mathrm{c}\mathrm{o}\mathrm{s}(3\theta \left)\right)/3\right],if\cos \left(3\theta \right)\ge 0\\ c_1\cos \left[{\mathrm{\pi }-\cos }^{-1}(c_2\mathrm{c}\mathrm{o}\mathrm{s}(3\theta \left)\right)/3\right],if\cos \left(3\theta \right)<0\end{array}\right.$ (9)

Still in the equation (9), the angle involved in this criteria is given by:

$\mathrm{\theta }=\frac{1}{3}{\cos }^{-1}\left(\frac{3\sqrt{J_3}}{2J_2^{\frac{3}{2}}}\right)$ (10)

There are several propositions for the determination of the four There are several propositions for the determination of the four parameters (α, β, c₁ e c₂) of the Ottosen model. The INSANE system has implemented the calibration proposals of OTTOSEN [⁹], CEB-FIP Model Code [¹⁰] and DAHL [¹¹].

According to OTTOSEN [⁹], 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) $\left(\frac{\epsilon }{f_c},\frac{\rho }{f_c}\right)=(-5,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.

Table 1
Parameters of the OTTOSEN [⁹]

k	α	β	c1	c2
0.08	1.8076	4.0962	14.4863	0.9914
0.10	1.2759	3.1962	11.7365	0.9801
0.12	0.9218	2.5969	9.9110	0.9647

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 [¹⁰], which also make use of the relation k = f_tm / f_cm.

$\alpha =\frac{1}{9k^{\mathrm{1,4}}}$ (11a)

$\beta =\frac{1}{\mathrm{3,7}{k}^{\mathrm{1,1}}}$ (11b)

$c_1=\frac{1}{\mathrm{0,7}{k}^{\mathrm{0,9}}}$ (11c)

$\begin{array}{c}\\ c_2=1-\mathrm{6,8}{\left(k-\mathrm{0,07}\right)}^2\end{array}$ (11d)

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 [¹¹] 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).

$k=\frac{f_{cm}}{100\left[MPa\right]}$ (12a)

$\alpha =-\mathrm{1,6}{k}^2+\mathrm{3,9}k+\mathrm{0,73}$ (12b)

$\beta =-\mathrm{0,19}{k}^2+-\mathrm{0,41}k+\mathrm{3,13}$ (12c)

$c_1=\mathrm{0,46}{k}^2-\mathrm{0,97}k+\mathrm{11,89}$ (12d)

$c_2=-\mathrm{0,02}{k}^2+\mathrm{0,04}k+\mathrm{0,974}$ (12e)

2.3 Hardening laws

NETO [¹²] 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 [¹³]).

Figure 4
a) Elastic-perfectly plastic; b) Elastic-linear work-hardening model; c) Nonlinear model (MALAVOLTA [¹³])

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:

${\mathcal{H}}_{\left(k\right)}=\mathcal{}\mathcal{}\mathcal{H}\mathcal{}$ (13)

The potential law, proposed by BOUCHARD et al [¹⁴] and VAZ JR; ROJAS [¹⁵], is characterized by:

${\mathcal{H}}_{\left(k\right)}=\frac{b.a.n}{{\left[\mathrm{K}+\left(\mathrm{n}.\mathrm{k}\right)\right]}^{1-\mathrm{b}}}$ (14)

where a e b are material parameters and k is the internal variable related to the hardening or softening.

2.4 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 [¹⁶]).

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., $f_{n+1}^{trial}>0$, such that λ > 0. Define the residual plastic flow R_n+1 and the yield condition.

$f_{n+1}^{\left(k\right)}≔f({\mathit{\sigma }}_{n+1}^{\left(k\right)},{\mathit{q}}_{n+1}^{\left(k\right)})$ (15)

(2) Get increments of consistent parameters.

${\mathrm{\Delta }}^2{\mathrm{\gamma }}_{\mathrm{n}+1}^{\left(k\right)}=\frac{f_{n+1}^{\left(k\right)}}{\left[{\partial }_{\sigma }{f}_{n+1}:C_{n+1}^{\left(k\right)}:{\partial }_{\sigma }{f}_{n+1}^{\left(k\right)}+{\partial }_{\sigma }{f}_{n+1}^{\left(k\right)}.n_{n+1}^{\left(k\right)}\right]}$ (16)

(3) Update the state variables and the consistent parameters

${\epsilon }_{n+1}^{p^{(k+1)}}={\epsilon }_n^{p^{\left(k\right)}}+{\mathrm{\Delta }}^2{\mathrm{\gamma }}_{\mathrm{n}+1}^{\left(k\right)}{\partial }_{\sigma }{f}_{n+1}^{\left(k\right)}$ (17a)

$q_{n+1}^{(k+1)}=q_n^{\left(k\right)}+{\mathrm{\Delta }}^2{\mathrm{\gamma }}_{\mathrm{n}+1}^{\left(k\right)}{n}_{n+1}^{\left(k\right)}$ (17b)

$\mathrm{\Delta }{\gamma }_{n+1}^{\left(k+1\right)}={\mathrm{\Delta }\gamma }_{n+1}^{\left(k\right)}+{\mathrm{\Delta }}^2{\mathrm{\gamma }}_{\mathrm{n}+1}^{\left(k\right)}$ (17c)

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.

3. 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 [¹⁷] and GORI et al. [¹⁸] presents an expansion of the theoretical framework proposed by CAROL et al. [¹⁹], 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.

3.1 Drucker-Prager Model

The mathematical representation of the Drucker-Prager criterion is given by the function:

$\mathrm{F}\left(\mathrm{\sigma },\mathrm{\alpha }\right)=\sqrt{J_2\left(s\left(\sigma \right)\right)}+\eta p\left(\sigma \right)-k$ (18)

For the case of isotropic hardening, the function σ (α), used for the determination of parameter k (given by [⁷] equation) is given by:

$\mathrm{\sigma }\left(\mathrm{\alpha }\right)={\sigma }_y+\mathcal{H}\alpha$ (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, ψ.

$\mathrm{Q}\left(\mathrm{\sigma },{\mathrm{\sigma }}_{\mathrm{y}}\right)=\sqrt{J_2\left(s\left(\sigma \right)\right)}+\stackrel{-}{\eta }p\left(\sigma \right)$ (20)

where ψ < ϕ and $\stackrel{-}{\eta }$ is an additional constant of material.

Therefore, F and Q gradients represented by tensors m and n are given by:

$m_{ij}=\frac{\partial Q}{\partial {\sigma }_{ij}}=\frac{\sqrt{2/3}}{3}\frac{1}{\sqrt{s:s}}\left[\begin{array}{ccc}\sigma & 0& 0\\ 0& -\sigma & 0\\ 0& 0& -\sigma \end{array}\right]+\frac{\stackrel{-}{\eta }}{3}\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]$ (21)

$n_{ij}=\frac{\partial F}{\partial {\sigma }_{ij}}=\frac{\sqrt{2/3}}{3}\frac{1}{\sqrt{s:s}}\left[\begin{array}{ccc}2\sigma & 0& 0\\ 0& -\sigma & 0\\ 0& 0& -\sigma \end{array}\right]+\frac{\eta }{3}\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]$ (22)

The inelastic modulus, obtained from a law associated with hardening or softening phenomenon, is given by:

$H=\left(\frac{\eta }{3}+\frac{1}{\sqrt{3}}\right)\mathcal{H}\mathcal{}\mathcal{}\mathcal{}\mathcal{}\mathcal{}\mathcal{}$ (23)

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 [²⁰] and SOUZA NETO et al. [³]. When the return occurs at the vertex, the yielding function (equation [¹⁸]) and plastic flow potential (equation [²⁰]) should be changed to (SZABÓ; KOSSA [²¹]):

$\mathrm{F}\left(\mathrm{\sigma },\mathrm{\alpha }\right)=\eta p\left(\sigma \right)-k$ (24)

$\mathrm{Q}\left(\mathrm{\sigma },\mathrm{\alpha }\right)=\stackrel{-}{\eta }p\left(\sigma \right)$ (25)

3.2 Ottosen model

The Ottosen criteria is given by the following yielding function

$\mathrm{F}\left(\mathrm{\sigma },\mathrm{\alpha }\right)={\left[\alpha J_2+{\sigma }_c\left(\lambda \sqrt{J_2}+\beta I_1\right)\right]}^{\frac{1}{2}}-{\sigma }_c$ (26)

where in σ_c = (σ_y + Hα).

The derivative of F for an isotropic material may be obtained by chain rule to:

$n_{ij}=\frac{\partial F}{\partial {\sigma }_{ij}}=\frac{\partial F}{\partial I_1}\frac{\partial I_1}{\partial \sigma }+\frac{\partial F}{\partial J_2}\frac{\partial J_2}{\partial \sigma }+\frac{\partial F}{\partial J_3}\frac{\partial J_3}{\partial \sigma }$ (27)

in which the invariant stress derivatives are:

$\frac{\partial I_1}{\partial \sigma }={\delta }_{ij}$ (28a)

$\frac{\partial J_2}{\partial \sigma }=s_{ij}$ (28b)

$\frac{\partial {\mathrm{J}}_3}{\partial \sigma }=t_{ij}=s_{ik}{s}_{kj}-\frac{1}{3}{J}_2{\delta }_{ij}$ (28c)

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:

$\frac{\partial F}{\partial I_1}=\frac{{\sigma }_c\beta }{h}$ (29a)

$\frac{\partial F}{\partial J_2}=\left[\alpha +{\sigma }_c\left(\frac{\lambda }{2\sqrt{J_2}}+\frac{\partial \lambda }{\partial \theta }\frac{\partial \theta }{\partial J_2}\right)\right]\frac{1}{h}$ (29b)

$\frac{\partial F}{\partial {\mathrm{J}}_3}=\frac{\partial \lambda }{\partial \theta }\frac{\partial \theta }{\partial J_3}\frac{\sqrt{J_2}}{h}$ (29c)

where

$h=2{\left[\alpha J_2+{\sigma }_c\left(\lambda \sqrt{J_2}+\beta I_1\right)\right]}^{\frac{1}{2}}$ (30)

$\frac{\partial \lambda }{\partial \theta }=\left\{\begin{array}{c}\frac{-c_1{c}_{2}sen\left(3\theta \right)sen\left(\right({\cos }^{-1}\left(c_2\cos \left(3\theta \right)\right))/3)}{sen({\cos }^{-1}\left(c_2\cos \left(3\theta \right)\right))},if\cos \left(3\theta \right)\ge 0\\ \frac{-c_1{c}_{2}sen\left(3\theta \right)sen\left(\right({\mathrm{\pi cos}}^{-1}\left(c_2\cos \left(3\theta \right)\right))/3)}{sen({\cos }^{-1}\left(c_2\cos \left(3\theta \right)\right))},if\cos \left(3\theta \right)<0\end{array}\right.$ (31)

$\frac{\partial \theta }{\partial J_2}=\frac{3\sqrt{3}{J}_3}{4J_2^{\frac{5}{2}}sen\left(3\theta \right)}$ (32)

$\frac{\partial \theta }{\partial J_3}=-\frac{\sqrt{3}}{2J_2^{\frac{3}{2}}sen\left(3\theta \right)}$ (33)

the model adopted is associative (n_ij = m_ij).

The inelastic module associated with hardening or softening phenomenon, is given by:

$H=\left(1-\frac{\lambda \sqrt{J_2}+\beta I_1}{h}\right)\mathcal{H}$ (34)

4. Numerical examples

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

Figure 5
Geometry and finite element mesh for the diametral compression test

The material parameters are shown in Table [2] and were based on the study made by CECÍLIO [²²].

Table 2
Material parameters

Young’s Modulus (E)	Poisson ratio (ν)	Cohesion (c)	Inelastic modulus (H)	Frictional angle (ϕ)
20,000.0 MPa	0.20	0.9218	1,000 MPa	20°

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:

${\sigma }_t=\frac{2P}{\mathrm{\pi }\mathrm{D}\mathrm{L}}$ (35)

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.

Figure 6
Start of plastic strain

The stress distribution is shown in Figure [7], where the condition of the indirect tensile test is clearly observed.

Figure 7
Distribution of stress σ_x

The tension stress limit in the elastic range, according to equation [35] on the presented conditions, is:

${\sigma }_t=\frac{2\times 60000}{\mathrm{\pi }\times 300\times 150}=0.848MPa$ (36)

The maximum value obtained for the tension shown in Figure [7] is 0.894 MPa, showing excellent agreement with the analytical solution. Figures [8] and [9] show the distributions of normalized stresses along the vertical axis of the model. The graphics of the Figures [8] and [9] emphasize the tension dominant state in the center of the specimen.

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

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

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

Figure 10
Dogbone-shaped panel - Finite element mesh, geometry and image of supports used during the test performed by PEREIRA et al [²³]. Adapted from PEREIRA et al. [⁸]

The material is considered a cement matrix composite reinforced with fibers with 2% of vinyl polyvinyl acetate (PVA), according to PEREIRA et al. [²³]. The material parameters of the experiments are given in Table [3].

Table 3
Material parameter of Dogbone-shape panel

Young’s modulus (E)	Poisson’s ratio (ν)	Compressive strength (σc)	Initial cracking tensile strength (σt)	Tensile strength (σu)	Hardening modulus (H)
20,000.0 MPa	0.20	60 MPa	3.5 MPa	4.0 MPa	1,000 MPa

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 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. [²³].

Figure 11
Numerical results obtained with Ottosen Model compared with experimental results obtained by PEREIRA et al. [⁸]

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-b]. The results are in excellent agreement with the experimental and numerical results presented by PEREIRA et al. [⁸] (Figure [12-c]).

Figure 12
Dogbone-shaped panel - (a) Deformed mesh, (b) displacement field obtained with Ottosen Model and (c) numerical results of PEREIRA et. al [⁸]

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 13
Dogbone-shaped panel - FEM meshes with: (a) 03, (b) 12, (c) 24 and (d) 128 quadrilateral eight-nodes elements

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.

Figure 14
Dogbone-shaped panel - load factor-displacement curves

4.3 Reinforced concrete beam

In this simulation, the results obtained with the elastoplastic Ottosen’s model are compared to experimental tests for a reinforced concrete beam. The experiments were performed by MAZARS and PIJAUDIER-CABOT [25] in a reinforced concrete beam subjected to a three-point bending, showed in Figure [ 15].

Figure 15
Reinforced concrete beam

The Ottosen’s criteria with potential hardening law for concrete was adopted. The materials parameters are given in Table [4].

Table 4
Reinforced concrete beam: material parameters for a concrete

Young’s modulus (E) - Concrete	Poisson’s ratio (ν) - Concrete	Compressive strength (σc)	Initial cracking tensile strength (σt)	Tensile strength (σu)	a	b	k0	n
30,000.0 MPa	0.20	20 MPa	2.0 MPa	2.0 MPa	0.96	2,850.0	n	1.0

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 one-dimensional three-nodes elements to represent the reinforcement.

Figure 16
Finite element mesh

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.

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 [²⁴].

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

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.

5. 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 [¹⁷] and GORI et al. [¹⁸]. 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:

1. Models showed appropriate behavior, and the responses were consistent;

2. Models showed appropriate behavior using the Cutting Plane return mapping algorithm;

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

6. Acknowledgements

The authors acknowledge the financial support in the form of scholarship granted by CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico) - Grants No. 309515/2017-3 and 311663/2017-6 and FAPEMIG (Fundação de Amparo à Pesquisa do Estado de Minas Gerais - Grants No: PPM-00747-18).

7. References

• ^[1]
  CHEN, W.F. and HAN, D.J. Plasticity for Structural Engineers. Ross Publishing Classics, J. Ross Pub, 2007.
• ^[2]
  LUBLINER, J. Plasticity Theory, EUA, Macmillan Publishing Company, 1990.
• ^[3]
  SOUZA NETO, E.A., PERIC, D. and OWEN, D.R.J. Computational Methods for Plasticity: Theory and Application, EUA, Wiley, 2008.
• ^[4]
  CHEN, W.F. and HAN, D.J. Plasticity in Reinforced Concrete. EUA, MacGraw-Hill, 1982.
• ^[5]
  PARK, H and KIM, J.Y. Plasticity Model Using Multiple Failure Criteria for Concrete in Compression. International Journal of Solids and Structures, vol. 42, 2005, p. 2303-2322.
• ^[6]
  GRASSL, P. and JIRÁSEK, M. Damage-plastic Model for Concrete Failure. International Journal of Solids and Structures, vol. 43(22-23), 2006, p. 7166-7196.
• ^[7]
  OTTOSEN, N.S. Nonlinear Finite Element Analysis of Concrete Structures. Technical Report, Riso National Laboratory, 1980.
• ^[8]
  PEREIRA, E.B., BARROS, J.A.O. and FISHER, G. 3-D 4-Parameter Elastoplasticity Model for Strain Hardening Cementitous Composites, Technical Report, 2012.
• ^[9]
  OTTOSEN, N.S. A Failure Criterion of Concrete. Journal of the Engineering Mechanics Division, vol.103, 1977, p. 527-535.
• ^[10]
  CEB. Ceb-fip Model Code 2010 - Final Draft. Bulletin D’Information, 2010.
• ^[11]
  DAHL, K.K.B. A Failure Criterion for Normal and High Strength Concrete. 1992.
• ^[12]
  NETO, J.M. Um Estudo da Formulação de Modelos Constitutivos Viscoelásticos e Elasto-Viscoplásticos e do Emprego de Algoritmos Implícitos e Explícitos Para a Sua Integração Numérica, São Carlos, 1998, Tese (Doutorado), Universidade de São Paulo.
• ^[13]
  MALAVOLTA, A.T. Metodologia Para Determinação de Parâmetros Utilizados em Uma Nova Superfície de Escoamento Anisotrópica para Processos de Conformação de Chapas Metálicas, São Carlos, 2008, Tese (Doutorado), Universidade de São Paulo.
• ^[14]
  BOUCHARD, P. O, BOURGEON, L., FAYOLLE, S. and MOCELLIN, K. An enhanced Lemaitre model formulation for materials processing damage computation. Int. J. Mater Form, vol 4, 2011, p. 299-315.
• ^[15]
  VAZ JR, M. and ROJAS, P. M. Damage evolution and thermal coupled effects in inelastic solids. International Journal of Mechanical Sciences, vol. 53, 2011, p. 387-398.
• ^[16]
  TAQIEDDIN, Z.N. Elasto-Plastic and Damage Modeling of Reinforced Concrete. Louisiana, 2008, Tese (Doutorado), Louisiana State University.
• ^[17]
  PENNA, S.S. Formulação Multipotencial para Modelos de Degradação Elástica: Unificação Teórica, Proposta de Novo Modelo, Implementação Computacional e Modelagem e Estruturas de Concreto. Belo Horizonte, 2011, Tese (Doutorado), Universidade Federal de Minas Gerais.
• ^[18]
  GORI, L., PENNA, S. S., PITANGUEIRA, R. L. S. A computational framework for constitutive modelling. Computer& Structures, vol. 187, 2017, p. 1-23.
• ^[19]
  CAROL, I., RIZZI, E. and WILLIAM, K. A Unified Theory of Elastic Degradation an Damage Based on a Loading Surface. International Journal of Solid and Structures, vol. 31(20), 1994, p. 2835-2865.
• ^[20]
  SIMO, J.C. and HUGHES, T.J.R. Computational Inelasticity, EUA, Springer-Verlag, 1998.
• ^[21]
  SZABÓ, L. and KOSSA, A. A New Exact Integration Method for the Drucker-Prager Elastoplastic Model with Linear Isotropic Hardening. International Journal of Solids and Structures, vol. 49, 2012, p.170-190.
• ^[22]
  CECILIO, D.L., Modelagem e Simulação Elastoplástica em Elementos Finitos, Campinas, 2011, Dissertação (Mestrado), Universidade Estadual de Campinas.
• ^[23]
  PEREIRA, E.B., BARROS, J.A.O. and FISHER, G. Direct Assessment of Tensile Stress Crack Opening Behavior of Strain Hardening Cementitious Composites (shcc). Cement and Concrete Research, vol. 42, 2012, p. 834-846.
• ^[24]
  MAZARS, J. and PIJAUDIER-CABOT, G. Continuum Damage Theory - Application to Concrete. Journal of Engineering Mechanins ASCE, vol. 115, 1989, p. 345-365.

Publication Dates

• Publication in this collection
  27 Mar 2020
• Date of issue
  Jan-Feb 2020

History

• Received
  10 Dec 2018
• Accepted
  11 Sept 2019