Customização de um programa de elementos finitos para análise de estruturas de concreto: efeitos de longa duração

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

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

1. Introdução

Desde que foi descoberto por Hatt [⁴], em 1907, diversos pesquisadores da ciência dos materiais e do cimento vêm estudando o fenômeno do comportamento diferido do concreto. Atualmente esse comportamento é agrupado em fluência (creep) e retração (shrinkage), sendo que uma das principais diferenças entre ambos é que o primeiro depende do carregamento, enquanto que o segundo não.

A fluência é caracterizada pelo acréscimo contínuo e gradual das deformações estando o material sob tensão. Essa deformação (excluído o efeito de Poisson) possui o mesmo sentido e direção da carga e é dividida, por conveniência analítica, em fluência básica e por secagem. A primeira ocorre sem troca de água com o meio externo e pode ser medida quando a peça está imersa em um ambiente com 100% de umidade. Por outro lado, a fluência por secagem é quando ocorre a troca de água com o meio ambiente e depende, portanto, dentre outros fatores, da umidade relativa do ar e da exposição da peça ao ambiente.

Já a retração é a redução do volume do material devido à perda gradual de água independente da tensão. É também dividida, por conveniência, em retração autógena e retração por secagem. A primeira ocorre sem perda de água para o meio ambiente e é consequência da remoção da água dos poros capilares pelas reações de hidratação do cimento, enquanto que a segunda, também chamada de retração hidráulica, ocorre devido à troca de água com o meio ambiente. Como a retração está associada à redução de volume, quando a peça está impedida de reduzir, esse fenômeno induz tensões de tração internas que podem ocasionar fissuras na peça. Nesse aspecto, a minimização da retração é um dos principais motivos para o procedimento de cura do concreto, que busca reduzir a perda de água nas primeiras horas que seguem a concretagem, uma vez que nessas idades o concreto está recém desenvolvendo sua resistência.

Essas deformações diferidas do concreto possuem importante influência no comportamento estrutural, uma vez que essas deformações podem atingir de 2 a 3 vezes o valor da deformação instantânea, conforme RILEM Thecnical Committees [⁵]. Como consequência, várias formulações semi-empíricas para previsão da deformação por fluência e retração encontram-se disponíveis na literatura e nas normas de projeto. Contudo, essas formulações, em geral, apresentam certa aproximação quando comparadas a ensaios experimentais. Conforme constatado por Fanourakis e Ballim [⁶], o erro pode variar de 25,9% a 67,4%, dependendo da formulação e ensaio utilizado na comparação.

O presente artigo se baseia nos trabalhos de Dias [⁷], Dias et al. [⁸], Quevedo [⁹] e Schmitz [¹⁰], que utilizam o modelo diferido do concreto proposto pelo Comité Euro-International du Béton (CEB-FIP MC90) [³], com a formulação da fluência adaptada à Teoria da Solidificação de Bazant e Prasannan ([¹], [²]) em seus estudos envolvendo vigas mistas, túneis e pontes, respectivamente. A formulação de fluência do CEB-FIP MC90 possui a vantagem de se encaixar perfeitamente na teoria da Solidificação de Bazant e Prasannan. Isso é possível, pois essa formulação separa o coeficiente de fluência em um fator que depende da idade do concreto e outro que depende exclusivamente da idade da carga. Tal separação, por exemplo, não ocorre em formulações mais atuais, como por exemplo, na Féderation Internationale du Béton [¹¹]. Em vista dessa decomposição, o algoritmo incremental desenvolvido por Bazant e Prasannan [²] permite modelar eficientemente os casos em que há um histórico de cargas (ou tensões) variáveis. O uso da formulação do CEB-FIP MC90 sem o algoritmo, em uma análise não linear, requer que sejam guardados todos históricos de tensões para posteriormente aplicar o princípio da superposição. Esse procedimento tornaria inviável a análise em situações com muitos elementos finitos.

Portanto, nesse artigo será apresentada a formulação da fluência considerando essa adaptação e alguns comentários sobre a customização do ANSYS e os elementos finitos utilizados. Ao final serão apresentados os resultados utilizados na validação dessa modelo.

2. Modelo numérico

Neste item é apresentado o modelo para o concreto considerando o comportamento diferido no tempo, especialmente a fluência e a adaptação da Teoria da Solidificação ao formulário do CEB-FIP MC90.

2.1 Adaptação da formulação da fluência do CEB-FIP MC90 à Teoria da Solidificação

