Customization of a software of finite elements to analysis of concrete structures: long-term effects

Abstract

Concrete is a material that presents time-delayed effects associated with creep and shrinkage phenomena. The numerical modeling of the creep is not trivial because it depends on the age of the concrete at the time of loading, which makes necessary to store the history of loads to apply the principle of superposition. However, a problem formulated in finite elements has many points of integration, making this storage impossible. The Theory of Solidification, together with incremental algorithm developed by Bazant e Prasannan ([¹], [²]), solves this problem. Therefore, this paper presents the adaptation of this algorithm to a finite element commercial software (ANSYS) considering the creep, according to the CEB-FIP Model Code 1990 [³]. The results demonstrate that the deformations provided by this adaptation are in accordance with the analytical solution given by the CEB-FIP MC90, including cases where the loads are applied at different ages of the concrete.

Keywords:
concrete long-term effects; creep and shrinkage on concrete; finite element method

Resumo

O concreto é um material que apresenta efeitos diferidos no tempo, que estão associados aos fenômenos de fluência e retração. A modelagem numérica da fluência não é trivial, pois depende da idade do concreto no instante de aplicação do carregamento, o que torna necessário armazenar a história de carga para posterior aplicação do princípio da superposição dos efeitos. Contudo, um problema formulado em elementos finitos possui muitos pontos de integração, tornando esse armazenamento impossível. A Teoria da Solidificação, juntamente com o algoritmo incremental desenvolvido por Bazant e Prasannan ([¹], [²]), permitem resolver esse problema. Esse artigo apresenta a adaptação desse algoritmo em um programa de elementos finitos (ANSYS), considerando a fluência e, também a retração, conforme formulado pelo Código Modelo CEB-FIP 1990 [³]. Os resultados comprovam que as deformações previstas por esta adaptação estão em conformidade com a solução analítica dada pelo CEB-FIP MC90, incluindo situações em que as cargas são aplicadas ao longo de diferentes idades do concreto.

Palavras-chave:
efeitos diferidos do concreto; fluência e retração do concreto; método dos elementos finitos

1. Introduction

Since it was discovered by Hatt [⁴] in 1907, several materials researchers have been studying long-term effects of concrete. Currently, this behavior is separated in creep and shrinkage, and the main difference between them is that the first one depends on the loading, while the second does not.

Creep is characterized by continuous and gradual advancement of deformations under constant stress. This deformation (excluding the Poisson effect) has the same direction of the loading and is divided, for analytical convenience, into basic and drying creep. The first occurs without water exchange with the environment and can be measured when the specimen is immersed in a 100% humidity ambience. On the other hand, the second is when water exchanged occurs with the environment and therefore depends on the relative humidity of the air and the exposure of the specimen to the ambience.

The shrinkage is the reduction of the volume due to the gradual loss of water independent of the tension. It is also divided, for convenience, into autogenous and drying shrinkage. The first occurs without loss of water to the environment and is a consequence of the removal of water from the capillary pores by the hydration reactions of the cement, while the second one, also called hydraulic shrinkage, occurs due to the loss of water to the environment. As the shrinkage is associated with the volume reduction, this phenomenon induces internal tensile stresses that can cause fissures on the concrete. In this aspect, the minimization of shrinkage is one of the main reasons for the concrete curing procedure, which seeks to reduce water loss in the early ages as the concrete is just developing its resistance.

These deformations of the concrete have important influence on the structural behavior, since these deformations can reach 2 to 3 times the value of the instantaneous deformation, according to RILEM Technical Committees [⁵]. Several semi-empirical formulations for predicting creep and shrinkage strain are available in the literature and in design standards. However, these formulations, in general, present a certain approximation when compared to experimental tests. As noted by Fanourakis e Ballim [⁶], the error can range from 25.9% until 67.4%, depending on the formulation and test used in the comparison.

This paper is based on the works of Dias [⁷], Dias et al. [⁸], Quevedo [⁹] and Schmitz [¹⁰], who used the long-term concrete formulation proposed by the Comité Euro-International du Béton (CEB-FIP MC90) [³] adapted to the Solidification theory of the Bazant and Prasannan ([¹], [²]) for their studies involving composite beams, tunnels and bridges, respectively. The creep formulation of the CEB-FIP MC90 has the advantage of fitting perfectly with the Bazant and Prasannan theory. This is possible because this formulation separates the coefficient of creep in a factor that depends on the age of the concrete and another one that depends exclusively on the age of the load. Such separation does not occur in current formulations, for example, in the Féderation Internationale du Béton [¹¹]. In view of this decomposition, the incremental algorithm developed by Bazant and Prasannan [²] allows to model efficiently the cases in which there are histories of variable loads (or stresses). The use of the CEB-FIP MC90 formulation without this algorithm, in a non-linear analysis, requires that all increase and decrease of stresses need stored to apply the superposition principle. This procedure would make it impossible to analyze in situations with many finite elements.

Therefore, this paper will present this adaptation with some comments on the customization and use of this model in ANSYS. At the end, some results are presented to validate this approach.

2. Numerical model

In this item, the model for the concrete is presented considering the long-term behavior, especially the creep and the adaptation of the CEB-FIP MC90 expressions to the Solidification theory.

