Estratégia h-adaptativa para a análise elasto-plástica bidimensional

Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería. Vol. 11,4, 519-528(1995) R.C.C. SILVA L. LANDAU e F.L.B. RIBEIRO Programa de Engenharia Civil COPPE/UFR J Caixa Postal 68506 - 21945 - 970 Rio de Janeiro, Brasil Neste trabalho apresenta-se uma estratégia h-adaptativa para problemas de elasticidade bidimensional, discretizados por elementos finitos triangulares lineares, considerando a nao linearidade do material. Esta estratégia é baseada na técnica de remalhamento (remeshing) onde, através de estimativa de erro a posteriori, se prediz o tamanho do elemento necessário para atingir um nível de precisa0 pré-especificado. O comportamento nao linear do material é tratado através do modelo de subcamadas, nao se limitando portanto a modelos bilineares. SUMMARY In this work an h-adaptive strategy for 2D elasticity is presented, using linear triangular elements and considering material non-linearity. This strategy is based on a remeshing technique where the size of the required element for a specific accuracy is provided by an a posteriori error estimation. The non-linear material behaviour is considered using a layered model and is not limited to bilinear models. Apesar do método de elementos finitos ser a técnica de aproximacáo de maior impacto na área de métodos numéricos, aplicados a problemas de engenharia, suas análises convencionais exigem grande esforco e experiencia por parte dos analistas, principalmente em aplicacóes industriais envolvendo análises nao lineares. Sendo assim, a automatizacáo de procedimentos, visando A análise de erro e refinamento (automático) adaptativo, torna-se uma questáo fundamental para o aperfeicoamento e consolidacáo do método. Recibido: Enero 1995 OUniversitat Politecnica de Catalunya (España) ISSN 0213-1315 520 R.C.C. SILVA, L. LANDAU E F.L.B. RIBEIRO A estratégia h-adaptativa apresentada é utilizada ao longo do processo incremental de problemas de elasticidade, considerando a náo linearidade física do material. A malha final, discretizada por elementos triangulares lineares, deverá atingir um nível de erro predeterminado pelo usuário. Duas alternativas de implementaciio do método h sáo freqüentemente utilizadas. Na primeira, a malha é sucessivamente refinada até a obtenciio do nível de precisa0 desejado. Neste caso, o elemento a ser refinado é identificado através de um critério apropriado. A outra prediz o tamanho do elemento necessário para que o nível de precisáo preestabelecido seja atingido. No caso náo linear optou-se por esta última técnica, pois é, a princípio, computacionalmente mais economica, já que evita a análise durante os sucessivos passos necessários para a obtencáo da malha final. Utilizou-se uma estimativa de erro a posteriori baseada nos conceitos propostos por Zienkiewicz e Zhu, que implica na recuperacao de um campo contínuo de tensóes suavizadas para substituir o campo exato de tensóes náo conhecido. Neste trabalho, adotou-se a técnica conhecida como superconvergent patch recovery (SPR) para a recuperacáo do campo de tensóes. Finalmente, o comportamento do material foi tratado através do modelo de subcamadas. As equacóes básicas necessárias para a resolucáo de um problema estático bidimensional náo linear podem ser deduzidas através do princípio dos trabalhos virtuais. Considera-se que para um corpo em equilíbrio (domínio O), sob a acao de forcas de volume b, e forqas de superfície i em Fa, o trabalho virtual das forcas externas é igual ao trabalho virtual das forcas internas, para todo campo de deslocamento virtual Su que se anule no contorno I',, onde o deslocamento u é prescrito (u = u em FU). Entáo, .adO=i SuT.bdCl+lo SuT.EdI' Em (1) a representa as deformacóes e a as tensóes. MODELO DISCRETO Utilizando os procedimentos normais de discretizacáo do método de elementos finitos foram obtidas as seguintes expressóes para os deslocamentos e deformacóes no elemento: onde a matriz N contém as funcóes de interpolacáo, u sao os valores nodais da solucáo aproximada e L é o operador diferencial linear. Substituindo-se (2) e (3) em (l), ESTRATÉGIA H-ADAPTATIVA PARA A ANÁLISE ELASTO-PLÁSTICA BIDIMENSIONAL 521 considerando que a expressáo assim obtida é válida para qualquer valor arbitrário de Su, tem-se Devido a nao linearidade do material existirá0 forcas residuais m, tais que m=P-F#O (5) onde P and F representam os vetores de forcas internas e os de forqas externas determinados, respectivamente, pelas parcelas do lado esquerdo e do direito de equacáo (4). Neste caso, para que se obtenha o equilíbrio estático, a análise deve ser realizada até que as forcas residuais se anulem. Técnicas de soluqáo O método utilizado para a obtencáo da resposta náo linear de um problema estático é essencialmente incremental. Assume-se que para uma determinada intensidade de carga (am), a solucáo (u,) é conhecida e procura-se a solu~áo (u, + Au,) para o nível de carga posterior (am + Aam). A solucáo através de método puramente incremental está sujeita a erros que váo se acumulando a cada passo da análise. Para contornar este inconveniente, pode-se utilizar um método iterativo tangente. Este método consiste em se obter uma melhor aproximacáo de u, utilizando uma expansáo de P através de séries de Taylor. Obtémse, entáo, uma formulacáo incremental-iterativa para solucionar problemas náo lineares na forma onde m e n representam, respectivamente, um contador de incrementos de carga e um contador de iteracoes e Kt é a matriz de rigidez tangente. Os novos deslocamentos sáo obtidos por Un+l - m+l - ~k+i + Aurn+1 a (7) A análise é realizada para este nível de carga até que se obtenha a convergencia para uma tolerancia pré-especificada. Critério de escoamento O critério de escoamento de Von Mises (Owen, Hinton3) é dado por f = ae - a0 = \/5. J:'~ - a. (8) onde a, é a tensáo efetiva; ao, a tensáo de escoamento da superfície e J2, o segundo invariante do estado desviatório. Foi considerado um endurecimento linear para o material, isto implica que, após o escoamento, as superfícies subseqüentes possuem uma tensáo de escoamento maior 522 R.C.C. SILVA, L. LANDAU E F.L.B. RIBEIRO do que a da superfície anterior. Além disso, considerou-se que o material com carregamento cíclico, possui um comportamento cinemático, isto é, as superfícies de escoamento preservam sua forma e orientaciio, mas transladam no espaco das tensóes como um corpo rígido. Para que o fluxo plástico ocorra (AA # O), as tensóes devem permanecer na superfície de escoamento (A f = O). Pode-se considerar, utilizando uma lei associativa, que as deformacoes plásticas (Aep = AA . a) ocorrem em uma direcáo normal A superfície de escoamento (a = S f /So). Encontra-se demostrado (Crisfieldl) que para as consideracoes acima é válida a expressiio Af =aT.&-H'.AA=O (9) com AA=aT.~.As/(aT.~.a+~') =aT.~.Ae/(3p+~') (10) onde D é a matriz constitutiva, Ae é o incremento total de deformacáo, H' é o módulo plástico e p = E/2(1 + u) é a constante de Lamé, onde E representa o módulo de elasticidade e u o coeficiente de Poisson. Modelo elasto-plástico de subcamadas O comportamento cinemático pode ser simulado através do modelo de subcamadas. Este modelo consiste em assumir que o sólido seja composto de várias subcamadas, sendo cada uma constituída de um material elasto-plástico perfeito. Estas subcamadas possuem propriedades diferentes (espessura e tensiio de escoamento), porém estáo sujeitas a uma mesma deformacáo (Owen et al2). Optou-se pela utilizacáo da curva característica de um material bilinear para exemplificar a deducáo das espessuras das subcamadas, embora o modelo permita a análise de materiais com curva tensiiodeformacao complexa. Neste caso, um estado de tensáo (ap) superior ao estado limite de escoamento (al) pode ser escrito na forma a, = al + Aal = al + D . (As - AA . a) (11) Considerando Ae = E, - €1, onde ep e el siio as deformacoes correspondentes As tens8es ap e al, e substituindo (10) em (ll), tem-se a, = D . el + [H1/(3p + H')] . D . (ep - el) (12) que pode ser reescrito de uma forma genérica como a, = 01 . ti + a2 . t2 (13) onde t e a para cada subcamada representam, respectivamente, a espessura e o estado de tensiio. ti = 1 - [H1/(3p + H')] e t2 = [H1/(3p + H')] (14) UI=D.E~ e u~=D.E, (15) -- - ESTRATÉGIA H-ADAPTATIVA PARA A ANÁLISE ELASTO-PLÁSTICA BIDIMENSIONAL 523 Considerando um estado de deformacáo total e qualquer, obtém-se um estado de tensa0 teste (at) com at = D E. A tensáo total (a) é obtida somando as diferentes contribuicoes de cada subcamada. Podem ocorrer duas situacoes: se a tensáo efetiva teste a: (8) for menor do que a tensáo de escoamento (a,,), entáo o material tem um comportamento elástico e a contribuiqiio para a subcamada i é at; se a: for maior do que a tensáo de escoamento (a,,), entáo houve plastificacáo e o estado de tensáo deve retornar i, superfície de escoamento da subcamada i. O retorno i, superficie de escoamento foi realizado através da técnica de retorno radial (Crisfieldl). Com o intuito de aprimorar a qualidade da solucáo numérica, utilizou-se uma estratégia h-adaptativa via remeshing, baseada em estimativas de erro a posteriori. Esta estratégia consiste em predizer o tamanho necessário de alguns elementos, identificados através de um critério apropriado, de tal forma que a malha final atinja o nível de erro predeterminado. Estimativa de erro A estimativa de erro desenvolvida por Zienkiewicz e Zhu e conhecida como z2 consiste em substituir o campo de tensoes exatas (a), desconhecido, por um campo suavizado contínuo (a*) obtido segundo uma interpolacáo de seus valores nodais (a*), aproximacáo esta igual i, realizada para os deslocamentos em (2). O erro aproximado, associado a solu~áo numérica, pode ser definido como sendo a funqáo que determina a diferenca entre o campo contínuo e o campo de tensoes aproximadas (ah). Utilizando a norma de L2 para determinar a magnitude do erro aproximado, tem-se ah)T . (a* - ah) dR (16) Sendo (le11 um escalar positivo, ele pode ser avaliado em todo o domínio ou em subdomínios (elementos) do problema. Pode-se definir um erro médio aproximado como sendo onde o índice i refere-se aos n elementos da malha. Como a ordem de grandeza da norma do erro depende das unidades empregadas no problema, procurou-se utilizar o conceito de erro relativo R.C.C. SILVA, L. LANDAU E F.L.B. RIBEIRO Suavizaciío das tensoes A estimativa de erro adotada depende da estratégia utilizada para a obtenciio do campo contínuo de tensoes. Optou-se pela técnica conhecida como superconvergent patch recovery (Zienkiewicz, Zhu6). O campo de tensoes suavizadas, no interior de cada patch (conjunto de elementos que circundam um nó), é determinado através de uma expansiio polinomial completa de grau p onde P é uma base do conjunto de todos os polinomios de grau p e c contém os parametros a serem determinados para o problema em questiio. Para funcoes lineares (p = 1) tem-se Os parametros c siio determinados através de um procedimento de ajuste, pelo processo dos mínimos quadrados, visando minimizar o quadrado das diferencas entre a expansiio polinomial e o campo de tensoes aproximadas nos pontos da amostragem (pontos de superconvergencia). Para o elemento linear triangular, o ponto de superconvergencia coincide com o baricentro do mesmo, ou seja, a tensiio obtida através da análise numérica é mais precisa no baricentro do que em qualquer outro ponto deste elemento. O processo dos mínimos quadrados envolve a avaliacao das derivadas parciais em relaciio a cada uma de suas variáveis, impondo-se a nulidade das mesmas. Obtém-se assim um sistema de equacoes que pode ser escrito de uma forma matricial onde o índice i refere-se aos n pontos de superconverg~ncia do patch considerado para o nó j. O valor da tensa0 nodal do campo contínuo para este nó pode ser obtido diretamente por Para a determinaciio dos coeficientes, é necessário que, no mínimo, tres elementos concorram para o nó analisado. Para os nós de contorno, esta condiciio nem sempre é satisfeita, neste caso, utilizou-se a tensa0 média nodal. ESTRATÉGIA H-ADAPTATIVA PARA A ANÁLISE ELASTO-PLÁSTICA BIDIMENSIONAL 525 Refinamento Deseja-se obter uma malha 'h-adaptativa vza remeshzbg que atenda a um nível de erro relativo preestabelecido (q). Determinou-se entáo um parametro (t,), por elemento, que é capaz de: identificar os elementos que devem ser subdivididos; prever o número de subdivisoes necessárias (Ribeiro, Landaus). Sabe-se que o tamanho (h) do elemento decresce na taxa de 2nS, onde ns é o número de subdivisoes, e que, para elementos lineares o erro é proporcional ao tamanho. Ou seja, onde h, e E, representam, respectivamente, o tamanho e o erro médio de uma malha que satisfaz o erro relativo predeterminado V. Obtendo o valor de E, através de (18) e substituindo em (26) temos O número de subdivisoes por elemento foi determinado por: Caso a malha nao atinja o erro predeterminado, uma nova análise e, conseqüentemente, novas divisoes deveráo ser efetuadas. Para a análise numérica realizada, foi considerada uma viga curta em balanqo (Figura l), no estado plano de deformaqáo, discretizada através de elementos triangulares lineares. A curva tensáo-deforma~áo do material, aproximada por seis trechos lineares, encontra-se representada na Figura 2, onde também foram indicadas a espessura (t,) e a tensáo de escoamento (a,,) de cada subcamada. As análises adaptativas (estimativa de erro e refinamento) foram realizadas ao longo de dez incrementos de carga, sendo cada incremento igual a 0.1 (Aa = 0.1). Estas análises foram executadas em cada passo de carga (aF,), com exceqao do primeiro, a partir da malha final obtida no passo anterior. Para o primeiro passo utilizou-se a malha inicial indicada na Figura 4. Se a malha inicial for pouco refinada, provavelmente apenas uma análise nao será suficiente para que a malha atinja o erro relativo predeterminado. Isto porque o critério de refinamento utilizado limita o número de subdivisoes por análise. Observou-se que para o material elasto-plástico o erro relativo náo se mantém constante durante todos os R.C.C. SILVA, L. LANDAU E F.L.B. RIBEIRO passos de carga. De uma forma geral, este erro aumenta com a evolucao da plastificaqiio. Sendo assim, em algumas etapas da análise incremental o erro relativo preestabelecido foi ultrapassado, acarretando uma nova subdivisiio da malha. Estado plano de deformaqáo 1.0 m carregamento- Fo =l. 0 N/m2 módulo de elasticidade E =1.0 N/m coeficiente de Poisson- V =0.3 Figura 1. Viga curta em balanco Overlays 1 2 3 4 5 6 ti 0.5357 0.2862 0.0604 0.0594 0.0293 0.0290 o,i 1.2 3.6 8.6 14.6 23.6 75.0 I O 10 20 30 40 50 60 70 80 deforrnag30 Figura 2. Curva tensao-deforrnacáo aproximada por seis overlays Figura 3. Curvas carga-deslocamento Estao indicados na Tabela 1, para todos os passos de carga, os erros relativos obtidos durante cada análise realizada, objetivando satisfazer o erro relativo predeterminado de 10%. O número de nós (NN) e o número de elementos (NE) para as malhas finais (última análise) encontram-se na mesma tabela. ESTRATÉGIA H-ADAPTATIVA PARA A ANÁLISE ELASTO-PLÁSTICA BIDIMENSIONAL 527 Outras duas análises foram realizadas com o objetivo de atingir erros relativos de 35% e 25%. Encontram-se representadas na Figura 3 as curvas carga-deslocamento para as malhas adaptativas finais (Figura 4) e para a malha inicial. Foi considerado o deslocamento (S) na diresilo vertical do nó situado no canto superior direito da viga curta. TABELA 1. Erros relativos (%) e malhas obtidas durante a análise adaptativa com erro relativo de 10% Malha inicial- NN=9 NE=8 Malha 02- q=25% NN=94 NE=158 Figura 4. Malha inicial e malhas adaptativas finais