Considerando que o concreto envelhece e a fluência é a parcela do comportamento diferido relacionada diretamente com as cargas aplicadas na estrutura, a idade do concreto nos instantes de aplicação das cargas influencia diretamente nas deformações por fluência. A Teoria da Solidificação proposta por Bazant e Prasannan ([¹], [²]) tem o intuito de facilitar a solução numérica desse comportamento diferido do concreto. Esses autores propuseram um modelo que considera o envelhecimento como um fator isolado, relacionado apenas com o volume de concreto solidificado, , ao longo do tempo.

A Figura 1 apresenta de forma esquematizada o modelo proposto por Bazant e Prasannan para o comportamento do concreto. A deformação total é dada pela soma de quatro parcelas: elástica, viscoelástica, viscosa e uma parcela que é independente das tensões, correspondente às deformações por retração, térmica e fissuração. Como pode ser visto a deformação por fluência está dividida em duas parcelas: uma chamada de “viscoelástica” e outra chamada de “viscosa”. Essa divisão tem o intuito de tratar a forma da curva de fluência em relação a idade do carregamento. Para idades de carregamento pequenas a função tem forma de potência (fluência primária), e para idades de carregamentos maiores, a fluência assume a forma predominantemente logarítmica (o que muitas vezes é convenientemente representado como uma reta na escala semi-log do tempo) (Bazant e Prasannan, [¹]).

Figura 1
Modelo do concreto (adaptado de Bazant e Prasannan [¹], [²])

Conforme indica a Figura 1, a parcela dita viscoelástica ε_v está relacionada com a fração de volume de concreto já solidificada v(t) e com um coeficiente de fluência dependente da idade do carregamento Ф(t - t₀), em que o tempo t é a idade atual do concreto e t₀ é a idade na qual foi aplicada a carga. De forma análoga, a parcela dita viscosa ε_f depende da fração de cimento hidratado h(t) e do coeficiente ψ(t - t₀) também dependente da idade do carregamento (Bazant e Prasannan, [¹]). Ainda na Figura 1, chama-se atenção para a parcela Ф(t - t₀), que está representada através de um modelo reológico de Kelvin-Generalizado (ou também conhecido como cadeia de Kelvin). Como essa função depende apenas da idade da carga, os parâmetros da cadeia de Kelvin se tornam constantes, facilitando o desenvolvimento do algoritmo incremental, posteriormente descrito.

A partir dessas considerações a função de fluência , conforme a Teoria da Solidificação é dada pela expressão (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)

Sendo E_0 o módulo de elasticidade referente aos agregados e partículas microscópicas da pasta de cimento, γ(t - t₀) a micro deformação viscoelástica, e η(t) a viscosidade macroscópica aparente dependente da idade do concreto. Já no CEB-FIP MC90 [³] a função de fluência é expressa por (2):

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

Sendo E_c (t₀) o módulo de elasticidade tangente em t₀, E_ci o módulo de elasticidade tangente aos 28 dias e Ф(t,t_0 ) o coeficiente de fluência. Em [³] o coeficiente de fluência é dividido em duas parcelas que se multiplicam: ϕ₀ que está relacionado com a idade de aplicação de carga e β_c (t - t₀) que descreve a forma da curva de fluência com a idade da carga. Comparando as expressões (1) e (2) tem-se que a primeira etapa da adaptação da Teoria da Solidificação com o CEB-FIP MC90 é dada pelas expressões (3), (4), (5) e (6).

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

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

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

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

A igualdade (6) ocorre, pois Bazant e Prasannan ([¹], [²]) tratam a curva de fluência separando-a em duas parcelas: uma correspondente às primeiras idades (viscoelástica) e outra se refere a idades posteriores (viscosa), contudo, a formulação do CEB-FIP MC90 não faz essa distinção representando o comportamento da curva de fluência diretamente pela função β_c (t - t₀).

Apesar da Teoria da Solidificação de Bazant e Prasannan ([¹], [²]) considerar um fator para incluir a não-linearidade da fluência em relação ao nível de tensão F[σ(t)], o modelo apresentado no CEB-FIP MC90 não considera esse comportamento e, portanto, esse fator não aparece nessa adaptação. Contudo, devido a essa simplificação é imposto o limite de tensão de 40% da resistência média à compressão. Esse limite é fundamental para se manter o concreto em regime viscoelástico, onde é válido o princípio da superposição dos efeitos. Portanto, durante as análises numéricas, considerando essa adaptação, deve-se sempre respeitar esse limite a fim de não subestimar as deformações.

2.2 Ajuste do fator do coeficiente de fluência que independe da idade do concreto