2.1 Adaptation of the CEB-FIP MC90 creep formulation to the Solidification theory

Considering that the concrete is an aging material and creep is the portion of the long-term behavior directly related to the loads applied to the structure, the age of the concrete at the times of application of the loads directly influences the creep deformations. The theory of Solidification proposed by Bazant and Prasannan ([¹], [²]) aims to facilitate the numerical solution of this particularity between creep and the age of loading. These authors proposed a model that considers aging as an isolated factor, associated only to the volume of solidified concrete, overtime.

Figure 1 shows in a schematic form the model proposed by Bazant and Prasannan for the concrete behavior. The total deformation is given by the sum of four parts: elastic, viscoelastic, viscous and a part that is independent of the stress, corresponding to the deformations by shrinkage, thermal and cracking. The creep strain is divided into two parts: one called “viscoelastic” and another called “viscous”. This division intends to treat the shape of the curve creep in relation to the loading age. For small loading ages, the function takes the form of power (primary creep), and for larger loading ages, the curve takes the predominantly logarithmic form (which is often conveniently represented as a straight line on the semi-log scale of time) (Bazant e Prasannan, [¹]).

Figure 1
Concret model (adapted from Bazant e Prasannan [¹], [²])

As shows in Figure 1, the viscoelastic portion ε_v is associated to the volume fraction of concrete already solidified v(t) and with a loading age depending of the creep coefficient Ф(t - t₀), where time t is the current age of the concrete and t₀ is the age at which the load was applied. Similarly, the viscous portion ε_f depends on the hydrated cement fraction h(t) and the coefficient ψ(t - t₀) dependent on the loading age (Bazant e Prasannan, [¹]). In Figure 1, the factor Ф(t - t₀) is represented by a Kelvin-Generalized rheological model (or also known as the Kelvin chain). Since this function depends only on the age of the load, the parameters of the Kelvin chain become constant, facilitating the development of the incremental algorithm described later.

From these considerations the creep function J(t, t₀), according to the theory of Solidification is given by expression (1):

$J\left(t,t_0\right)=\frac{1}{E_0}+\frac{\gamma \left(t-t_0\right)}{v\left(t\right)}+\frac{1}{\eta \left(t\right)}$ (1)

Where E₀ is the modulus of elasticity for the aggregates and microscopic particles of the cement paste, γ(t - t₀) is the microviscoelastic deformation, and η(t) is the apparent macroscopic viscosity dependent on the age of the concrete. In the CEB-FIP MC90 [³] the creep function is given by (2):

$J\left(t,t_0\right)=\frac{1}{E_c\left(t_0\right)}+\frac{\text{Ф}\left(t,t_0\right)}{E_{ci}}$ (2)

Where E_c (t₀) is the tangent modulus of elasticity at t₀, E_ci is the tangent modulus of elasticity at 28 days and Ф(t, t₀) is the creep coefficient. In [³] the coefficient of creep is divided into two factors: ϕ₀ which is related to the age of load application and β_c (t - t₀) which describes the shape of the creep curve with the age of the load. Comparing (1) and (2) we have expressions (3), (4), (5) and (6) to adapt both formulations.

$E_0=E_c\left(t_0\right)$ (3)

$\gamma \left(t-t_0\right)={\beta }_c\left(t-t_0\right)$ (4)

$\frac{1}{v\left(t\right)}=\frac{{\text{Ф}}_0}{E_c\left(t_0\right)}$ (5)

$\frac{1}{\eta \left(t\right)}=0$ (6)

Equality (6) occurs because Bazant and Prasannan ([¹], [²]) treat the creep curve separating it into two parts: one corresponding to the first ages (viscoelastic) and another referring to posterior ages (viscous). However, the CEB-FIP MC90 formulation does not make this distinction by representing the behavior of the creep curve only by the function β_c (t - t₀).

Although the Bazant and Prasannan Solidification theory ([¹], [²]) consider a factor to the stress level F[σ(t)], the model presented in CEB-FIP MC90 does not consider this behavior and, therefore, this factor does not appear in this adaptation. However, due to this simplification, the stress limit of 40% of the average compressive strength is imposed. This limit is fundamental to keep the concrete in viscoelastic regime, where the principle of superposition of the effects is valid. Therefore, during the numerical analysis, considering this adaptation, one must respect this limit in order do not underestimate the deformations.

2.2 Curve fitting to the coefficient of creep factor that is independent of the concrete age

As previously mentioned a Kelvin chain represents the coefficient of creep that is independent of the concrete age in the incremental algorithm of Bazant and Prasannan [¹]. Since this factor is independent of the concrete age t₀, the parameters of this chain (modulus of elasticity of the springs and coefficient of viscosity of the dashpots) also become constant and this facilitates the deduction of the incremental algorithm, which uses a numerical solution for the integral of the strain rate.

The equations that represent the behavior of the Kelvin chain are presented in expressions (7). Due to the serial configuration of the μ-th units of Kelvin, for each unit, the applied tension σ_μ is the same as the total σ and the total strain of the chain γ is given by the sum of the deformation of the units.

