Embed Size (px)
Citation preview
ANÁLISE DE MODELOS REOLÓGICOS VISCOELÁSTICOS ATRAVÉS DE
FORMULAÇÕES MISTAS EM ELEMENTOS FINITOS
João Paulo Lima Santos
DISSERTAÇÃO SUBMETIDA AO CORPO DOCENTE DA COORDENAÇÃO DOS
PROGRAMAS DE PÓS-GRADUAÇÃO DE ENGENHARIA DA UNIVERSIDADE
FEDERAL DO RIO DE JANEIRO COMO PARTE DOS REQUISITOS
NECESSÁRIOS PARA A OBTENÇÃO DO GRAU DE MESTRE EM CIÊNCIAS EM
ENGENHARIA CIVIL.
Aprovada por:
____________________________________________
Prof. Webe João Mansur, Ph.D.
____________________________________________ Prof. Eduardo Gomes Dutra do Carmo, D.Sc.
____________________________________________ Prof. José Antônio Fontes Santiago, D.Sc.
____________________________________________ Prof. Luiz Fernando Taborda Garcia, D.Sc.
____________________________________________ Prof. Marco Aurélio Chaves Ferro, D.Sc.
RIO DE JANEIRO, RJ – BRASIL
FEVEREIRO DE 2008
ii
SANTOS, JOAO PAULO LIMA
Análise de modelos reológicos viscoelásticos
através de formulações mistas em elementos finitos
[Rio de Janeiro] 2008
XII, 123 p. 29,7 cm (COPPE/UFRJ, M.Sc.,
Engenharia Civil, 2008)
Dissertação - Universidade Federal do Rio de
Janeiro, COPPE
1. Viscoelasticidade
2. Modelos reológicos
3. Elementos finitos
4. Métodos estabilizados
I. COPPE/UFRJ II. Título ( série )
iii
À minha família.
iv
AGRADECIMENTOS
A todos que colaboraram ao longo de minha trajetória e nas experiências vividas
desde a minha querida cidade do sertão de Alagoas até o Rio de Janeiro, na conclusão
desta etapa na vida.
Aos meus pais Francisco e Lourdes, pelo imensurável apoio e por ensinar que
sempre devemos ir em busca de nossos objetivos, com muito trabalho e perseverança.
Aos meus irmãos Fábio, Cássia e familiares pelo incentivo em encarar esse
desafio.
À minha namorada Jordana pelo apoio e compreensão constantes, ajudando a
tornar amenos os momentos mais temerosos.
Aos amigos de república, por tornar a estadia no Rio de Janeiro mais agradável,
em especial ao Daniel pelo auxílio na revisão de texto desta dissertação.
Aos colegas do LAMEC: Ivone, Fernanda B., Fernanda M., Tilene, Flávio e
demais colegas. Agradecimento especial ao capitão Vasconcellos e ao Felipe pelas
contribuições na concretização deste trabalho.
Aos professores do PEC pelo aprendizado, especialmente ao professor Marco
Ferro pelo pontapé inicial em um dos temas da dissertação.
Ao professor orientador Eduardo pela paciência e colaboração nos fundamentos
matemáticos utilizados no desenvolvimento do trabalho.
Ao professor orientador Webe Mansur pelas discussões temáticas essenciais no
desenvolvimento do trabalho e pela extrema confiança e conselhos prestados ao longo
do curso.
Ao CNPq pelo apoio financeiro.
v
Resumo da Dissertação apresentada à COPPE/UFRJ como parte dos requisitos
necessários para a obtenção do grau de Mestre em Ciências (M.Sc.)
ANÁLISE DE MODELOS REOLÓGICOS VISCOELÁSTICOS ATRAVÉS DE
FORMULAÇÕES MISTAS EM ELEMENTOS FINITOS
João Paulo Lima Santos
Fevereiro/2008
Orientadores: Webe João Mansur
Eduardo Dutra Gomes do Carmos
Programa: Engenharia Civil
Os materiais viscoelásticos vêm gradativamente atraindo interesse industrial e
adquirindo maior espaço entre aplicações práticas do cotidiano, o que tem estimulado o
interesse no desenvolvimento de estratégias de modelagem computacional desse tipo de
material. Nesse trabalho, aborda-se uma estratégia numérica para representação de
modelos reológicos para sólidos e fluidos viscoelásticos lineares, baseada no método
dos elementos finitos. Na estratégia empregada, toma-se como base um modelo integral
para representação do sólido viscoelástico padrão e o modelo de Maxwell,
transformando-os em um esquema recursivo passo-a-passo, onde a história temporal é
obtida por uma seqüência de solução de modelos lineares atualizados a cada etapa pelas
variáveis controladoras da viscoelasticidade. Faz-se uso de formulações mistas em
elementos finitos para o tratamento de problemas incompressíveis, empregando-se o
método Galerkin Least Square (GLS) e o método de estabilização direta de pressão
(Bochev). Propõe-se também uma estratégia para determinação do parâmetro de
estabilização para o modelo GLS, que representa atualmente uma das maiores
dificuldades no emprego do método. Para verificação dos modelos implementados,
diversas análises numéricas foram realizadas tanto para os modelos lineares quanto para
os modelos viscoelásticos.
vi
Abstract of Dissertation presented to COPPE/UFRJ as a partial fulfillment of the
requirements for the degree of Master of Science (M.Sc.)
ANALYSIS OF VISCOELASTIC RHEOLOGICAL MODELS USING FINITE
ELEMENTS MIXED FORMULATIONS
João Paulo Lima Santos
February/2008
Advisors: Webe João Masur
Eduardo Gomes Dutra do Carmo
Department: Civil Engineering
Viscoelastic materials are attracting industrial interest gradually and acquiring
larger space among practical applications of daily engineering design, thus stimulating
the development of computational modeling strategies considering this kind of material.
In this work, a numerical strategy is discussed for representation of rheological models
for linear viscoelastic solids and fluids, based on the finite element method. The used
strategy is based on the Maxwell viscoelastic model represented in the integral form,
transformed into a step-by-step recursive scheme, where the time history is obtained
through a sequency of linear models updated in each stage by the controlling variables
of the viscoelasticity. Finite elements mixed formulations are used for the treatment of
incompressible problems, being employed: the Galerking Least Square method (GLS)
and the method of direct stabilization of pressure (Bochev). It also proposed a strategy
for determination of the parameter of stabilization for the GLS model, that represents
one of the greatest difficulties nowadays in the use of the method. For verification of the
implemented models, several numerical analyses were accomplished so much for the
linear models as for the viscoelastic models.
vii
Índice 1. Introdução................................................................................................................1
2. Modelos reológicos viscoelásticos...........................................................................5
2.1. O comportamento elástico ................................................................................ 5
2.2. O comportamento viscoso ................................................................................ 6
2.3. O fenômeno de fluência (creep) ....................................................................... 7
2.4. O fenômeno de relaxação ................................................................................. 9
2.5. A viscoelasticidade e o fenômeno de fluência e relaxação ............................ 10
2.6. Modelo molecular........................................................................................... 14
2.7. Modelos físicos/matemáticos ......................................................................... 16
2.7.1 Modelo de Maxwell................................................................................ 17
2.7.2 Modelo de Kelvin ................................................................................... 21
2.7.3 Modelo de Burgers ................................................................................. 24
2.7.4 Modelo do sólido linear padrão.............................................................. 25
2.7.5 Generalização dos modelos básicos ....................................................... 27
2.8 Representação por integrais hereditárias ........................................................ 30
3. O método dos elementos finitos ............................................................................34
3.1 Conceitos básicos ........................................................................................... 34
3.2 Formulação mista e incompressibilidade ....................................................... 41
3.3 Condições de estabilidade na formulação (u-p) para a equação de Stokes .... 49
3.4 Métodos estabilizados .................................................................................... 53
3.4.1 Galerkin Least Square (GLS) ................................................................. 54
3.4.2 Estabilização direta de pressão (Bochev) ............................................... 57
4. IGLS: Uma estratégia para determinação do parâmetro de estabilização α do
método GLS ...................................................................................................................61
5. Modelo viscoelástico incompressível....................................................................66
5.1. Extensão dos modelos viscoelásticos unidimensionais .................................. 66
5.2. Formulação em elementos finitos para a viscoelasticidade linear
incompressível ............................................................................................................ 69
5.3. Estratégia numérica ........................................................................................ 75
6. Aplicações e discussões..........................................................................................81
6.1 Exemplo 1: Viga submetida à ação de momento de flexão............................ 81
viii
6.2 Exemplo 2: Escoamento de Stokes em regime estacionário .......................... 91
6.3 Exemplo 3: Viga viscoelástica submetida à ação de momento de flexão ..... 98
6.4 Exemplo 4: Membrana de Cook................................................................... 100
6.5 Exemplo 5: Análise de reversibilidade da fluência viscoelástica................. 107
7. Conclusões e recomendações ..............................................................................117
8. Referências bibliográficas...................................................................................120
ix
Lista de figuras Figura 2. 1 – Comportamento elástico linear ....................................................................6
Figura 2. 2 – Comportamento viscoso linear.....................................................................6
Figura 2. 3 – Curva típica experimental de fluência .........................................................8
Figura 2. 4 – Aspecto da curva de fluência em escala logarítmica ...................................8
Figura 2. 5 – Teste de relaxação........................................................................................9
Figura 2. 6 - Fluência em material viscoelástico: (A) sólido; (B) fluido .......................10
Figura 2. 7 – Correção da curva de fluência....................................................................12
Figura 2. 8 – Superposição tempo-temperatura...............................................................13
Figura 2. 9 – Relação entre a taxa de conformação (rate) e o tempo..............................16
Figura 2. 10 – Modelo de Maxwell .................................................................................17
Figura 2. 11 – Resposta ao teste de fluência para o elemento de Maxwell .....................19
Figura 2. 12 – Configuração típica de um corpo submetido a um teste de relaxação .....21
Figura 2. 13 – Representação do modelo de Kelvin........................................................22
Figura 2. 14 – Resposta ao teste de fluência para o elemento de Kelvin ........................23
Figura 2. 15 – Modelo reológico de Burgers...................................................................24
Figura 2. 16 – Elemento sólido viscoelástico linear........................................................25
Figura 2. 17 – Elemento de Maxwell generalizado em série ..........................................28
Figura 2. 18 – Elemento de Maxwell generalizado em paralelo (Maxwell-Weichert) ...28
Figura 2. 19 – Elemento de generalizado Kelvin Voigt em paralelo ..............................29
Figura 2. 20 – Modelo generalizado de Kelvin-Voigt em série. ....................................29
Figura 2. 21 – Variação da deformação no tempo...........................................................31
Figura 3. 1 – Caracterização do domínio e contorno.......................................................36
Figura 3. 2 – Travamento volumétrico da malha. ...........................................................42
Figura 3. 3 – Elemento Q1-P0 – adaptado de BATHE [32] ............................................51
Figura 3. 4- Elemento Q2r-P0 (8/1 em 2D e 20/1 em 3D): adaptado de BATHE [32] ...52
Figura 3. 5 - Elemento Q2-P1 (9/3 em 2D e 27/4 em 3D): adaptado de BATHE [32] ....52
Figura 3. 6 - Elemento P2-P0 (6/1 em 2D e 10/1 em 3D): adaptado de BATHE [32] .....52
Figura 3. 7 - Elemento P2+-P1 (7/3 em 2D e 11/4 em 3D): adaptado de BATHE [32] ...53
Figura 3. 8 – Elemento Q2-Q1 (9/4-c em 2D e 27/8-C em 3D): adaptado de
BATHE [32] ................................................................................................................... 53
x
Figura 3. 9 - Elemento P2-P1 (6/3-c em 2D e 10/4-c em 3D): adaptado de
BATHE [32] ................................................................................................................... 53
Figura 5. 1 – Modelo viscoelástico linear padrão............................................................68
Figura 5. 2 – Modelo sólido viscoelástico linear generalizado .......................................68
Figura 5. 3 – Elemento quadrilátero bilinear...................................................................76
Figura 5. 4 – Função de forma bilinear ...........................................................................77
Figura 6. 1 – Viga submetida à ação de um momento de flexão (adaptada de REDDY et.
al.[51]) ............................................................................................................................ 82
Figura 6. 2 – Domínio discretizado em 512 elementos (32x16) .....................................82
Figura 6. 3 – Convergência do deslocamento horizontal para ν = 0,499999999. ...........83
Figura 6. 4 – Convergência do deslocamento vertical para ν = 0,499999999.................83
Figura 6. 5 – Convergência do deslocamento horizontal para ν = 0,1. ...........................84
Figura 6. 6 – Convergência do deslocamento vertical para ν = 0,1.................................84
Figura 6. 7 – Deslocamento horizontal em função do Poisson (512 elementos) ............85
Figura 6. 8 – Deslocamento vertical em função do Poisson (512 elementos).................85
Figura 6. 9 – Deslocamento resultante para ν = 0,499999999. .......................................87
Figura 6. 10 – Deslocamento resultante para ν = 0,1. .....................................................87
Figura 6. 11 – Distribuição de pressão para ν = 0,499999999. .......................................87
Figura 6. 12 – Convergência para o deslocamento horizontal (ν = 0,499999999)..........88
Figura 6. 13 – Convergência para o deslocamento vertical (ν = 0,499999999)..............88
Figura 6.14 – Modo pressão espúria para (ν = 0,499999999): (a) 46 elementos;
(b) 85 elementos; (c) 191 elementos; (d) 466 elementos................................................ 90
Figura 6. 15 – Pressão estabilizado pela formulação IGLS (ν = 0,499999999)..............90
Figura 6. 16 – Caracterização do problema de escoamento de Stokes............................91
Figura 6. 17 – Variação da pressão na seção de corte BB...............................................92
Figura 6. 18 – Distribuição espacial da pressão pelo método GLS para α=0.02.............93
Figura 6. 19 – Distribuição espacial da pressão pelo método GLS para α=1.0...............93
Figura 6. 20 – Distribuição espacial da pressão pelo método GLS para α=20................93
Figura 6. 21 – Distribuição espacial da pressão pelo método Bochev para α=2.............94
Figura 6. 22 – Variação da velocidade vertical na seção de corte BB ............................95
Figura 6. 23 – Distribuição espacial da velocidade vertical obtido com o modelo GLS
para: (a) α=0.1; (b) α=0.5; (c) α=1.0; (d) α=20.0; .......................................................... 95
Figura 6. 24 – Variação da velocidade horizontal na seção de corte BB ........................96
xi
Figura 6. 25 – Distribuição espacial da velocidade horizontal obtido com o modelo GLS
para: (a) α=0.1; (b) α=0.5; (c) α=1.0; (d) α=20.0; ...........................................................96
Figura 6. 26 – Distribuição espacial do módulo de velocidade obtido com o modelo
GLS para: (a) α=0.1; (b) α=0.5; (c) α=1.0; (d) α=20.0;.................................................. 97
Figura 6. 27 – Discretização espacial em 128 elementos finitos.....................................98
Figura 6. 28 – Evolução temporal do deslocamento .......................................................99
Figura 6. 29 – Geometria e discretização espacial em 100 elementos (10 x 10) ..........100
Figura 6. 30 – Deslocamento vertical para a malha 10 x 10 (100 elementos). .............101
Figura 6. 31 – Deslocamento vertical para a malha 50 x 50 (2500 elementos). ...........101
Figura 6. 32 –(a) Malha de 900 elementos (b) malha de 2500 elementos.....................102
Figura 6. 33 – Distribuição da pressão para (a) α = 0 (b) α = 0,5. ...............................102
Figura 6. 34 – Estabilização da pressão para: (a) α = 0,1 ; (b) IGLS ............................103
Figura 6. 35 – Evolução temporal do deslocamento vertical (malha 30 x 30)..............104
Figura 6. 36 – Evolução temporal do deslocamento vertical (malha 50 x 50)..............105
Figura 6. 37 – Deslocamento resultante para t = 0. .......................................................106
Figura 6. 38 – Deslocamento resultante para t = 20. .....................................................106
Figura 6. 39 – Geometria e discretização da malha.......................................................107
Figura 6. 40 – Deslocamento horizontal para o modelo sólido viscoelástico padrão. ..108
Figura 6. 41 – Pressão para o modelo sólido viscoelástico padrão ...............................109
Figura 6. 42 – (a) Distribuição espacial da pressão para t = 0; (b) Evolução temporal da
pressão .......................................................................................................................... 110
Figura 6. 43 – Distribuição espacial da Pressão (estabilizada) .....................................110
Figura 6. 44 – Deslocamento resultante para t = 2,4 para: (a) ν = 0,5; (b) ν = 0,3. .....111
Figura 6. 45 – Deslocamento horizontal para o modelo de Maxwell............................112
Figura 6. 46 – Pressão- Modelo de Maxwell.................................................................112
Figura 6. 47 – Corpo submetido a estado de cisalhamento puro...................................113
Figura 6. 48 – Modo pressão espúria no estado de cisalhamento puro. ........................113
Figura 6. 49 – Deslocamento resultante para o GLS α = 0. ..........................................114
Figura 6. 50 – Evolução temporal do deslocamento vertical (Modelo de Maxwell). ...114
Figura 6. 51 – Evolução temporal da pressão (Modelo de Maxwell). .........................115
Figura 6. 52 – Evolução temporal do deslocamento vertical e pressão (Modelo do sólido
viscoelástico). ................................................................................................................115
Figura 6. 53 – Deslocamento resultando para o modelo de Maxwell em t = 100:........116
(a) ν = 0,2; (b) ν = 0,5....................................................................................................116
xii
Lista de tabelas Tabela 6. 1 – Resultado para o modelo de Bochev ....................................................... 86
Tabela 6. 2 – Resultado para o modelo GLS ................................................................ 86
Tabela 6. 3 – Resultado para o FEAP............................................................................ 86
Tabela 6. 4 – Resultados para a malha de 46 elementos (ν = 0,499999999) ................ 89
Tabela 6. 5 – Resultados para a malha de 85 elementos (ν = 0,499999999) ................ 89
Tabela 6. 6 – Resultados para a malha de 466 elementos (ν = 0,499999999) .............. 89
1
1. Introdução
Os estudos do comportamento reológico dos materiais remotam aos primórdios
da humanidade. Para o filósofo grego Heráclito de Éfeso (540-476 a.C.), a lei
fundamental do universo é dividir, o que significa contínuas transformações. Para ele,
todos os seres e a matéria estão em constante movimento, afirmando: “Tudo flui, nada
neste nosso mundo é permanente. Tudo está mudando o tempo todo. A mudança é a lei
da vida e do universo”.
O filósofo parecia observar que boa parte dos materiais quando submetidos a
estágios de solicitações distintas também seguiam os processos regidos por sua lei
fundamental, e já parecia impulsionar a origem dos estudos reológicos.
Nos séculos seguintes, a evolução dos estudos do comportamento reológico foi
descrito pelas “relações constitutivas”, que relacionavam a tensão e deformações
sofridas pelos meios em estudo. Dos trabalhos pioneiros sobre os materiais ideais
podem-se destacar: Archimedes e Newton (1687) sobre corpos perfeitamente rígidos;
Boyle (1600), Hook(1678), Young(1807) e Cauchy(1827) sobre materiais perfeitamente
elásticos; Pascal (1963), Bernoulli(1738) e Euller(1755) sobre fluidos imiscíveis;
Newton(1687), Navier(1823) e Stokes (1845) sobre a dinâmica dos fluidos
newtonianos. O termo reologia foi usado pela primeira vez em 1929, definida então
como o estudo do fluxo e deformação de todas as formas de matéria, apesar de alguns
autores recentemente limitarem o termo reologia ao estudo de materiais que apresentem
comportamento diferente dos sólidos perfeitamente elásticos e fluidos ideais. Maiores
detalhes sobre a origem e evolução da reologia são apresentados de forma resumida por
DORAISWAMY [1].
As relações constitutivas desenvolvidas até meados do século XVIII não
conseguiam representar com precisão o comportamento de determinados materiais, em
especial alguns que apresentam um comportamento intermediário entre o sólido
perfeitamente elástico e o fluido newtoniano, que ficaram conhecidos posteriormente
como materiais viscoelásticos.
Dessa forma, pode-se definir a viscoelasticidade como a propriedade de um
material exibir tanto comportamento elástico quanto viscoso. Os materiais
viscoelásticos têm características tanto de sólido (elasticidade, resistência ao fluxo e
2
estabilidade da forma) quanto de fluido, tais como a dependência do fluxo com o tempo,
temperatura e tensão aplicada.
Os materiais viscoelásticos são encontrados na indústria alimentícia, em
variedades de polímeros, em geomateriais, em materiais biológicos, em compósitos,
dentre outras diversas variedades. MACKERLE [2] afirma que as áreas de estudos e de
aplicações da viscoelasticidade são as mais variadas possíveis, a exemplo da engenharia
aeroespacial, engenharia biomédica, engenharia civil, engenharia de alimentos e
engenharia mecânica e de materiais.
De acordo com DORAISWAMY [1], os estudos iniciais dos materiais
viscoelásticos foram impulsionados pelo crescimento da indústria dos polímeros.
Webber (1835) foi um dos pioneiros nesse novo campo da reologia, estudando o
comportamento dos fios de seda, na época empregados como componentes de alguns
instrumentos eletromagnéticos.
Webber observou que após tensionar axialmente o fio de seda, observa-se não só
uma deformação instantânea elástica, mas também uma deformação que progredia em
função do tempo. Liberando-o em seguida do esforço de tração, havia uma contração
imediata do fio, que atingia o estado da deformação elástica, seguida por uma contração
gradual até atingir o estado de pré-carregamento inicial. Webber havia então
experimentado o fenômeno que ficou posteriormente conhecido como relaxação de
tensão, que ele chamou de “the after effect”, caracterizando dessa forma a dependência
temporal dessa até então nova classe de materiais.
Nos anos seguintes, cientistas e engenheiros dedicaram estudos experimentais e
matemáticos para obter uma representação das relações constitutivas para a
viscoelasticidade. Destacam-se os trabalhos desenvolvidos por Lord Kelvin e James
Clerk Maxwell. Demais contribuições podem ser consultadas nas referências
DOZDROV [3], FLÜGGE [4], BLEND [5] e CHRISTENSEN [6].
Os estudos numéricos aprofundados dos materiais viscoelásticos foram
desenvolvidos a partir da década de 60, sobretudo impulsionados pelo crescente
interesse industrial sobre os materiais poliméricos. Já nesse período, os métodos de
aproximação das equações diferenciais vinham sendo largamente empregados para a
representação de fenômenos inerentes à reologia, especialmente o método dos
elementos finitos (MEF). Do ponto de vista numérico, a linha de trabalho empregando o
método dos elementos finitos (MEF) na viscoelasticidade seguiu dois caminhos
distintos:
3
No primeiro, que inclui um maior número de trabalhos desenvolvidos, o
problema é resolvido diretamente no domínio do tempo usando um processo
incremental passo-a-passo. Faz-se uso de relações constitutivas nas formas integrais e
diferenciais, a exemplo dos trabalhos desenvolvidos por ZIENKIEWICZ et. al. [7] e
SRIBATHA e LEWIS [8].
No segundo tipo aplica-se o chamado princípio da correspondência elástico-
viscoelástico, discutido inicialmente por ALFREY [9], onde os problemas
viscoelásticos são transformados em elásticos no domínio de Laplace. Em seguida, o
problema elástico é resolvido via método dos elementos finitos. Depois disso, aplicando
a técnica de transformada inversa de Laplace, obtém-se a solução do problema
viscoelástico inicial. Empregando esse método, ADEY e BREBIA [10] analisaram um
cilindro espesso oco de aço sob pressão interna. Alguns autores, a exemplo de
CARPENTER [11] comentam que o emprego da técnica baseada no princípio da
correspondência pode tornar-se inviável em casos de configurações geométricas
complexas e em casos não lineares.
Um dos trabalhos pioneiros nessa linha foi desenvolvido por WEBBER [12],
que aplicou o princípio da correspondência para problemas termo-viscoelásticos
baseados no modelo de Maxwell, empregando o MEF na formulação clássica, ou seja,
empregando a sentença variacional apenas em termos de deslocamento/velocidade.
Os trabalhos pioneiros que empregavam diretamente a solução no domínio do
tempo eram bastante simplificados. CARPENTER [11] discorre que no trabalho
desenvolvido por WHITE [13] assumia-se o módulo de elasticidade volumétrica
constante no tempo e o material homogêneo e isotrópico.
Fortemente influenciados pela temperatura e pela alteração das propriedades do
material ao longo do tempo, observou-se a necessidade de se incorporar formulações
estáveis que permitissem análises de meios próximos à condição de incompressibilidade
(elevados valores do coeficiente de Poisson), uma vez que a formulação de elementos
finitos baseada apenas na sentença variacional do deslocamento (ou velocidade) pode
induzir a severas oscilações na aproximação da pressão e no campo da
tensão/deformação, além de sofrer o efeito de “travamento” volumétrico da malha. Um
dos trabalhos pioneiros foi proposto por YADAGIRI e REDDI [14], que fizeram uso de
elementos isoparamétricos associados a esquemas de integração reduzida. Nos anos
seguintes, algumas formulações mistas foram sendo empregadas e ainda hoje diversas
4
estratégias são testadas, sobretudo para um maior ganho de precisão atrelado à busca
por formulações com menor custo computacional.
O presente trabalho objetiva a análise e implementação de modelos reológicos
viscoelásticos que permitam a representação de meios compressíveis e incompressíveis.
Para tanto, investigam-se algumas formulações estabilizadas em elementos finitos,
sobretudo para a garantia de estabilidade numérica em situações de incompressibilidade.
Uma breve revisão bibliográfica relacionada aos modelos viscoelásticos é
apresentada no capítulo 2, destacando-se parte dos modelos conceituais mais difundidos
na literatura.
A abordagem de formulações em elementos finitos para modelos
incompressíveis é apresentada no capítulo 3, enfatizando principalmente o uso de
formulações mistas. Já no capítulo 4 é discutida uma nova estratégia para determinação
do parâmetro de estabilização para o método de aproximação misto em elementos
finitos GLS (Galerkin Least Square).
No capítulo 5 são abordados alguns aspectos relacionados à implementação
computacional e à estratégia numérica para representação de modelos viscoelásticos
incompressíveis.
Algumas aplicações são apresentadas no capítulo 6, sobretudo para análise da
acurácia dos modelos implementados. As considerações finais e sugestões para
trabalhos futuros são apresentadas no capítulo 7.
5
2. Modelos reológicos viscoelásticos
Os materiais viscoelásticos são caracterizados pela observância de um
comportamento intermediário entre o sólido elástico e o fluido newtoniano. DROZDOV
[3] afirma que para um material ser considerado viscoelástico, é necessário que o
mesmo experimente os fenômenos de fluência e relaxação de tensões. Esses fenômenos
estão associados à combinação dos efeitos elásticos e viscosos, podendo ser
representados por diversos modelos reológicos cujo comportamento peculiar é função
do arranjo em associação dos elementos básicos: mola (representando a parcela elástica)
e amortecedor (representado a parcela viscosa), que podem caracterizar sólidos
viscoelásticos ou fluidos viscoelásticos.
Os estudos discutidos nesse trabalho serão restritos a pequenas deformações,
caracterizando apenas os aspectos da viscoelasticidade linear.
Para uma melhor caracterização dos fenômenos atrelados ao comportamento
viscoelástico, será apresentado nesse capítulo uma breve descrição dos fenômenos
inerentes à viscoelasticidade.
2.1. O comportamento elástico
Os materiais são classificados em elásticos quando após aplicação de um
carregamento interino e conseqüente deformação observa-se a recuperação de sua forma
original, apresentando assim uma deformação reversível, havendo então conservação de
sua energia interna. Nesse material, o estado de deformação é atingido imediatamente
após a aplicação da carga.
Um material básico elástico é geralmente representado por uma mola, conforme
mostrado na figura (2.1a). Aplicando-se uma tensão constante 0)( σσ =t (fig 2.1b), é
dito que o corpo apresenta comportamento elástico linear caso se observe uma
deformação constante e proporcional à tensão aplicada (fig 2.1c).
6
Figura 2. 1 – Comportamento elástico linear
Os materiais elásticos lineares obedecem à lei de Hooke, sendo a constante de
proporcionalidade E uma característica intrínseca do material, conhecida como módulo
de elasticidade. A relação constitutiva para o material elástico linear é então definida
pela seguinte relação:
( ) . ( )et E tσ ε= (2. 1)
2.2. O comportamento viscoso
A viscosidade de um corpo descreve a resistência que um material oferece ao
escoamento. O efeito da viscosidade é mais evidente nos fluidos, que por definição são
substâncias que se deformam continuamente sob ação de qualquer força tangencial.
Nesta definição, não é levada em conta a estrutura molecular do fluido, que é composto
de diversas moléculas em movimento.
Figura 2. 2 – Comportamento viscoso linear
Um material viscoso é usualmente representado pelo elemento da figura (2.2a).
Aplicando-se nesse elemento uma tensão constante 0)( σσ =t (fig 2.2b), é dito que o
elemento se comporta como fluido newtoniano se sua taxa de deformação é diretamente
proporcional à tensão aplicada:
7
( ) . ( )it tσ η ε•
= (2. 2)
Nesse caso, a constante de proporcionalidade η é conhecida como viscosidade, e
seu valor irá determinar o grau de resistência ao cisalhamento do fluido, sendo um
parâmetro característico para a descrição da curva de deformação apresentada na figura
(2.2c).
2.3. O fenômeno de fluência (creep)
A fluência, também conhecida na literatura por sua terminologia inglesa (creep),
está relacionada à tendência das partículas constituintes dos materiais sofrerem
movimentos conseqüentes da aplicação continuada de cargas de intensidade constante.
Para um melhor entendimento do fenômeno, discute-se o experimento oriundo do teste
de fluência, apresentado por DROZDOV [3]. Nesse teste, idealiza-se um corpo
inicialmente em seu estado natural (sem carregamento). No instante t = 0 é aplicada
uma força de intensidade P em uma seção de área S obtendo-se uma tensão 0σ . Em
resposta à aplicação da carga, uma deformação 0ε é observada de forma instantânea, e
no caso dos materiais perfeitamente elásticos essa deformação é mantida constante para
t > 0. Entretanto, para materiais viscoelásticos, observa-se uma deformação incremental
)(tε , que em geral apresenta tendência crescente ao longo do tempo, sendo determinada
por uma taxa de variação 0)( ≥•
tε . Esse fenômeno de alteração no valor da deformação
sob condições de carregamento uniforme no tempo é chamado de fluência.
Segundo DROZDOV [3], o processo de fluência pode ser caracterizado pela
deformação de fluência )(tcε , que é a diferença entre a deformação no instante )(tε e a
deformação oriunda do processo elástico 0ε . Entretanto, a representação mais comum é
feita através da chamada função de fluência, que corresponde a uma relação entre a
deformação de fluência e a tensão aplicada 0σ , conforme relação apresentada a seguir:
0
0 0
( ) ( )( ) c t tJ t ε ε εσ σ
−= =
(2. 3)
Na figura (2.3) é apresentada uma curva típica de fluência para um sólido
viscoelástico. Em função de resultados experimentais, em geral ajustam-se as
8
expressões analíticas para a curva de fluência, usualmente sob forma exponencial, onde
J corresponde ao módulo de fluência.
Figura 2. 3 – Curva típica experimental de fluência
(adaptado de DROZDOV [3])
De acordo com ROYLANCE [15], a curva de fluência quando plotada em escala
logarítmica apresenta um ponto de inflexão, que caracteriza o “tempo de relaxação” τ
do processo de fluência, conforme apresentado na figura (2.4). Os valores Jg e Jr
representam, respectivamente, a fluência do material na “fase vítrea (glassy)” e “fase
elástica (ruberry)”, que serão discutidas na seção 2.5.
Figura 2. 4 – Aspecto da curva de fluência em escala logarítmica
(adaptado de ROYALANCE [15])
9
2.4. O fenômeno de relaxação
Define-se como relaxação de tensões o fenômeno de alívio nos esforços internos
que ocorre em alguns materiais quando sujeitos a uma deformação constante. Para
melhor entendimento do fenômeno, será apresentada a experiência oriunda do teste de
relaxação. Nesse experimento, um corpo é subitamente posto em situação de
tracionamento axial e conseqüentemente é observado um incremento em seu
comprimento, induzindo a um estado de deformação axial que é mantido constante para
efeitos do teste ao longo do tempo. O comportamento dos valores associados à tensão e
deformação são ilustrados na figura (2.5).
Figura 2. 5 – Teste de relaxação
Para um material perfeitamente elástico a tensão também se manteria constante
ao longo da análise. Entretanto, para um material viscoelástico esse valor diminui
monotonicamente a partir do tempo t >0, tendendo a um valor assintótico limσ . A curva
do comportamento da tensão no tempo apresentada pela figura (2.5) é conhecida como
curva de relaxação, e é geralmente obtida através de ensaios experimentais, onde em
seguida ajustam-se funções para caracterização do comportamento, sendo estas
conhecidas como funções de relaxação. Um parâmetro importante associado aos ensaios
de relaxação é, portanto, a função de relaxação, que determina o módulo de relaxação
em cada instante t, definido pela relação:
0
)()(εσ ttG =
(2. 4)
De acordo com ROYLANCE [15], a fluência e relaxação são resultantes do
mesmo mecanismo molecular, e por essa razão o valor associado à fluência J(t) e o
módulo de relaxação )(tG na maioria das situações podem ser correlacionados. No
caso de um material viscoelástico linear, esses valores são inversos entre si. Entretanto,
nem sempre essa relação é fácil de ser obtida nos demais materiais viscoelásticos, em
particular, porque a resposta ao teste de relaxação é obtida de forma mais rápida do que
o equilíbrio de forças obtido no teste de fluência.
10
2.5. A viscoelasticidade e o fenômeno de fluência e relaxação
Na literatura, a dependência temporal do campo de tensão e deformação advinda
dos fenômenos de fluência e relaxação confere aos materiais viscoelásticos a
denominação de materiais com memória, que podem ser subclassificados em fluidos
viscoelásticos e sólidos viscoelásticos.
ROYLANCE [15] apresenta alguns elementos diferenciadores do
comportamento de sólidos e fluidos viscoelásticos. De acordo com sua definição, um
fluido viscoelástico submetido a um estado de deformação constante (teste de
relaxação) irá produzir um estado de tensão que irá decair a zero. Contrariamente, os
sólidos viscoelásticos quando submetidos a um estado de deformação constante (teste
de relaxação) apresentarão uma componente de tensão que remanesce não-nula.
DROZDOV [3] afirma que para um meio sólido viscoelástico a deformação ( )tε tende
a um valor limite ε∞ e sua derivada ( )tε•
tende a zero quanto t →∞ . Entretanto, para
um fluido viscoelástico a taxa de deformação ( )tε•
tende para um valor limite ε•
∞ , que
corresponde a um fluxo newtoniano. Em outras palavras, um sólido viscoelástico
quando submetido ao teste de fluência apresentará uma deformação limitada, enquanto
que o fluido viscoelástico tenderá a uma deformação ilimitada cuja taxa de deformação
se aproximará à de um fluído newtoniano. Esse comportamento é apresentado na figura
(2.6).
Figura 2. 6 - Fluência em material viscoelástico: (A) sólido; (B) fluido
(adaptado de DROZDOV [3])
11
A distinção entre o comportamento de um fluido e de um sólido pode ser feito
através do número de Deborah. De acordo com LEAL [16], o número de Deborah, De, é
uma grandeza adimensional que permite distinguir, em função da sua magnitude, se um
determinado material tem um comportamento mais próximo de um fluido ou de um
sólido. Essa grandeza pode ser obtida pela relação tempo de relaxação/tempo de
observação. Pequenos valores do número de Deborah correspondem a materiais que têm
tempo suficiente para atingir um estado de relaxação e, por isso, comportam-se como
fluidos. Contrariamente, os materiais que apresentam um valor elevado do número de
Deborah estão mais próximos do comportamento de um material sólido.
DROZDOV [3] apresenta uma expressão para a deformação total ( )tε baseada
nas componentes da deformação de fluência, instantânea (elástica) e newtoniana
(viscosa). Para tal, considera-se um elemento viscoelástico submetido a uma tensão
uniforme 0σ , ao qual se pode associar uma deformação newtoniana nε proporcional ao
tempo (eq. 2.5), uma deformação inicial 0ε instantânea e uma componente deformação
de fluência )(~ tcε , que geralmente não é tão simples de ser determinada, correspondendo
a uma parcela variável ao longo do tempo. A deformação total deve corresponder à
soma de uma eventual deformação instantânea e de uma deformação de fluência,
resultante da soma das parcelas newtoniana e )(~ tcε variável, conforme equação (2.6).
0( )n t tσεη
= (2. 5)
)()( 0 tt cεεε +=
onde: )(~)()( ttt cnc εεε +=
(2. 6)
As propriedades de resposta inicial ao estímulo de forças também não é
uniforme nos materiais viscoelásticos. Baseado no trabalho de Giesekus (1995),
DROZDOV [3] apresenta uma classificação de sólidos e fluidos viscoelásticos em
função do comportamento mecânico inicial, conforme apresentado em seguida:
12
• Sólidos com viscosidade instantânea ( 0 0ε = ,ε∞ < ∞ , 0=∞
•
ε );
• Sólidos com elasticidade instantânea ( 0 0ε ≠ ,ε∞ < ∞ , 0=∞
•
ε );
• Fluidos com viscosidade instantânea ( 0 0ε = ,ε∞ = ∞ , ε•
∞ ≠ ∞ );
• Fluidos com elasticidade instantânea ( 0 0ε ≠ ,ε∞ = ∞ , ε•
∞ ≠ ∞ );
(2. 7)
Alguns polímeros apresentam modelo mecânico com comportamento diferente
dos modelos de fluência convencional dos sólidos e fluidos viscoelásticos. Segundo
DROZDOV [3], em alguns casos a deformação ( )tε em sólidos poliméricos aumenta
monotonicamente e não tende a um valor limitado finito, enquanto que sua taxa de
deformação tende a zero com o aumento do tempo.
De acordo com ROYLANCE [15], a temperatura influencia fortemente as
propriedades do material viscoelástico. Essa influência é observada em um
deslocamento relativo da curva de fluência. Com base nisso, pode-se fazer uma analogia
tempo x temperatura, efetuando-se uma correção no tempo de relaxação do material a
partir de um fator de troca Ta (eq. 2.8). Tal procedimento ocorre através da diferença
dos valores apresentados na curva de fluência, conforme figura (2.7), onde
refτ corresponde a temperatura de referência e τ a temperatura para o qual se deseja
efetuar a correção.
log log logref Taτ τ= + (2. 8)
Figura 2. 7 – Correção da curva de fluência.
(adaptado de ROYLANCE [15]).
13
Para determinar o fator de troca Ta , parte-se de um gráfico de referência que
apresente diferentes valores do módulo de fluência de um mesmo material para distintas
temperaturas. Deve-se então construir uma curva para a temperatura a ser corrigida
através de uma interpolação dos valores das curvas de temperaturas existentes (fig.
2.8a). Com base nos incrementos necessários para o ajuste, retorna-se à curva principal
corrigindo-se os valores deslocados no eixo ordenado, conforme procedimento
apresentado na figura (2.8).
Figura 2. 8 – Superposição tempo-temperatura
(adaptado de ROYLANCE [15).
O procedimento acima é válido apenas para materiais que possuam tempo de
relaxação simples. Para materiais viscoelásticos complexos, que podem apresentar
tempo múltiplo de relaxação (mais de um ponto de inflexão) devem-se efetuar correções
através das relações de Arhenius, conforme discutido por ROYLANCE [15].
Um outro fator que tem influência significativa no fenômeno da
viscoelasticidade é o efeito do envelhecimento do material, que pode trazer mudanças
nas propriedades viscoelásticas com o tempo. Normalmente o efeito do envelhecimento
induz ao aumento do módulo de elasticidade e à diminuição das deformações de
fluência. DROZDOV [3] argumenta que é conveniente distinguir os efeitos físicos
(reversíveis) e químicos (irreversíveis) do envelhecimento. Os efeitos das alterações
químicas, que são geralmente irreversíveis se processam através de reações químicas
sofridas pelos materiais constituintes e irão induzir a modificações mecânicas com o
tempo, a exemplo do processo de hidratação do concreto. Já os processos físicos estarão
associados ao equilíbrio termodinâmico, e tendem, por exemplo, a reduzir a energia de
entropia e por conseqüência as taxas de fluência e relaxação. Convém relembrar que a
entropia está relacionada à quantidade de energia no corpo que não pode ser
14
transformada em trabalho. Esses processos ocorrem de maneira peculiar a cada classe
de material viscoelástico, cabendo então uma análise pontual para sua caracterização.
2.6. Modelo molecular
Os principais mecanismos controladores da viscoelasticidade estão relacionados
à reacomodação de partículas da estrutura molecular constituinte. Esse comportamento
é explicado pela teoria molecular. Os mecanismos apresentados se baseiam no trabalho
“Engineering Viscoelasticity”, apresentada por ROYLANCE [15] para polímeros
viscoelásticos.
Segundo essa teoria, quando um material viscoelástico é submetido a um estado
de tensão, os comprimentos e ângulos das ligações moleculares são alterados, e os
átomos são movidos para uma posição de maior energia interna. Caso o material tenha
mobilidade molecular suficiente, rearranjos moleculares de larga-escala podem ser
observados, podendo haver grandes mudanças na conformação molecular.
A combinação da primeira e segunda leis da termodinâmica, encarada a partir de
um incremento no trabalho mecânico fdx , faz interpretar que, para o equilíbrio do
sistema, deverá ocorrer um aumento na energia interna dU ou uma diminuição na
energia de entropia TdS, regida pela seguinte relação:
fdx dU TdS= − (2. 9)
De acordo com a formulação apresentada, a entropia dS é fortemente
influenciada pela temperatura T, e pode ser uma forma conveniente de representar
experimentalmente a relação entre as propriedades físicas do material e a energia de
entropia. A entropia está relacionada com o número de configurações (ou arranjos) de
mesma energia que um dado sistema pode assumir. A interpretação molecular da
entropia sugere que, em uma situação puramente geométrica, quanto maior o número de
configurações, maior a entropia [17].
Nesse contexto, nota-se que a força de retração interna necessária para um
eleastômero (borracha) alongado recuperar sua forma original deverá aumentar em
função da elevação de temperatura, caracterizando uma alta parcela de energia
entrópica. Por outro lado, no caso do aço (metal), essa força de retração irá diminuir
com a elevação da temperatura, ou seja, a expansão térmica atuará aliviando a tensão
15
interna, apresentando dessa forma uma baixa parcela de energia entrópica. A diferença
nesse comportamento é devida à mobilidade molecular, que está relacionada com a
energia entrópica.
Um mesmo material pode apresentar diferentes valores em sua força interna, que
dependerá da temperatura e conseqüentemente do grau de entropia, caracterizando fases
físicas diferenciadas. Na fase vítrea, atingida em baixas temperaturas, o deslocamento
molecular é relativamente pequeno em relação a sua configuração de origem e, por isso,
na realização do trabalho mecânico fdx observa-se um maior armazenamento da energia
elástica, refletido pela elevada rigidez do material e conseqüente menor deformação. Já
em temperaturas mais elevadas (fase elástica), as cadeias moleculares possuem maior
mobilidade, associadas a um maior valor da conformação molecular.
A mobilidade molecular é influenciada por fatores físicos e químicos, tal qual a
arquitetura molecular, a temperatura ou a presença de fluidos de absorção. O grau de
mobilidade molecular é medido pela taxa de conformação (rate), que pode ser
quantificada pela expressão de Arrhenius:
'exp( )ErateRT−
∝ (2. 10)
Nessa expressão, E’ representa a energia de ativação aparente e R é a constante
universal dos gases perfeitos, R=8,314J/mol-K. Em uma zona de transição do
comportamento vítreo-elástico, define-se uma temperatura de transição Tg, tal que para
T>>Tg (fase elástica) o valor da taxa de conformação (rate) é elevado, permitindo aos
materiais atingir reversibilidade completa de deformação. Para T<<Tg (fase vítrea) os
valores da taxa de conformação são pequenos e a resposta ao estímulo de forças é mais
lenta e, por essa razão, a recuperação de deformação é mais lenta. A região
intermediária é conhecida como “leathery” ou viscoelástica, pois sua resposta é
combinação dos efeitos viscosos das fases vítrea e elástica, conforme apresentado na
figura (2.9).
16
Figura 2. 9 – Relação entre a taxa de conformação (rate) e o tempo
Para EIRICH [18], o comportamento do sólido viscoelástico resulta do rearranjo
molecular que pode ser reversível ou irreversível, sugerindo que os principais
mecanismos controladores da fluência em sólidos viscoelásticos são: a impulsão
mecânica oriunda de movimentos diferenciais, o escorregamento das partículas e a
difusão molecular. No caso dos fluidos viscoelásticos, a viscoelasticidade é controlada
pelo entrelaçamento de cadeias moleculares com as cadeias adjacentes, sendo que a
dinâmica dessas interações controlam a função de fluência do material, podendo
inclusive não ser constante.
Em função dos efeitos viscosos atuantes nos materiais viscoelásticos, o efeito da
fluência introduz modificações nas propriedades intrísecas do material permitindo uma
nova acomodação molecular, que pode induzir a um comportamento de
incompressibilidade.
2.7. Modelos físicos/matemáticos
Para representar os efeitos viscosos e elásticos presentes nos materiais
viscoelásticos diversos modelos físico-matemáticos são encontrados na literatura e
geralmente são constituídos por diferentes arranjos geométricos de molas
(representando a parcela elástica) e amortecedores (representando a parcela viscosa).
Os principais modelos observados na literatura são apresentados em síntese nas
subseções seguintes.
17
2.7.1 Modelo de Maxwell
O modelo de Maxwell é constituído pela associação de um elemento elástico
ligado em série a um elemento viscoso, conforme ilustrado na figura (2.10). Esse
modelo foi proposto inicialmente por James Clerk Maxwell, em analogia aos modelos
elétricos.
Figura 2. 10 – Modelo de Maxwell
Aplicado-se uma tensão σ nas extremidades do elemento viscoelástico
apresentado na figura (2.10), observa-se que a deformação total corresponde à soma das
parcelas elástica e viscosa (eq. 2.11), representadas pelos índices e e i respectivamente.
Por outro lado, a tensão se distribui de forma constante nas parcelas elástica e inelástica,
conforme expressão (2.12).
ie ttt )()()( εεε += (2. 11)
ie ttt )()()( σσσ == (2. 12)
As relações constitutivas dos elementos componentes da parcela elástica e
viscosa assumem a seguinte forma:
( ) . ( )e et E tσ ε= (2. 13)
( ) . ( )i it tσ η ε•
= (2. 14)
Combinando as expressões anteriores chega-se a uma expressão diferencial para
o modelo de Maxwell, conforme apresentado na equação (2.15):
ie ttt )()()(•••
+= εεε ησσε )()()( t
Ett +=
••
(2. 15)
18
A partir da equação (2.15) pode-se estudar o comportamento dos materiais
viscoelásticos de Maxwell. Primeiramente, será efetuada uma análise a partir do teste de
fluência. Para tal, aplica-se uma tensão 0σ que se mantém constante para t >0,
conforme expressão:
)(.)( 0 tut σσ = (2. 16)
Nessa expressão acima, u(t) representa a função passo unitário (Heaviside) e δ a
função delta de Dirac. Admite-se inicialmente que o material não está sob ação de
carregamento, conforme condições iniciais apresentadas na expressão (2.17).
0
0
0
0,0;( ) 1,0;
0,5;
t tu t t t
t t
<⎧⎪= >⎨⎪ =⎩
; 0( ) 0; 0
( ) 1
tt t
t dt
δ
δ∞
−∞
⎧ ⎫⎪ ⎪∞ =⎪ ⎪⎪ ⎪= ≠⎨ ⎬⎪ ⎪⎪ ⎪=⎪ ⎪⎩ ⎭∫
0)0( =−ε 0)0( =−σ (2. 17)
Aplicando o operador transformada de Laplace na equação (2.15) obtém-se:
⎭⎬⎫
⎩⎨⎧
+⎪⎭
⎪⎬⎫
⎪⎩
⎪⎨⎧
=⎭⎬⎫
⎩⎨⎧
••
ησσε )()()( tL
EtLtL
[ ] ησσσεε )()(.)0()(.)0( tL
EtLstLs +
+−=+−
−−
(2. 18)
Introduzindo as condições iniciais (eq. 2.17) na equação (2.18), resulta:
( ) ( )
. ( )L t L t
s L t sEσ σ
εη
⎡ ⎤⎣ ⎦= +
(2. 19)
Após a aplicação da transformada Laplaciana na equação (2.16), tem-se:
19
s
tuLtL 00 )()(
σσσ ==
(2. 20)
Rearranjando a equação (2.19) e com base no resultado de (2.20) temos:
200
..)(
ssEtL
ησσ
ε +=
(2. 21)
Reescrevendo (2.21) na forma do operador Laplaciano e aplicando a função
passo unitário chega-se facilmente a uma função que irá fornecer a deformação sofrida
pelo material no instante de tempo t:
⎭⎬⎫
⎩⎨⎧
+⎭⎬⎫
⎩⎨⎧= )(.)()( 00 tutLtu
ELtL
ησσ
ε
)(.)( 00 tutE
t ⎟⎟⎠
⎞⎜⎜⎝
⎛+=ησσ
ε
(2. 22)
A resposta ao teste de fluência para o elemento de Maxwell é apresentada na
figura (2.11), destacando-se a variação temporal da tensão e deformação ao longo do
teste de fluência.
Figura 2. 11 – Resposta ao teste de fluência para o elemento de Maxwell
Com base no gráfico da figura (2.11), observa-se que no instante de aplicação da
carga ( 0 0t t= = ) o elemento de Maxwell responde com deformação equivalente à
20
deformação elástica. Nos instantes de aplicação de carga seguintes ( 10 ttt << ),
observa-se um incremento na deformação advinda do comportamento viscoso,
conforme descrito na equação (2.22). No instante 1tt = , após a remoção da carga
aplicada, nota-se o fenômeno de reversibilidade de fluência. Dessa forma, o material
recupera o estado de deformação sofrido pela parcela elástica, atingindo um estado de
deformação menor que será equivalente ao grau de deformação perdida pela parcela
viscosa.
Através da equação (2.22), nota-se que se a mola for infinitamente rígida (E =∞)
o modelo se reduz ao fluido newtoniano. Da mesma forma se o amortecedor se tornar
infinitamente viscoso (η =∞), o modelo se reduz a uma mola (elástico).
A análise dos limites da equação (2.22) permite extrair as seguintes conclusões:
E0)0(
σε =+
(2. 23)
∞→∞)(ε (2. 24)
De acordo com CHRISTENSEN [6], o modelo de Maxwell prevê que há um
aumento ilimitado da deformação. Isto é uma característica de muitos fluidos, e por isso
os materiais que podem ser descritos pela equação (2.22) são conhecidos como fluidos
de Maxwell.
Analisa-se agora o elemento de Maxwell submetido ao teste de relaxação,
através da aplicação de uma deformação constante:
)().()( 0 tutt εε = (2. 25)
)(1)(.)(. tLtLEstLs σ
ησε +=
)(1. 0 tLEs
ss σ
ηε
⎟⎟⎠
⎞⎜⎜⎝
⎛+=
sE
EtL+
=
η
εσ 1..)( 0
⎪⎭
⎪⎬⎫
⎪⎩
⎪⎨⎧
=− tE
eLEtL.
0 ..)( ηεσ
21
tE
etEt ηεσ−
= ).(.)( 0 (2. 26)
A curva de relaxação para um elemento de Maxwell é apresentada na figura
(2.12). Através da análise dos limites da expressão (2.26), pode-se observar que
inicialmente a tensão no elemento é equivalente a tensão em regime elástico (eq. 2.27),
que decai exponencialmente até um valor que tenderá a zero (eq. 2.28).
0.)0( εσ E=+ (2. 27)
0)( =∞σ (2. 28)
Figura 2. 12 – Configuração típica de um corpo submetido a um teste de relaxação
O comportamento do elemento de Maxwell é estritamente peculiar ao de um
fluido viscoelástico, tratando-se dessa forma de um fluido não-newtoniano. Na
literatura, observa-se uma grande quantidade de materiais que podem ser enquadrados
como fluidos viscoelásticos, que são principalmente observados em alimentos e
componente biológicos.
NAGAYAMA et. al [19] utilizaram um modelo de Maxwell composto por
dois elementos ligados em paralelo para modelar o comportamento de células do tecido
muscular de artérias submetidas a pressões do fluxo sanguíneo.
No artigo intitulado “Rheology for the food industry”, MUNIZAGA e
CÁNOVAS [20] descrevem uma série de alimentos que podem ser representados pelo
modelo de Maxwell, a exemplo da maionese e ketchup.
2.7.2 Modelo de Kelvin
De acordo com o modelo reológico proposto por Lord Kelvin, o efeito da
viscoelasticidade pode ser representado pela combinação de um elemento viscoso e um
elástico, ligados em paralelo, conforme figura (2.13).
22
Figura 2. 13 – Representação do modelo de Kelvin
A aplicação de uma tensão σ nas extremidades do elemento de Kelvin provocará
uma distribuição de tensões nas parcelas elástica e viscosa, de forma que a tensão total
atuante no elemento será equivalente a soma de cada parcela (eq. 2.29). Por outro lado,
como os elementos estão interconectados em suas extremidades, as deformações
sofridas pela parcela elástica e viscosa deverão apresentar igual valor, conforme
equação (2.30), onde e representa a parcela elástica e i representa a parcela viscosa.
ie ttt )()()( σσσ += (2. 29)
ie ttt )()()( εεε == (2. 30)
Combinando as expressões (2.29) e (2.30) às relações constitutivas elásticas e
viscosas, eq. (2.13) e eq. (2.14), obtém-se a equação diferencial do modelo de Kelvin,
apresentada a seguir:
ie ttt )()()( σσσ += ησε
ηε )()()()( tttEt =+•
(2. 31)
Para analisar o comportamento do modelo de Kelvin quando submetido ao teste
de fluência, a equação diferencial (2.31) será resolvida através da aplicação do operador
Laplaciano, conforme procedimento apresentando anteriormente para o modelo de
Maxwell. Chega-se então a expressão (2.32), onde u(t) corresponde a função
Heavisidade. Para a tensão 0σ atuando no intervalo de (0, t1), o comportamento da
história de deformação )(tε é o indicado na figura 2.14.
)(.1)(.
0 tueE
ttE
⎥⎥⎦
⎤
⎢⎢⎣
⎡−=
−ησ
ε (2. 32)
23
Figura 2. 14 – Resposta ao teste de fluência para o elemento de Kelvin
Através da análise de limites da equação (2.32) pode-se provar que
imediatamente após a aplicação da carga, o elemento de Kelvin submetido ao teste de
fluência apresenta uma deformação instantânea nula (eq. 2.33). Por outro lado, se o
carregamento for mantido constante no tempo, o valor da deformação tende para um
valor assintótico apresentado na equação (2.34). Entretanto, após a remoção da carga, a
deformação decai exponencialmente até o valor de deformação nula.
0)0( =+ε (2. 33)
E0)(
σε =∞
(2. 34)
Submetendo o elemento de Kelvin ao teste de relaxação, através da análise da
equação (2.32), sob a aplicação súbita de uma deformação 0ε mantida constante, irá se
obter:
ση
εη
ε LLEL .1. =+•
ση
εη
Ls
Es .11. 0 =⎟⎟⎠
⎞⎜⎜⎝
⎛+
)(..)(..)(..)(.... 000000 tuEtLtuLEtLsEL εδεηεδεηεεησ +=+=+=
[ ])(.)(.)( 0 tuEtt += ηδεσ
(2. 35)
Com a análise dos limites da expressão (2.35) pode-se observar que, após a
aplicação da carga, a tensão se mantém constante, conforme a equação (2.37 ).
24
∞→+ )0(σ (2. 36)
0.)( εσ E=∞ (2. 37)
O comportamento do elemento de Kelvin é típico de um sólido viscoelástico. Na
literatura, podem-se destacar aplicações do modelo de Kelvin no campo da engenharia
geotécnica, geologia e indústria de alimentos. MAKRIS [21] apresenta um modelo
viscoelástico de Kelvin, combinado a partir de 2 elementos, para análise de fundações
frente a resposta sob o efeito de terremotos. Como aplicação na indústria de alimentos,
KUO e WANG [22] utilizaram um esquema composto de 6 elementos de Kelvin para
modelar o comportamento de diversas variedades de queijo, objetivando o
desenvolvimento de uma metodologia de rotulação de qualidade através de uma análise
de fluência viscoelástica.
2.7.3 Modelo de Burgers
O modelo proposto por Burger consiste na junção de um elemento de Kelvin
associado a um elemento viscoso e outro elástico em série, conforme figura (2.15).
Figura 2. 15 – Modelo reológico de Burgers
Delimitando 3 zonas de deformação, em que a zona (1) corresponde à parcela
viscosa, a zona (2) ao elemento de Kelvin e a zona (3) à parcela elástica, pode-se
afirmar que a deformação total no elemento é dada por:
)()()()( 321 tttt εεεε ++= (2. 38)
Redefinindo a expressão (2.38) em termos de taxa de deformação, e combinando
com as expressões (2.13), (2.14) e (2.31), chega-se a formulação diferencial para o
modelo de Burges apresentada na equação (2.39).
25
1 2 1 2 1 2 11
2 2 1 2 1 2
E E E E E E EEσ σ σ ε εη η η η η η
•• • •• •⎛ ⎞+ + + + = +⎜ ⎟⎝ ⎠
(2. 39)
O modelo de Burgers tem sido um dos modelos reológicos mais utilizados na
descrição do comportamento viscoelástico de misturas betuminosas, a exemplo dos
estudos realizados por VALE e MELO [23].
2.7.4 Modelo do sólido linear padrão
Aparentemente parece não haver uma padronização na definição do chamado
modelo do sólido viscoelástico padrão. Entretanto, os dois modelos apresentados na
literatura estão relacionados à associação de dois elementos elásticos e um elemento
viscoso.
O modelo mais difundido associa um elemento de Maxwell a um elemento
elástico conectados em paralelo (fig. 2.16). Esse modelo é discutido por
CHRISTENSEN [6], SIMO e HUGHES [24], DROZDOV [3], FINDLEY et. al. [25] e
VINOGRADOV e MALKIN [26]. A outra representação não será abordada nesse
trabalho.
Figura 2. 16 – Elemento sólido viscoelástico linear
Sendo eσ e mσ as tensões na parcela elástica e no elemento de Maxwell
respectivamente, representados pelos índices 1 e 2, a distribuição de tensões no
elemento é dada por:
)()()()( 221 tttt iee σσσσ +== (2. 40)
Rearranjado a expressão (2.40) e aplicando as equações (2.13) e (2.14) chega-se
à equação diferencial para o modelo do sólido padrão apresentado na equação (2.41).
26
)()()1()()(
2
1
2
1
22
tEtEEt
Et ε
ηε
ησσ
++=+•
•
(2. 41)
O sólido viscoelástico linear ao ser submetido ao ensaio de relaxação, fazendo-
se 0)( εε =t e considerando a equação (2.41), passará a ser representado pela equação
diferencial (2.42).
02
1
22
)()( εηη
σσ EtE
t+=+
•
(2. 42)
A solução geral da equação (2.42) é dada pela expressão (2.43), onde )0( +σ
representa a tensão no corpo viscoelástico imediatamente após a aplicação da carga.
[ ])0()exp()( 012
201
++−⎥⎦
⎤⎢⎣
⎡−+= σεη
εσ EtE
Et (2. 43)
Integrando a equação (2.41) com relação ao tempo, objetivando a determinação
de )0( +σ a partir do instante imediatamente anterior à aplicação da deformação no
tempo t(0-), obtém-se:
∫∫∫∫−−−−
++=+•
•nnnn tttt
dttEdttEEdttdt
Et
02
1
02
1
0 20 2
)()()1()()( εη
εησσ
(2. 44)
Admitindo que no instante inicial da aplicação da deformação (teste de
relaxação) instante t=t(0-), o corpo apresenta estado de tensão e deformação nula, o
resultado da equação (2.44) no instante t=tn se apresenta na forma:
)(1)(
2
1
2n
n tEE
Et
εσ
⎟⎟⎠
⎞⎜⎜⎝
⎛+=
(2. 45)
Fazendo tn=0+ na equação (2.45), observa-se que imediatamente após a
aplicação da deformação, o modelo do sólido viscoelástico linear apresenta uma
27
resposta em tensão total que corresponde a um comportamento elástico com rigidez
equivalente a soma dos elementos elásticos associados em paralelo, ou seja:
( ) 00)0( εσ E=+ ,
onde 210 EEE +=
(2. 46)
Substituindo o valor da tensão inicial, advindo da equação (2.46), na expressão
(2.43), pode-se estabelecer o valor da tensão para cada instante de tempo )(tσ na
forma:
⎥⎦
⎤⎢⎣
⎡−⎟⎟
⎠
⎞⎜⎜⎝
⎛−+= )exp(1)(
2
2
0
1
0
10 tE
EE
EEt
ηεσ
(2. 47)
Pela definição da eq. (2.04), a razão entre a tensão total e a deformação
constante é definida como função de relaxação )(tG , de forma que pode-se reescrever a
equação (2.47) na forma:
⎥⎦
⎤⎢⎣
⎡−⎟⎟
⎠
⎞⎜⎜⎝
⎛−+== )exp(1)()(
2
2
0
1
0
1
0
tE
EE
EEttG
ηεσ
(2. 48)
O modelo viscoelástico do sólido linear padrão tem sido amplamente utilizado
em modelagens de aplicações na engenharia estrutural, mecânica e geotécnica.
Tipicamente, aplicações podem ser feitas em materiais que apresentem fluência
viscoelástica, a exemplo do concreto, algumas variedades de madeira, materiais
poliméricos, dentre outros.
2.7.5 Generalização dos modelos básicos
Os modelos mais simples podem ser generalizados por um número finito de
combinações dos elementos básicos. Esse processo permite particularizar alguns
fenômenos viscoelásticos que não conseguem ser representados com precisão pelos
modelos tradicionais.
Para o modelo de Maxwell, as combinações podem ser feitas através da
interligação de elementos básicos em série ou em paralelo. Na combinação em série
28
apresentada na figura (2.17), o comportamento é semelhante ao modelo básico, sendo
que as parcelas elástica e viscosa resultantes correspondem ao somatório da
contribuição de cada elemento básico, conforme equação (2.49).
Figura 2. 17 – Elemento de Maxwell generalizado em série
∑∑==
••
+=n
i i
i
i i
tE
tt11
1)(1)()(η
σσε (2. 49)
Já na combinação dos elementos de Maxwell em paralelo, chamado de modelo
de Maxwell-Weichert, apresentado na figura (2.18), a influência de cada elemento de
Maxwell produz alteração no tempo de relaxação. Nesse modelo, a tensão e deformação
total podem ser obtidas através da análise individual de cada elemento.
Figura 2. 18 – Elemento de Maxwell generalizado em paralelo (Maxwell-Weichert)
De acordo com FINDLEY et. al [25], cada equação diferencial do elemento de
Maxwell pode ser reescrita sob a forma de um operador diferencial em relação ao tempo
D = d/dt. No modelo de Maxwell-Weichert, a tensão total corresponde à soma das
tensões em cada elemento de Maxwel que compõe o novo sistema, conforme a seguinte
equação:
1 1 1
n n
ii i
i i
DDE
σ σ ε
η= =
⎛ ⎞⎜ ⎟⎜ ⎟= =⎜ ⎟+⎜ ⎟⎝ ⎠
∑ ∑
(2. 50)
Eliminando o operador D do denominador da equação (2.50) através de uma
manipulação algébrica, chega-se facilmente à equação (2.51), que representa uma
expressão para o modelo generalizado de Maxwell-Weichert.
29
εηηηη
ση ⎥
⎥⎦
⎤
⎢⎢⎣
⎡+⎟⎟
⎠
⎞⎜⎜⎝
⎛+⎟⎟
⎠
⎞⎜⎜⎝
⎛++⎟⎟
⎠
⎞⎜⎜⎝
⎛+⎟⎟
⎠
⎞⎜⎜⎝
⎛+=
⎥⎥⎦
⎤
⎢⎢⎣
⎡⎟⎟⎠
⎞⎜⎜⎝
⎛+∏
=
......1.1...1.11
331133221 ED
ED
ED
ED
EDn
i ii
(2.51)
Já na associação de um número finito de elementos de Kelvin em paralelo,
dando origem ao modelo generalizado de Kelvin-Voigt (fig. 2.19), o sistema irá se
comportar como um elemento com rigidez e viscosidade equivalentes ao somatório de
cada propriedade associada ao n-ésimo elemento de Kelvin, conforme equação (2.52).
Figura 2. 19 – Elemento de generalizado Kelvin Voigt em paralelo
∑∑=
•
=
+=n
ii
n
iiE
11ηεεσ
(2. 52)
Na combinação de elementos de Kelvin em série, observa-se que a deformação
total corresponde à soma das deformações nos elementos básicos (fig. 2.20), conforme
equação (2.53).
Figura 2. 20 – Modelo generalizado de Kelvin-Voigt em série.
⎟⎟⎠
⎞⎜⎜⎝
⎛+
== ∑∑==
n
i ii
n
ii ED11
1η
σεε (2. 53)
Para a eliminação do operador diferencial D da equação (2.53), faz-se
novamente uma operação algébrica, cuja expressão final é dada pela equação (2.54).
30
( ) ( )( ) ( )( )[ ]......... 331133221
++++++=⎟⎟⎠
⎞⎜⎜⎝
⎛+
=
EDEDEDEDEDn
iii ηηηησηε
(2.54)
2.8 Representação por integrais hereditárias
Até a presente seção, os modelos viscoelásticos descritos foram apresentados na
forma diferencial. Uma outra maneira bastante difundida na literatura de se representar
o comportamento viscoelástico é através da chamada forma integral.
Esse procedimento é possível, no caso da viscoelasticidade linear, graças ao
emprego do princípio de superposição de Boltzmann (BERGER [27]). Para descrever a
formulação integral da viscoelasticidade linear adotam-se novamente das hipóteses as
seguintes considerações (DROZDOV [3]):
Idealiza-se um corpo viscoelástico inicialmente em seu estado natural (livre de
tensões). No instante 0τ ≥ uma força de valor unitário é aplicada, e a deformação
longitudinal e(t,τ) em cada instante t τ≥ é apresentado pela forma:
0 1( , ) ( ) ( , )e t e e tτ τ τ= + (2. 55)
Na formulação apresentada em (2.55), 0e representa a deformação elástica
instantânea e0, ao passo que 1e traduz a deformação adicional advinda dos efeitos
viscoelásticos. De acordo com o princípio de superposição de Boltzmann, a deformação
( )tε resultante da história de tensão ( ) | (0 )t tσ τ≤ ≤ é igual à soma das deformações
causadas pelas variações elementares de tensões:
0
( )( ) ( , )t dt e t d
dσ τε τ ττ
= ∫ (2. 56)
Para tal, o comportamento da tensão no tempo é idealizado a partir da curva
apresentada na figura (2.21), onde se podem observar os elementos infinitesimais e o
tempo 0τ ≥ de aplicação da força.
31
Figura 2. 21 – Variação da deformação no tempo
Integrando por partes a equação (2.56), obtém-se:
τσττσττσ
ττττσε τ
τ dttYtEtdtetet
ttt )(),(
)()()(),(|),()()(
000 ∫∫ ∂
∂+=
∂∂
−= ==
⎥⎦
⎤⎢⎣
⎡+−=Υ ),(
)(1),( ττ
τ tCE
t
(2. 57)
Como a tensão aplicada apresenta valor unitário, observa-se que na expressão
(2.57) a função ),( τtΥ corresponde à deformação total advinda da parcela 10
−= Ee e
1),( etC =τ . Naturalmente, se E for constante, a derivada da parcela ),( τtΥ e,
conseqüentemente, o respectivo integrando da expressão (2.57) é nulo.
A partir do operador K(t,τ) aplicado a ),( τtΥ , definido por DROSDOV [3]
como “creep kernel”, a equação (2.57) pode ser representada sob a forma da equação
(2.58), onde ),( τtC contabiliza uma medida dos incrementos de deformação
viscoelásticos:
⎥⎦
⎤⎢⎣
⎡+
∂∂
−= ),()(
1)(),( τττ
τ tCE
tEtK
⎥⎦⎤
⎢⎣⎡ += ∫
tdtKt
tEt
0)(),()(
)(1)( ττστσε
(2. 58)
32
Na prática, o operador “creep kernel” é um controlador do parâmetro de
fluência, e a parcela ),( τtC irá representar a influência do envelhecimento no material
viscoelástico. Idealiza-se que o módulo de elasticidade se mantém constante, de forma
que EtE =)( . Pode-se também caracterizar o operador K “creep kernel” unicamente
em função da diferença de tempo ( τ−t ), de forma que a equação (2.59) pode ser
reescrita sob a forma de uma integral de convolução:
⎥⎦⎤
⎢⎣⎡ −+= ∫
tdtKt
tEt
0)()()(
)(1)( ττστσε
⎥⎦⎤
⎢⎣⎡ −+= ∫
•tdtCt
tEt
0)()()(
)(1)( ττστσε
onde: •
= 0K(t) C e )(0 tECC =
(2. 59)
A equação (2.59) é conhecida como equação de fluência. De forma similar,
podemos estabelecer uma equação integral que correlacione a tensão viscoelástica em
função da deformação, também conhecida como equação de relaxação.
ττεττχεσ dtttEt
t)(),()()()(
0∫ ∂∂
−=
Onde ),()(),( τττχ tQEt +=
(2. 60)
A função ),( τtQ quantifica uma medida do parâmetro de relaxação do material
viscoelástico, que no instante τ = 0 corresponde a 0),( =τtQ . Com base no funcional
),( τtR , definido por DROSDOV [3], a equação (2.60) pode ser reescrita sob a forma:
⎥⎦⎤
⎢⎣⎡ −= ∫ ττετεσ dtRttEt
t)(),()()()(
0
onde:
1 1( , ) ( , )( ) ( )
R t Q tE t t E
τ ττ
⎡ ⎤∂= +⎢ ⎥∂ ⎣ ⎦
(2. 61)
De maneira similar às simplificações adotadas na equação de fluência, a equação
de relaxação também pode ser reescrita unicamente em função do tempo ( τ−t ), de
forma que a equação (2.61) pode ser reapresentada como:
33
⎥⎦⎤
⎢⎣⎡ −−=⎥⎦
⎤⎢⎣⎡ −−= ∫∫
•
ττετεττετεσ dtQttEdtRttEttt
)()()()()()()()()(0 00
onde:
)()( 0 tQtR•
−= e EtQtQ )()(0 =
(2. 62)
Em geral, as equações integrais são expressas em termos das funções de fluência
e de relaxação, conforme observado em CHRISTENSEN [6], BLAND [5], BERGEN
[27] e VINOGRADOV e MALKIN [26]. Dessa forma, DOZDROV [3] correlaciona a
função Q(t) com a função de relaxação )(tG , de forma que:
1)()( 0 += tQtG (2. 63)
Com o uso da expressão (2.63) e integração por partes, a equação (2.62) pode ser
apresentada na forma (2.64), obtendo desta maneira uma equação integral para a tensão
em função da taxa de deformação e da função de relaxação )(tG .
=⎥⎦
⎤⎢⎣
⎡−+−−= ∫
•== ττεττετεσ τ
τ dtGtGttEttt )()(|)()()()()(0 00
⎥⎦⎤
⎢⎣⎡ −=
•
∫ ττετσ dtGEtt
)()()(0
(2. 64)
No capítulo 4, serão discutidas as técnicas empregadas para a solução das
equações integrais apresentadas, baseadas principalmente no método dos elementos
finitos, que será focado no próximo capítulo.
34
3. O método dos elementos finitos
Em paralelo aos avanços tecnológicos, a engenharia tem aprimorado seus
métodos e processos em busca da obtenção de melhores resultados. Os métodos
numéricos têm acompanhado essa tendência e a cada dia busca-se por formulações
numéricas mais precisas e computacionalmente eficientes.
O método dos elementos finitos já é bastante difundido na solução de problemas
que envolvem equações diferenciais. Teve suas origens na década de 50, sobretudo na
aplicação à indústria aeroespacial e à engenharia estrutural, sendo posteriormente
estendido na década de 70 à mecânica dos fluidos e outras áreas a partir dos princípios
da análise funcional. Formulações distintas foram desenvolvidas nas décadas seguintes,
principalmente em função das oscilações numéricas observadas nos problemas da
mecânica dos fluidos. Uma das linhas de maior destaque desenvolvida nesse período
envolve a chamada formulação mista em elementos finitos. Maiores detalhes sobre o
desenvolvimento das formulações em elementos finitos ao longo da história podem ser
observados em LOGAN [28] e HUEBNER e THORNTON [29].
Serão apresentados na seção (3.1) os conceitos elementares da formulação
clássica do método dos elementos finitos. Em seguida, discute-se a formulação mista do
MEF, destacando-se os critérios de estabilidade numérica. As formulações estabilizadas
são discutidas na seção (3.4), enfatizando a formulação Galerkin Least Square,
doravante chamada de (GLS), e a de estabilização direta da pressão, que será
referenciada como (Bochev), formulações estas empregadas no desenvolvimento deste
trabalho.
3.1 Conceitos básicos
O foco principal do método dos elementos finitos consiste na transformação da
formulação diferencial do problema em uma forma integral variacional, também
conhecida como formulação variacional ou formulação fraca. A partir da formulação
variacional e da discretização espacial e temporal (problemas transientes) da região em
análise (domínio), a solução que em geral pertence a um espaço de dimensão infinita
(solução contínua) é aproximada por elementos de espaço de dimensão finita (solução
discreta). Esta solução aproximada é chamada de aproximação de elementos finitos. O
problema então recai na resolução de um sistema de equações algébricas cuja solução
corresponde aos valores aproximados nos nós dos elementos.
35
Para tal procedimento, é necessário ter o conhecimento da equação diferencial
regente do problema e das restrições às quais a mesmo está sujeita, caracterizando dessa
forma um problema de contorno e de valor inicial.
De forma geral, o objetivo é obter a solução u de um problema que atenda ao
conjunto de equações diferenciais no domínio Ω definidas pelo operador diferencial
matricial A(u) (3.1), respeitando as condições no contorno Γ regidas pelo operador
diferencial matricial B(u) (3.2).
1
2
( )( )
( )n
AA
A
⎧ ⎫⎪ ⎪⎪ ⎪= =⎨ ⎬⎪ ⎪⎪ ⎪⎩ ⎭
uu
Α(u) b
u
em Ω
(3.1)
g
u
uu
B(u) =
⎪⎪⎭
⎪⎪⎬
⎫
⎪⎪⎩
⎪⎪⎨
⎧
=
)(
)()(
2
1
nB
BB
em Γ
(3. 2)
Através da discretização espacial 1,..., e
hnΠ = Ω Ω , a geometria da região será
decomposta em partições do domínio, denominadas de malha de elementos finitos,
composta por ne elementos e nn pontos nodais, satisfazendo: 'e eΩ ∩Ω ≠∅ se 'e e≠
(fig. 3.1). Dessa forma, eΩ pode ser mapeado no correspondente elemento padrão
regular através de mapeamentos isoparamétricos, sendo o domínio de modelagem dado
pela união de todos os elementos finitos, tal que: 1 1
e ecn n
e ece ec= =
Ω∪Γ = Ω ∪ Γ∪ ∪ (fig. 3.1),
onde ecΓ denota o contorno dos elementos em contato com eΓ . A forma geométrica
dos elementos da malha pode assumir diversas configurações, sendo as mais comuns
para problemas bidimensionais as geometrias triangulares e quadriláteras.
36
Figura 3. 1 – Caracterização do domínio e contorno
No domínio da região ilustrada na figura (3.1), a solução aproximada em cada
elemento da malha (um) é obtida através da interpolação dos valores nodais ami de cada
i-ésimo ponto nodal do elemento (constituído por r pontos nodais ) da malha m, através
das funções Ni, conforme equação (3.3). Convém ressaltar que a aproximação é feita
para todos os graus de liberdade nodal (componentes do vetor solução).
1
rh
m m mi ii
u u a N=
≈ =∑ com m = 1, ..., ne (3. 3)
De acordo com ZIENKIEWICZ e TAYLOR [30], os valores de ami da expressão
(3.3) são parâmetros a serem determinados. As funções Ni utilizadas para aproximar a
solução nodal no interior dos elementos são conhecidas como funções de interpolação
(funções de forma), que devem ser linearmente independentes. MANSUR [31] discorre
que a escolha das funções de interpolação deve ser feita de tal maneira que, quando o
número de partições no domínio tenda a um valor infinito, a solução aproximada deve
tender ao valor exato, sendo este critério conhecido por completividade.
Pode-se também representar a solução aproximada em qualquer ponto do
domínio u = u1, u2, ..., unglT (vetor contendo aproximações de distintos graus de
liberdade) através de uma relação global entre os parâmetros a serem determinados em
cada elemento da malha (arranjados em uma forma matricial global) e suas respectivas
funções de forma (interpolação) (eq. 3.4).
37
h≈ =u u aN
Onde: 11 1
1
r
m mr
a a
a a
⎡ ⎤⎢ ⎥= ⎢ ⎥⎢ ⎥⎣ ⎦
a…
1j
jm
⎡ ⎤⎢ ⎥= ⎢ ⎥⎢ ⎥⎣ ⎦
NN
N e 1 2[ ]j T
kN N N=N
(3. 4)
A partir da discretização do domínio em elementos finitos, as funções de forma
Ni serão definidas localmente em cada elemento, tendo seus valores atribuídos a partir
de seus nós. É comum na literatura rotular a função de interpolação através de sua
classe. Uma função de interpolação pertence à classe Cn-1 se ela é de ordem n, ou seja,
se sua derivada de ordem (n-1) existe (é finita). Neste caso, tanto a função como suas
derivadas até ordem (n-1) devem ser contínuas.
Conforme discutido anteriormente, o objetivo é transformar a forma diferencial
em uma forma integral variacional, que pode ser de maneira geral representada pela
equação (3.5). A partir dessa formulação será permitida a obtenção dos parâmetros
desconhecidos ami da expressão (3.3), calculados com uma malha de ne elementos:
1
e
e e
n
j j j j je
d d d d dΩ
=Ω Γ Ω Γ ΩΓ
⎛ ⎞Ω+ Γ = Ω+ Γ = Ω⎜ ⎟⎜ ⎟
⎝ ⎠∑∫ ∫ ∫ ∫ ∫G g G g F
(3. 5)
A forma integral da expressão (3.5) para cada problema específico pode ser
obtida através da aplicação do método dos resíduos ponderados. Uma outra maneira de
se obter a forma (3.5) é através da determinação do funcional variacional.
No método dos resíduos ponderados, os resíduos ( ) ( )h −A u A u e ( ) ( )h −B u B u
são multiplicados por funções vetoriais ( )1
T
j jw w=w com j = (1,...,r),
linearmente independentes (funções de ponderação), e em seguida são integrados em Ω
e Γ , respectivamente na seguinte forma:
( ) ( ) 0TT
jj d dΩ Γ
Ω+ Γ =∫ ∫w A aN w B aN Para j = 1 a ne (3. 6)
38
Na prática, é possível reduzir a ordem dos operadores A(u) e B(u) da expressão
(3.6) através da integração por partes, para casos unidimensionais, e do teorema de
Gauss para casos bidimensionais, obtendo-se assim a formulação variacional associada.
Maiores detalhes sobre o procedimento podem ser observados em MANSUR [31],
ZIENKIEWICZ e TAYLOR [30] e BATHE [32].
A forma de definição da função de ponderação wj irá determinar o método de
aproximação em elementos finitos. Os mais difundidos na literatura são:
• Método de Galerkin: No método de Galerkin, adota-se como função de
ponderação a mesma função utilizada como função de forma
(interpolação). Por produzir matrizes simétricas, esse método é bastante
empregado na literatura.
• Ponderação por mínimos quadrados (least square): No método dos
mínimos quadrados, adota-se como função de ponderação o próprio
resíduo A(u) derivado do operador diferencial.
De acordo com ZIENKIEWICZ e TAYLOR [30], na literatura, a denominação
de método de Petrov-Garlerkin é sempre associado às estratégias em que o uso da
função de ponderação é tal que j jw N≠ .
As aplicações dos modelos empregados nesse trabalho atendem à equação de
Stokes para a elasticidade linear compressível e incompressível, o escoamento de
fluidos incompressíveis em regime estacionário (fluxo de Stokes) e para a
viscoelasticidade linear, sendo regidas pelas equações (3.7) e (3.8), onde div
corresponde ao operador divergente. Além disso, b está assocido à ação de forças de
massa, p é a pressão e K o módulo de elasticidade volumétrica (Bulk Module).
div= − =A(u) σ b (3.7)
div pK
= − =A(u) u (3. 8)
A equação (3.7) é na verdade uma simplificação da equação geral de
conservação de momento de Navier-Stokes, obtida através da análise de equilíbrio de
um volume de controle elementar, desprezando-se a variação temporal de u. Além
disso, a expressão (3.8) representa a relação constitutiva volumétrica, relacionada à
39
conservação da massa. Maiores detalhes podem ser observados em FOX &
MCDONALD [33].
As componentes de uT=[u,v,w] na mecânica dos sólidos representam as
componentes dos deslocamentos, enquanto que na mecânica dos fluidos o vetor u
corresponde as componentes da velocidade [30]. Em ambos os casos, as componentes
do tensor de tensão σ serão de maneira geral dependentes das componentes de u através
das relações constitutivas, discutidas previamente no capítulo 2.
O tensor σ é simétrico, a menos que se considere efeitos de rotação na equação
de equilíbrio, e pode ser decomposto na forma expressa em (3.9a) (FOX &
MCDONALD [33]). Em função da simetria, pode-se também representar σ em uma
forma vetorial, expressa em (3.9b).
0 0
0 0
0 0
xx xy xz xx xy xz
yx yy yz yx yy yz
zx zy zz zx zy zz
p p
p p
p p
σ τ τ σ τ τ
τ σ τ τ σ τ
τ τ σ τ τ σ
⎡ ⎤ ⎡ ⎤+⎛ ⎞⎢ ⎥ ⎢ ⎥⎜ ⎟⎢ ⎥ ⎢ ⎥⎜ ⎟= = − + +⎢ ⎥ ⎢ ⎥⎜ ⎟⎢ ⎥ ⎢ ⎥⎜ ⎟ +⎢ ⎥ ⎢ ⎥⎝ ⎠⎣ ⎦ ⎣ ⎦
σ
T
xx yy zz xy xz yzσ σ σ τ τ τ⎡ ⎤= ⎣ ⎦σ
(3. 9a) (3. 9b)
Na equação (3.9a), as componentes iiσ se associam à tensão normal, τ à tensão
tangencial (tensão de cisalhamento), e p representa a pressão estática associada a parte
isotrópica do tensor de tensões.
A formulação aqui apresentada é complementada pelas condições de contorno
especificadas em (3.10) e (3.11). A equação (3.10) faz referência à condição de
contorno de Dirichlet, utilizada quando o valor da velocidade (ou deslocamento) é
conhecido em uma região do contorno uΓ . No caso em que se conhecem as
componentes das forças atuantes por unidade de área da superficie em tΓ , faz-se uso da
equação (3.11), conhecida como condição de Neumann.
gu = em uΓ (3.10)
tσ.n = em tΓ (3.11)
onde u tΓ = Γ ∪Γ
40
De posse da equação diferencial regente do problema e das condições de
contorno associadas, pode-se finalmente estabelecer uma forma fraca da equação de
Stokes para a aproximação em elementos finitos. A seguir, será apresentada uma síntese
do processo utilizado por ZIENKIEWICZ e TAYLOR [30], através do estabelecimento
de uma expressão funcional via princípio dos trabalhos virtuais. Para tanto, a equação
de equilíbrio (3.7) pode ser agora reescrita em função de suas componentes, de maneira
que:
⎪⎭
⎪⎬
⎫
⎪⎩
⎪⎨
⎧=
⎪⎭
⎪⎬
⎫
⎪⎩
⎪⎨
⎧+
⎪⎪⎪
⎭
⎪⎪⎪
⎬
⎫
⎪⎪⎪
⎩
⎪⎪⎪
⎨
⎧
∂
∂+
∂∂
+∂
∂∂
∂+
∂
∂+
∂
∂∂∂
+∂
∂+
∂∂
000
z
y
x
yzxzzz
yzyxyy
xzxyxx
bbb
yxz
zxy
zyx
ττσ
ττσ
ττσ
(3. 12)
Partindo da expressão geral (3.6), a variação [ ], , Tu v wδ δ δ δ=u corresponde à
função de ponderação (forma vetorial) equivalente ao método dos resíduos ponderados:
xyxx xzx
yy yx yzTy
V V
yzxzzzz
u bx y z
dV v b dVy x z
w bz x y
τσ τδ
σ τ τδ δ
ττσδ
⎡ ⎤∂⎛ ⎞∂ ∂− + + + +⎢ ⎥⎜ ⎟∂ ∂ ∂⎝ ⎠⎢ ⎥⎢ ⎥∂ ∂ ∂⎛ ⎞⎢ ⎥= − + + + +⎜ ⎟∂ ∂ ∂⎢ ⎥⎝ ⎠⎢ ⎥
∂⎛ ⎞∂∂⎢ ⎥− + + +⎜ ⎟⎢ ⎥∂ ∂ ∂⎝ ⎠⎣ ⎦
∫ ∫u A(u)
(3. 13)
Através da integração por partes e da realização de algumas operações algébricas
chega-se à forma variacional dada pela expressão (3.14), que pode ser reescrita
utilizando a notação de resíduos ponderados, bastando substituir Tuδ =u w .
0=Γ−Ω−Ω ∫∫∫Γ
ddd T
V
T
V
T tubuσε δδδ
0Tu u u
V V
d d dΓ
∇ Ω− Ω− Γ =∫ ∫ ∫w σ w b w t
(3. 14)
Na formulação acima, δ δ=ε S u , onde ∇ corresponde ao operador gradiente:
41
0 0 0
0 0 0
0 0 0
x y z
y x z
z y x
⎡ ⎤∂ ∂ ∂⎢ ⎥∂ ∂ ∂⎢ ⎥⎢ ⎥
∂ ∂ ∂⎢ ⎥∇ = ⎢ ⎥∂ ∂ ∂⎢ ⎥⎢ ⎥∂ ∂ ∂⎢ ⎥⎢ ⎥∂ ∂ ∂⎣ ⎦
e εδ é um operador linear definido por:
( )
( )
( )
( ) ( )
( ) ( )
( ) ( )
ux
vy
wz
u vy x
v wz y
w ux z
δ
δ
δ
δ δ δ
δ δ
δ δ
δ δ
∂⎧ ⎫⎪ ⎪∂⎪ ⎪⎪ ⎪∂⎪ ⎪
∂⎪ ⎪⎪ ⎪⎪ ⎪∂⎪ ⎪
∂⎪ ⎪= = = ∇⎨ ⎬
∂ ∂⎪ ⎪+⎪ ⎪∂ ∂⎪ ⎪⎪ ⎪∂ ∂⎪ ⎪+⎪ ⎪∂ ∂⎪ ⎪⎪ ⎪∂ ∂⎪ ⎪+∂ ∂⎩ ⎭
ε S u u
(3.15)
Após a discretização, a equação (3.14) passa a ser representada por um conjunto de
equações algébricas e é geralmente caracterizada na forma matricial Ku=f.
3.2 Formulação mista e incompressibilidade
O comportamento incompressível ou quase incompressível dos materiais pode
induzir a algumas dificuldades na simulação numérica da formulação clássica do MEF,
tais como o “travamento” volumétrico (locking), com conseqüente imprecisão de
solução e oscilações na distribuição de pressão. Esse efeito em geral ocorre na
aproximação de fluidos incompressíveis, na elasticidade linear incompressível, em
42
materiais hiper-elásticos e em matérias que exibem comportamento viscoso, a exemplo
da viscoelasticidade e viscoplasticidade.
De acordo com BATHE [32], o termo “travamento” (locking) é tradicionalmente
usado quando o deslocamento (ou velocidade) aproximado pela solução do MEF é
menor do que o valor real. Isso irá induzir a uma propagação do erro ainda maior na
aproximação das tensões/deformações, que dependem das derivadas dos
deslocamentos/velocidades.
Uma interpretação física do fenômeno de “travamento” é discutida por
PANNACHET e BOONPICHETVONG [34], que afirmam que o efeito de
“travamento” ocorre quando a malha em elementos finitos não é capaz de descrever, de
forma contínua, a preservação de volume. No caso da elasticidade, esse efeito se dá
quando o material está próximo à situação de incompressibilidade (elevados valores do
coeficiente de Poisson: 5,0→ν ), implicando que a deformação volumétrica tende a ser
nula. Para ilustrar o fenômeno, apresentam-se os elementos da figura (3.2). Nota-se que
o elemento superior da região destacada pode apenas deslocar verticalmente seu nó
livre. Por outro lado, o único nó livre do elemento inferior pode apenas ser deslocado
horizontalmente. Como o único nó livre é comum a ambos os elementos, a única forma
de atender a ambas as situações de forma simultânea é não realizar movimento,
travando dessa forma o deslocamento no elemento.
Figura 3. 2 – Travamento volumétrico da malha.
(Adaptado de PANNACHET e BOONPICHETVONG [34]] )
Convém ressaltar que o efeito de “travamento” ocorre somente em problemas
em que há restrição à deformação. Por exemplo, esse fenômeno não está presente em
problemas em que se admite um estado plano de tensão. Entretanto, valores elevados do
coeficiente de Poisson ( 5,0→ν ) podem induzir ao modo de pressão espúria,
43
comprometendo dessa forma a aproximação das tensões e da pressão. Esse fenômeno
será descrito com mais detalhes na seção 3.3.
Baseado no trabalho de Nagtegall (1974), BITEENCOUR [35] discute o efeito
do refinamento da malha no fenômeno de “travamento”. Segundo o mesmo, o efeito do
“travamento” ocorre quando a razão entre o número de graus de liberdade na malha e o
número de restrições é igual ou inferior ao valor unitário. Para contornar o problema,
poder-se-ia aumentar o número de graus de liberdade através do refinamento de malha,
de forma que quando o número de elementos tende a um valor infinito, a formulação
clássica de Galerkin se aproxima da solução exata, removendo dessa forma o efeito do
“travamento” de malha. Entretanto, a aproximação da pressão ainda continua
comprometida.
Uma interpretação em termos de conservação de energia sobre a origem do
fenômeno de “travamento” é discutida por SORIANO [36]. Segundo o mesmo, a
inabilidade de certos elementos finitos em representar quantidades críticas de
determinados tipos de energia, com conseqüente grande rigidez espúria, é o que
denomina de “travamento”. Há então inconsistência de alguns campos de variáveis
quando um ou mais tipos de energia tende a ser anulada, o que geralmente é requerido
pelas condições físicas limites. Por exemplo, no caso da elasticidade incompressível, a
energia potencial correspondente à variação de volume tende a ser anulada à medida
que o coeficiente de Poisson tende a 0,5.
Uma das alternativas pioneiras empregadas para remover o efeito de
“travamento” consiste no uso do modelo de deslocamentos/velocidades com integração
reduzida ou integração reduzida seletiva.
A integração reduzida seletiva substitui a parcela relacionada à deformação
volumétrica em cada ponto de integração por uma deformação volumétrica média no
elemento. Porém, esse método não é capaz de prever o “travamento cisalhante” em
problemas predominados por flexão [37].
A integração reduzida uniforme, por possuir apenas 1 ponto de integração, é
mais eficiente do que o método de integração reduzida seletiva. Contudo, a energia
artificial introduzida para controlar o “travamento cisalhante” pode afetar
negativamente a precisão da solução e requerer um refinamento de malha significativo
[37].
Com os avanços na fundamentação matemática do método dos elementos finitos,
sobretudo baseada na análise funcional, começaram a ser desenvolvidas formulações
44
alternativas à formulação clássica (deslocamentos/velocidades), principalmente para
remover o efeito do “travamento” característico da formulação da Galerkin em situações
próximas a incompressibilidade.
Uma das vertentes se baseia na chamada formulação mista, cuja idéia central é
utilizar como variável de estimação não apenas o campo de velocidade/deslocamento,
mas também outras variáveis de interesse de aproximação. De acordo com BATHE
[32], os trabalhos pioneiros que empregavam formulações mistas se basearam no
funcional de Hu-Washizu, sendo posteriormente substituído pelo uso do funcional de
Hellinger-Reissner.
Dentro do contexto lingüístico da formulação mista, o número de variáveis
utilizadas na aproximação sub-definirá a classe do problema em uma terminologia
denominada de campo. Na literatura, são bastante difundidas formulações de 2 campos
(u-σ e u-p), 3 campos (u-σ-ε e u-p-ε) e 4 campos (u-σ-ε-T). Esse trabalho se restringirá
ao estudo da formulação mista de 2 campos (u-p). Maiores detalhes sobre as demais
formulações são disponibilizados em BATHE [32] e ZIENKIEWICZ e TAYLOR [30].
De acordo com ZIENKIEWICZ e TAYLOR [30], o principal problema no uso
da formulação clássica baseada apenas no deslocamento para a classe de problemas
incompressíveis ou quase-incompressíveis reside na determinação da tensão média ou
pressão. Por essa razão, é conveniente efetuar a separação do tensor de tensões em suas
componentes desviatória e hidrostática, de maneira similar à expressão (3.9). A pressão
hidrostática corresponde à média das tensões nas direções ortogonais, de forma que:
1 1( ) ( )3 3
Txx yy zzp trσ σ σ= + + = =σ σ m
(3.16)
onde:
000111[=Tm ] (3.17)
Admitindo-se que o material seja isotrópico, a deformação volumétrica εv se
relaciona com a pressão hidrostática através do módulo de elasticidade volumétrica K
(3.18):
vpK
ε =
onde:
(3. 18)
45
Tv xx yy zzε ε ε ε= + + = m ε
T
xx yy zz xy xz yzε ε ε ε ε ε⎡ ⎤= ⎣ ⎦ε
)21(3 ν−=
EK (3.19)
De maneira similar à componente do tensor de tensão, pode-se definir uma
componente de deformação desviatória, expressa pela equação (3.20):
1 1 1( )3 3 3
T T Td v dε ⎛ ⎞= − ≡ − = − ⊗ =⎜ ⎟
⎝ ⎠ε ε m I mm ε I m m ε I ε
(3. 20)
Ainda de acordo com a notação de ZIENKIEWICZ e TAYLOR [30], a tensão
desviatória pode ser reescrita na forma: p= −s σ m (3. 21)
Com base nas prerrogativas discutidas acima, pode-se representar a tensão total
pela soma da parcela desviatória s e da pressão hidrostática p.
0123
Tuu upp G p p⎛ ⎞= + = − + ≡ +⎜ ⎟
⎝ ⎠σ s m I mm ε m D ε D
(3. 22)
onde:
0
2 02
2112
10 1
⎡ ⎤⎢ ⎥⎢ ⎥⎢ ⎥
= ⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎣ ⎦
I
2ij
ij =ε
ε se i j≠
(3. 23)
)1(2 ν+=
EG (3. 24)
Para se obter a forma fraca, aplica-se o método dos resíduos ponderados,
conforme discutido na seção 3.1. Os espaços correspondentes às funções de forma e
46
função de ponderação devem estar restritos às definições de (3.25) a (3.27). Maiores
detalhes podem ser observados em GUERMOND [38] e ADAMS [39]. 1 2 ( ) | S H u g= ∈ Ω =u em uΓ (3. 25)
* 2 ( )P q L= ∈ Ω se não existe incompressibilidade total
* 2 ( ); 0P q L qdΩ
= ∈ Ω Ω =∫ no caso de incompressibilidade total
(3. 26)
1 2 ( ) ; 0w H w= ∈ Ω =w em uΓ (3. 27)
Para uma melhor compreensão das definições supracitadas, será apresentado
resumidamente cada espaço associado às funções de interpolação e de ponderação para
a formulação mista. Detalhes adicionais podem ser consultados nas referências [36],
[37] e [38]. Convém ressaltar que as derivadas são no sentido fraco:
• 2 2( ) : ; L dϕ ϕΩ
Ω = Ω Ω < ∞∫ ;
• 1 2 2( ) : ; ( ), ( )i
H L Lxϕϕ ϕ ∂
Ω = Ω ∈ Ω ∈ Ω∂
;
• 1 11,...,( ) ( ); ( )n
n iH Hϕ ϕ ϕ ϕΩ = = ∈ Ω com 1n ≥ ;
O objetivo então é encontrar o par de soluções (u, p) ∈ *S P× que satisfaça à
sentença de resíduos ponderados para a equação de conservação de momentum e
conservação de massa. Para a conservação de momentum, o tensor de tensões definido
na expressão (3.14) pode ser agora aplicado na forma mista, através da substituição do
tensor de tensões (3.24) em suas parcelas desviatória e hidrostática:
: ( ) d du uu up u up dΩ Ω Γ
∇ + Ω = Ω+ Γ∫ ∫ ∫w D ε D w b w t (3. 28)
Na equação (3.28), u∇w denota o gradiente da função de ponderação uw . O
deslocamento (ou velocidade) e a pressão são finalmente calculados fazendo-se a
aproximação (3.31) (forma matricial), escolhendo-se como funções de ponderação para
o deslocamento e pressão as mesmas adotadas como funções de interpolação (3.29 e
3.30), caracterizando assim o uso do método de Galerkin. Na notação utilizada, uδd e
pδP correspondem respectivamente aos valores infinitesimais (discretos) do campo de
deslocamento e pressão, conforme notação apresentada em [44]:
47
( ) ( )T Tu u u=w δd N (3. 29)
( ) ( )T T
p pq = δP N (3. 30)
hu≈ =u u N u
hpp p≈ = N P
(3.31)
A sentença variacional resultante da expressão (3.28) assume a forma da
equação (3.32) (forma fraca):
( )
( )
( ) ( ) : ( )d ( ) d
( ) d
TT T Tuu u up p u u
TTu u
Ω Ω
Γ
∇ ∇ + Ω = ⋅ Ω+
+ ⋅ Γ
∫ ∫
∫
u uδd N D N u D N P δd N b
δd N t
(3.32)
Além disso, para a sentença volumétrica (conservação da massa), a forma fraca é
obtida de acordo com a expressão (3.33), onde divu corresponde ao divergente aplicado
ao campo u.
( )div ( ) divT T
p pp pq d dK KΩ Ω
⎛ ⎞ ⎛ ⎞− Ω = − Ω⎜ ⎟ ⎜ ⎟⎝ ⎠ ⎝ ⎠∫ ∫u δP N u *q P∈
(3.33)
Matricialmente, a sentença resultante pode ser representada pela expressão
(3.34), conforme formulação apresentada por ZIENKIEWICZ e TAYLOR [30] e
BATHE [32]:
uu up u
up pp p
⎡ ⎤ ⎧ ⎫⎧ ⎫⎪ ⎪ ⎪ ⎪⎢ ⎥ =⎨ ⎬ ⎨ ⎬⎢ ⎥− ⎪ ⎪ ⎪ ⎪⎩ ⎭⎣ ⎦ ⎩ ⎭
T
K K Ru
K K RP
sendo:
1
ene
uu uue=
=K K∪
1
ene
up upe=
=K K∪
1
ene
pp ppe=
=K K∪
(3.34)
48
1
ene
u ue=
=R R∪
1
ene
p pe=
=R R∪
onde:
de Tuu uu
Ω
= Ω∫K B D B
de Tup up p
Ω
= Ω∫K B D N
1 de Tpp p pKΩ
= Ω∫K N N
d de T Tu u u
Ω Γ
= Ω + Γ∫ ∫R N b N t
0ep =R
(3. 35)
s u= ∇ =ε N u Bu (3. 36)
1[ ]es u nB B= ∇ =B N … (3. 37)
0123
Tuu G ⎛ ⎞= −⎜ ⎟
⎝ ⎠D I mm
(3. 38)
up =TD m (3. 39)
De acordo com BATHE [32], uma atenção especial deve ser dada às condições
de contorno no caso da incompressibilidade total. Segundo o mesmo, quando todo o
contorno estiver restrito à condição de Dirichlet (apenas deslocamento/ velocidade
prescrito), o deslocamento/velocidade u deve ser compatível com a deformação
volumétrica do domínio. Convém ressaltar que esse caso particular de condição de
contorno é mais comum em fluidos. Essa observação física é expressa por:
0t
vd dtεΩ Γ
Ω = ⋅ =∫ ∫ u n (3. 40)
Na expressão (3.40), n é o vetor unitário normal à superfície do domínio. A
relação é atendida caso o deslocamento/velocidade prescrito na direção normal a
superfície do domínio seja tal que o volume do corpo seja conservado. Uma segunda
condição é que a pressão deve ser precrita em algum ponto do domínio, visto que a
relação constitutiva é função tanto do deslocamento/velocidade quanto da pressão, e a
49
não prescrição de um valor pontual poderá torná-la indeterminada. Em todo caso, pode-
se prescrever qualquer valor para a pressão, uma vez que a pressão não irá conduzir a
nenhuma deformação. Detalhes adicionais podem ser consultados na referência [36].
Observa-se na formulação apresentada em (3.35) que a matriz Kpp assume
valores próximos a zero para elevados valores de Poisson ( 5,0→ν ), uma vez que o
módulo K tende ao infinito. Esse fato induzirá a alguns critérios para garantia de
estabilidade.
O estudo da estabilidade nas formulações mistas (u-p) são geralmente associadas
à condição de elipsidade e inf sup, relacionada às condições de (BNB) Babuska-Necas-
Brezzi (GUERMOND [38]), discutidas na próxima seção.
3.3 Condições de estabilidade na formulação (u-p) para a equação de Stokes
Conforme discutido por DE [40], a formulação em elementos finitos baseada
unicamente na sentença variacional do deslocamento/velocidade (formulação clássica) é
sempre estável, mas é sujeita ao fenômeno “travamento” em situações de
incompressibilidade. Por outro lado, a formulação mista consegue remover o efeito do
“travamento”, mas nem sempre é estável para todas as classes de funções de
interpolação da pressão e deslocamento/velocidade.
De acordo com DE [40], a estabilidade da formulação mista (u-p) é garantida
caso os critérios de elipisidade e condições inf-sup sejam atendidos. A condição de
elipsidade é automaticamente atendida para problemas incompressíveis, desde que se
faça uso da integração total, não usando dessa forma esquemas de integração reduzida.
Para a discussão do critério inf sup, consideram-se espaços das funções de
ponderação wu como sendo 1
0HVh ⊂ . Toma-se também 2 ( )hQ L⊂ Ω como o espaço das
funções de ponderação wp.
De acordo com a condição inf sup, 0>∃β tal que (GUERMOND [38], BATHE
[32]):
*
0 0 0 1
( )(div )inf sup
u hq wu
q P w V
q d
qβ
≠ ≠
Ω
∈ ∈
Ω≥∫ u
u
w
w
(3. 41)
onde:
50
( ) Ω⋅=⋅ ∫Ω
d22
0 (3. 42)
22 2
2 2
1 01 1
( )i
i j jx= =Ω
⎛ ⎞∂ ⋅⋅ = ⋅ + ⎜ ⎟⎜ ⎟∂⎝ ⎠
∑∑∫ (3. 43)
De acordo com GUERMOND [38], quando a condição (3.41) é atendida, a
seguinte estimativa pode ser obtida:
1 0( , ) l
h hu u p p C p h− + − ≤ u com 1l ≥ (3. 44)
Em (3.44) h é o maior diâmetro avaliado para os elementos da malha. Maiores
detalhes sobre as condições inf-sup podem ser obtidos nas referências [38] e [42].
De forma bastante simplificada, quando o resultado da expressão (3.41) conduz
a β=0 é dito que a formulação encontra-se no “modo de pressão espúria”, visto que a
condição inf-sup não é satisfeita. Na prática, isso ocorre em situações próximas a
incompressibilidade total, quando o módulo K tende a um valor infinito e,
conseqüentemente, a matriz Kpp tende a ser anulada. Dessa forma, a matriz resultante
(3.34) irá se tornar singular, e o sistema poderá dessa forma retornar múltiplas soluções,
não garantindo a acurácia da solução. Os fundamentos matemáticos dessa interpretação
podem ser observados em [38] e [40].
De acordo com ZIENKIEWICZ e TAYLOR [30], o resultado mais importante
das restrições da inf-sup refere-se ao fato de que a ordem da função de interpolação
utilizada para aproximação do deslocamento/velocidade deve ser maior do que a
utilizada na pressão (3.45). Em função disso, a formulação mista (u-p) foi largamente
empregada através do uso de elementos finitos isoparamétricos com elevados números
de nós.
pu ηη > (3. 45)
CHAPELLE e BATHE [41] propuseram uma estratégia numérica para
verificação das condições inf-sup. Nessa estratégia, o parâmetro β é calculado para uma
seqüência de malhas discretizadas do problema. É dito que a formulação mista (u-p)
passa no teste inf-sup numérico se o parâmetro β se aproxima assintoticamente para um
51
valor positivo à medida que a malha vai sendo refinada. Maiores detalhes são
apresentados nas referências [40] e [41].
Em geral, observam-se na literatura duas estratégias distintas para o tratamento
da pressão no uso da formulação mista (u-p). Na primeira, a pressão não é considerada
contínua entre os pontos nodais adjacentes dos elementos. De acordo com BATHE [32],
nessa estratégia a pressão pode ser estaticamente condensada, fazendo com que a
solução seja aproximada apenas em termos do campo de deslocamento (ou velocidade).
Segundo o mesmo, essa estratégia pode ser eficiente para situações próximas à
incompressibilidade somente quando a malha for suficientemente refinada. Um
arquétipo dessa estratégia é o chamado método B-bar. Segundo ZIENKIEWICZ e
TAYLOR [30], esse método é frequentemente utilizado em problemas não-lineares, tais
como plasticidade e problemas com não linearidade geométrica, principalmente pela
simplicidade da formulação e respectiva facilidade no acoplamento da formulação mista
a problemas mais complexos. Entretanto, a condensação estática inviabiliza o uso da
formulação para incompressibilidade total. Convêm ressaltar que o uso dessa
formulação está restrito às condições de estabilidade discutidas previamente.
Uma outra maneira de tratar a pressão nos problemas mistos consiste no
emprego de formulações que garantam a continuidade da pressão nos nós adjacentes aos
elementos finitos. Nesse caso, a pressão não pode ser condensada, sendo então obtida
diretamente através da aproximação nodal.
Com base nos critérios de estabilidade, BATHE [32] apresenta em síntese as
seguintes considerações para classes de elementos finitos submetidas ao uso da
formulação mista (u-p):
• Elemento Q1-P0 (4/1 em 2D e 8/1 em 3D): O elemento aproxima razoavelmente
o deslocamento, mas a aproximação da tensão pode não ser acurada, visto que o
tratamento da pressão constante no elemento pode induzir a possíveis flutuações
da pressão. O elemento não satisfaz à condição inf-sup.
Figura 3. 3 – Elemento Q1-P0 – adaptado de BATHE [32]
52
• Elemento Q2r-P0 (8/1 em 2D e 20/1 em 3D): O elemento satisfaz à condição inf-
sup, mas o fato de se assumir a pressão constante pode requerer uma malha
bastante refinada para a aproximação da tensão com acurácia.
Figura 3. 4- Elemento Q2
r-P0 (8/1 em 2D e 20/1 em 3D): adaptado de BATHE [32]
• Elemento Q2-P1 (9/3 em 2D e 27/4 em 3D): Corresponde ao primeiro membro
da família de elementos da forma Qn-Pn-1, com n≥2. Essa família de elementos
satisfaz à condição inf-sup.
Figura 3. 5 - Elemento Q2-P1 (9/3 em 2D e 27/4 em 3D): adaptado de BATHE [32]
• Elemento P2-P0 (6/1 em 2D e 10/1 em 3D): O elemento satisfaz à condição inf-
sup. Entretanto, o fato de se assumir a pressão constante pode requerer uma
malha bastante refinada para a aproximação da tensão com acurácia.
Figura 3. 6 Elemento P2-P0 (6/1 em 2D e 10/1 em 3D): adaptado de BATHE [32]
• Elemento P2+-P1 (7/3 em 2D e 11/4 em 3D): Corresponde ao primeiro membro
da família de elementos Pn+-Pn-1, com n≥2. Essa família de elementos satisfaz à
condição inf-sup.
53
Figura 3. 7 - Elemento P2
+-P1 (7/3 em 2D e 11/4 em 3D): adaptado de BATHE [32]
• Elemento Q2-Q1 (9/4-c em 2D e 27/8-C em 3D): Corresponde à família de
elementos Qn-Qn-1 , com n≥2. Essa família de elementos satisfaz à condição inf-
sup. Os nós destacados são pontos de aproximação do deslocamento (ou
velocidade) e da pressão, que agora é contínua entre os elementos.
Figura 3. 8 – Elemento Q2-Q1 (9/4-c em 2D e 27/8-C em 3D): adaptado de BATHE [32]
• Elemento P2-P1 (6/3-c em 2D e 10/4-c em 3D): Corresponde à família de
elementos Pn-Pn-1 , com n≥2. Essa família de elementos satisfaz à condição inf-
sup. Os nós destacados são pontos de aproximação do deslocamento (ou
velocidade) e da pressão, que agora é contínua entre os elementos.
Figura 3. 9 - Elemento P2-P1 (6/3-c em 2D e 10/4-c em 3D): adaptado de BATHE [32]
3.4 Métodos estabilizados
As aplicações envolvendo a formulação mista (u-p) foram bastante difundidas
para a elasticidade incompressível e para o problema de Stokes através do uso de
elementos finitos isoparamétricos com funções de interpolação de ordem superior.
54
Entretanto, o custo computacional necessário para resolver sistemas de equações
oriundos dos elementos finitos isoparamétricos com elevados números de pontos nodais
é relativamente significativo.
Por essa razão, novas formulações continuaram sendo desenvolvidas para que
fosse possível a aproximação da formulação mista (u-p) em qualquer classe de
elementos e, conseqüentemente, qualquer ordem de função de interpolação para o
deslocamento/velocidade e para a pressão. Entretanto, essas formulações devem de
alguma maneira atender ou burlar os critérios de estabilidade de Babuska-Necas-Brezzi
[38], principalmente em relação à condição inf-sup.
De acordo com FRANCA [42], os métodos estabilizados são construídos a partir
da modificação da forma variacional do problema, de tal maneira que o ganho na
estabilidade numérica não comprometa a consistência da aproximação.
No presente trabalho, optou-se pelo emprego das formulações mistas (u-p) GLS
e Bochev, que serão descritas nas seções seguintes. 1
3.4.1 Galerkin Least Square (GLS)
O método Garlerkin least square (GLS) foi inicialmente proposto por HUGHES et
al.[43], e consiste em adicionar na sentença variacional de Garlekin a parcela
correspondente a ponderação por mínimos quadrados, conforme (3.46).
*[ ] [ ] 0
TT d dδ δΩ Ω
− Ω+ − Ω =∫ ∫u A(u) b A(u) τ A(u) b (3. 46)
A primeira parcela da sentença (3.46) é idêntica à formulação (3.28). Deve-se,
agora, estabelecer a sentença variacional para a parcela de mínimos quadrados,
conforme expressão (3.47):
*| | [div ( , )] [div ( , ) ]T Td p dδ δ
Ω
− Ω = + Ω∫ ∫Ω
δA(u) τ A(u) b σ u τ σ u p b (3. 47)
Nessa formulação, τ corresponde a uma matriz de estabilização, definida por
HUGHES et al. [43] na forma: * valor em módulo
55
IτGh
2
2α−=
(3. 48)
Na expressão (3.48), α corresponde a um parâmetro que deve ser ajustado para a
boa performance do método, sendo recomendado ser de ordem O(1) para elementos
triangulares e quadriláteros. De acordo com HUGHES et al. [43], h é um parâmetro
característico dos elementos da malha, e pode ser representado como: 1/dim( )e
h dΩ
= Ω∫ ,
onde dim corresponde à dimensão do domínio.
Para a aproximação por mínimos quadrados, a função de ponderaçãoδA(u)
pode ser decomposta nas parcelas desviatória e esférica, na forma da equação (3.49). A
formulação empregada nessa seção se baseia na referência [44].
[div ( , )] [div( )]T u T Tuuq p= = +δA(u) σ w D ε I (3. 49)
onde: T
uT
u )()( Nδdε ∇= (3. 50)
Tp
Tpp )()( NδP= (3. 51)
Utilizando as expressões (3.50) e (3.51) em (3.49) obtém-se:
[div ( , )] ( ) [div( ) ] ( ) div( )T u T T T T Tu uu u p pq ⎡ ⎤= = ∇ + ⎣ ⎦δA(u) σ w δd D N δP mN (3. 52)
De maneira similar, pode-se reescrever a sentença residual na forma forte A(u)
de acordo com a expressão:
div ( ) div( ) div( )Tuu u p uu upp= + = ∇ + + = + +A(u) σ u, b D N u mN P b L u L P b (3. 53)
onde:
div( )Tuu uu u= ∇L D N (3. 54)
div( )up p=L mN (3. 55)
Finalmente, pode-se escrever a sentença da ponderação por mínimos quadrados
na forma da equação (3.56).
56
[( ) ( ) ( ) ( ) ] [ ]T T T T Tu uu P up uu upd d
Ω Ω
Ω = + + + Ω∫ ∫δA(u) τA(u) δd L δP L τ L u L P b
[( ) ( ) [ ]
[( ) ( ) [ ]
T Tu uu uu up
T Tp up uu up
d
dΩ
Ω
= + + Ω+
+ + + Ω
∫
∫
δd L τ L u L P b
δP L τ L u L P b
(3. 56)
A formulação GLS é então obtida somando a parcela de ponderação do resíduo
(sentença variacional) dos método de Galerkin (3.32) e de mínimos quadrados (3.56):
( ) ( ) ( ) ( )T T T Tuu up p uu uu uu up
T T Tu u uu
d d d d
d d dΩ Ω Ω Ω
Ω Ω Ω
Ω + Ω + Ω + Ω
= ⋅ Ω+ ⋅ Ω − Ω
∫ ∫ ∫ ∫
∫ ∫ ∫
B D B u B D N P L τL u L τL P
(N ) b (N ) t (L ) τb
(3. 57)
A parcela corresponde à equação de conservação de massa (3.8) para o GLS
assume a forma:
( ) div 0T T T
p up uu upp d dK
δ δΩ Ω
⎛ ⎞− Ω+ + + Ω =⎜ ⎟⎝ ⎠∫ ∫TP (N ) u ( P) (L ) τ[L u L P b]
(3. 58)
Combinando as expressões (3.57) e (3.58) pode-se chegar a uma representação
matricial para o modelo GLS, conforme apresentado pelo conjunto de equações (3.59).
T
⎡ ⎤ ⎡ ⎤+ + +⎧ ⎫=⎢ ⎥ ⎢ ⎥⎨ ⎬+ +⎢ ⎥ ⎢ ⎥⎩ ⎭⎣ ⎦ ⎣ ⎦
G G Guu uu up up u u
G G Gup up pp pp p
K K K K R Ru(K K ) K K RP
(3. 59)
As matrizes correspondentes à expressão (3.59) podem ser calculadas a partir do
conjunto de equações (3.60):
e
dΩΩ
⎡ ⎤ =⎣ ⎦ ∫G Tuu uu uuK (L ) τL
edΩ
Ω
⎡ ⎤ =⎣ ⎦ ∫G Tup uu upK (L ) τL
57
edΩ
Ω
⎡ ⎤ =⎣ ⎦ ∫G Tpu up uuK (L ) τL
èdΩ
Ω
⎡ ⎤ =⎣ ⎦ ∫G Tpp up upK (L ) τL
ed
Ω
⎡ ⎤ = − Ω⎣ ⎦ ∫TG
u uuR (L ) τb
edΩ
Ω
⎡ ⎤ = −⎣ ⎦ ∫TG
p upR (L ) τb
Onde:
* *1
ene
e=
⎡ ⎤= ⎣ ⎦G GK K∪
** **1
en e
e=
⎡ ⎤= ⎣ ⎦G GR R∪
(3. 60)
O método GLS é bastante difundido na literatura, principalmente em aplicações
da dinâmica dos fluidos computacional (CFD). A estabilidade do método é garantida
pela ponderação entre a soma dos resíduos dos métodos de Galerkin e de mínimos
quadrados. De acordo com BARTH et al. [45], essa formulação tem obtido sucesso por
ser uma das únicas formulações mistas a garantir a conexão entre as modificações da
sentença variacional por penalização a um problema de otimização do funcional
Lagrangiano na forma mista, o que aproxima a solução a um ponto de sela (saddle-
point) (u-p), garantindo dessa forma acuradas aproximações para a
velocidade/deslocamento e a pressão.
Um fator complicador nessa formulação está relacionado à definição do
parâmetro de estabilização α , que possui um valor ótimo para cada classe de problema
e requer testes numéricos para sua calibração, sendo definido globalmente para todos os
elementos da malha. Testes numéricos revelam uma influência significativa do
parâmetro de calibração na precisão da aproximação, sendo a interferência ainda maior
em malhas não estruturadas. Para contornar essa dificuldade, no presente trabalho se
investiga uma estratégia para determinação do parâmetro de estabilização de forma
automática a nível de elemento e com base em parâmetros geométricos da malha,
conforme será discutido no capítulo 4.
3.4.2 Estabilização direta de pressão (Bochev)
58
O método proposto por DOHRMANN e BOCHEV [46] é obtido através da
modificação da sentença variacional mista através do uso da projeção da pressão local
em L2. Nesse método, a estabilidade resulta de uma relação especial entre os espaços da
velocidade e pressão, que assegura que o espaço da pressão coincida com o intervalo do
operador do espaço divergente da velocidade. De acordo com DOHRMANN e
BOCHEV [46], essa relação é equivalente à condição inf-sup.
A estabilização do método é obtida sem o uso da sentença variacional da
equação de momentum, usando apenas a equação de conservação de massa. Dessa
forma, não há necessidade de cálculo de derivadas de ordens superiores para a
implementação desse esquema. Detalhes do fundamento teórico podem ser consultados
na referência [46].
De acordo com ZIENKIEWICZ e TAYLOR [30], esse método parece estar
totalmente correto e é obtido sem ignorar nenhum termo da sentença variacional global.
A sentença variacional mista da equação de momento é mantida intacta (3.61). Para
garantir a estabilidade, DOHRMANN e BOCHEV [46] introduzem na sentença
variacional da parcela de conservação de massa uma expressão dependente da pressão p
(aproximada por uma função de forma contínua de grau de polinômio k ), e da projeção
da pressão p em um espaço polinomial de ordem (k-1), conforme expressão (3.62).
0=Γ−Ω−Ω+Ω ∫∫∫∫ΓΩΩΩ
ddpdd tδubδumδεεDδε TTTd
T (3. 61)
( ) ( ) Ω−−−Ω⎥⎦⎤
⎢⎣⎡ − ∑∫∫ Ω
Ω
dppG
ppdpK
pe e
αδδδ 1mε (3. 62)
Na expressão (3.62), α corresponde a um parâmetro de estabilização a ser
ajustado para boa performance da formulação. Nota-se claramente que esse parâmetro
não é dependente de parâmetros da malha. DOHRMANN e BOCHEV [46] sugerem o
uso de α = 1. Já ZIENKIEWICZ e TAYLOR [30], com base em exemplos numéricos,
sugerem o valor α = 2 para ótima performance da aproximação.
De acordo com ZIENKIEWICZ e TAYLOR [30], à parcela correspondente a
estabilização de pressão, a nível de elemento, é garantida pela expressão (3.63), e
apresenta baixo custo computacional, pois não envolve derivadas de ordem superiores.
59
( ) 0=Ω−∫Ω dppp
e
δδδ (3. 63)
Na formulação (3.63) a pressão é aproximada pela expressão (3.64), onde Na
contém o conjunto de polinômios de ordem k:
pN~~ˆ ==≈ ∑
aaa pNpp (3. 64)
Por outro lado, a pressão projetada p é calculada com base no conjunto de
polinômios hb(x) de ordem (k-1) , conforme expressão (3.65).
ββ h(x)=≈ ∑b
bb xhp )( (3. 65)
Com base nas equações (3.64) e (3.65), pode-se estabelecer uma relação entre a
pressão p e p , de forma que:
pNhhh TT ~∫∫
ΩΩ
Ω=Ωee
dd β (3. 66)
pGH ~=β (3. 67)
Finalmente, a projeção da pressão pode ser calculada pela expressão (3.68):
pGh(x)H 1 ~−=p (3. 68)
Com base na expressão (3.68) e no desenvolvimento das sentenças variacionais
(3.61) e (3.62), chega-se à aproximação para o deslocamento/velocidade e pressão,
conforme expressões matriciais (3.69) e (3.70) [30].
⎭⎬⎫
⎩⎨⎧
=⎭⎬⎫
⎩⎨⎧⎥⎦
⎤⎢⎣
⎡− 0
fpu
VCCK
Td
~~
(3. 69)
e d
Ω
= Ω∫ Td dK B D B
e dΩ
= Ω∫ TC B mN
(3. 70)
60
1e dK G G
α α −
Ω
⎡ ⎤= + Ω−⎢ ⎥⎣ ⎦∫ T 1V N N GH G
onde:
1
ene
e=
=d dK K∪
1
ene
e=
=C C∪
1
ene
e=
=V V∪
ZIENKIEWICZ e TAYLOR [30] discorrem que para a matriz V, quando G/α é
muito menor que K, o efeito da estabilização direta da pressão corresponde a uma forma
de penalidade equivalente à diferença entre a pressão interpolada e a projetada. Além
disso, para as condições de existência e estabilidade, nota-se que o fato de se utilizar a
pressão projetada de ordem (k-1) por si só equivale à condição inf-sup.
61
4. IGLS: Uma estratégia para determinação do parâmetro de estabilização α do método GLS
A grande dificuldade nas aplicações do método GLS proposto por HUGHES et
al. [43] está relacionada à definição do parâmetro de estabilização α, que assume um
valor ótimo para cada classe de problema.
Essa dificuldade passa a ser mais evidente em problemas que envolvem malhas
não estruturadas, principalmente em casos onde haja grande contraste entre as áreas
características dos elementos e em casos de distorções excessivas da malha.
Para contornar essas dificuldades, neste capítulo investiga-se uma formulação
para a determinação do parâmetro α a nível de elemento, tomando por base as
propriedades físicas do problema e a característica geométrica de cada elemento da
malha. A formulação aqui apresentada foi desenvolvida em CARMO et al. [47].
Para facilitar a compreensão do desenvolvimento teórico, reapresenta-se a
formulação do problema de Stokes, cujo objetivo é encontrar o par de solução do campo
de deslocamento e pressão (u, p), pertencente ao espaço 1 2( ) ( )nH LΩ × Ω satisfazendo:
ii
i fxpuG =
∂∂
+∇⋅∇− )( em Ω para i = 1, n*2 (Conservação de momento) (4. 1)
0=⋅∇+ uBp em Ω (Conservação da massa) (4. 2)
ii gu = em gΓ para i = 1, n* (Condição de Dirichlet) (4. 3)
rnIu =−∇ )( pG em qΓ (Condição de Newman) (4. 4)
Na formulação acima, G corresponde ao módulo elástico cisalhante (para o caso
da elasticidade) ou a viscosidade (para o caso de fluidos) (4.5), e junto com λ (4.6)
define as chamadas constantes de Lamé. A partir das mesmas pode-se definir a relação
(4.7) para B:
)1(2 ν+=
EG (4. 5)
)21)(1( νννλ−+
=E
(4. 6)
λ+=
GB 1
(4. 7)
* n = 2 para problemas bidimensionais; n = 3 para problemas tridimensionais
62
Nota-se claramente que quando 5,0→ν o valor de ∞→λ e,
conseqüentemente, 0→B .
Conforme discutido no capítulo 3, o método GLS consiste em adicionar na
formulação de ponderação de Galerkin por uma sentença de ponderação em mínimos
quadrados, assumindo a forma da expressão (4.8):
Ω
⎥⎥⎥⎥
⎦
⎤
⎢⎢⎢⎢
⎣
⎡∂∂
+∇⋅∇−
⎥⎦
⎤⎢⎣
⎡∂∂
+∇⋅∇−∫ ∑Ω =
dG
xpuG
xpuGh i
in
i iie
e
e
)()()(
1
2α
(4. 8)
Pode-se reescrever em função da parcela da pressão a identidade (4.9), pela
adição e subtração do termo difusivo:
)()( ii
ii
uGxpuG
xp
∇⋅∇+∂∂
+∇⋅−∇=∂∂ para i = 1, n
(4. 9)
Elevando-se os lados esquerdo e direito da expressão (4.9) ao quadrado,
observa-se que a relação (4.10) é verdadeira:
( ) ( )
( )⎥⎥⎦
⎤
⎢⎢⎣
⎡∇⋅∇+⎟⎟
⎠
⎞⎜⎜⎝
⎛∂∂
+∇⋅∇−
≤∇⋅∇+⎟⎟⎠
⎞⎜⎜⎝
⎛∂∂
+∇⋅∇−+∇⋅∇−=⎟⎟⎠
⎞⎜⎜⎝
⎛∂∂
22
222
)()(2
)()(2)(
ii
i
ii
iii
uGxpuG
uGxpuGuG
xp
(4. 10)
Realizando a multiplicação de cada termo da expressão (4.10) por uma parcela
constante, a nível de elemento, chega-se à seguinte relação:
( ) ( )22
2222
22
)()(2
)()(2
)(2
ie
e
iie
e
ie
e
uGhGx
puGhG
xph
G
∇⋅∇+⎟⎟
⎠
⎞
⎜⎜
⎝
⎛⎟⎟⎠
⎞⎜⎜⎝
⎛∂∂
+∇⋅∇−
≤⎟⎟⎠
⎞⎜⎜⎝
⎛∂∂
αα
α
(4. 11)
A partir da desigualdade (4.11), são estabelecidas as seguintes conclusões:
63
• O primeiro termo da desigualdade, 2
2)(2 ⎟⎟
⎠
⎞⎜⎜⎝
⎛∂∂
ie
e
xph
Gα , corresponde à
parcela de estabilização da pressão;
• A parcela ( )22
22 )()(2 ⎟
⎟
⎠
⎞
⎜⎜
⎝
⎛⎟⎟⎠
⎞⎜⎜⎝
⎛∂∂
+∇⋅∇−i
ie
e
xpuGh
Gα corresponde ao termo
fornecido pelo GLS;
• A parcela ( )22 )()(2 ie
e
uGhG
∇⋅∇α corresponde a um termo não fornecido
pela sentença variacional do GLS.
No caso de incompressibilidade total ( 5,0=ν ), nota-se que a parcela
relacionada ao divergente do campo de velocidade/deslocamento tende a ser anulado na
equação de conservação de massa. Conseqüentemente, 0)( =∇⋅∇ iuG , e, portanto,
0eα > é exigido apenas para melhoria da aproximação, mas não para a estabilização da
pressão, que é estabilizada para qualquer 0≠eα .
Para um eα pequeno o suficiente, a conseqüência para os espaços de polinômios
a nível de elemento na sentença variacional de Galerkin pode ser reescrita na seguinte
desigualdade:
( ) ( ) Ω∇=Ω∇⋅∇<Ω∇⋅∇ ∫∫∫ ΩΩΩduGduuGduG
Gh
eeeiiii
ee
222
)()(
ββα
(4. 12)
Como 10 << β , é observada uma diminuição na parcela correspondente à parte
difusiva na sentença variacional de Galerkin. Dessa forma, a estabilidade do GLS é
conseguida por um empréstimo da estabilidade da parcela difusiva para o campo de
pressão.
O que deve ser feito então é determinar um eα que consiga ponderar de forma
eficiente essa distribuição numérica da parcela que irá garantir a estabilidade da pressão
no método GLS. Esse procedimento pode ser realizado pela estratégia descrita a seguir:
Sendo e1η , e
2η , ... , enpelη as funções de forma nos npel pontos nodais do
elemento eΩ , podem-se definir as funções:
64
enpele
npel
eie
iei d
d
e
e ηη
ηηϕ
Ω
Ω−=∫∫Ω
Ω ,
com i = 1, 2, ..., (npel-1)
(4. 13)
Definindo o espaço ∑ −
=ℜ∈=ΨΨ=Ω
1
1*, ;;)( npel
i ieiie
k AAP ϕ , nota-se que
eλ corresponde ao máximo da seguinte relação limitada superiormente:
( ) [ ]
∫∫
Ω
Ω
≠Ω∈ ΩΨ∇
ΩΨ∇⋅∇
=
e
e
ek dG
dG
Ghe
P
e2
22
0),(
)(
Sup*, ψψ
λ
(4. 14)
Se a relação (4.14) for verdadeira, a seguinte sentença deve ser obedecida:
( ) [ ]Ω
Ψ∇⋅∇≥ΩΨ∇ ∫∫ Ω
Ω
dG
GhdG
ee
ee22
2 )(λ *, ( )n
eP∀Ψ∈ Ω
(4. 15)
Multiplicando ambos os termos da expressão (4.15) por 0eα > , obtém-se:
( ) [ ] ΩΨ∇⋅∇≥ΩΨ∇ ∫∫ ΩΩ
dGGh
dGe
e
ee
ee 22
2 )(α
λα *, ( )neP∀Ψ∈ Ω
(4. 16)
Para que ocorra um maior equilíbrio entre a estabilização do campo u e p adota-
se:
21
=eeλα
(4. 17)
Portanto, pode-se finalmente estabelecer uma relação para eα nos moldes da
expressão (4.18) em função de eλ , que traduzirá informações relativas ao grau de
distorção do elemento.
0,5eα = se 810−≤eλ
1min 0,5;2( 0,25)
eeα
λ⎧ ⎫
= ⎨ ⎬+⎩ ⎭ se 810−>eλ
(4. 18)
65
É necessário, portanto, determinar eλ para obter informações da geometria do
elemento e, conseqüentemente, estabelecer o valor para o parâmetro de estabilização eα . Isso é feito pela definição do funcional (4.19):
( ) [ ]
∫∫
Ω
Ω
ΩΨ∇
ΩΨ∇⋅∇
=Ψ
e
e
dG
dG
Ghe
2
22 )(
)(λ *, ( )neP∀Ψ∈ Ω
(4. 19)
Pode-se então estabelecer um valor ótimo para eλ através da maximização do
funcional (4.19):
0)(=
∂Ψ∂
iAλ com i = 1, 2, ..., (npel-1)
(4. 20)
Para facilitar a notação resultante, definem-se as matrizes gradeij
,M e diveij,M ,
expressas respectivamente pelas equações (4.21) e (4.22)
Ω∇⋅∇= ∫Ω
dGM ej
ei
gradeij
e
ϕϕ, com i = 1, 2, ..., (npel-1) (4. 21)
( ) Ω∇⋅∇∇⋅∇= ∫Ω
dGGhM ej
eie
diveij
e
)]()][([2, ϕϕ com i = 1, 2, ..., (npel-1) (4. 22)
Do resultado da maximização (4.20) do funcional (4.14), conclui-se que eλ é o
maior auto-valor da seguinte relação:
( ) ( ) λXXMM =− dive
ijgrade
ij,1,
Com X = (X1, X2 , ..., Xnpe-1)
(4. 23)
66
5. Modelo viscoelástico incompressível
A propriedade intrínseca do material viscoelástico de ir adquirindo um
comportamento incompressível ao longo tempo, sobretudo pelos efeitos da fluência e
relaxação, discutidos no capítulo 2, tem despertado ao longo dos anos o interesse por
formulações numéricas que garantam consistência e estabilidade numérica para a
modelagem dos materiais viscoelásticos.
Conforme discutido no capítulo 3, a técnica de aproximação numérica utilizada
nesse trabalho é baseada no método dos elementos finitos (MEF). Será apresentada
nesse capítulo uma extensão dos conceitos de viscoelasticidade apresentados no capítulo
2. Em seguida, na seção 5.2, é discutida a formulação em elementos finitos para a
viscoelasticidade linear. Na seção seguinte, são abordadas algumas considerações sobre
a estratégia numérica adotada para a implementação computacional para a
viscoelasticidade incompressível.
5.1. Extensão dos modelos viscoelásticos unidimensionais
Os modelos reológicos unidimensionais apresentados no capítulo 2 podem ser
estendidos facilmente para as formulações bidimensional e tridimensional. Dessa forma,
os valores da tensão total σ, deformação total ε, tensão desviatória s e deformação
desviatória e, até então idealizados como escalares, passam a ser interpretados como
tensores, expresso num arranjo vetorial, conforme representado pela expressão (5.1)
para 3 dimensões e (5.2) para 2 dimensões.
[ ]Txzzxzyyzyxxyzyx σσσσσσσσσ=σ
[ ]Txzzxzyyzyxxyzyx εεεεεεεεε=ε
[ ]Txzzxzyyzyxxyzyx sssssssss=s
[ ]Txzzxzyyzyxxyzyx eeeeeeeee=e
(5.1)
67
[ ]Txyyx σσσ=σ
[ ]Txyyx εεε=ε
[ ]Txyyx sss=s
[ ]Txyyx eee=e
(5.2)
De acordo com SIMO e HUGHES [24], os efeitos da viscoelasticidade são
desprezíveis na parcela hidrostática do campo de tensão e de deformação, sendo estes
admitidos como puramente elásticos. Dessa forma, o efeito da viscoelasticidade
formulado neste trabalho se restringirá à parcela desviatória. Partindo da equação
integral do modelo viscoelástico linear apresentado na equação (2.57), pode-se
reescrever a expressão para a parcela desviatória da tensão em função da deformação
desviatória e na forma da equação (5.3), conforme discutido por ZIENKIEWICZ e
TAYLOR [48].
⎥⎦⎤
⎢⎣⎡
∂∂
−= ∫ ∞−τ
ττ dtGt
t es )(2)( (5.3)
As análises aqui desenvolvidas serão restritas ao modelo do sólido viscoelástico
linear padrão (fig 5.1). Através dessa formulação, será possível particularizar a
representação final tanto para um fluido viscoelástico (modelo de Maxwell) quanto para
um sólido viscoelástico (sólido viscoelástico padrão). Partindo da forma integral (5.3),
nota-se que é necessário obter uma função de relaxação G(t-τ), de tal maneira que a
expressão diferencial do modelo do sólido viscoelástico linear apresentado na equação
(2.41) seja atendida. Uma representação da função de relaxação que atende à equação
diferencial resultante é apresentada na expressão (5.4), expressa agora em termos do
módulo elástico cisalhante G0. Note que a expressão assume a mesma forma da função
de relaxação apresentada na equação (2.48).
[ ]0 0 1( ) exp( / )G t G tμ μ λ= + −
onde:
210 GGG +=
210 EEE +=
(5.4)
68
0
1
0
10 G
GEE
==μ
0
2
0
21 G
GEE
==μ
2
2
Eηλ =
Figura 5. 1 – Modelo viscoelástico linear padrão
Generalizando o modelo apresentado na figura (5.1), através de uma associação
de N elementos de Maxwell ligados paralelamente a um elemento elástico (fig. 5.2),
observa-se que a tensão resultante nesse elemento assume a forma da expressão (5.5),
onde qi representa a deformação parcial advinda de cada elemento de Maxwell.
∑=
−=N
i
E1
0 iqεσ (5.5)
Figura 5. 2 – Modelo sólido viscoelástico linear generalizado
Tomando-se como referência a expressão (5.4) e a equação (2.41), redefine-se
uma expressão diferencial para o sólido viscoelástico padrão generalizado. Nesse caso,
com base na solução da equação diferencial resultante, conclui-se que a função de
relaxação deve assumir a seguinte forma:
69
0 01
( ) exp( / )N
N Ni
G t G tμ μ λ=
⎡ ⎤= + −⎢ ⎥
⎣ ⎦∑ , onde:
01
1N
iiμ μ
=
= −∑
00 G
G∞=μ
0GGN
N =μ
N
NN E
ηλ =
(5.6)
De acordo com ZIENKIEWICZ e TAYLOR [48], a formulação (5.6) equivale a
uma expressão da série de Dirichelet-Prony. Dessa forma, as formulações integral e
diferencial do modelo viscoelástico são equivalentes a partir do uso da série de Prony
como função de relaxação. Na prática, os parâmetros associados à série de Prony são
ajustados a partir do ensaio de relaxação.
A solução de modelos viscoelásticos em elementos finitos através da série de
Prony é bastante difundida, sendo inclusive incorporada em programas comerciais, a
exemplo do ABAQUS®, ANSYS® e FEAP®. Maiores detalhes sobre a estratégia
numérica utilizada nesses programas podem ser observados nas referências [37], [49] e
[50].
5.2. Formulação em elementos finitos para a viscoelasticidade linear
incompressível
A estratégia numérica utilizado no desenvolvimento deste trabalho se baseia na
formulação apresentada por ZIENKIEWICZ e TAYLOR [48] e SIMO e HUGHES [24].
A idéia central é transformar a integral de convolução (5.3) em esquema recursivo
passo-a-passo, envolvendo variáveis internas que permitam efetuar a integração
numérica a partir do histórico de tensões e deformações em cada incremento de tempo.
Dessa forma, em cada intervalo de tempo obtém-se uma configuração típica de uma
formulação linear, onde em cada passo subseqüente ocorre a atualização incremental
das variáveis intrínsecas ao processo viscoelástico.
70
Incorporando a expressão da função de relaxação (5.6) na equação integral (5.3)
é possível reescrever a componente da tensão desviatória em função dos parâmetros da
série de Prony.
ττ
λτμμ dtGtt M
nn∫ ∑∞−
= ∂∂
⎥⎦
⎤⎢⎣
⎡−−+=
es1
100 )/)(exp(2)( (5.7)
Para o modelo do sólido viscoelástico linear padrão, particularizado a partir da
expressão (5.7) fazendo-se M=1 tem-se:
ττ
λτμμ dtGtGtt
∫ ∞− ∂∂
−−+=eεs )/)(exp(2)(2)( 11000
(5.8)
A partir da introdução da variável interna h(t) a expressão (5.8) pode ser
reescrita na forma:
ττ
λτ dttt
∫ ∞− ∂∂
−−=eh )/)(exp()( 1
[ ]0 0 1( ) 2 ( ) ( )t G t tμ μ= +s ε h
(5.9)
De acordo com ZIENKIEWICZ e TAYLOR [48], nas aplicações dos modelos
viscoelásticos assume-se que inicialmente o material não está submetido à ação de
carregamento. No instante t0, o corpo viscoelástico sofre um valor deformacional
instantâneo, apresentando variação em função dos tempos subseqüentes. Para calcular a
deformação nos instantes seguintes à aplicação da carga, recorre-se ao artifício:
∫∫∫∫∫+
+
+
−
−+
+++=∞−∞−
11
(.)(.)(.)(.)(.)0
0
0
0 n
n
nn t
t
tt
ddddd τττττ (5.10)
Pelas condições impostas, nota-se que o primeiro termo do lado direito é nulo. O
segundo termo se relaciona à deformação instantânea (elástica), enquanto os últimos 2
termos cobrem a história das deformações. O resultado da expressão (5.10) aplicado a
71
(5.9) explicitado em termos do instante t(n+1) resulta em uma expressão que pode ser
representada pela combinação das equações (5.11) e (5.12).
hhh Δ+Δ−=+ nn t )/exp( 11 λ (5.11)
∫+
∂∂
−−=Δ +1
11 )(exp(n
n
t
t nn dtt ττ
λ eh (5.12)
Para obter uma solução numérica da expressão (5.12), pode-se adotar que a taxa
de formação em cada intervalo de tempo Δt é constante, sendo expressa em termos da
diferença finita progressiva (5.13), conforme discutido por ZIENKIEWICZ e TAYLOR
[48] e SIMO e HUGHES [24]. Aplicando (5.13) na equação (5.12) é possível
representar a variável interna hΔ na forma da equação (5.14).
tnn
Δ−
=∂∂ + eee 1
τ (5.13)
∫+
−−−Δ
=Δ +++1 ]][/)(exp[1
1111n
n
t
t nnnnn dttt
τλ eeh (5.14)
Integrando a expressão (5.14) em função do valor médio do intervalo [tn,tn+1],
obtém-se o valor da variável h(t) indicada na equação (5.9) para cada incremento de
tempo t(n+1) sob a forma da equação (5.15).
*
1 1 1 1exp( / ) ( )n n n n nt hλ+ + += −Δ + Δ −h .h e e
onde
)]/exp(1[ 11
1* λ
λ tt
h n Δ−−Δ
=Δ +
(5.15) (5.16)
ZIENKIEWICZ e TAYLOR [48] sugerem que para pequenos intervalos de
tempo Δt, a expressão (5.16) deve ser substituída pela solução da série apresentada na
equação (5.17), visto que para valores de Δt próximos a zero a função (5.16) vai
assumindo comportamento singular.
72
...!4
1!3
1211
3
1
2
11
1* +⎟⎟
⎠
⎞⎜⎜⎝
⎛ Δ−⎟⎟
⎠
⎞⎜⎜⎝
⎛ Δ+⎟⎟
⎠
⎞⎜⎜⎝
⎛ Δ−=Δ +
λλλttth n
(5.17)
A relação constitutiva viscoelástica apresentada na equação (5.9) pode então ser
discretizada para cada intervalo de tempo t(n+1) em função da variável interna (5.15),
de forma que a tensão em cada passo subseqüente é dada por:
⎟⎠
⎞⎜⎝
⎛+= ∑
=+++
M
m
mnmnn G
1
)(11001 ...2 hes μμ
(5.18)
Do ponto de vista computacional, a partir da formulação de elementos finitos
discutida no capítulo 3, é necessário calcular a matriz constitutiva Cn+1** para o modelo
viscoelástico em estudo, que pode ser obtida através da diferenciação do tensor de
tensões desviatório apresentado na equação (5.18) em relação à deformação total,
conforme estabelecido por ZIENKIEWICZ e TAYLOR [48] na forma da equação
(5.19). 3
dIes
εsC
1
1
1
11
+
+
+
++ ∂
∂=
∂∂
=n
n
n
nn
(5.19)
A derivada parcial da expressão da parcela desviatória do tensor de tensões
assume a forma apresentada na equação (5.20), onde I corresponde ao tensor identidade.
Para exemplificar o processo de obtenção analítica das expressões, restringe-se o
número de termos para M=1, correspondente a um modelo do sólido viscoelástico
padrão. Entretanto, para a extensão da formulação para M>1, deve-se calcular o valor
da variável interna hn+1 para cada M-ésimo elemento de Maxwell ligado ao modelo,
alocando a contribuição de cada parcela de forma linear.
⎥⎦⎤
⎢⎣⎡
∂∂
+=∂∂
ehI
es
1002 μμG (5.20)
Por outro lado, derivando a expressão (5.15) com relação a deformação
desviatória tem-se: *** Para manter a notação original de ZIENKIEWICZ e TAYLOR [48], a matriz constitutiva Duu passa agora a ser representada por Cn
73
Iheh )(1
1
1 tnn
n ΔΔ=∂∂
++
+ (5.21)
Combinando as equações (5.20) e (5.21), pode-se reescrever a expressão (5.19)
sob a forma:
[ ] dnn
nn tG Ih
εsC )(2 1100
1
11 ΔΔ+=
∂∂
= ++
++ μμ
(5.22)
ZIENKIEWICZ e TAYLOR [48] discorrem que a formulação para matriz
constitutiva apresentada na equação (5.22) difere da formulação elástica linear apenas
na parcela do módulo de cisalhamento, que agora passa a ser dependente dos parâmetros
viscoelásticos e da variável interna.
O tensor de tensões σ em cada passo de tempo tn+1 é atualizado a partir da
parcela desviatória, conforme apresentado pela equação (5.18). Dessa forma, é
necessário efetuar a atualização do vetor de forças internas, conforme formulação
apresentada por ZIENKIEWICZ e TAYLOR [48] e SIMO e HUGHES [24]:
Ω= ∫Ω ++ dn
Tnvisc .1
1 σBf (5.23)
Objetivando uma generalização dos conceitos discutidos, será apresentada uma
formulação alternativa para a viscoelasticidade, baseada na função de energia
armazenada Wo(ε), conforme apresentado por SIMO e HUGHES [24].
Essa função de energia Wo(ε) pode também ser decomposta em uma parcela
desviatória e volumétrica, de forma que:
)()(~)( 000 Θ+= UWW eε (5.24)
onde Θ corresponde ao traço do tensor de deformação ε.
A resposta elástica volumétrica é caracterizada pela parcela U0(Θ). Pode-se
definir também uma resposta elástica instantânea para o tensor de tensão, dado pela
relação:
εεεεσ ε ∂
∂≡∂=
)(~)(~)(
000 WW
(5.25)
74
Com base no elemento apresentado na figura (5.2), e de maneira similar à
obtenção da expressão (5.5), pode-se reescrever o tensor de tensão sob a forma:
∑=
−=N
iitt
1
0 )()( qσσ (5.26)
A equação constitutiva para o modelo viscoelástico discutido previamente
também pode ser reescrito com base na função da energia Wo, resultando na expressão
(5.27).
( )∫ ∞−∂−+Θ=
tdssW
dsdstGUt )]([~dev)()(')( 00 eIσ e
(5.27)
Similarmente ao procedimento realizado para a obtenção do esquema recursivo
passo-a-passo apresentado anteriormente, adotando-se a série de Prony como função de
relaxação, pode-se a partir da expressão (5.27) redefinir a variável interna h(t), agora
expressa em termos da função de energia:
( )∫ ∞−∂−−=
tdssW
dsdstt )]([~dev]/)(exp[)( 0 eh eτ
(5.28)
Dessa forma, em cada n-ésima iteração temporal a variável interna h(t) é
atualizada com base na equação (5.29), obtida de maneira similar à expressão discutida
em (3.14), através de um esquema de integração de ponto médio em [tn tn+1]. SIMO e
HUGLES [24] enfatizam que a formulação apresentada é incondicionalmente estável e
acurada em segunda-ordem. Entretanto, para elementos de baixa ordem a formulação
está sujeita ao fenômeno de “travamento” próximo ao limite de incompressibilidade.
Por essa razão, é extremamente importante fazer uso da formulação mista em elementos
finitos.
))(2/exp()/exp( 0011 nnnnnn tt sshh −Δ−+Δ−= ++ ττ (5.29)
onde:
)](~[dev 10
1 ++ ∂= nn W es e
(5.30)
][dev 11 ++ = nn εe (5.31)
75
Assumindo-se d/ds [s0(s)] constante no intervalo [tn tn+1], chega-se finalmente a
uma expressão equivalente à forma apresentada anteriormente na equação (5.15):
)(/
)/exp(1)/exp( 00
11 nnn
nnnn t
tt sshh −
ΔΔ−−
+Δ−= ++ ττ
τ (5.32)
Para a matriz constitutiva, a expressão resultante é obtida através da relação
(5.33):
11
11 1 +
+
++ +
∂=∂∂
= nn
nn n
σεσC ε
(5.33)
SIMO e HUGHES [24] apresentam a formulação da matriz constitutiva em
função de uma parcela 01+nC , obtida pela linearização do tensor 0
1+ns (eq. 5.34) e de
alguns termos da variável interna auxiliar hn+1, resultando na expressão (5.35).
( ) ( )[ ] ( )[ ]
( )[ ] IIII
IIIIC
ee
eeeeee
⊗∂+
+⊗∂−∂⊗−∂=
+
++++
:~:91
:~31:~
31~
12
12
12
120
1
n
nnnn
W
WWW
(5.34)
01
1
1''011 )( +
+
+++ Δ+
∂Θ∂
⊗= nn
nnn tGU C
εIC
(5.35)
Nota-se que se o comportamento inicial do material for admitido como elástico
linear, e que os efeitos da viscoelasticidade estão presentes apenas na parcela
desviatória, a expressão (5.35) é equivalente à particularização da formulação (5.22),
discutida anteriormente.
5.3.Estratégia numérica
Conforme discutido anteriormente, a idéia central da estratégia numérica
adotada para a solução do modelo viscoelástico linear consiste na transformação da
expressão (5.3) em um esquema incremental, onde em cada passo de tempo Δt é obtida
76
uma configuração de um modelo linear, que é atualizada a cada etapa de solução através
das variáveis controladores da viscoelasticidade.
Duas estratégias de solução para o modelo linear foram empregadas nesse
trabalho: A primeira se baseia na formulação mista proposta por HUGHES et. al. [43],
conhecida como GLS[25]. A segunda na formulação de estabilização direta da pressão,
proposta por DOHRMANN e BOCHEV [46] .
Adotou-se para implementação das formulações descritas no capítulo 3 o
elemento quadrilátero bilinear apresentado na figura (5.3). O deslocamento/velocidade e
a pressão são aproximados diretamente a partir dos 4 pontos nodais do elemento.
Convém ressaltar que o procedimento de cálculo agora passa a ser efetuado a nível de
elemento, em função das coordenadas locais ξ e η (bidimensional).
Figura 5. 3 – Elemento quadrilátero bilinear
A função de interpolação adotada assume a forma da expressão (5.36), que foi
empregada tanto para a aproximação do deslocamento/velocidade quanto para a
pressão.
)1)(1(41
iiiN ηηξξ ++= i = 1,2,3 e 4 (5.36)
Os valores de iξ e iη estão associados a cada índice i da expressão (5.36) e,
conseqüentemente, a cada ponto nodal do elemento (fig 5.3), na forma apresentada na
tabela seguinte:
77
Tabela 5. 1 – Parâmetros da função de interpolação bilinear i iξ iη1 -1 -1 2 +1 -1 3 +1 +1 4 -1 +1
Nota-se que a expressão (5.36) garante continuidade entre os nós adjacentes dos
elementos finitos tanto no deslocamento/velocidade quanto na pressão. O
comportamento gráfico da função de interpolação para i = 1 é ilustrado na figura (5.4).
Figura 5. 4 – Função de forma bilinear
De acordo com BATHE [32], de forma geral, nos elementos isoparamétricos as
coordenadas locais se relacionam com as coordenadas globais através da matriz
Jacobiana J, relação esta que pode ser facilmente obtida através da aplicação da regra
da cadeia, resultando na expressão (5.37), onde ς representa a coordenada natural na
direção z, nos casos tridimensionais.
⎥⎥⎥⎥⎥⎥⎥
⎦
⎤
⎢⎢⎢⎢⎢⎢⎢
⎣
⎡
∂∂∂∂∂∂
=
⎥⎥⎥⎥⎥⎥
⎦
⎤
⎢⎢⎢⎢⎢⎢
⎣
⎡
∂∂∂∂∂∂
−
ς
η
ξ
i
i
i
i
i
i
N
N
N
zNy
Nx
N
1J
(5.37)
onde:
78
⎥⎥⎥⎥⎥⎥⎥
⎦
⎤
⎢⎢⎢⎢⎢⎢⎢
⎣
⎡
∂∂
∂∂
∂∂
∂∂
∂∂
∂∂
∂∂
∂∂
∂∂
=
∑∑∑
∑∑∑
∑∑∑
ii
ii
ii
ii
ii
ii
ii
ii
ii
zN
yN
xN
zN
yN
xN
zN
yN
xN
ςςς
ηηη
ξξξ
J
(5.38)
Na implementação do método GLS, uma importante etapa está relacionada ao
cálculo das derivadas de segunda ordem das funções de interpolação. A estratégia
adotada baseia-se na formulação apresentada na referência [44], sob a forma:
⎟⎟⎟⎟⎟⎟⎟
⎠
⎞
⎜⎜⎜⎜⎜⎜⎜
⎝
⎛
⎪⎪⎭
⎪⎪⎬
⎫
⎪⎪⎩
⎪⎪⎨
⎧
∂∂∂∂
⎥⎥⎥⎥⎥⎥⎥
⎦
⎤
⎢⎢⎢⎢⎢⎢⎢
⎣
⎡
∂∂∂
∂∂∂
∂∂
∂∂
∂∂
∂∂
−
⎪⎪⎪
⎭
⎪⎪⎪
⎬
⎫
⎪⎪⎪
⎩
⎪⎪⎪
⎨
⎧
∂∂∂∂∂∂∂
⎥⎥⎥⎥⎥⎥⎥⎥
⎦
⎤
⎢⎢⎢⎢⎢⎢⎢⎢
⎣
⎡
∂∂
∂∂
+∂∂
∂∂
∂∂
∂∂
∂∂
∂∂
∂∂
∂∂
⎟⎟⎠
⎞⎜⎜⎝
⎛∂∂
⎟⎟⎠
⎞⎜⎜⎝
⎛∂∂
∂∂
∂∂
⎟⎟⎠
⎞⎜⎜⎝
⎛∂∂
⎟⎟⎠
⎞⎜⎜⎝
⎛∂∂
=
⎥⎥⎥⎥⎥⎥⎥
⎦
⎤
⎢⎢⎢⎢⎢⎢⎢
⎣
⎡
∂∂∂∂∂∂
−
yNx
N
yx
yx
yx
N
N
N
xyyxyyxx
yxyx
yxyx
xyN
yN
xN
i
i
ii
ii
ii
i
i
i
iiiiiiii
iiii
iiii
i
i
i
ηξηξ
ηη
ξξ
ηξ
η
ξ
ξηξηηξηξ
ηηηη
ξξξξ
22
2
2
2
2
2
2
2
2
2
2
2
2
2
1
22
22
2
2
2
2
2
.
.2
2
(5.39)
A integração numérica foi realizada através do método da quadratura de Gauss.
A expressão resultante para a integração tanto no modelo GLS quando no BOCHEV
requer apenas 2 pontos de quadratura de integração. De acordo com BATHE [32], o
número de pontos de integração é função do grau p do polinômio do integrando,
definido pela relação ⎟⎠⎞
⎜⎝⎛ +
21p , para p ímpar; para p par considera-se o número inteiro
imediatamente superior ao resultado da razão descrita.
A correlação, a nível de elemento, entre as integrações nos sistemas de
coordenadas global e local é estabelecida mediante a expressão (5.40), onde det J
corresponde ao determinante da matriz Jacobiana, npi o número de pontos de integração
e w o peso correspondente a cada quadrante de integração de Gauss. Maiores detalhes
podem ser consultados em BATHE [32] e ZIENKIEWICZ e TAYLOR [30].
79
ςηξςηξ
ςηξςηξ
dddwwwh
dddgedzyxf
ijknpii npij npik
iii
e
J
J
det),,(
det),,(),,(
,1 ,1 ,1
1
1
1
1
1
1
∑ ∑ ∑
∫ ∫ ∫∫
= = =
− − −Ω
=
==Ω
(5.40)
Com base nas prerrogativas enfatizadas nas seções preliminares, o modelo
viscoelástico incompressível é obtido por um esquema recursivo passo-a-passo,
recaindo em um modelo linear a cada iteração temporal n, conforme esquema
apresentado em (5.41).
( ) ),(),( 111 pupu nvis
nn +++ −= ffXK δ (5.41)
No caso da viscoelasticidade linear, nota-se que a matriz Kn+1 depende apenas do
intervalo de tempo Δt, conseqüente da dependência da matriz constitutiva Cn+1(ver
equação 5.22).
Em função das hipóteses discutidas nos capítulos preliminares, observa-se que,
para t = 0, a formulação corresponde à solução de um modelo elástico linear. Para tanto,
atribui-se:
00 11 =→= +n
visfh (5.42)
Nas etapas seguintes, a força interna viscoelástica 1+nvisf é calculada com base na
expressão (5.23), que é dependente das tensões e deformações deviatórias e por
conseqüência da variável interna hn+1.
O deslocamento e a pressão em cada etapa n são atualizados através da solução
incremental da expressão (5.41), conforme expressões (5.43) e (5.44). nnn uuu += ++ 11 δ (5.43) nnn ppp += ++ 11 δ (5.44)
Para melhor compreensão dos procedimentos descritos, o algoritmo com os
procedimentos de cálculo é apresentado na figura seguinte:
80
Figura 5. 5 – Algoritmo para o modelo viscoelástico linear
Início das iterações no tempo n
tn+1=tn+Δt
Início do modelo elástico incompressível [GLS ou Bochev]
Cn+1 Motangem da matriz constitutiva (eq. 5.22)
Iteração nos elementos
-Montagem da matriz Kn+1 [GLS ou bochev]
nnn ppp += ++ 11 δ (nodal) – Atualização da pressão nodal
Iteração nos pontos de integração de Gauss
εn+1 = εn + Δε - (Atualização da deformação total)
en+1=Id εn+1 (Deformação desviatória)
hn+1=exp(-Δt /λ).hn + Δh(en-en+1) (Variável interna)
⎟⎠
⎞⎜⎝
⎛+= ∑
=+++
M
m
mnmnn G
1
)(11001 ...2 hes μμ (Tensão desviatória)
pn
np Np 11
++ = (Pressão no ponto de integração)
111 +++ += nnn pmsσ (Tensão total)
Ω= ∫Ω ++ dn
Tn .11
int σBf (Força interna viscoelástica)
Fim da iteração nos pontos de integração de Gauss
Fim das iterações nos elementos
-Computar K e 11 ++ −= nvis
nr fff global;
-Computar solução 1),( += nrspu fKX 1+nuδ e 1+npδ (eq. 5.41)
- nnn uuu += ++ 11 δ (atualizar deslocamento)
Fim das iterações no tempo
81
6. Aplicações e discussões
Este capítulo é voltado à realização de análises das estratégias e métodos
desenvolvidos no trabalho. Nas primeiras aplicações (exemplos 6.1 e 6.2), avaliam-se as
formulações em elementos finitos implementadas frente a situações de quase-
incompressibilidade e incompressibilidade total (fluidos), sobretudo para verificar a
acurácea da solução dos modelos lineares GLS e Bochev, necessários para o processo
de solução do modelo viscoelástico.
Nos exemplos seguintes, investiga-se o comportamento reológico de sólidos e
fluidos viscoelásticos, efetuando a análise comparativa através do software FEAP,
desenvolvido por TAYLOR [50]. Na última aplicação, avalia-se o comportamento
numérico frente à realização de um ensaio de fluência. Em situações em que há
presença de malhas não estruturadas, investigam-se também os resultados obtidos com a
formulação IGLS, proposta no capítulo 4.
Convém ressaltar que a formulação utilizada para o modelo misto no FEAP
corresponde à formulação de 3 campos (u-p-θ) com pressão descontínua entre os nós
adjacentes dos elementos na malha, podendo eventualmente, em função da malha,
apresentar um pequeno desvio quando comparada às formulações de 2 campos (u-p) de
pressão contínua avaliados nesse trabalho.
Todas as unidades físicas utilizadas nos exemplos são do Sistema Internacional
de unidades (SI), e essas unidades não serão explicitadas nos exemplos apresentados.
6.1 Exemplo 1: Viga submetida à ação de momento de flexão
O objetivo dessa primeira análise é avaliar globalmente os modelos lineares
implementados. Para tanto, os modelos serão aplicados para obter a aproximação das
componentes do deslocamento de uma viga submetida à ação de um momento de flexão
em estado plano de deformação, conforme apresentado na figura (6.1). O momento
aplicado pode ser substituído pelo sistema binário de força equivalente (carregamento)
apresentado na figura (6.1).
82
Figura 6. 1 – Viga submetida à ação de um momento de flexão (adaptada de REDDY et.
al.[51])
Esse problema é discutido por REDDY e DJOCO [51], que apresentam a
solução analítica para os deslocamentos vertical e horizontal, explicitados
respectivamente pelas equações (6.1) e (6.2).
⎟⎠⎞
⎜⎝⎛ −
−= ylx
ELfyxu
2)1(2),(
2ν (6. 1 )
( )⎥⎦⎤
⎢⎣⎡ −
−+
−= lyyx
ELfyxv
ννν
1)1(),( 2
2
(6. 2)
Primeiramente, será analisada a convergência dos modelos GLS e Bochev frente
ao grau de refinamento da malha. Adotam-se como parâmetros os seguintes valores: L =
10; l = 2; f = 100; E = 1500.
Os resultados foram extraídos a partir dos valores da solução aproximada e exata
no ponto nodal da extremidade superior da face direita (fig. 6.1). As malhas utilizadas
nas análises apresentam distribuição de refinamento em elementos retangulares, sendo o
domínio subdividido em: 8 elementos (4x2); 32 elementos (8x4); 128 elementos (16x8)
e 512 elementos (32x16). Para fins ilustrativos, a malha de 512 elementos é apresentada
na figura seguinte (fig. 6.2). Nessa primeira análise, optou-se apenas pelo uso de malhas
estruturadas.
Figura 6. 2 – Domínio discretizado em 512 elementos (32x16)
83
Cabe ressaltar que os algoritmos implementados devem apresentar boa
performance tanto para valores elevados de coeficiente de Poisson quanto para valores
mais baixos. Por essa razão, a análise de convergência foi realizada para os valores de
coeficiente de Poisson ν = 0,1 e ν = 0,499999999.
Adotou-se como parâmetro característico da malha a raiz quadrada do número
de elementos. A norma do erro absoluto foi utilizado como critério indicador da
qualidade da aproximação. Nas figuras seguintes são apresentadas as variações do erro
absoluto em função da raiz quadrada do número de elementos, que será um indicativo
da taxa de convergência dos algoritmos em função da variação dos parâmetros de
estabilização.
0.00
0.20
0.40
0.60
0.80
1.00
1.20
1.40
0.0 5.0 10.0 15.0 20.0 25.0
Raiz quadrada do número de elementos
|Err
o ab
solu
to|
Bochev a = 1
Bochev a = 2
GLS a = 0.1
GLS a = 0.5
GLS a = 1.
Figura 6. 3 – Convergência do deslocamento horizontal para ν = 0,499999999.
0.00
1.00
2.00
3.00
4.00
5.00
6.00
7.00
0.0 5.0 10.0 15.0 20.0 25.0
Raiz quadrada do número de elementos
|Err
o ab
solu
to|
Bochev a = 1
Bochev a = 2
GLS a = 0.1
GLS a = 0.5
GLS a = 1.
Figura 6. 4 – Convergência do deslocamento vertical para ν = 0,499999999.
84
0.00
0.10
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.0 5.0 10.0 15.0 20.0 25.0
Raiz quadrada do número de elementos
|Err
o ab
solu
to|
Bochev a = 1
Bochev a = 2
GLS a = 0.1
GLS a = 0.5
GLS a = 1.
Figura 6. 5 – Convergência do deslocamento horizontal para ν = 0,1.
0.00
0.50
1.00
1.50
2.00
2.50
3.00
3.50
4.00
0.0 5.0 10.0 15.0 20.0 25.0
Raiz quadrada do número de elementos
|Err
o ab
solu
to|
Bochev a = 1
Bochev a = 2
GLS a = 0.1
GLS a = 0.5
GLS a = 1.
Figura 6. 6 – Convergência do deslocamento vertical para ν = 0,1.
Observa-se que para um número reduzido de elementos a formulação mista não
apresenta boa performance, principalmente para elevados valores do coeficiente de
Poisson (fig 6.3 e 6.4). Entretanto, aumentando-se o número de elementos (de 8 para
32) já se observa uma melhora significativa na aproximação.
Os testes mostraram que os algoritmos implementados conseguiram evitar o
efeito de travamento característico da formulação clássica de Garlekin para elevados
valores de Poisson em estado plano de deformação. Esse comportamento é apresentado
na figura (6.7) para o deslocamento horizontal e (6.8) para o deslocamento vertical.
Os resultados foram condizentes com a resposta analítica. Para facilitar a
comparação e avaliação de qualidade dos métodos implementados, em função do
parâmetro de estabilização, são apresentadas as tabelas (6.1) e (6.2), que trazem a
informação do erro relativo comparado com a solução analítica. Na tabela (6.3)
85
apresentam-se os resultados obtidos através do FEAP®, sobretudo para que se possa
utilizá-lo como parâmetro de análise de qualidade nos exemplos seguintes.
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5-0.7
-0.6
-0.5
-0.4
-0.3
-0.2
-0.1
0
Coeficiente de poisson
u - d
eslo
cam
ento
hor
izon
tal
AnaliticaBochev a=1Bochev a=2GLS a=0.1GLS a=0.5FEAP-F.MistaFEAP-F.Desl.
Figura 6. 7 – Deslocamento horizontal em função do Poisson (512 elementos)
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.50
0.5
1
1.5
2
2.5
3
3.5
Coeficiente de poisson
v - d
eslo
cam
ento
ver
tical
AnaliticaBochev a=1Bochev a=2GLS a=0.1GLS a=0.5FEAP-F.MistaFEAP-F.Desl.
Figura 6. 8 – Deslocamento vertical em função do Poisson (512 elementos)
86
Tabela 6. 1– Resultado para o modelo de Bochev Poisson Bochev α =1 Bochev α = 2
u v Erro u Erro v u v Erro u Erro v 0 -0,6641 3,3215 0,39% 0,36% -0,6647 3,3243 0,30% 0,27%
0,1 -0,6585 3,2932 0,23% 0,21% -0,6593 3,2969 0,11% 0,09%0,2 -0,6396 3,1983 0,06% 0,05% -0,6406 3,2031 0,10% 0,10%0,3 -0,6073 3,0364 0,11% 0,10% -0,6086 3,0425 0,32% 0,30%0,4 -0,5616 2,8073 0,29% 0,26% -0,5633 2,8150 0,58% 0,53%0,49 -0,5089 2,5435 0,46% 0,41% -0,5109 2,5527 0,84% 0,77%
0,4999 -0,5025 2,5111 0,47% 0,43% -0,5045 2,5205 0,88% 0,80%0,499999 -0,5024 2,5108 0,47% 0,43% -0,5044 2,5202 0,88% 0,80%
0,499999999 -0,5024 2,5108 0,47% 0,43% -0,5044 2,5202 0,88% 0,80%
Tabela 6. 2– Resultado para o modelo GLS Poisson GLS α = 0.1 GLS α = 0.5
u v Erro u Erro v u v Erro u Erro v 0 -0,6639 3,3206 0,41% 0,38% -0,6655 3,3282 0,17% 0,15%
0,1 -0,6582 3,2917 0,27% 0,25% -0,6601 3,3006 0,02% 0,02%0,2 -0,6391 3,1960 0,14% 0,13% -0,6414 3,2064 0,21% 0,20%0,3 -0,6066 3,0332 0,01% 0,01% -0,6092 3,0452 0,42% 0,39%0,4 -0,5606 2,8030 0,12% 0,11% -0,5637 2,8167 0,65% 0,59%0,49 -0,5077 2,5379 0,21% 0,19% -0,5111 2,5534 0,88% 0,80%
0,4999 -0,5012 2,5054 0,22% 0,20% -0,5046 2,5211 0,91% 0,82%0,499999 -0,5011 2,5050 0,22% 0,20% -0,5046 2,5207 0,91% 0,82%
0,499999999 -0,5011 2,5050 0,22% 0,20% -0,5046 2,5207 0,91% 0,82%
Tabela 6. 3– Resultado para FEAP Poisson FEAP – Formulação mista FEAP – Formulação u
u v Erro u Erro v u v Erro u Erro v 0 -0,6618 3,3156 0,73% 0,53% -0,6603 3,3103 0,97% 0,70%
0,1 -0,6573 3,2891 0,41% 0,33% -0,6554 3,2826 0,70% 0,53%0,2 -0,6394 3,1960 0,09% 0,13% -0,6370 3,1879 0,46% 0,38%0,3 -0,6080 3,0358 0,21% 0,08% -0,6049 3,0249 0,30% 0,28%0,4 -0,5628 2,8081 0,49% 0,29% -0,5578 2,7896 0,39% 0,37%0,49 -0,5096 2,5441 0,58% 0,44% -0,4847 2,4343 4,53% 4,05%
0,4999 -0,5027 2,5113 0,52% 0,44% -0,0852 0,4274 486,7% 485,1%0,499999 -0,5026 2,5110 0,52% 0,44% -0,0013 0,0053 3,8.104% 4,7.104%
0,499999999 -0,5023 2,5085 0,45% 0,34% -1,4,10-6 5,3,10-6 3,6.107% 4,7.107%
Com base nos resultados apresentados, observa-se que os modelos Bochev e
GLS removeram eficientemente o efeito do “travamento” característico da formulação
clássica de Galerkin em estado plano de deformação, obtidos através do FEAP®. Nota-
se também que o grau de aproximação depende essencialmente do parâmetro α adotado
nas análises.
O deslocamento resultante (módulo do deslocamento) e a pressão são
apresentados nas figuras a seguir:
87
Figura 6. 9 – Deslocamento resultante para ν = 0,499999999.
Figura 6. 10 – Deslocamento resultante para ν = 0,1.
Figura 6. 11 – Distribuição de pressão para ν = 0,499999999.
Analisa-se agora o comportamento das formulações lineares frente à
aproximação em malhas não estruturadas e com grau de distorção entre os elementos,
sobretudo para avaliar o desempenho das formulações GLS, Bochev e IGLS. Os estudos
foram realizados com malhas subdivididas em 46, 85, 191 e 466 elementos. A relação
entre a raiz quadrada do número de elementos e o erro absoluto obtido na aproximação
dos deslocamentos horizontal e vertical é apresentada, respectivamente, nas figuras
(6.12) e (6.13). Os estudos foram realizados adotando-se ν = 0,499999999.
88
6 8 10 12 14 16 18 20 220
0.02
0.04
0.06
0.08
0.1
0.12
Raiz quadrada do número de elementos
|Erro
abs
olut
o|
Convergencia - Deslocamento horizontal
IGLSBochev a=1.bochev a=2.GLS a=0.1
Figura 6. 12 – Convergência para o deslocamento horizontal (ν = 0,499999999).
6 8 10 12 14 16 18 20 220
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5Convergencia - Deslocamento vertical
Raiz quadrada do número de elementos
|Erro
abs
olut
o|
IGLSBochev a=1.bochev a=2.GLS a=0.1
Figura 6. 13 – Convergência para o deslocamento vertical (ν = 0,499999999).
Os resultados mostraram um grande aumento na sensibilidade do parâmetro de
calibração α para malhas irregulares, principalmente no modelo GLS. Nota-se que a
formulação proposta IGLS apresentou os melhores resultados, tanto em termos de
aproximação quanto em termos de taxa de convergência.
Para facilitar as análises comparativas, são apresentados nas tabelas de (6.4) a
(6.6) os resultados numéricos para a aproximação dos deslocamentos vertical e
89
horizontal oriundos das diversas formulações implementas nesse trabalho,
apresentando-se também o erro relativo obtido em cada aproximação.
Tabela 6. 4– Resultados para a malha de 46 elementos (ν = 0,499999999) Formulação u v erro u erro v
IGLS -5.20E-01 2.59E+00 4.06% 3.72% Bochev α = 1 -5.52E-01 2.76E+00 10.42% 10.25% Bochev α = 2 -6.01E-01 2.99E+00 20.18% 19.72% GLS α = 0.1 -5.71E-01 2.83E+00 14.17% 13.24% GLS α = 0.5 -1.40E+00 6.51E+00 179.81% 160.26% GLS α =0. -4.71E-01 2.33E+00 5.85% 6.85%
Tabela 6. 5– Resultados para a malha de 85 elementos (ν = 0,499999999) Formulação u v erro u erro v
IGLS -5.10E-01 2.55E+00 1.92% 1.96% Bochev α = 1 -5.20E-01 2.60E+00 3.99% 4.14% Bochev α = 2 -4.83E-01 2.43E+00 3.42% 2.67% GLS α = 0.1 -5.22E-01 2.61E+00 4.43% 4.31% GLS α = 0.5 -6.32E-01 3.16E+00 26.41% 26.31% GLS α =0. -4.83E-01 2.43E+00 3.42% 2.67%
Tabela 6. 6– Resultados para a malha de 466 elementos (ν = 0,499999999) Formulação u v erro u erro v
IGLS -5.01E-01 2.50E+00 0.26% 0.11% Bochev α = 1 -5.06E-01 2.53E+00 1.11% 1.05% Bochev α = 2 -5.12E-01 2.56E+00 2.32% 2.27% GLS α = 0.1 -5.05E-01 2.52E+00 1.07% 0.99% GLS α = 0.5 -5.35E-01 2.69E+00 6.94% 7.52% GLS α =0. -4.97E-01 2.48E+00 0.63% 0.63%
Um fato que vale a pena observar está relacionado ao excelente grau de
aproximação do campo de deslocamento obtido pela formulação GLS com α = 0, que
corresponde ao uso apenas da parcela de Galerkin. Entretanto, conforme discutido no
capítulo 3, em situações próximas à incompressibilidade o uso da formulação mista de
Galerkin é sujeita ao chamado modo “pressão espúria”. Esse comportamento ocorre no
exemplo analisado. A distribuição da pressão para GLS, α = 0, é apresentada na figura
seguinte:
90
Figura 6. 14 – Modo pressão espúria para (ν = 0,499999999): (a) 46 elementos; (b) 85
elementos; (c) 191 elementos; (d) 466 elementos.
A formulação IGLS conseguiu contornar com eficiência o modo de pressão
espúria (instabilidades no campo de pressão), bem como obteve grande ganho de
precisão no campo de deslocamento. A distribuição de pressão obtida pelo IGLS para a
malha de 466 elementos é apresentada na figura (6.15).
Figura 6. 15 – Pressão estabilizado pela formulação IGLS (ν = 0,499999999)
91
6.2 Exemplo 2: Escoamento de Stokes em regime estacionário
Nesta análise procura-se avaliar a performance das formulações mistas
implementadas frente à situação de incompressibilidade total. Tradicionalmente,
emprega-se o problema de fluxo de Stokes em regime estacionário, que coincide com o
problema de incompressibilidade da elasticidade linear (ν = 0,5). Esse problema é
descrito por ZIENKIEWICKZ e TAYLOR [30] e BATHE [32].
O problema consiste na representação de um domínio bidimensional em forma
de um quadrado de dimensão unitária. Idealiza-se que todos os lados do contorno
apresentam velocidade nula nas direções vertical e horizontal, exceto no contorno
superior que apresenta velocidade unitária na direção horizontal em todos os nós, com
exceção dos nós das extremidades. Admite-se que o fluido possui viscosidade unitária, e
através da relação de Lamé facilmente pode ser obtido o valor do módulo de
elasticidade correspondente.
Como o problema analisado não possui solução analítica, serão adotados os
mesmos parâmetros utilizados por ZIENKIEWICKZ e TAYLOR [30], sobretudo para
que seja possível efetuar uma análise comparativa.
Tratando-se de um problema totalmente incompressível, e que dispõe somente
de velocidades prescritas no contorno, é necessário especificar o valor da pressão em
um ponto no domínio, conforme discutido no capítulo 3. A geometria e condições de
contorno são apresentadas no modelo conceitual, mostrado na figura (6.16).
O domínio foi discretizado em elementos quadrangulares de comprimento de
lado 0.1, obtendo-se dessa forma uma malha em elementos finitos composta por 100
elementos.
Figura 6. 16 – Caracterização do problema de escoamento de Stokes
92
ZIENKIEWICKZ e TAYLOR [30] apresentam a solução numérica do problema
obtida por diversas técnicas distintas através da aproximação de uma malha de 200
elementos triangulares, cuja geometria se sobrepõe à obtida pela divisão de cada
elemento quadrangular da figura (6.16) em 2 elementos triangulares a partir de sua
diagonal.
Os resultados aqui obtidos foram semelhantes aos obtidos por ZIENKIEWICKZ
e TAYLOR [30]. Na figura (6.17), a solução numérica da pressão em uma seção de
corte BB (fig. 6.16) é apresentada. Foram avaliados distintos parâmetros de calibração
para os modelos GLS e Bochev. Uma solução obtida por ZIENKIEWICKZ e TAYLOR
[30] para elementos triangulares é afixada no gráfico para fins comparativos.
Figura 6. 17 – Variação da pressão na seção de corte BB
De acordo com o gráfico apresentado na figura (6.17), observa-se que o cálculo
da pressão é extremente sensível e dependente dos parâmetros de penalização utilizados
nos modelos, principalmente pelo fato de se considerar apenas velocidades prescritas.
Para melhor compreensão do fenômeno, será apresentado o comportamento da
distribuição espacial da pressão para o modelo GLS em função do parâmetro de
estabilização α (figuras 6.18, 6.19 e 6.20).
93
Figura 6. 18 – Distribuição espacial da pressão pelo método GLS para α=0.02
Figura 6. 19 – Distribuição espacial da pressão pelo método GLS para α=1.0
Figura 6. 20 – Distribuição espacial da pressão pelo método GLS para α=20
Observa-se que para valores muito baixos de α a pressão permanece estável (fig.
6.18), mas apresenta gradientes acentuados e altos valores numéricos nas regiões
94
próximas à extremidade superior do domínio. À medida que o parâmetro α vai
assumindo maiores valores os gradientes nas extremidades vão diminuindo (fig. 6.19).
Entretanto, para valores muito elevados de α observa-se uma forte dissipação nos
valores da pressão (fig. 6.20). Especificamente nesse problema ZIENKIEWICKZ e
TAYLOR [30] sugerem um valor para o parâmetro de estabilização GLS compreendido
entre 15.0 ≤≤ α .
Nota-se na figura (6.17) que o comportamento do parâmetro α=0.5 no modelo
GLS foi similar ao modelo Bochev com α=2. A figura (6.21) mostra a distribuição
espacial da pressão para o modelo do Bochev.
Figura 6. 21 – Distribuição espacial da pressão pelo método Bochev para α=2
Esse valor condiz com os resultados analisados por DORHMANN e BOCHEV
[46], que sugerem o uso do parâmetro estabilizador para esse modelo de projeção de
pressão (Bochev) entre 2.1 ≤≤ α .
De maneira geral, o campo de velocidade nos problemas incompressíveis é bem
mais comportado quando comparado ao de pressão. Na figura (6.22) é apresentado o
perfil de velocidade vertical (v) na seção de corte BB. A distribuição espacial do campo
de velocidade obtido através do método GLS é apresentada na figura (6.23). Já o perfil
de velocidade horizontal (u) é representado na figura (6.24), tomado como base a seção
de corte AA. Em seqüência, o campo de velocidade horizontal é mostrado na figura
(6.25).
Convém ressaltar que para elevados valores do parâmetro de estabilização α
observa-se um comportamento anômalo do campo de velocidade (fig. 6.23d e 6.25d).
95
-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5-0.2
-0.15
-0.1
-0.05
0
0.05
0.1
0.15
0.2
Coordenada x
Vel
ocid
ade
verti
cal (
v)
Velocidade no corte BB
GLS 0.1GLS 0.5GLS 1.0Bochev 1.Bochev 2.
Figura 6. 22 – Variação da velocidade vertical na seção de corte BB
Figura 6. 23 – Distribuição espacial da velocidade vertical obtido com o modelo GLS
para: (a) α=0.1; (b) α=0.5; (c) α=1.0; (d) α=20.0;
96
-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
Coordenada x
Vel
ocid
ade
horiz
onta
l (u)
Velocidade no corte AA
GLS 0.1GLS 0.5GLS 1.0Bochev 1.Bochev 2.
Figura 6. 24 – Variação da velocidade horizontal na seção de corte BB
Figura 6. 25 – Distribuição espacial da velocidade horizontal obtido com o modelo GLS
para: (a) α=0.1; (b) α=0.5; (c) α=1.0; (d) α=20.0;
97
Observa-se que as oscilações no campo de velocidade em função do parâmetro
estabilizador são mínimas. No entanto, a combinação dos efeitos oscilatórios na parcela
da velocidade vertical e/ou horizontal pode comprometer a velocidade resultante
(módulo da velocidade) das partículas fluidas e, conseqüentemente, a direção de fluxo.
Objetivando efetuar essa análise, são apresentadas as componentes gráficas da figura
(6.26), que mostra a variação do módulo de velocidade em função do parâmetro
estabilizador para o modelo GLS.
Figura 6. 26 – Distribuição espacial do módulo de velocidade obtido com o modelo
GLS para: (a) α=0.1; (b) α=0.5; (c) α=1.0; (d) α=20.0;
Conforme já esperado, as interferências das oscilações no módulo da velocidade
das partículas fluidas são praticamente inexistentes no intervalo 15.0 ≤≤ α para o
modelo GLS. De maneira similar ao campo de pressão, o campo do módulo de
velocidade para o modelo Bochev é idêntico ao campo de velocidade do modelo GLS
com α=1.0, conforme previamente apresentado na figura (6.26c).
98
6.3 Exemplo 3: Viga viscoelástica submetida à ação de momento de flexão
Para analisar o efeito temporal da viscoelasticidade, reapresenta-se o problema
discutido no item 6.1. A análise foi realizada com a malha de 128 elementos, ilustrada
na figura seguinte:
Figura 6. 27 – Discretização espacial em 128 elementos finitos.
Para a análise, adota-se um Poisson ν = 0,3 e módulo de elasticidade E = 1000.
Os parâmetros viscoelásticos correspondem a μ1 = 0,99 e τ = 1, admitindo-se que a viga
apresenta comportamento típico de um sólido viscoelástico linear padrão. Adota-se nas
análises um intervalo de tempo Δt = 0,1.
Nota-se que a função de relaxação para o problema descrito corresponde à
expressão (6.3):
[ ]1000( ) 0,01 0,99exp( )2,6
G t t= + − (6.3)
Conforme discutido no capítulo 2, o rearranjo molecular advindo dos efeitos da
viscoelasticidade pode induzir a situação de incompressibilidade. Esse efeito pode ser
avaliado através de um coeficiente de Poisson efetivo, que pode ser calculado pela razão
entre o módulo de elasticidade volumétrica K e o módulo de cisalhamento ( )G t ,
resultando na seguinte expressão:
3 2 ( )( )6 2 ( )ef
K G ttK G t
ν −=
+
(6.4)
No gráfico (6.28), é apresentada a evolução temporal do deslocamento vertical
na extremidade superior da face direita (ver ponto de análise na figura 6.1). Os
resultados obtidos através de diferentes parâmetros de estabilização para o modelo
viscoelástico através dos métodos GLS e Bochev são comparados aos resultados obtidos
através do software FEAP®. Convém ressaltar que a análise no FEAP® foi realizada a
partir da formulação baseada apenas no método dos deslocamentos, principalmente para
avaliar o efeito do “travamento” característico dessa formulação para elevados valores
de Poisson.
99
0 5 10 15 20 25 30 35 40 45 500
10
20
30
40
50
60
70
80
90
Tempo
Des
loca
men
to v
ertic
al -
v
Feap - desl.Bochev a = 1Bochev a = 2GLS a = 0.1GLS a = 1.0
Figura 6. 28 – Evolução temporal do deslocamento
A partir da equação (6.4), pode-se observar que o Poisson efetivo para ∞→t
corresponde a νef = 0.49769585253456. Um valor significativo do coeficiente de
Poisson efetivo é alcançado em t = 10, correspondendo a νef (10) = 0.49768551230893.
A partir dos instantes seguintes, nota-se um crescente aumento do erro absoluto entre a
formulação MEF baseada no método dos deslocamentos, obtidos com o FEAP®, e as
aproximações dos demais modelos.
Nos resultados apresentados na tabela (1.3), nota-se que essa formulação do
FEAP® apresenta erro numérico de 4,05% para um ν = 0,49, sendo que os erros
crescem consideravelmente quando 5,0→ν . Entretanto, para valores mais baixos de ν
os resultados são aceitáveis. Dessa forma, pode-se concluir que os resultados obtidos
com os modelos GLS e Bochev foram condizentes com os resultados obtidos com o
FEAP® dentro da margem de validade da formulação nele utilizada para análise. Por
outro lado, para tempos elevados, onde o Poisson efetivo se aproxima de 0,5,
aparentemente a formulação em deslocamento (relativa ao FEAP®) sofreu o efeito de
“travamento”.
100
6.4 Exemplo 4: Membrana de Cook
Nesse exemplo, apresenta-se um problema em estado plano de deformação
tradicionalmente utilizado para avaliação de resultados das formulações mistas. O
problema de Cook foi introduzido por SIMO e RAFAI [52] e consiste em uma placa
engastada na extremidade esquerda e submetida à ação de forças tangenciais na face da
direita. Nas demais interfaces a placa é livre e não há aplicação de forças. Tal
configuração torna-se interessante para análise, uma vez que se observa o aparecimento
de esforços normais e cortantes. A geometria e condições de contorno podem ser
visualizadas na figura (6.29). Adota-se um módulo de elasticidade E = 250. Observa-se
que os elementos da malha apresentam um certo grau de distorção, o que pode
complicar o desempenho da formulação GLS.
Figura 6. 29 – Geometria e discretização espacial em 100 elementos (10 x 10)
Primeiramente, é feita uma análise da estabilidade e convergência dos modelos
lineares frente a situações de quase-incompressibilidade. Os estudos foram realizados
com a região discretizada em 100 elementos (10 x 10), 400 elementos (20 x 20), 900
elementos (30 x 30) e 2500 elementos (50 x 50).
Os resultados foram comparados aos do programa FEAP®. Nas figuras (6.30) e
(6.31) são apresentados a variação do deslocamento vertical no ponto de referência
(extremidade superior da face da direita) em função do coeficiente de Poisson, para as
malhas de 100 (10x10) e 2500 (50x50) elementos respectivamente.
101
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.50
1
2
3
4
5
6
7
8
9
10
Coeficiente de Poisson
Des
loca
men
to v
ertic
al -
v
GLS a = 0.05GLS a = 0.0Feap - mist.Feap - desl.Bochev a = 1Bochev a = 2
Figura 6. 30 – Deslocamento vertical para a malha 10 x 10 (100 elementos).
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.51
2
3
4
5
6
7
8
9
10
Coeficiente de Poisson
Des
loca
men
to v
ertic
al -
v
GLS a = 0.05Feap - mist.Feap - desl.Bochev a = 1Bochev a = 2
Figura 6. 31 – Deslocamento vertical para a malha 50 x 50 (2500 elementos).
Através das figuras (6.30) e (6.31) nota-se claramente que à medida que a malha
atinge um certo grau de refinamento, a aproximação dos modelos implementados (GLS
e Bochev) são equivalentes à aproximação obtida com a formulação mista do FEAP®.
102
Independentemente do grau de refinamento da malha, a formulação baseada apenas na
sentença variacional do deslocamento (u) apresenta “travamento” para elevados valores
do coeficiente de Poisson. Apenas para efeito ilustrativo, na figura seguinte pode-se
observar a distribuição de elementos para a malha (30 x 30) e (50x50):
Figura 6. 32 –(a) Malha de 900 elementos (b) malha de 2500 elementos
Por outro lado, independentemente do grau de refinamento da malha, a
distribuição de pressão mostrou-se bastante sensível em função dos parâmetros de
estabilização. No modelo GLS, para α = 0 se observou novamente o modo pressão
espúria. Já para α = 0,5 as pressões foram fortemente amortecidas, gerando também
oscilações numéricas, conforme apresentado na figura (6.33).
(a) (b)
Figura 6. 33 – Distribuição da pressão para (a) α = 0 (b) α = 0,5.
103
Os melhores resultados relativos a estabilização da pressão foram obtidos com
valores α<0,1. Já a formulação IGLS novamente obteve grande êxito na aproximação,
apresentando ganho de precisão numérica quando comparado a α = 0,1, conforme
apresentado na figura (6.34).
(a) (b)
Figura 6. 34 – Estabilização da pressão para: (a) α = 0,1 ; (b) IGLS
Admite-se para as análises seguintes que o material em estudo apresenta um
comportamento reológico típico de um sólido viscoelástico padrão (linear), com os
seguintes parâmetros: ν = 0,49; μ1 = 0,5 ; τ = 1. A função de relaxação assume então a
forma:
[ ]250( ) 0,5 0,5exp( )2,98
G t t= + − (6. 5)
Através da relação (6.4) pode-se observar que o material apresenta um
coeficiente de Poisson efetivo de 0,495 para ∞→t .
A evolução temporal do deslocamento é apresentada na figura seguinte. Os
resultados foram obtidos com o uso da malha de 900 elementos (30x30).
104
0 2 4 6 8 10 12 14 16 18 207
8
9
10
11
12
13
14
15
16
Tempo
Des
loca
men
to v
ertic
al -
v
GLS a = 0.05Feap - mist.Feap - desl.Bochev a = 1Bochev a = 2
Figura 6. 35 – Evolução temporal do deslocamento vertical (malha 30 x 30)
A diferença percentual máxima entre os resultados obtidos com a formulação
mista do FEAP® foi de 2,8% para o GLS, de 1,6% para o Bochev α = 1 e de 1,7 %
para o Bochev α = 2. Nota-se que a formulação baseada na sentença variacional do
deslocamento (Feap – desl.) apresentou valores inferiores às demais formulações para
t >1,5, aparentemente conseqüente do efeito de “travamento” para elevados valores do
coeficiente de Poisson.
A diferença percentual entre os resultados obtidos através da formulação mista
do FEAP® e dos modelos desenvolvidos tende a diminuir em função do grau de
refinamento da malha. Na figura seguinte, é apresenta a evolução temporal do
deslocamento obtida com uma malha de 2500 elementos (50x50). Os resultados foram
obtidos utilizando um coeficiente de Poisson inicial de ν=0,499.
105
0 2 4 6 8 10 12 14 16 18 207
8
9
10
11
12
13
14
15
16
Tempo
Des
loca
men
to v
ertic
al -
v
GLS a = 0.05Feap - mist.Bochev a = 1Bochev a = 2
Figura 6. 36 – Evolução temporal do deslocamento vertical (malha 50 x 50)
Nessa análise, a diferença percentual máxima entre a formulação mista do FEAP
e dos modelos implementados corresponde a 1,1% para o GLS, 0,9% para o Bochev
α = 1, e 1,2% para o Bochev α = 2.
Nas figuras seguintes, são apresentados o deslocamento resultante (módulo do
deslocamento) nos instantes t = 0 (deformação elástica instantânea) e t = 10 (advinda
dos efeitos da viscoelasticidade). Nota-se que o efeito da viscoelasticidade incorporou
um incremento significativo do deslocamento resultante.
106
Figura 6. 37 – Deslocamento resultante para t = 0.
Figura 6. 38 – Deslocamento resultante para t = 20.
107
6.5 Exemplo 5: Análise de reversibilidade da fluência viscoelástica
Nesta aplicação, analisa-se o comportamento numérico em analogia ao
comportamento físico de um corpo viscoelástico posto em ensaio de fluência (carga
constante) e submetido à remoção da carga q em um instante t’>0.
Para tanto, idealiza-se inicialmente um corpo engastado submetido à ação de
força de tração originada por um carregamento q, cuja configuração geométrica e
discretização espacial é apresentada na figura (6.39), onde pode-se observar a
subdivisão do domínio em 64 elementos. Adota-se: a =b =2; q=5; t’=2,5;
Figura 6. 39 – Geometria e discretização da malha.
Primeiramente, considera-se um material com propriedades características de um
sólido viscoelástico padrão, adotado como incompressível (ν = 0,5) e com viscosidade
G = 3,333 ( equivalente a um E = 100). Adota-se como tempo de relaxação o valor
unitário (τ = 1) e μ1 = 0,6.
O comportamento temporal do deslocamento horizontal no ponto de análise (ver
figura 6.39) é apresentado na figura (6.40).
108
0 1 2 3 4 5 6 7 8 9 100
0.05
0.1
0.15
0.2
0.25
0.3
0.35
Tempo
Des
loca
men
to h
oriz
onta
l - u
Bochev a = 2GLS a = 0.5GLS a = 0.1
Figura 6. 40 – Deslocamento horizontal para o modelo sólido viscoelástico padrão.
Observa-se no instante inicial a resposta instantânea correspondente ao
deslocamento de um modelo elástico linear. Nos instantes seguintes de aplicação da
carga constante q, o deslocamento aumenta continuamente, apresentando tendência a
assumir um valor assintótico característico dos exemplos apresentados anteriormente.
No instante t’ = 2,6, após a remoção da carga, há uma recuperação imediata equivalente
a resposta instantânea inicial. Em seguida, observa-se a recuperação da parcela viscosa,
que deve ser integralmente recuperada em função do arranjo do modelo.
A interpretação do comportamento ilustrado na figura (6.40) é facilmente
compreendida se o analisarmos a partir do modelo físico apresentado na figura (2.16).
Inicialmente, a resposta instantânea deve corresponder ao deslocamento advindo da
parcela elástica do modelo, que corresponde a soma dos elementos elásticos que estão
ligados em paralelo. Após a remoção da carga, o elemento elástico que está
paralelamente ligado ao elemento de Maxwell forçará a recuperação total do
deslocamento, condizendo com o comportamento obtido pelo MEF mostrado na figura
(6.40).
109
Pela própria definição do ensaio de fluência, a pressão deve ser mantida
constante no elemento ao longo do teste. Esse comportamento é compatível com os
resultados obtidos, onde nos instantes de aplicação da carga t<2,6 a pressão é mantida
constante e proporcional ao valor de carregamento aplicado. Após a retirada da carga, é
esperado que pressão seja nula. Dessa forma, os resultados numéricos apresentados na
figura (6.41) estão condizentes com o comportamento físico do ensaio de fluência.
0 1 2 3 4 5 6 7 8 9 10-1
0
1
2
3
4
5
6
Tempo
Pre
ssão
Bochev a = 2GLS a = 0.5
Figura 6. 41 – Pressão para o modelo sólido viscoelástico padrão
Um fato importante a destacar está relacionado aos resultados obtidos através do
modelo GLS com α = 0, que equivale ao uso somente da formulação mista de Galerkin.
Conforme discutido no capítulo 3, para elevados valores do coeficiente de Poisson essa
formulação está sujeita ao modo “pressão espúria”. Tem-se aqui novamente um caso
típico, onde os deslocamentos são razoavelmente aproximados (fig 6.40) mas a pressão
está totalmente comprometida, conforme ilustrado em (6.42).
110
00.5
11.5
2
0
0.5
1
1.5
2
-100
-50
0
50
100
Eixo X
Pressão (p)
Eixo Y 0 1 2 3 4 5 6 7 8 9 10-50
0
50
100
150
200
250
300
350
400
Tempo
Pre
ssão
GLS a = 0.
Figura 6. 42 – (a) Distribuição espacial da pressão para t = 0; (b) Evolução temporal da
pressão
Dessa forma, as formulações implementadas conseguiram efetivamente remover
o efeito do modo “pressão espúria” e o efeito de “travamento”, garantindo estabilidade
para o modelo viscoelástico. A pressão estabilizada é apresentada na figura (6.43),
obtida pelo modelo Bochev.
0
0.5
1
1.5
2
0
0.5
1
1.5
2
0
5
10
15
20
Eixo X
Pressão (p)
Eixo Y
Pressao
4 5 6 7 8 9 10 Figura 6. 43 – Distribuição espacial da Pressão (estabilizada)
Para avaliar os efeitos temporais da viscoelasticidade, consideraram-se também
análises com diferentes valores do coeficiente de Poisson, obtendo-se resultados
similares aos apresentados. Entretanto, para a configuração em estudo, nota-se
claramente a dependência dos deslocamentos com o valor do coeficiente de Poisson. Na
figura seguinte, os resultados em t = 2,4 para coeficientes de Poisson 0,5 e 0,2 são
apresentados. Destaca-se o comportamento próximo às extremidades do engaste, onde
observa-se a conservação de volume nos elementos da malha para ν =0,5.
111
Figura 6. 44 – Deslocamento resultante para t = 2,4 para: (a) ν = 0,5; (b) ν = 0,3. (escala
alterada para visualização – 2x)
Considera-se agora o material com propriedades de um fluido viscoelástico de
Maxwell. Nesse caso, têm-se μ1 = 1,0 e demais propriedades idênticas às utilizadas na
análise anterior (ν = 0,5).
O comportamento temporal do deslocamento horizontal no ponto de análise é
apresentado na figura (6.45). Conforme esperado, nota-se uma resposta elástica inicial
seguida por um aumento do deslocamento no tempo, que de acordo com as análises
teóricas apresentadas no capítulo 2 crescerá indefinidamente, cuja tendência da taxa de
deformação (coeficiente angular) se aproximará a de um fluido newtoniano. Nota-se que
a partir do instante de retirada da carga (t>2,6), há recuperação imediata da parcela
elástica, sendo a parcela viscosa irrecuperável, conforme discutido previamente.
Observa-se uma pequena oscilação numérica no instante imediatamente posterior a
retirada da carga, originada pela descontinuidade do carregamento. Convém relembrar
que a taxa de deformação é aproximada através de um esquema numérico em diferença
finita progressiva.
112
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 50
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
Tempo
Des
loca
men
to h
oriz
onta
l - u
Bochev a = 2GLS a = 0.5GLS a = 0.1
Figura 6. 45 – Deslocamento horizontal para o modelo de Maxwell
O comportamento da pressão aproximada é mostrado no gráfico seguinte, e
conforme discussão prévia está de acordo com a interpretação física do problema.
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5-1
0
1
2
3
4
5
6
Tempo
Pre
ssao
Bochev a = 2GLS a = 0.5
Figura 6. 46 – Pressão- Modelo de Maxwell
113
Agora, faz-se uma análise de um corpo submetido a um estado de cisalhamento
puro. Nesse caso, os efeitos da viscoelasticidade estão presentes apenas na parcela
desviatória, uma vez que a pressão deve ser nula. A distribuição da carga q é
apresentada na figura seguinte:
Figura 6. 47 – Corpo submetido a estado de cisalhamento puro
Novamente, para elevados valores de coeficiente de Poisson, o uso da
formulação mista de Galerkin revelou o modo “pressão espúria”, conforme apresentado
na figura (6.48).
Pressao
-3000 -2000 -1000 0 1000 2000 0
0.511.52
0
0.5
1
1.5
2
-4000
-2000
0
2000
4000 Pressão (p)
Eixo X
Eixo Y
Figura 6. 48 – Modo pressão espúria no estado de cisalhamento puro.
114
Entretanto, os valores aproximados para o deslocamento foram novamente
razoáveis, evitando dessa forma o efeito do “travamento”. A configuração do
deslocamento resultante é apresentada na figura seguinte:
0 0.5 1 1.5 20
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2Deslocamento resultante
Figura 6. 49 – Deslocamento resultante para o GLS α = 0.
Os resultados obtidos para o modelo de Maxwell com ν = 0,5 são apresentados
na figura (6.50), onde se observa um comportamento condizente com os resultados
teóricos previamente discutidos.
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 50
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
Tempo
Des
loca
men
to h
oriz
onta
l - u
Bochev a = 2GLS a = 0.5GLS a = 0.1
Figura 6. 50 – Evolução temporal do deslocamento vertical (Modelo de Maxwell).
115
De maneira esperada, a pressão aproximada toma valores nulos ao longo de toda
análise, conforme apresentado na figura (6.51).
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 50
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1x 10
-8
Tempo
Pre
ssao
Bochev a = 2GLS a = 0.5
Figura 6. 51 – Evolução temporal da pressão (Modelo de Maxwell).
Em outra análise, adotou-se o material com comportamento típico de um sólido
viscoelástico padrão, com μ1=0,6. Os resultados obtidos para o deslocamento e pressão
foram novamente condizentes com o comportamento físico esperado, conforme
apresentado na figura (6.52).
0 1 2 3 4 5 6 7 8 9 100
0.2
0.4
0.6
0.8
1
1.2
1.4
Tempo
Des
loca
men
to h
oriz
onta
l - u
Bochev a = 2GLS a = 0.5GLS a = 0.1
0 1 2 3 4 5 6 7 8 9 100
1
2
3
4
5
6
7
8
9
10x 10-9
Tempo
Pre
ssao
Bochev a = 2GLS a = 0.5
Figura 6. 52 – Evolução temporal do deslocamento vertical e pressão (Modelo do sólido
viscoelástico).
116
Ambas as análises foram realizadas com diversos valores do coeficiente de
Poisson. Para ilustrar a influência desse parâmetro nas análises, apresenta-se na figura
seguinte o deslocamento resultante (módulo do deslocamento) para o modelo de
Maxwell realizado com Poissons ν = 0,2 e ν = 0,5 em t = 100.
Figura 6. 53 – Deslocamento resultando para o modelo de Maxwell em t = 100:
(a) ν = 0,2; (b) ν = 0,5.
a b
117
7. Conclusões e recomendações
Com base nos resultados apresentados, observa-se que a aplicação da
formulação mista em elementos finitos é altamente eficaz no que diz respeito à
eliminação do fenômeno de “travamento” de malha, que ocorre com o emprego da
formulação clássica de Galerkin em situações próximas à incompressibilidade.
Entretanto, é necessário o uso de formulações estabilizadas para a garantia da
aproximação do campo de pressão com acurácia, evitando dessa forma as instabilidades
numéricas advindas do chamado “modo pressão espúria”. Assim, pretende-se enfatizar
a idéia de que nem sempre os resultados com bom grau de aproximação do campo de
deslocamento/velocidade induzirá a boas aproximações do campo de pressão (tensão e
deformação). Por essa razão, os métodos estabilizados são altamente recomendados no
emprego das formulações mistas em elementos finitos.
O método Galerkin Least Square (GLS) [43] mostrou-se eficiente na remoção
do fenômeno de “travamento” de malha e do modo “pressão espúria”. Entretanto, o fato
de a formulação depender de um parâmetro de calibração α, que detém um valor ótimo
para cada classe de problemas, passa a dificultar o uso do modelo. Nos exemplos
analisados, pode-se observar que em valores de α próximos a zero a pressão pode não se
comportar de maneira estabilizada, visto que a parcela relacionada à ponderação por
mínimos quadrados tende a ser anulada, predominando dessa forma a parcela de
ponderação de Galerkin. Entretanto, para valores elevados de α nota-se um
amortecimento artificial na aproximação da pressão, que pode dessa forma apresentar
valores não compatíveis com a solução exata.
Já a formulação linear de estabilização direta de pressão Bochev [46] mostrou
ser menos sensível à variação do parâmetro de estabilização, pois as modificações do
parâmetro α alteravam de maneira pouco significativa os resultados aproximados para a
pressão e deslocamentos/velocidades.
É notório, de acordo com os resultados apresentados no exemplo 4, que em
função do grau de refinamento da malha, o parâmetro de estabilização passa a não ter
influência acentuada na aproximação dos resultados.
Observou-se que a sensibilidade do parâmetro de estabilização é mais crítica em
situações onde todo o contorno do domínio apresenta a condição de Dirichlet
(velocidades/deslocamento prescritos). Nesse caso, é extremamente necessária a
118
prescrição do valor da pressão em algum ponto do domínio, conforme discutido no
capítulo 3.
No tocante ao comportamento dos modelos lineares em situações em que a
malha não seja estruturada e contenha elementos com alto grau de distorção, foi
perceptível uma perda na precisão da aproximação numérica (exemplo 1). Esse efeito é
mais evidente no modelo GLS [43], onde para os exemplos numéricos analisados foi
observado que quanto maior é a parcela atribuída à sentença de ponderação em mínimos
quadrados, maiores são os desvios numéricos. Aparentemente isso ocorre ao se fazer
atribuições ao parâmetro de ponderação α maiores do que o valor necessário,
provocando dissipações numéricas. Em função da capacidade de determinar o valor de α
para cada elemento, o método IGLS proposto obteve grande ganho na precisão do
campo de deslocamento/velocidade, garantindo também a estabilidade do campo de
pressão. Trata-se de uma estratégia interessante, visto que ainda hoje a formulação GLS
é alvo de enumeras discussões em função da necessidade de determinação do parâmetro
de calibração de forma automática.
Ambos os modelos lineares apresentaram significativas taxas de convergência
entretanto, conforme discutido no exemplo 1, a formulação mista não tem boa
aproximação para malhas com o grau de refinamento ínfimo (pobre). Observou-se em
todos os exemplos consistência, convergência e conseqüentemente estabilidade
numérica das formulações.
Em termos de tempo de processamento, a formulação proposta por Bochev
mostrou-se computacionalmente mais eficiente, uma vez que não são necessários
cálculos de derivadas de ordens superiores, tal qual ocorre na formulação GLS. Nos
exemplos executados, o tempo de processamento na formulação GLS foi em média 9%
mais lento quando comparado à formulação do Bochev. Entretanto, para tratamento de
problemas que apresentem distorções na malha, observou-se numericamente que os
resultados obtidos através da estratégia proposta (IGLS) foram mais precisos que os
advindos da formulação do Bochev. Por essa razão, para se o obter o mesmo grau de
precisão da formulação IGLS seria necessário o refinamento de malha para o modelo do
Bochev, o que aumentaria o custo de processamento numérico.
No que diz respeito aos resultados advindos do modelo viscoelástico, pode-se
destacar a influência do grau de aproximação da pressão na evolução temporal dos
descolamentos\velocidades, uma vez que é necessária a aproximação da pressão em
cada etapa do incremento temporal para o cálculo da tensão total e conseqüente
119
atualização do vetor de forças internas resultante. Como discutido no capítulo 2, a
reacomodação molecular oriunda da combinação dos efeitos viscosos e elásticos do
material viscoelástico pode induzir a um comportamento de incompressibilidade. Assim
sendo, a formulação do modelo linear deve responder de maneira eficiente à remoção
dos efeitos de “travamento de malha” e do modo “pressão espúria”, de maneira similar
aos resultados dos exemplos (6.3) e (6.4). Com base nas análises realizadas, é possível
afirmar que a consistência e convergência numérica dos modelos lineares (Bochev, GLS
e IGLS) garantem estabilidade para o modelo viscoelástico incompressível.
Com base nas análises comparativas realizadas através do software FEAP®,
observou-se que o modelo viscoelástico linear se comportou de maneira compatível às
respostas obtidas com o FEAP®, o que de certa forma valida os resultados numéricos
obtidos.
No caso da última análise, realizada a partir da verificação da reversibilidade da
fluência, constatou-se que os resultados obtidos são compatíveis com o comportamento
físico previsto.
Em função dos incrementos temporais dos deslocamentos (ou velocidades) e
pressão, observa-se em alguns casos, principalmente no modelo de Maxwell, a presença
de grandes deslocamentos resultantes. Dessa forma, sugere-se incorporar estratégias
para o tratamento dos efeitos da não linearidade geométrica em trabalhos futuros.
Pretende-se em trabalhos futuros incorporar estratégias para a representação de
fenômenos inerentes à poromecânica e plasticidade, o que após a integração de módulos
com as estratégias já desenvolvidas permitirá o acoplamento de modelos
poroviscoelásticos e poro-viscoelastoplásticos.
Em todo caso, os algoritmos desenvolvidos nesse trabalho apresentaram
robustez e consistência numérica, podendo ser aplicados em análises reológicas
envolvendo a elasticidade linear compressível e incompressível e na modelagem de
escoamento de fluidos em regime permanente (equação de Stokes). Permitem também
aplicações na modelagem de sólidos viscoelásticos (modelo do sólido viscoelástico
linear padrão) e fluidos viscoelásticos (modelo de Maxwell).
120
8. Referências bibliográficas [1] DORAISWAMY, D. “The origens of rheology: A short historical excursion”,
Rheology Bulletin, v. 71, Jan. 2002.
[2] MACKERLE, J. “Finite element/boundary element analysis of viscoeleastic and
viscopleastic problems - A bibliography (1994-1996)”, Finite elements in
analysis and design, v. 27, n. 3, pp 273-287, Nov 1997.
[3] DOZDROV, A. D. Finite eleasticity and viscoeleasticity: a course in the nonlinear
mechanics of solids, Massachusettes, Word scientific publishing Co. Pte.
Ltd, 1996.
[4] FLÜGGE, W. Viscoeleasticity. Stanford, Blaisdell Publishing Company, 1967.
[5] BLAND, D. R. Theory of linear viscoeleasticity, Oxford, Pergamon Press, 1960.
[6] CRISTENSEN, R. M. Theory of viscoeleasticity, 2 ed. New York, Dover
Publication Inc., 2003.
[7] ZIENKIEWICZ, O. C; WATTSON, M.; KING, I.P. “A numerical method of visco-
eleastic stress analysis”. Int. J. Mech. Sci., vol. 10, pp. 807-827, 1968.
[8] SRINATHE, H. R; LEWIS, R.W. “A finite element method for thermoviscoeleastic
analysis of plane problems”. Comput. Meth. Appl. Mech. Eng, vol. 25, pp.
21-23, 1981.
[9] ALFREY, T. “Methods of representation the properties of viscoeleastic materials”.
Quarterly of applied mathematics, vol. 3, pp 143-150, 1945.
[10] ADEY, R. A.; BREBBIA, C. A. “Efficient method for solution of viscoeleastic
problems”. J. Engng Mech, vol. 99-6, pp. 119-1127, 1973.
[11] CARPENTER, W. C. “Viscoeleastic stress analysis”. Int. J. Num. Meth. Ingin., vol.
4, pp. 357-366, 1972.
[12] WHITE, J. L. “Finite Elements in linear viscoeleasticity”. Proc. Conf. Struct.
Mech. Dayton, Ohio, 1958.
[13] WEBBER, J. P. H. “Stress analysis in viscoeleastic bodies using finite elements
and a correspondece rule with eleasticity”. J. Strain Analysis, vol. 4, pp.
236-243, 1969
[14] YADAGIRI, S. REDDY, C.P. “Viscoeleastic analysis of nearly incompressible
solids”. Comp. & Struc., vol. 20, n. 5, pp. 817-825, 1985.
[15] ROYLANCE, D. “Engineering Viscoeleasticity. Overview of nonlinear
viscoeleastic theory, see for instance W.N. Findley et al., Creep and
121
Relaxation of Nonlinear Viscoeleastic Materials”. Dover Publications, New
York, 1989. Disponível em
http://web.mit.edu/course/3/3.11/www/modules/visco.pdf. Acessado em
02/03/2007.
[16] LEAL, C.R. Breve introdução à reologia. Instituto superior de engenharia de
Lisboa. Disponível em: http://sites.isel.ipl.pt/fisica/lr_f/main.html. Acessado
em 10 de novembro de 2007.
[17] WINKIPEDIA. Definição: entropia. Disponível em:
http://pt.wikipedia.org/wiki/Entropia. Acessado em 12 de novembro de
2007.
[18] EIRICH, F. R. Rheology : Theory and Applications, Vol. 1, New York, Academic
Press, 1956.
[19] NAGAYAMA, K.; YANAGIHARA, S; MATSUMOTO, T. “Actin Filaments
Affect on Not Only Eleasticity But Also Late Viscous Response in Stress
Relaxation of Single Isolated Aortic Smooth Muscle Cells”. Journal of
Biomechanical Science and Engineering. vol 2, n.3, pp.93-104, 2007.
[20] MUNIZAGA, G.T; CÁNOVAS, G.V.B. “Rheology for the food industry”.
Journal of Food Engineering , vol. 67, Issues 1-2, pp. 147-156, 2005.
[21] MAKRIS, N. “Complex-parameter kelvin model for eleastic foundations”.
Earthquake Engineering & Structural Dynamics, vol. 23, Issue 3 , pp. 251 –
264, 2007.
[22] KUO, M.I; WANG Y.C; GUNASEKARAN, S. “A Viscoeleasticity Index for
Cheese Meltability Evaluation”. J Dairy Sci, vol. 83, pp. 412–417., 2000.
[23] VALE e MELO. “Análise do comportamento de pavimentos rodoviários flexíveis:
os modelos materiais das misturas betuminosas”. In: Anais da V jornadas
luso-brasileiras de pavimentos: políticas e tecnologias. Recife, Brasil, julho
de 2006.
[24] SIMO, J.C; HUGHES, T.J.R. Computational ineleasticity: Interdisciplinary
Applied Mathematics. Berlin, Springer-Verlag, 1998.
[25] FINDLEY, W.N. Creep and relaxation of nonlinear viscoeleastic materials;with
an introduction to linear viscoeleasticity. Amsterdam, North-Holland, 1976.
[26] VINOGRADOV, G.V; MALKIN, A.Y. Rheology of polymer, Moscow, Mir
Publisher, 1980.
122
[27] BERGER, J.T. Viscoeleasticity: Phenomenological aspects, New York, Academic
Press Inc, 1906.
[28] LOGAN, D.L. A first course in the finite element method. Boston, PWS Publishing
Company, 1992.
[29] HUEBNER, K.H; THORNTON, E.A. The finite element methods for engineers.
New York, A wiley-intersciense publication. 1942.
[30] ZIENKIEWICZ, O.C.; TAYLOR, R.L. Finite Element Method: Its Basis and
Fundamentals, ed.6. Londres, Published Elsevier, 2005.
[31] MANSUR, W.J. Métodos numéricos em engenharia. Notas de Aula – Elementos
finitos, Universidade Federal do Rio de Janeiro, 2006.
[32] BATHE , K. J. Finite Element Procedures in Engineering Analysis. New Jersey,
Prentice-Hall, Inc., , 1982.
[33] FOX R.W; MCDONALD, A.T. Introduction to Fluid Mechanics, ed. 4. Wiley &
Sons, 1992.
[34] PANNACHET, T; BOONPICHETVONG, M. “An application of hierarchical
finite elements for overcoming volumetric locking in plasticity”. The
Enginnering Institute of Thailand. Disponível em:
http://www.eit.or.th/article/data/01080002.pdf. Acessado em setembro de
2007.
[35] BITTENCOURT, E. Mecânica da fratura e do dano – Notas de aula, cap. 8.
Universidade Federal do Rio Grande do Sul. Porto Alegre, 2005. Disponível
em: http://www.cpgec.ufrgs.br/bittenco/. Acessado em outubro de 2007.
[36] SORIANO, H.L. Método de elementos finitos em análise de estruturas. São Paulo,
Editora da Universidade de São Paulo, 2002.
[37] ANSYS, Inc. Theory reference. Version 10. Ansys Company, 2000.
[38] GUERMOND, J.L; ERN, A. Theory and practice of finite elements. New York,
Springer-Verlag, 2004.
[39] ADAMS, R.A. Sobolev Spaces, Academic Press, New York, 1975.
[40] DE, S. On the development of an efficient truly meshless discretization procedure
in computational mechanics. D.Sc dissertation, Massachusetts Institute of
Technology, Massachusetts, 2001.
[41] CHAPELLE, D. BATHE, K.J. “The inf-sup test”. Computers & Structures, vol. 47,
pp. 537-545, 1993.
123
[42] FRANCA, L.P. Advances in stabilized methods in computational mechanics-
Preface. Comput. Methods Appl. Mech. Engrg, vol. 166, pp 1-2, 1998.
[43] HUGHES, T.J.R; FRANCA, L.P; HULBERT, G.M. “A new finite element
formulation for computational fluid dynamics. The Galerkin/least-squares
method for advective-difusive equations”. Comp. Meth. Appl. Mech. Eng.,
vol. 73, pp. 173-189, 1989.
[44] XIA, K. YAO, H. “A Galerkin/least-square finite element formulation for nearly
incompressible eleasticity/stokes flow”. Applied Mathematical Modelling,
vol. 31, pp. 513–529, 2007.
[45] BARTH, T.; BOCHEV, P.; SHADID, J.; GUNZBURGER, M. “A taxonomy of
consistently stabilized finite element methods for the Stokes problem”. SIAM
J. Sci. Comp., vol. 25, pp. 1585-1607, 2004.
[46] DOHRMANN, C.R; BOCHEV, P.B. “A stabilized finite element method for the
Stokes problem based on polynomial pressure projections”. International
Journal for Numerical Methods in Fluids, vol. 46/2 , pp. 183-201, 2004.
[47] CARMO, E.G.D; SANTOS, J.P.L; MANSUR, W.J. An improvement in the
Galerkin Least Square method (GLS) for Stokes and incompressible
elasticity problems. To appear.
[48] ZIENKIEWICZ, O.C.; TAYLOR, R.L. The Finite Element Method for Solid and
Structural Mechanics, ed.6. Londres, Published Elsevier, 2005.
[49] ABAQUS, Inc. Tutorials. Version 6.6, Abaqus documentation. Disponível em
http://lcadm.rrzn.uni-hannover.de:2080/v6.6/index.html. Acessado em
06/08/2007.
[50] TAYLOR, R.L. FEAP-A Finite Element Analysis Program. Version 8, Theory
Manual. Berkeley, University of California at Berkeley, 2007.
[51] REDDY, B.D; DJOKO, J.K. “An extended Hu–Washizu formulation for
eleasticity”. Comput. Methods Appl. Mech. Engrg, vol. 195, pp. 6330–6346,
2006.
[52] SIMO, J.C.; RIFAI, M.S. “A class of mixed assumed strain methods and the method
of incompatible modes”. Int. J. Num. Meth. Eng., 29, 1595-638, 1990.