Como visto anteriormente, o fator do coeficiente de fluência, que independe da idade do concreto, no algoritmo incremental de Bazant e Prasannan [¹], é representado através de uma cadeia de Kelvin. Uma vez que esse fator é independente da idade do concreto, os parâmetros dessa cadeia (módulo de elasticidade das molas e coeficiente de viscosidade dos amortecedores) também passam a ser constantes e isso facilita a dedução do algoritmo incremental, que utiliza uma solução numérica para a integral da taxa de deformação.

As equações que representam o comportamento da cadeia de Kelvin estão apresentadas nas expressões (7). Devido à configuração em série das μ-ésimas unidades de Kelvin, para cada unidade, a tensão aplicada σ_μ é a mesma que a total σ e a deformação total da cadeia γ é dada pela soma da deformação das unidades.

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

Sendo E_μ o módulo de elasticidade, η_μ a viscosidade e γ_μ a deformação da μ-ésima unidade da cadeia. Solucionando a equação diferencial apresentada em (7), obtém-se a expressão (8) que é conhecida como série de Dirichlet (ou Prony). O termo τ_μ é denominado tempo de retardamento da μ-iésima unidade da cadeia.

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

A fim de utilizar a Teoria da Solidificação, adaptada ao CEB-FIP MC90, a expressão (8) deve representar a função β_c (t - t₀), conforme apresentado em (4). Isso pode ser feito discretizando essa função em L pontos e ajustando-a à (8) através do método dos mínimos quadrados. Conforme Dias [⁷] isso pode ser feito através da montagem e resolução do sistema [A]{X}=[B], cujos elementos são dados por (9), (10) e (11).

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

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

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

Através de (11) obtêm-se os módulos de elasticidade da cadeia de Kelvin. É importante destacar que o ajuste do modelo pode levar a valores negativos de módulo de elasticidade em determinados elementos, o que não têm correspondência física com o modelo mecânico, se tratando apenas de um ajuste matemático. Ademais, a resolução desse sistema linear necessita da escolha dos tempos de retardamento dos elementos da cadeia τ_μ. Para que a matriz [A] seja bem condicionada, Bazant e Prasannan [²] propõem que os tempos de retardamento sejam escolhidos como um incremento logarítmico de base 10 conforme (12).

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

Segundo Dias [⁷], os tempos de retardamento devem abranger metade do período de análise t_máx e τ₁ deve ser pequeno o suficiente para representar o comportamento da curva de fluência em pequenas idades, tal como expresso em (13).

${\tau }_1=0.01t_0;\mathrm{}{\tau }_N=\mathrm{0,5}{t}_{m\mathrm{á}x}$ (13)

De acordo com Dias [⁷], a quantidade N de unidades de Kelvin foi determinada como sendo a máxima possível até que se tenha um τ_μ superior ou igual a τ_N. Porém, não mais do que 6 unidades de Kelvin.

Além dos tempos de retardamento, para se resolver o sistema, também é necessário escolher os L pontos da função β_c (t - t₀) que serão utilizados para o ajuste. Como essa função está relacionada com a idade do carregamento, na verdade, a questão está em definir os intervalos de idades (t_k - t₀), em que t₀ é a idade do concreto no instante da primeira aplicação de carga e k=1,…,L pontos. Considerando que o comportamento da função é mais acentuado nas primeiras idades de carregamentos, uma escolha eficaz é manter os intervalos de tempo constantes em escala logarítmica. Dessa forma, os pontos para o ajuste podem ser obtidos através da fórmula de recorrência (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);\mathrm{}\mathrm{}\mathrm{}k=2,\dots ,\mathrm{}L$ (14)

Em que o parâmetro m é o número de passos por década, sendo este valor, segundo [⁷], em torno de 10 para uma boa precisão. Além disso, o primeiro intervalo de tempo escolhido corresponde a um décimo da primeira idade de carregamento da estrutura [⁷], conforme a expressão (15).

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

E a quantidade de pontos L pode ser determinada como sendo a necessária para atender o tempo em estudo, ou seja, até que t_k+1 - t₀ seja maior do que t_máx - t₀.

2.3 Algoritmo incremental

Uma vez que se tenham os parâmetros (E_μ, τ_μ) da cadeia de Kelvin ajustados pode-se utilizá-los no algoritmo incremental. Expressando a deformação total em termos de taxas tem-se (16).

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

Sendo $\dot{{\epsilon }_e}=\dot{\sigma }\left(t\right)/E\left(t\right)$, a taxa de deformação elástica, $\dot{{\epsilon }_v}=\dot{\gamma }(t-t_0)/v\left(t\right)$, a taxa de deformação viscosa e $\dot{{\epsilon }_0}$ o somatório de todas as taxas de deformação independentes da tensão (retração, deformação térmica e fissuração).