$\sigma =E_{\mu }{\gamma }_{\mu }+{\eta }_{\mu }{\dot{\gamma }}_{\mu }\gamma ={\displaystyle \sum }_{\mu =1}^N{\gamma }_{\mu }$ (7)

Where E_μ is the modulus of elasticity, η_μ is the viscosity and γ_μ is the deformation of μ-th unit of the chain. Solving the differential equation presented in (7) we get the expression (8) which is known as the Dirichlet (or Prony) series. The term τ_μ is called the retardation time of the μ-th unit of the chain.

$\gamma \left(t-t_0\right)={\displaystyle \sum }_{\mu =1}^N\frac{1}{E_{\mu }}\left(1-e^{-\left(t-t_0\right)/{\tau }_{\mu }}\right);{\tau }_{\mu }=\frac{{\eta }_{\mu }}{E_{\mu }}$ (8)

In order to use the Solidification theory, adapted to CEB-FIP MC90, the expression (8) must represent the function β_c (t - t₀), as shown in (4). This can be done by discretizing this function into L points and fitting it to (8) through the least squares method. According to Dias [⁷] this can be done by assembling and solving the system [A]{X}=[B], whose elements are given by (9), (10) and (11).

$A\left(i,j\right)={\displaystyle \sum }_{k=1}^L\left[1-e^{-\left(t_k-t_0\right)/{\tau }_i}\right]\left[1-e^{-\left(t_k-t_0\right)/{\tau }_j}\right]$ (9)

$B\left(i\right)={\displaystyle \sum }_{k=1}^L{\beta }_c\left(t_k-t_0\right)\left[1-e^{-\left(t_k-t_0\right)/{\tau }_i}\right]$ (10)

$X\left(j\right)=1/E_j$ (11)

Through (11) the moduli of elasticity of the Kelvin chain are obtained. It is important to emphasize that the fitting of the model can lead to negative values of modulus of elasticity in certain elements, which do not have physical correspondence with the mechanical model, because it is only a mathematical adjustment. In addition, the resolution of this linear system requires the choice of the retardation times of the elements of the chain τ_μ. In order for the matrix [A] to be well conditioned, Bazant and Prasannan [²] propose that the retardation times are chosen as a logarithmic increase of base 10 according to (12).

${\tau }_{\mu }={\tau }_1\cdot {10}^{\mu -1};\mu =\mathrm{1,2,}\dots ,N$ (12)

According to Dias [⁷], the retardation times should cover half the analysis period t_max and τ₁ be small enough to represent the behavior of the creep curve at early ages, as expressed in (13).

${\tau }_1=0.01t_0;{\tau }_N=0.5t_{m\text{a}x}$ (13)

According to Dias [⁷], the quantity N of Kelvin units was determined to be the maximum possible until one has τ_μ greater than or equal to τ_N. However, no more than six Kelvin units.

In addition to the retardation times, for the system solution, it is also necessary to choose the L points of the function β_c (t - t₀) for the adjustment. Since this function is related to the age of the load, in fact, the question is to define the age intervals (t_k - t₀), where t₀ is the age of the concrete at the time of the first load application and k=1,…,L are the points. Considering that the behavior of the function is more pronounced in the first loading ages, an effective choice is to keep the time intervals constant in a logarithmic scale. Thus, the points for adjustment can be obtained through the recurrence formula (14) (Bazant e Prasannan [¹], [²]).

$\left(t_k-t_0\right)={10}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$m$}\right.}\left(t_{k-1}-t_0\right);k=\mathrm{2,}\dots ,L$ (14)

So in that, the parameter m is the number of steps per decade, according to [⁷], and around 10 for a good precision. In addition, the first time interval chosen corresponds to one tenth of the first loading age of the structure [⁷], according to expression (15).

$\left(t_1-t_0\right)=0.1t_0$ (15)

And the number of points L can be determined to be that necessary to cover study time, in other words, until t_(k+1) - t₀ is greater than t_max - t₀.

2.3 Incremental algorithm

Then, the adjusted parameters (E_μ, τ_μ) of the Kelvin chain can be used in the incremental algorithm. Expressing the total deformation in terms of rates, we have (16).

$\dot{\epsilon }=\dot{{\epsilon }_e}+\dot{{\epsilon }_v}+\dot{{\epsilon }_0}$ (16)

Where $\dot{{\epsilon }_e}=\dot{\sigma }\left(t\right)/E\left(t\right)$ is the elastic deformation rate, $\dot{{\epsilon }_v}=\dot{\gamma }(t-t_0)/v\left(t\right)$ is the viscous strain rate and $\dot{{\epsilon }_0}$ is the sum of all strain rates independent of stress (shrinkage, thermal and cracking).

Integrating expression (16) through the trapezoid rule, for an integration interval from t_i to t_(i+1) we obtain (17), where E_(i+1/2) e v_(i+1/2) are the values corresponding to the middle of the interval between t_i and t_(i+1).

$\Delta \epsilon =\frac{\Delta \sigma }{E_{i+1/2}}+\frac{\Delta \gamma }{v_{i+1/2}}+\Delta {\epsilon }_0$ (17)

If the increment of elastic deformation is separated from the remainder of the deformation, we can write (18), the part ∆ε* being all deformation which is non-elastic (creep, shrinkage, thermal and cracking) and E* the equivalent modulus of elasticity, both to be deducted.