Integrando a expressão (16) através da regra dos trapézios, para um intervalo de integração de t_i até t_i+1 obtém-se (17), em que E_i+1/2 e v_i+1/2 são os valores correspondentes ao meio do intervalo entre t_i e t_i+1.

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

Se o incremento de deformação elástica for separado do restante da deformação, pode-se escrever (18), sendo a parcela ∆ε* toda a deformação que não é elástica (fluência, retração, térmica e fissuração) e E* o módulo de elasticidade equivalente, ambos a serem deduzidos.

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

A deformação por fluência é obtida assumindo-se que a tensão varia linearmente entre dois passos de tempo consecutivos, conforme (19). Essa suposição, como será visto, impõe a necessidade de um incremento de tempo pequeno.

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

Introduzindo (19) em (7), tem-se (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)

Sendo ${\gamma }_{{\mu }_{i+1}}$ a deformação da μ-ésima unidade de Kelvin no passo de tempo seguinte (i+1) e ${\gamma }_{{\mu }_i}$ é a deformação μ-ésima unidade de Kelvin para o passo de tempo atual i. Da mesma forma σ_i é a tensão atual e ∆σ é o incremento de tensão considerando o intervalo de tempo entre t_i e t_i+1. Assim o incremento de deformações pode ser calculado como ${∆\gamma }_{\mu }={\gamma }_{{\mu }_{i+1}}-{\gamma }_{{\mu }_i}$, resultando em (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)

Sendo o primeiro fator do primeiro termo à direita da igualdade a deformação viscosa da μ-ésima unidade de Kelvin do passo i, expressa por (22).

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

E de forma análoga a expressão (22), em (23) tem-se a deformação viscosa da unidade de Kelvin do passo seguinte i+1.

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

Conforme deduzido por Dias [⁷], substituindo ${\gamma }_{{\mu }_i+1}$ dado pela equação (20), na expressão (23), obtém-se a expressão (24), que é a deformação viscosa ainda desconsiderando o envelhecimento. O envelhecimento é considerado incluindo o volume de concreto solidificado na expressão (24), de forma que a deformação viscosa é dada pela equação (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)

Para obter o incremento de deformação viscosa, substitui-se a expressão (21) em (17) obtendo-se, primeiramente o incremento de deformação total, dado pela expressão (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)

Considerando a separação da parcela de deformação elástica proposta pela equação (18) podem-se definir as expressões para ∆ε* e E*, que ficam conforme as equações (27) e (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)

Generalizando as expressões para um domínio tridimensional, considerando um material isotrópico tem-se que (27) torna-se uma equação matricial, conforme (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)

Sendo a parcela $\left\{{\epsilon }_{{\mu }_i}^{*}\right\}$ calculada pela expressão (30), generalizada de (25), considerando que [D_μ]^-1 é o inverso da matriz constitutiva isotrópica linear, sendo o módulo de elasticidade igual a $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)

Na expressão (30) pode-se verificar que é necessário guardar apenas a deformação viscosa correspondente ao passo anterior $\left\{{\epsilon }_{{\mu }_{i-1}}^{\mathrm{*}}\right\}$. Evitando-se assim guardar todo o histórico de carregamento.

Durante o processo de solução iterativo, além da determinação do incremento de deformação não elástica, é necessário também atualizar a matriz constitutiva e corrigir o incremento de deformação total no cálculo do incremento de tensões. Generalizando para um domínio tridimensional a expressão (18) e isolando o incremento de tensão obtém-se (31).

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

Sendo [D*] a matriz constitutiva isotrópica linear considerando o módulo de elasticidade E* calculado pela expressão (28) e o coeficiente de Poisson do concreto.

3. Customizando o ANSYS

3.1 Funcionamento do recurso programável Usermat

O ANSYS possui uma série de recursos programáveis (User Programmable Features - UPF) que permitem ao usuário customizar aspectos do modelo, através de sub-rotinas escritas em Fortran77. Uma vez alteradas, essas sub-rotinas são então compiladas e associadas ao programa principal do ANSYS. Uma dessas sub-rotinas, a UserMat, é específica para a customização do comportamento do material e, portanto, ideal para introduzir o modelo diferido do concreto.

O processo iterativo utilizado pelo ANSYS para resolver problemas não-lineares consiste em aplicar o método de Newton-Raphson. Para cada incremento n de tempo ou carga (também chamado de subpasso) é executado i-ésimas iterações de equilíbrio até que o critério de convergência seja satisfeito. As expressões (32), (33) e a Figura 2 (adaptados de ANSYS [¹³]) exemplificam esse processo.

${\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)

Figura 2
Ilustração do método de Newton-Raphson com iterações de equilíbrio i à cada subpasso n (ANSYS, [¹³])

Sendo {u}_n,i+1 o vetor de deslocamentos nodais do subpasso n atualizado, {u}_n,i vetor de deslocamentos nodais da iteração i, {∆u}_n,i incremento de deslocamentos nodais da iteração i, [K^T]_n,i matriz de rigidez tangente da iteração i (ou matriz Jacobiana do método Newton-Raphson), {F^a}_n+1 vetor de forças atuantes, {F^nr}_n,i vetor de forças internas da iteração i, {R}_n,i vetor de desequilíbrio ou resíduo da iteração i utilizado no critério de convergência. Com os deslocamentos nodais, para cada iteração de equilíbrio i, são obtidas as deformações e tensões em cada ponto de integração através das relações de deslocamento-deformação e deformação-tensão formuladas através das funções de interpolação do método dos elementos finitos. Com isso, para cada ponto de integração, o programa principal entrega para a UserMat as tensões totais, deformações totais e o incremento de deformação total (vindo do processo iterativo de Newton-Raphson) no corrente incremento de tempo. As instruções programadas dentro da sub-rotina ficam responsáveis por coletar as constantes do material (definidas pelo usuário), atualizar as tensões e a matriz Jacobiana utilizada no processo iterativo. A Figura 3 apresenta um fluxograma resumindo o funcionamento da UserMat dentro da análise não linear.

Figura 3
Funcionamento da UserMat no processo de solução não linear

3.2 Implementação do modelo viscoelástico na Usermat

O ANSYS disponibiliza uma UserMat com uma relação tensão-deformação bilinear isotrópica com critério de plasticidade de von Mises. Cabe ao usuário entender essa sub-rotina e modificar conforme sua necessidade. É muito importante compreender todos os argumentos de entrada e saída da sub-rotina para criar e utilizar o menor número de variáveis possíveis no momento da customização do comportamento do material. Atenção especial deve ser dada às variáveis que servem de entrada e saída da UserMat, já que estas não devem ser apagadas. Já as variáveis locais, que em grande parte estão referenciadas ao comportamento plástico de von Mises, podem ser eliminadas sem comprometerem o funcionamento da sub-rotina.

A Usermat possui quatro sub-rotinas internas dependendo do número de componentes de deformação (ou tensão) e das dimensões envolvidas no problema (1D, 2D ou 3D). No presente artigo foi customizado a sub-rotina mais geral das quatro, a UserMat3D, que corresponde a problemas que possuem mais de 3 componentes, tal como problemas tridimensionais (6 componentes e 3 direções), estado plano de deformação e axissimétrico (ambos com 4 componentes de tensão e 2 direções).

O usuário é livre para criar quantas variáveis locais forem necessárias. Contudo, variáveis locais perdem seu valor quando a sub-rotina encerra durante as iterações de equilíbrio. Quando é necessário manter o valor de alguma variável entre subpassos é necessário guardá-lo no array de variáveis de estado ustatev (Figura 4-(d)). Inclusive, pode-se utilizar diretamente esse array sem a necessidade de criar uma variável. Contudo, para que o array ustatev esteja dimensionado é necessário declarar seu tamanho no script através do comando TB,STATE durante a atribuição do material (Figura 4-(a)). Não é aconselhável guardar o valor das variáveis em um bloco COMMON, exceto as que mantêm os seus valores constantes ao longo da análise. Esse tipo de declaração traz problemas junto à paralelização do ANSYS, pois vários processos paralelos podem acessar e sobrescrever esse espaço na memória simultaneamente. É possível também visualizar os valores do array ustatev no pós-processamento através do comando PLESOL,SVAR,[posição da variável no array], desde que sejam salvos através do comando OUTRES,SVAR,ALL.

Figura 4
Comandos APDL para (a) definição do material, (b) definição do tempo e do incremento de tempo da análise. No interior da Usermat3D: (c) atribuição das variáveis referente ao material e (d) guardando variáveis de tensões principais no array ustatev

O array de entrada de dados prop contém as constantes do material que são preenchidas através do script durante a definição do material. Esse array é dimensionado através do comando TB,USER (que define a quantidade nprop de propriedades) e os valores das propriedades são atribuídos através do comando TB,DATA (Figura 4-(a), Figura 4-(c)).

Outro aspecto importante da utilização da UserMat em viscoelasticidade do concreto é a definição dos tempos e idades do concreto durante a análise. O ANSYS possui uma variável time e dtime, compartilhada entre o script e a UserMat que representa o tempo corrente da análise e o incremento de tempo, respectivamente (Figura 4-(b)).