$∆\epsilon =\frac{∆\sigma }{E^{\mathrm{*}}}+{∆\epsilon }^{\mathrm{*}}$ (18)

The creep strain is obtained by assuming that the stress varies linearly between two consecutive time steps, according (19). This assumption, as will be seen later, imposes the need for a small increment of time.

$\sigma \left(t\right)={\sigma }_i+\left(t-t_i\right)∆\sigma /∆t$ (19)

By introducing (19) into (7), we have (20).

$\begin{array}{l}{\gamma }_{{\mu }_{i+1}}={\gamma }_{{\mu }_i}{e}^{-\Delta y_{\mu }}+\frac{{\sigma }_i}{E_{\mu }}\cdot \left(1-e^{-\Delta y_{\mu }}\right)+\frac{1-{\lambda }_{\mu }}{E_{\mu }}\cdot \Delta \sigma \\ \Delta y_{\mu }=\frac{\Delta t}{{\tau }_{\mu }};{\lambda }_{\mu }=\frac{\left(1-e^{-\Delta y_{\mu }}\right)}{\Delta y_{\mu }}\end{array}$ (20)

Since ${\gamma }_{{\mu }_{i+1}}$ is the strain of the μ-th Kelvin unit in the next time step (i + 1) and ${\gamma }_{{\mu }_i}$ is the μ-th Kelvin unit strain for the current time step i. In the same way σ_i is the current stress and ∆σ is the stress increment considering the time interval between t_i and t_i+1. This increase of deformations can be calculated as ${∆\gamma }_{\mu }={\gamma }_{{\mu }_{i+1}}-{\gamma }_{{\mu }_i}$, resulting in (21).

${∆\gamma }_{\mu }=\left(\frac{{\sigma }_i}{E_{\mu }}-{\gamma }_{{\mu }_i}\right)\bullet \left(1-e^{-{∆y}_{\mu }}\right)+\frac{1-{\lambda }_{\mu }}{E_{\mu }}\bullet ∆\sigma$ (21)

The first factor of the first term to the right side of the equality in (21) is the viscous deformation of the μ-th Kelvin unit of step i, expressed by (22).

${\epsilon }_{{\mu }_i}^{\mathrm{*}}=\frac{{\sigma }_i}{E_{\mu }}-{\gamma }_{{\mu }_i}$ (22)

In the same way as expression (22), in (23) we have the viscous deformation of the Kelvin unit of the next step i +1.

${\epsilon }_{{\mu }_{i+1}}^{\mathrm{*}}=\frac{{\sigma }_{i+1}}{E_{\mu }}-{\gamma }_{{\mu }_i+1}$ (23)

As deduced by Dias [⁷], by substituting ${\gamma }_{{\mu }_i+1}$ given by equation (20), in the expression (23) results the expression (24), which is still ignoring the viscous deformation aging. The aging is considered including the volume of solidified concrete in expression (24), so that the viscous deformation is given by equation (25).

${\epsilon }_{{\mu }_{i+1}}^{\mathrm{*}}=\frac{{\lambda }_{\mu }}{E_{\mu }}\bullet ∆\sigma +{\epsilon }_{{\mu }_i}^{\mathrm{*}}\bullet e^{-{∆y}_{\mu }}$ (24)

${\epsilon }_{{\mu }_{i+1}}^{\mathrm{*}}=\frac{{\lambda }_{\mu }}{E_{\mu }{v}_{i+1/2}}\bullet ∆\sigma +{\epsilon }_{{\mu }_i}^{\mathrm{*}}\bullet e^{-{∆y}_{\mu }}$ (25)

For the increment of viscous deformation, first replace the expression (21) in (17) having the increment of total strain, given by the expression (26).

$∆\epsilon =∆\sigma \bullet \left(\frac{1}{E_{i+1/2}}+\sum _{\mu =1}^N\frac{1-{\lambda }_{\mu }}{E_{\mu }\bullet V_{i+1/2}}\right)+\left[\sum _{\mu =1}^N{\epsilon }_{{\mu }_i}^{\mathrm{*}}\bullet \left(1-e^{-{∆y}_{\mu }}\right)\right]+{∆\epsilon }_0$ (26)

Considering the separation of the elastic deformation part proposed by equation (18) we can define the expression for ∆ε* and E*, which are according to equations (27) and (28).

${∆\epsilon }^{\mathrm{*}}=\left[\sum _{\mu =1}^N{\epsilon }_{{\mu }_i}^{\mathrm{*}}\bullet \left(1-e^{-{∆y}_{\mu }}\right)\right]+{∆\epsilon }_0$ (27)

$\frac{1}{E^{\mathrm{*}}}=\left(\frac{1}{E_{i+1/2}}+{\sum }_{\mu =1}^N\frac{1-{\lambda }_{\mu }}{E_{\mu }\bullet v_{i+1/2}}\right)$ (28)

Generalizing the expressions for a three-dimensional domain, considering an isotropic material, we have (27) that become a matrix equation, according (29).