De uma forma geral a UserMat3D customizada segue o seguinte roteiro:

• 1) cálculo do incremento de deformação não elástica {∆ε*};

• 2) cálculo do incremento de deformação elástica: {∆ε_e} = {∆ε} - {∆ε*}, sendo {∆ε} o incremento de deformação total informado pelo ANSYS (vindo do processo iterativo de Newton-Raphson);

• 3) cálculo da matriz constitutiva [D];

• 4) cálculo do incremento de tensão {∆σ} = [D] {∆ε_e};

• 5) atualização da matriz Jacobiana (∂∆σ_ij) ⁄ (∂ε_ij) = [D];

• 6) cálculo da tensão atualizada {σ_(i+1)} = {σ_i} + {∆σ}.

Ao se retirar o incremento de deformação não elástica do incremento de deformação total (passo 2), o método de Newton-Raphson compensará aumentando o incremento de deformação total durante as iterações de equilíbrio até que a convergência seja satisfeita. Dessa forma é introduzida a deformação viscosa do concreto no problema.

Na UserMat customizada, o concreto possui um instante inicial t_i, opção que pode ser útil quando se utiliza o recurso de ativação/desativação dos elementos para simular o processo de construção de uma ponte ou túnel, por exemplo. A idade do material é dada por t_mat = time - ti. Porém, os fenômenos de retração e fluência, possuem, respectivamente, o tempo de secagem t_s e o tempo em que a carga é aplicada t₀. Dessa forma, o fenômeno da retração inicia quando time > t_i, ou seja, no instante em que o concreto foi iniciado, porém considerando t_s, nas expressões do CEB-FIP MC90. Já o fenômeno da fluência inicia quando time > t₀, ou seja, a partir do momento em que a peça entra em carga, sendo a idade da carga dada por t_mat - t₀.

Em Quevedo [⁹], tem-se a sub-rotina UserMat3D para o concreto viscoelástico customizada e acompanhada das as instruções para adicioná-la ao programa principal do ANSYS.

3.3 Elementos finitos compatíveis

A UserMat3D pode ser utilizada em problemas tridimensionais genéricos, estado plano de deformação ou então axissimétricos. Para o primeiro caso podem ser utilizados os elementos SOLID185 e SOLID186 já para as duas outras situações existem os elementos PLANE182 e PLANE183.

Os elementos PLANE182 e PLANE183, apresentados na Figura 5, são elementos bastante semelhantes, ambos são bidimensionais quadriláteros e possuem dois graus de liberdade por nó (translação nas duas direções do plano em que estão inseridos). A principal diferença entre eles é que o PLANE182 possui 4 nós e, portanto, funções de interpolações de deslocamentos lineares enquanto que o PLANE183 possui 8 nós e, portanto, funções quadráticas de interpolação. Além da forma quadrilátera, ambos os elementos englobam também a forma triangular, sendo que essa forma não é aconselhável para o PLANE182 por apresentar pouca precisão de interpolação.

Figura 5
Elementos planos: (a) PLANE182 e (b) PLANE183 (ANSYS, [¹⁴])

Os elementos SOLID185 e SOLID186, ilustrados na Figura 6, são elementos sólidos, sendo assim, seus nós possuem três graus de liberdade de deslocamentos nas direções x, y e z. Novamente a principal diferença entre os dois tipos de elementos é o grau das funções de interpolação, o SOLID185 possui 8 nós e, portanto, funções lineares enquanto que o SOLID186 possui 20 nós e, portanto, funções de interpolação quadráticas. Ambos os elementos englobam também a forma tetraédrica e piramidal.

Figura 6
Elementos tridimensionais: (a) SOLID185 e (b) SOLID186 (ANSYS, [¹⁴])

A escolha pelo uso de elementos quadráticos deve levar em conta, além da quantidade de elementos da malha, o tempo de análise dos problemas, pois apesar de ser possível o uso de uma malha menos refinada, o tempo de processamento costuma ser muito maior em comparação aos elementos lineares devido à quantidade de nós. Além disso, por se tratar de uma análise não-linear, em que o sistema de equações é resolvido diversas vezes em função do incremento de tempo, o tempo de processamento pode se tornar alto.

Para análises de concreto armado, ambas as classes de elementos (PLANE e SOLID) possuem o recurso de armadura incorporada, através da inclusão de elementos REINF263 (para os elementos planos) e REINF264 (para os elementos sólidos) após a geração da malha dos elementos de concreto. Esses elementos são ideais para representar armadura, pois possuem apenas rigidez axial e podem ser inseridos em qualquer posição e direção dentro dos elementos de concreto. Esse aspecto traz uma enorme vantagem em relação ao elemento de concreto do ANSYS (SOLID65) no qual a armadura não é incorporada exigindo uma malha muito mais refinada.