$\left\{{∆\epsilon }^{\mathrm{*}}\right\}=\sum _{\mu =1}^N\left\{{\epsilon }_{{\mu }_i}^{\mathrm{*}}\right\}\left(1-e^{{-∆y}_{\mu }}\right)+\left\{{∆\epsilon }_0\right\}$ (29)

Where $\left\{{\epsilon }_{{\mu }_i}^{*}\right\}$ is calculated by the expression (30), generalized of the (25), considering which [D_μ]^-1 is the inverse of the linear isotropic constitutive matrix, and the elastic modulus is equal to $E_{\mu }{v}_{i+1/2}/{\lambda }_{\mu }$.

$\left\{{\epsilon }_{{\mu }_i}^{\mathrm{*}}\right\}={\left[D_{\mu }\right]}^{-1}\bullet \left\{∆\sigma \right\}+\left\{{\epsilon }_{{\mu }_{i-1}}^{\mathrm{*}}\right\}\bullet e^{-{∆y}_{\mu }}$ (30)

In expression (30), it can be verified that it is necessary to store only the viscous deformation corresponding to the previous step $\left\{{\epsilon }_{{\mu }_{i-1}}^{\mathrm{*}}\right\}$, avoiding to store all the loading history.

During the iterative solution process, in addition to determining the increment of non-elastic deformation, it is also necessary to update the constitutive matrix and correct the increment of total deformation in the calculation of the stress increase. Generalizing to a three-dimensional domain the expression (18) and isolating the stress increment is obtained (31).

$\left\{∆\sigma \right\}=\left[D^{\mathrm{*}}\right](\left\{∆\epsilon \right\}-\left\{{∆\epsilon }^{\mathrm{*}}\right\})$ (31)

Where [D*] is a linear isotropic constitutive matrix considering the modulus of elasticity E* calculated by the expression (28) and Poisson’s ratio of the concrete.

3. ANSYS customization

3.1 Operation of programmable tool UserMat

ANSYS offers a series of programmable features (User Programmable Features - UPF) that allows the user customize models aspects through subroutines wrote in Fortran77. Once the change is made, the subroutines are compiled and associated to the main program ANSYS. One of these subroutines, the UserMat, is specific to customize the material laws; therefore, it is ideal to insert the concrete model.

The iterative process used by ANSYS to solve nonlinear problems consists to apply the Newton-Raphson procedure. For each time/load increment n (also called substep) it is calculated i equilibrium iterations until the convergence criterion is satisfied. The expressions (32), (33) and the Figure 2 (adapted of ANSYS [¹³]) explain this process.

Figure 2
Illustration of Newton-Raphson Method with equilibrium iterations i for each substep n (ANSYS, [¹³])

${\left\{\mathrm{u}\right\}}_{n,i+1}={\left\{\mathrm{u}\right\}}_{n,i}+{\left\{∆\mathrm{u}\right\}}_{n,i}$ (32)

${\left\{∆\mathrm{u}\right\}}_{n,i}={{\left[K^T\right]}_{n,i}}^{-1}\left({\left\{{\mathrm{F}}^a\right\}}_{n+1}-{\left\{{\mathrm{F}}^{nr}\right\}}_{n,i}\right)={{\left[K^T\right]}_{n,i}}^{-1}{\left\{R\right\}}_{n,i}$ (33)

Where {u}_n,i+1 is the actualized nodes displacement vector of the substep n, {u}_n,i is the nodes displacement vector of the iteration i, {∆u}_n,i is the nodes displacement increment vector of the iteration i, [K^T]_n,i is the tangent stiffness matrix (or Jacobian matrix of Newton-Raphson method), {F^a}_n+1 is the applied forces vector, {F^nr}_n,i is the internal forces vector of the iteration i, {R}_n,i is the unbalance vector or residual of the iteration i used in the convergence criterion. Using the nodes displacement, in each equilibrium iteration i, it is obtained the deformations and stress in each integration point through the relations displacement-strain and stress-strain formulated by shape functions of the Finite Element Method. Therefore, for each integration point, the main program gives to UserMat the total stress, total deformations and increment of total deformations (resultants of the Newton-Raphson iterative process) in the current time increment. The instructions programmed inside of the subroutines are responsible to collect the material constants (defined by the user), to update the stress and the Jacobian matrix used in the iterative process. The Figure 3 shows a flow chart that summarize the UserMat operation inside the nonlinear analysis.

Figure 3
UserMat operation in the nonlinear solution process

3.2 Implementation of viscoelastic model in UserMat

ANSYS provides one UserMat with a bilinear isotropic relation of stress-strain with the von Mises plasticity criterion. The user needs understand this subroutine and modify according the necessity. It is important to understand all input and output arguments of the subroutine for creating and using the less number of variables that is possible in the material behavior customization. Special attention should be given to UserMat input and output variables, since they must not be deleted. However, the local variables, which in the most part are referenced to the von Mises plastic behavior, can be excluded without compromise the subroutine operation.

The UserMat consists of four internal subroutines depending on the strain components number (or stress) and of the dimensions involved in the problem (1D, 2D, 3D). In the present paper was customized the more generic subroutine of the four, the UserMat3D, that corresponds to problems with more than 3 components, as three-dimensional problems (6 components, 3 directions), plane strain and axisymmetric (both with 4 stress components and 2 directions).

The user is free to create how many local variables are necessary. However, local variables lose their value when the subroutine finish during the equilibrium iterations. When it is necessary holding the value of some variables between the substeps it is necessary to store in array of state variables ustatev (Figure 4-(d)). This array can be used directly without creating a variable. However, in order that the ustatev array is dimensioned it is necessary to declare its size in the script through the command TB, STATE during the material assignment (Figure 4-(a)). It is not recommended save the variables values in a COMMON instruction, except the ones that maintain constant values during the analysis. This type of declaration brings problems in the ANSYS parallelization, because many parallel processes can access and overwrite this memory space simultaneously. It is also possible to visualize the ustatev array values in the post processor through the command PLSOL,SVAR, [variable position in the array], provided that they are saved through the command OUTRES,SVAR,ALL.

Figure 4
APDL commands for (a) material definition, (b) definition of time and analysis time increment. Whitin the Usermat3D: (c) assignment of the variables related to the material and (d) storing principal stress variables in the ustatev array

The input data array prop contain the material constants that are completed through the script during the material definitions. This array is dimensioned through the command TB, USER (that define the quantity nprop of proprieties) and the proprieties values are attributed through the command TB,DATA (Figure 4-(a), Figure 4-(c)).

Other important aspect about the use of UserMat in concrete viscoelasticity is the definition of times and concrete ages during the analysis. The ANSYS provides the variable time and dtime that are shared between script and the UserMat and represent the analysis current time and the time increment, respectively (Figure 4-(b)).

In general way the customized UserMat3D follows the steps:

• 1) calculation of the non-elastic strain increment {∆ε*};

• 2) calculation of the elastic strain increment: {∆ε_e} = {∆ε} - {∆ε*}, where {∆ε} is the total strain increment informed by ANSYS (comes from the Newton-Raphson iterative process);

• 3) calculation of the constitutive matrix [D];

• 4) calculation of the stress increment {∆σ} = [D] {∆ε_e};

• 5) update of the Jacobian matrix (∂∆σ_ij) ⁄ (∂ε_ij) = [D];

• 6) calculation of the update stress {σ_(i+1)} = {σ_i} + {∆σ}.

When it is decreased the non-elastic strain increment of the total strain increment (step 2), the Newton-Raphson procedure will compensate this increasing the total strain increment during the equilibrium iterations until the convergence will be satisfacted. In this way it is introduced the concrete viscous deformation in the problem.

In the customized UserMat, the concrete is defined with an initial instant t_i, option that can be useful when it is used the birth/death recurs of the elements to simulate the construction procedure of a bridge or a tunnel, for example. The material age is given by t_mat = time - t_i. However, the phenomena of shrinkage and creep depend of, respectively, the age of the beginning of shrinkage t_s and the time when the load is applied t₀. In this way, the phenomenon of shrinkage starts when time > t_i, that is at the moment of the concrete birth, but considering t_s, according to the CEB-FIP MC90 expressions. The creep phenomenon starts when time > t₀, that is at the moment that the piece gets loaded, where the load age is equal to t_mat - t₀.

In Quevedo [⁹], there is the customized UserMat subroutine for the viscoelastic concrete, with creep and shrinkage phenomena) and there are instructions to add it in the ANSYS main program.

3.3 Compatible finite elements

The UserMat3D can be used in three-dimensional generic problems, plain strain or axisymmetric. For the first case can be used the elements SOLID185 and SOLID186, for the others two situations there are the elements PLANE182 and PLANE183.

The elements PLANE182 and PLANE183, presented in Figure 5, are elements quite similar, both are quadrilaterals and they have with two degrees of freedom at each node (translations in the two plane directions where are inserted). The main difference between them is that the PLANE182 has four nodes and, therefore, linear shape functions while the PLANE183 has eight nodes and, therefore, quadratic shape functions. Besides the quadrilateral shape, both encompass the triangular shape, being this shape not advisable for PLANE182 because it presents little precision on the interpolation.

Figure 5
Plane elements: (a) PLANE182 and (b) PLANE183 (ANSYS, [¹⁴])

The elements SOLID185 and SOLID186, illustrated in Figure 6, are solid elements; therefore, their nodes have three degrees of freedom of displacements in the directions x, y e z. So again the main difference between the two elements is the shape functions degree, the SOLID185 has eight nodes and, therefore, linear shape functions while the SOLID186 has twenty nodes and, therefore, quadratic shape functions. Both elements encompass also the tetrahedral and pyramidal shape.

Figure 6
Three-dimensional elements: (a) SOLID185 and (b) SOLID186 (ANSYS, [¹⁴])

The choice of quadratic elements should consider the number of elements in the mesh. Due to a nonlinear analysis, where the system of equations is solved at each increment of time, the processing time may be large.

For reinforced concrete analysis, both classes of elements (PLANE and SOLID) have the incorporated reinforcement feature through the inclusion of the elements REINF263 (for plane elements) and REINF264 (for solid elements) after setting of the concrete mesh elements. These elements are ideals to represent reinforcement, because they have only axial stiffness and they can be inserted in any position and direction inside the concrete elements. This aspect brings a large advantage in comparison with the ANSYS concrete element (SOLID65) that the reinforcement is not incorporated and it needs a much more refined mesh.