4. Validação e exemplos numéricos

A fim de validar o modelo programado na UserMat, foram feitas algumas simulações numéricas, sendo que os primeiros dois testes serviram para comparar os resultados com a solução analítica obtida pelas formulações do CEB-FIP MC90. Esses testes consistiram na modelagem de um ensaio de compressão não confinado, tal como na Figura 7. As propriedades consideradas nestas simulações iniciais estão apresentadas na Tabela 1. O primeiro teste consistiu em aplicar incrementos de tensão de 5MPa de compressão aos 10, 50 e 75 dias, sendo a análise feita até os 100 dias. Já o segundo teste, consistiu na aplicação de 5MPa constantes até os 100 dias. As Figuras 8 e 9 apresentam os resultados obtidos decompostos em deformação total, elástica, fluência e retração. Nos dois testes foi verificada a concordância perfeita entre as simulações numéricas e o cálculo analítico. Esses testes também foram repetidos considerando um elemento e vários elementos.

Figura 7
Modelo para o teste da implementação do modelo viscoelástico

Figura 8
Teste com cubo sujeito a tensões variáveis utilizando o elemento SOLID185

Figura 9
Teste com cubo sujeito a tensões constante utilizando o elemento SOLID185

Tabela 1
Parâmetros para o concreto nos testes iniciais

Parâmetro	Símbolo	Unidade	Valor
Resistência característica à compressão	fck	MPa	40
Coeficiente que depende do tipo de cimento	s	adm	0,25
Coeficiente de Poisson	ν	adm	0,2
Umidade relativa do ar	RH	%	70
Espessura fictícia	hf	cm	54,54
Idade do concreto no final da cura	ts	dias	7
Coeficiente que depende do tipo de cimento - retração	βsc	adm	5
Temperatura	temp	°C	20
Coeficiente que depende do tipo de cimento	α	adm	1

O próximo teste consistiu na verificação da influência do incremento de tempo utilizado nas análises. Até então o incremento de tempo considerado foi apenas de 1 dia. Porém, o algoritmo incremental da Teoria da Solidificação supõe uma variação linear da tensão entre dois intervalos de tempos consecutivos. Porém, essa variação não é de fato linear, e dessa forma, quanto maior o incremento de tempo e o tempo total de análise, maior será o erro no cálculo da deformação por fluência. Na Figura 10 estão apresentados os resultados dessa análise considerando apenas a deformação elástica mais a deformação por fluência para 5 incrementos de tempo: 1, 2, 5, 10 e 20 dias. Observa-se que com o incremento de 5 dias, o erro ainda é pequeno, cerca de 2,5%. Entretanto para o incremento de 10 dias, o erro atinge 4,8% e para 20 dias, é de 7,5%. Portanto, para melhor precisão nos resultados sugere-se utilizar incrementos de tempo menores do que 5 dias.

Figura 10
Influência do incremento de tempo

A fim de comparar os resultados com dados experimentais, foram simulados numericamente 5 ensaios experimentais apresentados por Ross [¹²]. Esses ensaios experimentais consistem de dois corpos de prova cilíndricos, com diâmetro de 117,5mm e 305mm de altura, um deles sujeito a um histórico de cargas (Tabela 2) e o outro sem carregamentos, para que pudesse ser medida a deformação gerada pela retração e posteriormente subtraída dos resultados do corpo de prova sujeito ao carregamento, tendo dessa forma apenas resposta elástica com fluência. As propriedades do concreto nos ensaios de Ross estão apresentadas na Tabela 3. Para as simulações numéricas foram testados os dois tipos de elementos, o SOLID185 e o PLANE182 em compressão axial, tal como na Figura 7. Os resultados obtidos encontram-se nas Figuras 11 a 15. Pode-se observar a perfeita concordância entre o modelo SOLID185 e a formulação do CEB-FIP MC90, sendo que as curvas ficam sobrepostas. Em relação ao elemento PLANE182, percebe-se que em todas as situações existe um pequeno erro. Esse erro é atribuído ao fato do ensaio não se tratar rigorosamente de um estado plano de deformações ou de axissimetria, que são o escopo desse elemento. Apesar disso, pode-se confirmar que a UserMat3D não está restrita apenas a problemas tridimensionais.

Figura 11
Comparação com ensaios: teste 1

Figura 12
Comparação com ensaios: teste 2

Figura 13
Comparação com ensaios: teste 3

Figura 14
Comparação com ensaios: teste 4