4. Validation and numerical examples

In order to validate the programmed model in UserMat, it was done some numerical simulations, being that the first two tests served to compare the results with the analytical solution obtained with CEB-FIP MC90 formulations. These tests consisted in modeling the compression experiment, as in Figure 7. The proprieties considered in these initial simulations are presented in Table 1. The first test consisted to apply compression stress increments equal to 5 MPa at 10, 50 and 75 days, being the analysis made until 100 days. The second test consisted to apply 5 MPa constant until 100 days. The Figures 8 and 9 present the obtained results decomposed in total strain, elastic, creep and shrinkage. On the two tests, it was verified perfect agreement with the numerical simulations and analytical calculation. These tests were also repeated considering one element and various elements.

Figure 7
Model to verify the implementation of the viscoelastic model

Figure 8
Test 1 with cube subjected to variable stress using the element SOLID185

Figure 9
Test 2 with cube subjected to variable stress using the element SOLID185

Table 1
Concrete parameters in the initial tests

Parameter	Symbol	Unit	Value
Characteristic compressive strength	fck	MPa	40
Coefficient which depends on concrete type	s	adm	0.25
Poisson’s ratio	ν	adm	0.2
Relative humidity of ambient environmental	RH	%	70
Fictitious thickness	hf	cm	54.54
Age at the beginning of shrinkage	ts	days	7
Coefficient which depends on concrete type - shrinkage	βsc	adm	5
Temperature	temp	°C	20
Coefficient which depends on concrete type	α	adm	1

The next test consisted in check the increment time influence used in the analysis. Until now, the time increment considered was only 1 day. However, the incremental algorithm of Solidification theory suppose a stress linear variation between two consecutive time intervals. However, in fact, this variation is not linear, and this way, how bigger the time increment and the time of whole analysis, bigger will be the error in the creep deformation calculation. In Figure 10 are presented the results of this analysis considering only the elastic strain added the creep strain for 5 time increments: 1, 2, 5, 10 and 20 days. It can be observed that with the increment of 5 days, the error is little, around 2.5%. However, for the increment of 10 days, the error reach 4.8% and for the 20 days, it is 7.5%. Therefore, to have better precision on the results it is suggested to use time increments minors than 5 days.

Figure 10
Influence of time increment

In order to compare the results with experimental data, it was numerically simulated 5 experimental tests presented by Ross [¹²]. These experimental tests consist in two cylindrical proof bodies, with diameter equal to 117.5mm and 305mm of height. One of them subjected to load history (Table 2) and the other without loads, to measure the deformation generated by shrinkage, and after to decrease it of the results of proof bodies subjected to load, having in this way only the elastic response with creep. The concrete proprieties in the Ross’ experiments are in Table 3. To the numerical simulations, it was tested the two elements, SOLID185 and PLANE182 in axial compression, like in Figure 7. The obtained results are in Figure 11 until Figure 15. It can be observed the perfect agreement between the model with SOLID185 and the CEB-FIP MC90 formulation, being that the two curves are overlapping. About the PLANE182, it realized that in all situations there is a little error. This error is attributed to the fact that tests are not in a plain strain or axisymmetric conditions, that are the scope of the element. Despite this, it can be confirmed that the UserMat3D is not restrict only to three-dimensional problems.

Figure 11
Comparison with experimental test 1

Figure 12
Comparison with experimental test 2

Figure 13
Comparison with experimental test 3

Figure 14
Comparison with experimental test 4

Figure 15
Comparison with experimental test 5

Table 2
Stress history in Ross’ tests [¹²]

Tests		t0	t1	t2	t3	t4	t5	tf
1	Age (days)	14	60	-	-	-	-	140
	Stress (MPa)	1.503	0	-	-	-	-	0
2	Age (days)	28	60	91	120	154	-	190
	Stress (MPa)	1.503	1.127	0.751	0.376	0	-	0
3	Age (days)	8	14	28	63	90	120	180
	Stress (MPa)	1.379	1.103	0.827	0.551	0.275	0	0
4	Age (days)	8	16	28	63	90	120	180
	Stress (MPa)	2.75	0.551	0.827	1.103	1.379	0	0
5	Age (days)	8	14	28	63	90	120	180
	Stress (MPa)	1.379	0.827	0.275	0.827	1.379	0	0

Table 3
Concrete parameters in Ross’ tests [¹²]

Parameter	Symbol	Unit	Value
Characteristic compressive strength	fck	MPa	38
Coefficient which depends on concrete type	s	adm	0.2
Poisson’s Ratio	ν	adm	0.15
Relative humidity of ambient environmental	RH	%	93
Fictitious thickness	hf	cm	3.939
Age at the beginning of shrinkage	ts	days	7
Coefficient which depends on concrete type - shrinkage	βsc	adm	8
Temperature	temp	°C	17
Coefficient which depends on concrete type	α	adm	1

In addition, the difference between the experimental results and the analytical curves is because of the approximation of the CEB-FIB MC90 formulation. This variation is expected, since the formulation is based on a series of experiments and considerations, including statistic, of many studies. These differences are common and inherent to the used formulation, conform it was verified by [⁶].