Figura 15
Comparação com ensaios: teste 5

Tabela 2
Histórico de tensão nos testes de Ross [¹²]

Testes		t0	t1	t2	t3	t4	t5	tf
1	Idade (dias)	14	60	-	-	-	-	140
	Tensão (MPa)	1,503	0	-	-	-	-	0
2	Idade (dias)	28	60	91	120	154	-	190
	Tensão (MPa)	1,503	1,127	0,751	0,376	0	-	0
3	Idade (dias)	8	14	28	63	90	120	180
	Tensão (MPa)	1,379	1,103	0,827	0,551	0,275	0	0
4	Idade (dias)	8	16	28	63	90	120	180
	Tensão (MPa)	2,75	0,551	0,827	1,103	1,379	0	0
5	Idade (dias)	8	14	28	63	90	120	180
	Tensão (MPa)	1,379	0,827	0,275	0,827	1,379	0	0

Tabela 3
Parâmetros para o concreto nos testes de Ross [¹²]

Parâmetro	Símbolo	Unidade	Valor
Resistência característica à compressão	fck	MPa	38
Coeficiente que depende do tipo de cimento	s	adm	0,2
Coeficiente de Poisson	ν	adm	0,15
Umidade relativa do ar	RH	%	93
Espessura fictícia	hf	cm	3,939
Idade do concreto no final da cura	ts	dias	7
Coeficiente que depende do tipo de cimento - retração	βsc	adm	8
Temperatura	temp	°C	17
Coeficiente que depende do tipo de cimento	α	adm	1

Ademais, a diferença entre os resultados experimentais e as curvas analíticas deve-se à aproximação da formulação do CEB-MC90. Essa variação é esperada, uma vez que a formulação é baseada em uma série de ensaios e considerações, inclusive estatísticas, de diversos estudos. Essas diferenças são comuns e inerentes à formulação utilizada, conforme já constatado por [⁶].

5. Conclusões

O presente trabalho propôs apresentar os resultados da implementação do comportamento diferido para o concreto, considerando os efeitos de fluência e retração, especialmente a fluência, através da teoria da Solidificação de Bazant e Prasannan ([¹], [²]) adaptada com as expressões do CEB-MC90 [³]. A adaptação foi baseada em Dias [⁷], entretanto sua programação foi feita dentro do recurso de customização do ANSYS. O comportamento do material programado foi testado e validado com a solução analítica do CEB-FIP MC90, além da comparação com alguns ensaios experimentais realizados por Ross [¹²]. Obteve-se ótima concordância entre os resultados numéricos e analíticos, sendo a diferença em relação aos resultados experimentais inerente ao modelo do CEB-FIP MC90.

Também foi constatada a importância do uso de incrementos de tempo pequenos (menor que 5 dias) durante a análise para obtenção de precisão na predição das deformações. O que ocorre devido à simplificação da variação linear das tensões na integração das taxas de deformação adotada durante a dedução do algoritmo incremental da Teoria da Solidificação.

Nessas simulações numéricas foram testados tanto elementos planos (PLANE182 e PLANE183) quanto tridimensionais (SOLID185 e SOLID186) em compressão não confinada, com um e vários elementos. Mostrando que é possível utilizar essa customização para problemas tridimensionais, estado plano de deformações e estado axissimétrico que envolvam análise diferida do concreto. Dessa forma, considera-se que o modelo para o comportamento do concreto está validado e pode ser aplicado em diversos tipos de estruturas, tendo a vantagem do uso do software ANSYS que possui diversos recursos de pré-processamento, processamento, pós-processamento, otimização e incorporação de armadura. São exemplos de trabalhos: Quevedo [⁹], que utilizou essa programação para análise de efeitos de longo prazo em túneis revestidos em concreto projetado considerando o processo de construção, e Schmitz [¹⁰] que fez análises numéricas de ensaios experimentais de vigas mistas de aço-concreto e utilizou essa customização para predições de deformações durante a vida útil de uma ponte rodoviária em estrutura mista.

Duas vantagens são evidentes dessa abordagem: eliminar a necessidade de salvar o histórico de tensões numa análise não linear com efeitos diferidos no concreto e o uso conjunto dessa abordagem com elementos que permitem armadura incorporada.

6. Agradecimentos

Os autores agradecem a CAPES e ao CNPq pelo apoio financeiro e também ao CEMACOM/UFRGS pela disponibilização de sua infraestrutura para realização deste trabalho.

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.

Datas de Publicação

• Publicação nesta coleção
  Ago 2018

Histórico

• Recebido
  23 Ago 2017
• Aceito
  25 Jan 2018