5. Conclusions

This paper proposes to present the implementation of the concrete long-term behavior, considering the effects of creep and shrinkage, especially the creep, through the Bazant and Prasannan ([¹], [²]) Solidification theory adapted with CEB-MC90 [³] expressions. The adaptation was based on Dias [⁷], however it was programmed in the customization feature of ANSYS. The behavior was tested and validated with analytical solution of CEB-FIP MC90 and with some experimental tests performed by Ross [¹²]. It was obtained an excellent agreement between the numerical and analytical results, being the difference in relation to the experimental results inherent to the CEB-FIP MC90 model.

It was verified the importance of using small time increments (minors than 5 days) during the analysis to obtain precision in the prediction of deformations. This is due to the simplification of the linear variation of the stresses in the integration of the strain rates adopted during the deduction of the Solidification theory incremental algorithm.

In these numerical simulations, plane (PLANE182 and PLANE183) and three-dimensional (SOLID185 and SOLID186) elements were tested in unconfined compression with one and several elements. It shows that it is possible to use this customization for three-dimensional, plane strain and axisymmetric conditions that involve long-term deformations concrete. In this way, it is considered that the model for the concrete behavior is validated and can be applied in several types of structures, having the advantage of using ANSYS software that has several preprocessing, processing, post-processing, optimization and incorporation of reinforcement. Studies that use this procedure are: Quevedo [⁹], who used this programming to analyze long-term effects in deep lined tunnels and Schmitz [¹⁰] who did numerical analysis of experiments with composite steel-concrete beams and used this customization to predict the deformations during the useful life of a highway bridge with composite structure.

Two advantages are evident from this approach: to eliminate the need to save stress history in a nonlinear analysis with long-term effects and use of this approach with elements that allow incorporate reinforcement.

6. Acknowledgments

The authors are grateful to CAPES and CNPq for the financial support and also the CEMACOM/UFRGS for providing its infrastructure for the development of this work.

7. References

• ^[1]
  BAZANT, Z.; P.; PRASANNAN, S. Solidification theory for aging creep. Cement and concrete research. USA, v. 18, n. 6, p. 923-932, 1988.
• ^[2]
  BAZANT, Z.; P.; PRASANNAN, S. Solidification theory for aging creep II: verification and application. Journal of Engineering Mechanics. USA, v. 115, n. 8, p. 1704-1725, 1989.
• ^[3]
  COMITÉ EURO-INTERNATIONAL DU BÉTON. CEB-FIP Model Code 1990. ThomasTelford: London, 1993. 437p.
• ^[4]
  HATT, W. K. Notes on the effect of time element in loading reinforced concrete beams. In: ASTM. Proceeding…, p.421-33, 1907.
• ^[5]
  RILEM TECHNICAL COMMITEES. Measurement of time-dependent strains of concrete. Materials and Structures/Matériaux et Contructions, v. 31, p. 507-512, 1998.
• ^[6]
  FANOURAKIS, G. C.; BALLIM, Y. Predicting creep deformation of concrete: a comparison of results from different investigations. In: 11th FIG SYMPOSIUM ON DEFORMATION MEASUREMENTS, 11, 2003. Santorini, Greece. Proceedings… Santorini, Greece: 2003.
• ^[7]
  DIAS, M. M. Análise numérica de vigas mistas aço-concreto pelo método dos elementos finitos: efeitos de longa duração, 2013. 177 f. Dissertação (Mestrado) - Universidade Federal do Rio Grande do Sul, Porto Alegre. 2013.
• ^[8]
  DIAS, M. M.; TAMAYO, J. L. P.; MORSCH, I. B.; AWRUCH, A. M. Time dependente finite element analysis of steels-concrete composite beams considering partial interaction. Computers and Concrete, v. 15, n. 4, p. 687-707, 2015.
• ^[9]
  QUEVEDO, F. P. M. Comportamento a longo prazo de túneis profundos revestidos com concreto: modelo em elementos finitos. 2017. 203f. Dissertação (Mestrado em Engenharia Civil) - Programa de Pós-Graduação em Engenharia Civil, UFRGS, Porto Alegre, 2017.
• ^[10]
  SCHMITZ, R. J. Estrutura mista aço-concreto: análise de uma ponte composta por vigas de alma cheia. 2017. 212 f. Dissertação (Mestrado Engenharia Civil) - Programa de Pós Graduação em Engenharia Civil, Universidade Federal do Rio Grande do Sul, Porto Alegre, 2017.
• ^[11]
  FÉDERATION INTERNATIONALE DU BÉTON. fib Model Code 2010. Final Draft. V. 1, (Bulletins 65), 2012.
• ^[12]
  ROSS, A. Creep of concrete under variable stress. Journal of ACI Proceedings, p. 739-758, 1958.
• ^[13]
  ANSYS, Inc. Ansys Mechanical APDL theory reference. Release 15.0, 2013.
• ^[14]
  ANSYS, Inc. Element Reference. Release 21.0, 2009.

Publication Dates

• Publication in this collection
  Aug 2018

History

• Received
  23 Aug 2017
• Accepted
  25 Jan 2018