HENRIQUE CAMPELO GOMES

Mtodo dos Elementos Finitos com Fronteiras Imersas aplicado a Problemas de Dinmica dos Fluidos e Interao Fluido-Estrutura

So Paulo 2013

HENRIQUE CAMPELO GOMES

Mtodo dos Elementos Finitos com Fronteiras Imersas aplicado a Problemas de Dinmica dos Fluidos e Interao Fluido-Estrutura

Tese apresentada Escola Politcnica da Universidade de So Paulo para obteno do ttulo de Doutor em Engenharia

So Paulo 2013

HENRIQUE CAMPELO GOMES

Mtodo dos Elementos Finitos com Fronteiras Imersas aplicado a Problemas de Dinmica dos Fluidos e Interao Fluido-Estrutura

Tese apresentada Escola Politcnica da Universidade de So Paulo para obteno do ttulo de Doutor em Engenharia

rea de Concentrao: Engenharia Civil Orientador: Prof. Dr. Paulo de Mattos Pimenta

So Paulo 2013

Este exemplar foi revisado e corrigido em relao verso original, sob responsabilidade nica do autor e com anuncia de seu orientador.

So Paulo, de abril de 2013. Assinatura do autor __________________________________ Assinatura do orientador _____________________________

FICHA CATALOGRFICA

Gomes, Henrique Campelo

Mtodo dos elementos finitos com fronteiras imersas aplica- do a problemas de dinmica dos fluidos e interao fluido-estru-tura / H.C. Gomes. -- verso corr. -- So Paulo, 2013.

98 p.

Tese (Doutorado) - Escola Politcnica da Universidade de So Paulo. Departamento de Engenharia de Estruturas e Geotc-nica.

1. Mtodo dos elementos finitos 2. Interao fluido-estrutura 3. Dinmica dos fluidos 4. Fronteiras imersas I. Universidade de So Paulo. Escola Politcnica. Departamento de Engenharia de Estruturas e Geotcnica II. t.

DEDICATRIA

minha esposa Mnica

s minhas filhas Amanda e Helena

Aos meus Pais Jonas e Zeza

AGRADECIMENTOS

Ao professor, mestre e amigo Paulo Pimenta pela confiana, orientao e ensinamentos acadmicos e extra-acadmicos.

Ao Concelho Nacional de Desenvolvimento Cientfico e Tecnolgico - CNPq, pelo suporte financeiro.

Ao Servio Alemo de Intercmbio Acadmico DAAD, pelo suporte financeiro durante os seis meses de doutorado sanduche na Universidade Tcnica de Munique (TUM).

John Argyris Foundation pelo prmio, em forma de ajuda financeira, que permitiu a realizao de cursos, participao em congressos e visitas a grupos de pesquisa no exterior.

Ao professor Wolfgang Wall, pela cooperao e pelas poucas, porm produtivas, conversas durante a estada no Lehrstuhl fr Numerische Mechanik (LNM) da TUM.

Aos colegas do LNM, em especial, ao Ismail, Andreas, Key Mller, Peter, Volker e Axel Gerstenberger.

turma do JAC, pela amizade e ajuda. Agradecimentos especiais ao Jorge, Eduardo, Campello, Marcelo, Evandro, Fernando, Leo, Alexandre e Paulo.

Aos casais Flvio e Cris, Lus e Gleicy e Tarsis e Darlene, por todos os bons momentos vividos durante o perodo em So Paulo.

Aos meus pais pelo suporte emocional e financeiro, sem os quais o incio, meio e fim deste trabalho no seriam possveis.

minha esposa Mnica pelo companheirismo durante todas as fases do doutorado e por me dar os dois presentes mximos da minha vida: Amanda e Helena.

Aos meus irmos Jonas e Felipe pelas conversas e conselhos sempre muito sbrios e teis.

SUMRIO

Resumo ....................................................................................................................................... i

Abstract ..................................................................................................................................... ii

Objetivo .................................................................................................................................... iii

1. Introduo ............................................................................................................................ 1

2. Reviso Bibliogrfica ........................................................................................................... 6

3. O Problema de Fluidos ...................................................................................................... 10 3.1. Equaes de Navier-Stokes ......................................................................................... 10 3.2. Forma fraca .................................................................................................................. 12 3.3. Tenses e pseudotenses no contorno ......................................................................... 14 3.4. Integrao no tempo..................................................................................................... 15

Mtodo de Newmark .................................................................................................... 15 3.5. Discretizao do espao por elementos finitos ............................................................ 17 3.6. Estabilizao da condio LBB ................................................................................... 22

Problema 3.1: Escoamento de Stokes com soluo analtica ..................................... 22 Estabilizao da condio LBB usando Galerkin/Least-Squares ............................... 27

3.7. Estabilizao da conveco ......................................................................................... 32 3.8. O problema tangente .................................................................................................... 34 3.9. Exemplos Numricos ................................................................................................... 37

Problema 3.2: Escoamento numa cavidade ................................................................ 37 Problema 3.3: Escoamento em torno de um cilindro .................................................. 41

3.10. Concluso Parcial ........................................................................................................ 47

4. O Problema Estrutural ...................................................................................................... 48 4.1. Problema esttico ......................................................................................................... 48 4.2. Modelo constitutivo Neo-Hookiano de Ciarlet-Simo .................................................. 50 4.3. Forma fraca .................................................................................................................. 52 4.4. Discretizao do espao por elementos finitos ............................................................ 53 4.5. Equilbrio ..................................................................................................................... 55 4.6. O problema tangente .................................................................................................... 56 4.7. Problema dinmico ...................................................................................................... 57 4.8. Forma fraca .................................................................................................................. 58 4.9. Discretizao por elementos finitos ............................................................................. 59 4.10. Integrao no tempo..................................................................................................... 60 4.11. O problema tangente .................................................................................................... 61 4.12. O problema plano......................................................................................................... 62

5. O Problema de Interao Fluido-Estrutura .................................................................... 64 5.1. Introduo .................................................................................................................... 64 5.2. Imposio das condies na interface .......................................................................... 66 5.3. O problema fluido com interfaces mveis ................................................................... 66

Discretizao do espao por elementos finitos............................................................ 68 Discretizao da interface ........................................................................................... 69 Problema matricial ...................................................................................................... 70

5.4. O problema acoplado ................................................................................................... 71 Algoritmo do problema acoplado ................................................................................ 71

5.5. Simulaes numricas.................................................................................................. 71 Problema 5.1: escoamento estacionrio em torno de um cilindro (Re=20) ............... 72 Problema 5.2: escoamento em torno de um cilindro em movimento (Re=100) .......... 76 Problema 5.3: escoamento estacionrio em um canal com obstculo elstico .......... 78 Problema 5.4: escoamento em torno de um cilindro com barra flexvel atrs ........... 79

6. Discusso e Concluses ...................................................................................................... 84

7. Propostas de Trabalhos Futuros ...................................................................................... 86

8. Referncias Bibliogrficas ................................................................................................. 87

9. Apndice ............................................................................................................................. 97 9.1. Matriz convectiva: ....................................................................................................... 97 9.2. Matriz viscosa: ............................................................................................................. 97 9.3. Operador gradiente: ..................................................................................................... 97 9.4. Vetor das foras de volume e tenso aplicada: ............................................................ 98 9.5. Matriz viscosa (3.35) ................................................................................................... 98

i

Resumo

Este trabalho pode ser dividido em trs etapas principais. Inicialmente proposta uma

formulao estabilizada do mtodo dos elementos finitos (MEF) para soluo de problemas

de escoamento incompressvel governado pela equao de Navier-Stokes. Esta formulao foi

implementada em um cdigo computacional e testada atravs de diversos exemplos

numricos. Alguns elementos finitos com diferentes pares de funo de interpolao da

velocidade e presso, consagrados na literatura, e tambm elementos finitos menos populares,

foram investigados e seus resultados e performance comparados. A segunda etapa consiste na

formulao do problema estrutural. Buscou-se por uma formulao dinmica, no linear,

capaz de simular movimentos complexos de estruturas sujeitas a grandes deslocamentos e

grandes deformaes durante longos intervalos de tempo. A etapa final deste trabalho a

proposio de um mtodo para soluo de problemas de Interao Fluido Estrutura (IFE) que

utiliza o conceito de fronteiras imersas como alternativa a abordagens ALE (Arbitrary

Lagrangian Eulerian) clssicas. Elementos Finitos Generalizados, juntamente com

Multiplicadores de Lagrange, so utilizados para prover descontinuidade nos campos de

velocidade e presso do fluido ao longo da interface com a estrutura. O acoplamento dos dois

problemas realizado utilizando um mtodo implcito e alternado (staggered scheme), que

possui a vantagem de permitir, facilmente, a implementao de cdigos computacionais

desenvolvidos para resolver isoladamente o problema fluido e/ou estrutural.

ii

Abstract

This work is divided in three parts. Initially, it is presented a stabilized Finite Element Method

formulation to solve fluid flow problems governed by the incompressible Navier-Stokes

Equations. This formulation was implemented in a computer code and validated throughout

several numeric simulations. Some well-known finite elements with different pairs of

velocity/pressure aproximations, as well as some other less popular elements, were

investigated and their performance compared. The second part describes the Structural

Problem formulation. This formulation is able to simulate nonlinear dynamic problems

involving large displacements and finite strains during long period of time. In the final part of

this work, it is proposed a Fluid-Structure Interaction method based on an immersed interface

approach in opposition to classical ALE (Arbitrary Lagrangian Eulerian) approaches.

Generalized Finite Elements, together with Lagrange Multipliers, are used to provide velocity

and pressure discontinuities on the fluid domain across the immersed interface. To couple

both fluid and structural problems, an implicit staggered scheme is adopted, which allows the

easy implementation of already developed black box computer codes.

iii

Objetivo

O objetivo deste trabalho o desenvolvimento de um cdigo computacional eficiente e

robusto de soluo de problemas da Mecnica dos Fluidos e Interao Fluido-Estrutura.

Optou-se por uma formulao de Fronteiras Imersas para permitir simulaes de problemas

que envolvem movimentos e transformaes complexas da estrutura. Em problemas com

essas caractersticas, abordagens ALE clssicas costumam perder robustez em razo da

necessidade de reconstruo da malha de fluidos para evitar distoro excessiva dos

elementos. Num contexto mais geral, este trabalho tambm objetiva a incluso do grupo de

pesquisa do Laboratrio de Mecnica Computacional (LMC) da Escola Politcnica da USP

nas reas de Mecnica dos Fluidos Computacional e IFE, como uma etapa inicial de uma

proposta mais ambiciosa, j em curso, de estender o espectro dos problemas aos quais se

aplicam as teorias e formulaes desenvolvidas neste laboratrio. Citam-se aqui como

exemplos de aplicaes futuras, problemas de biomecnica, aeroelasticidade de obras civis e

estruturas aeroespaciais e multifsica. Para este fim, usufruir-se- de trabalhos desenvolvidos

por este grupo de pesquisa, no campo da mecnica das estruturas computacional, que j se

encontram num estgio de maturidade bastante avanado.

1

1. Introduo

Problemas de Interao Fluido-Estrutura (IFE) so de grande importncia para a Engenharia e

esto presentes tanto na natureza quanto em obras e mquinas construdas pelo Homem.

Trata-se de um problema que envolve a Mecnica dos Fluidos e das Estruturas e no qual a

soluo de um domnio depende da soluo do outro, caracterizando assim um sistema

acoplado. Se resolver um problema de mecnica dos fluidos com condies de contorno

diversas j uma tarefa difcil, um problema de IFE acrescenta como desafio a soluo

simultnea do sistema acoplado no qual as condies de contorno na interface entre o fluido e

a estrutura so desconhecidas a priori, pois dependem da soluo do problema. Alm disso,

pelo fato de grande parte dos problemas de interesse para a engenharia envolver grandes

deslocamentos da estrutura e conveco do fluido, a IFE por natureza um problema

fortemente no linear.

Da descrio acima, o leitor j poder perceber que o apelo por mtodos numricos para se

resolver problemas de IFE inevitvel. O advento de computadores mais velozes a preos

cada vez mais acessveis e o uso de computao paralela so tambm fortes propulsores para

o avano de simulaes computacionais de problemas complexos com fidelidade crescente

realidade.

A natureza dos problemas de IFE pode ser muito variada. A Figura 1, por exemplo, ilustra

algumas situaes bastante distintas nas quais se observa este tipo de fenmeno acoplado.

Estruturas offshore como plataformas de explorao de petrleo, de interesse para Engenharia

Civil e de Petrleo, so exemplos tpicos que envolvem IFE (Figura 1 (a)). Nestes casos,

necessria uma formulao do problema capaz de representar a superfcie livre do fluido e,

em algumas situaes, interessante tambm conseguir simular o contato entre as estruturas

flutuantes. Essa linha de estudo atrai hoje diversos grupos de pesquisa no mundo

[1,2,3,4,5,6,7,8] e em especial no Brasil [9,10,11] em virtude da oportunidade de explorao

de petrleo em guas profundas.

2

(a) Plataformas offshore de explorao de petrleo.

(b) Simulao computacional de um paraquedas [12].

(c) Ponte de Tacoma momentos antes do colapso [13].

(d) Tenses na artria aorta de um paciente com aneurisma [14].

(e) Resultados da simulao computacional do escoamento atravs de uma turbina elica [15].

Figura 1. Exemplos de problemas de Interao Fluido-Estrutura.

No contexto da multidisciplinaridade, a bioengenharia ou, mais especificamente, a

biomecnica tem ganhado muito destaque na comunidade cientfica internacional. A

ilustrao da Figura 1 (d), por exemplo, mostra o mapa de tenses tpico de uma simulao

computacional da artria aorta abdominal de um paciente com aneurisma [14] e cujo

carregamento provm do bombeamento de sangue atravs da mesma. Nesta simulao de IFE,

o sangue corresponde ao domnio fluido do problema enquanto que o domnio estrutural

constitudo pelo tecido da artria. O aneurisma uma patologia cardiovascular causada pela

dilatao da artria e cuja ruptura leva ao bito do paciente na maioria dos casos. Quantificar

o risco de ruptura de um aneurisma, entretanto, uma tarefa muito difcil e ao mesmo tempo

de extrema importncia, pois a interveno cirrgica em muitos casos tambm envolve um

alto risco. Neste contexto, a geometria do aneurisma, que varia bastante caso a caso,

desempenha um papel crucial no tratamento desta enfermidade e a simulao numrica do

problema pode revelar informaes mais quantitativas e individualizadas e, portanto,

substanciar melhor a escolha do tratamento mdico mais adequado para cada paciente. Do

ponto de vista numrico, quando se tem fluido e estrutura com densidades semelhantes, como

3

o caso do sangue e artrias do corpo humano, o problema de IFE requer cuidados especiais

para controlar possveis instabilidades na resoluo do sistema acoplado [16,17]. Outras

dificuldades inerentes a este tipo de simulao so a aquisio das imagens do paciente para

gerar o modelo geomtrico, estabelecer a condio inicial do problema de IFE, definir o perfil

de velocidades na fronteira de entrada do escoamento, dentre outras. Para mais referncias de

trabalhos de biomecnica, o leitor ir encontrar em [14,18,19,20,21,22,23,24] trabalhos

tambm sobre aneurisma cerebral, influncia da presso arterial e da geometria do aneurisma

no risco de ruptura do mesmo, simulao de vlvulas cardacas e bombeamento artificial de

sangue em pacientes com deficincia cardaca.

Outro fenmeno no qual a IFE est presente o flutter ou drapejamento. Este fenmeno um

problema de instabilidade aeroelstica e de grande interesse da Engenharia Aeronutica,

porm, tambm algo que pode causar problemas em obras da Engenharia Civil.

Carregamentos no conservativos, como foras devidas ao vento, podem despertar a

instabilidade de uma estrutura uma vez que, em geral, realizam trabalho no nulo num ciclo

fechado e podem aumentar a energia do sistema sucessivamente. Em [25], Argyris promove

uma extensa discusso sobre instabilidade dinmica e esttica de uma estrutura simples

submetida a carregamentos no conservativos. O caso mais famoso de instabilidade

aeroelstica na Engenharia Civil provavelmente o da Ponte de Tacoma (Tacoma Narrows

Bridge [13]) sobre o Estreito de Tacoma no estado de Washington, Estados Unidos. Esta

ponte suspensa colapsou pouco mais de quatro meses aps sua inaugurao no ano de 1940 e

a foto mostrada na Figura 1 (c) ilustra a impressionante magnitude dos deslocamentos com

que a ponte oscilava instantes antes do colapso. Fenmenos de instabilidade aeroelstica

tambm podem ocorrer em superfcies aerodinmicas de avies, como asas, empenagens,

antenas e etc. Uma dificuldade em simular computacionalmente problemas desta natureza est

na alta velocidade do escoamento do fluido, que normalmente requer modelos complexos de

turbulncia [26,27]. Existe tambm o fato de que o tempo de simulao, em geral, precisa ser

grande, para despertar o fenmeno fsico, e a incrementos de tempo pequenos, exigncia de

qualquer simulao numrica de escoamento de fluidos. Esta combinao, naturalmente,

demanda computao de alto desempenho que, muitas vezes, favorece a utilizao de

mtodos para reduo da ordem do problema [28]. A simulao do escoamento do ar em

turbinas elicas ou atravs de um para-quedas so tambm exemplos de problemas de

aeroelasticidade. Na Figura 1 (e), so mostrados resultados da simulao computacional de

4

uma turbina elica com detalhes, direita, do mapa de tenses aerodinmicas (acima) e da

deformada de uma das ps (abaixo). Na Figura 1 (b), mostrado o campo de velocidades do

escoamento atravs de um paraquedas. O leitor interessado em problemas de aeroelasticidade

e sua soluo pelo Mtodo dos Elementos Finitos pode consultar as referncias [29,30]. Em

[31] e [32] so apresentados problemas interessantes de aeroelasticidade aplicada

otimizao de estruturas aeronuticas e controlabilidade de aeronaves, respectivamente.

Em razo da natureza variada dos problemas de IFE, difcil imaginar que uma determinada

abordagem ou mtodo seja superior em todas as classes de problemas. Na realidade, muitos

programas computacionais, comerciais ou acadmicos, perdem robustez em certos problemas

que envolvem transformaes complexas da estrutura, com grandes deslocamentos ou

rotaes, ou problemas cuja interface entre o fluido e a estrutura seja complicada, s para citar

alguns exemplos.

De uma maneira geral, os mtodos numricos existentes hoje para resoluo de problemas de

IFE podem ser classificados em dois grandes grupos dependendo da forma com que tratam a

interface entre o fluido e a estrutura: mtodos de fronteira coincidente e mtodos de fronteira

imersa. A Figura 2 ilustra a diferena entre as duas abordagens mostrando uma disposio

tpica da malha do fluido e do domnio da estrutura s em cada caso.

(a) fronteira coincidente.

(b) fronteira imersa

Figura 2. Classificao dos mtodos numricos de soluo de problemas de IFE quanto descrio da interface: (a) fronteira coincidente; (b) fronteira imersa.

Os mtodos de fronteira coincidente utilizam em sua formulao a abordagem ALE (sigla em

ingls para Arbitrary Lagrangian Eulerian), na qual o problema fluido formulado e

resolvido a partir de uma malha que se deforma para acompanhar o movimento da estrutura.

Isso se deve ao fato da malha do fluido ser construda limitando-se ao domnio fsico do

s s

5

fluido, conforme ilustra a Figura 2 (a). Dessa forma, a malha do fluido se deforma

acompanhando a deformao da estrutura na interface e esta deformao se propaga numa

regio predefinida do fluido. Como a malha do fluido est conectada superfcie molhada da

estrutura, no possvel preserv-la quando a estrutura for submetida a transformaes

complexas e, portanto, a gerao de uma nova malha inevitvel. Alternativamente, mtodos

que se valem do conceito de fronteiras imersas tratam o fluido e a estrutura a partir de dois

domnios superpostos. As malhas so construdas de forma independente e, portanto, o fluido

pode convenientemente ser resolvido utilizando uma malha fixa e indeformvel numa

abordagem Euleriana clssica, como ilustrado na Figura 2 (b). A estrutura, por sua vez,

resolvida de forma Lagrangiana como usual. Como as malhas do fluido e da estrutura so

construdas de forma independente, torna-se necessria alguma tcnica especial para impor o

efeito do movimento da superfcie molhada da estrutura no fluido, pois o movimento desta

interface atua como uma condio de contorno essencial ao fluido. Entretanto, em geral, no

existem ns da malha do fluido ao longo da interface para que esta condio de contorno seja

imposta diretamente, prtica usual do mtodo dos elementos finitos.

possvel tambm se pensar em um mtodo que trata o fluido e a estrutura a partir de uma

abordagem Lagrangiana unificada. Este mtodo, desenvolvido por Oate e Idelson

[33,34,35,36,37] e chamado de mtodo das partculas, mostra-se muito interessante em

aplicaes que envolvem escoamento de fluidos com superfcie livre e desprendimento de

parte do domnio fluido em subdomnios, fenmeno comum em situaes de impacto de

ondas em estruturas offshore ou quebra mares. A grande vantagem desta abordagem est na

eliminao do termo convectivo da equao que descreve o movimento do fluido, porm se

paga o preo de gerar uma nova malha para cada incremento de tempo da simulao.

6

2. Reviso Bibliogrfica

Conforme j mencionado no Resumo deste texto, uma etapa necessria a este trabalho foi o

desenvolvimento de um programa computacional para soluo de problemas de escoamento

incompressvel de fluidos utilizando o Mtodo dos Elementos Finitos (MEF). A literatura

sobre o assunto bastante vasta, sendo que os primeiros trabalhos remontam ao final da

dcada de 1960 e incio da dcada de 1970. Oden, Taylor e Cheng so alguns dos expoentes

precursores e suas valiosas contribuies podem ser observadas em [38,39,40]. Inicialmente

se achou que o mesmo sucesso que o Mtodo dos Elementos Finitos obteve em problemas da

mecnica dos slidos tambm pudesse se repetir em problemas da mecnica dos fluidos.

Entretanto, a existncia de um funcional quadrtico cuja minimizao equivale a satisfazer a

equao diferencial que rege o problema, caracterstica intrnseca dos problemas de mecnica

das estruturas, que so regidos por equaes elpticas auto-adjuntas, no se verifica em

problemas de fluidos. A natureza hiperblica da equao de Navier-Stokes, acentuada em

problemas nos quais a conveco dominante em relao aos efeitos viscosos, faz com que a

soluo numrica por elementos finitos usando projees de Galerkin usuais no seja tima,

alm de susceptvel a problemas de instabilidade numrica. Este fato contribuiu para

descentralizar o esforo da comunidade cientfica na dcada de 1970 na busca por mtodos

numricos alternativos ao MEF para soluo de problemas da mecnica dos fluidos. O

Mtodo das Diferenas Finitas e, principalmente, o Mtodo dos Volumes Finitos [41,42]

atraram a ateno de muitos pesquisadores e profissionais da rea e, como consequncia, este

ltimo ainda hoje utilizado de forma predominante por programas comerciais.

Em 1982, Brooks e Hughes [43] deram uma grande contribuio ao desenvolverem uma

formulao estabilizada de elementos finitos usando o conceito de SUPG (do ingls,

Streamline Upwind/Petrov-Galerkin) para soluo de problemas de escoamento com

conveco dominante. Este artigo foi durante muito tempo uma das referncias mais citadas

no assunto e parece ter dado novo nimo aos pesquisadores adeptos do MEF. Diversos

trabalhos foram publicados logo em seguida como, por exemplo, Argyris et.al. [44]. A

literatura sobre mtodos de estabilizao da conveco bastante ampla e o tpico continua

sendo objeto de pesquisa at os dias de hoje, como se pode observar em [45] e [46] recentes

trabalhos sobre estabilizao em malhas distorcidas e com incrementos de tempo muito

pequenos, respectivamente. Obviamente no se tem aqui a pretenso de explorar todas as

7

referncias no assunto, inclusive porque este no o foco principal do trabalho. Entretanto

alguns trabalhos no podem deixar de ser citados, pois foram, de uma forma ou de outra,

importantes para o desenvolvimento do programa computacional proposto aqui como um dos

objetivos. Dentre eles, citam-se [47,48,49] e em especial o trabalho de Tezduyar e Osawa [50]

que foi usado como base para formulao estabilizada do problema de Navier-Stokes a ser

apresentado no Captulo 3.

Outra dificuldade numrica inerente formulao pelo MEF do problema de Navier-Stokes

para escoamentos incompressveis a satisfao da condio de incompressibilidade. Esta

restrio ao campo de velocidades conduz normalmente a uma formulao mista de elementos

finitos que exige um cuidado especial na escolha dos subespaos de aproximao da

velocidade e da presso. Pensando localmente nos graus de liberdade do elemento finito,

alguns pares de velocidade e presso podem gerar modos esprios de presso caso no

satisfaam a condio de compatibilidade conhecida como inf-sup ou LBB [51], em

referncia a Ladyzhenskaya, Babuska e Brezzi por sua descoberta. Inicialmente, achou-se que

a escolha de pares de velocidade e presso deveria necessariamente satisfazer a condio

LBB, porm sucessivos esforos no sentido de contornar esta restrio levaram a formulaes

de elementos finitos estveis para qualquer escolha de subespaos de aproximao. Ver, por

exemplo, [52,50,53].

Apesar de ter recebido muita importncia nos ltimos anos, o assunto Interao Fluido-

Estrutura (IFE) j atrai pesquisadores h muito mais tempo. No domnio da mecnica

computacional, trabalhos de Peskin [54,55] da dcada de 70 j mostravam ser possvel a

simulao de modelos de vlvulas cardacas. Os modelos usados por Peskin nesta poca eram

bastante simplificados e as referncias citadas no utilizavam ainda em sua formulao o

Mtodo dos Elementos Finitos. Os primeiros trabalhos a utilizarem o MEF para resolver

problemas de IFE o fizeram a partir da abordagem ALE [56,57]. Esta abordagem utilizada

at os dias de hoje [58,59,60,61,62,63] e teve sua origem ainda no Mtodo das Diferenas

Finitas [64]. Entretanto, em situaes nas quais a estrutura est sujeita a transformaes

complexas, ou a prpria interface entre o fluido e a estrutura seja geometricamente

complicada, desejvel se dispor de mtodos alternativos que no utilizem a abordagem

ALE. Estes mtodos so apresentados com diferentes nomes na literatura: fronteira imersa

(immersed boundary), interface imersa (immersed interface), domnio imerso (embedded

domain) ou domnio fictcio (fictitious domain). Este ltimo tambm chamado de mtodos

8

dos volumes fictcios. interessante notar que as referncia [54,55], citadas como pioneiras

em IFE, so tambm os primeiros trabalhos a utilizarem o conceito de volumes fictcios e a

empregar a nomenclatura immersed boundary. Neste texto, ser utilizada a nomenclatura de

fronteira imersa para se referir a esses mtodos.

Se por um lado a utilizao do conceito de fronteira imersa possui como grande vantagem a

utilizao de malhas simples e regulares para resolver o problema fluido, a imposio da

condio de contorno na interface entre fluido e estrutura passa a ser no trivial. No problema

fluido essa condio se traduz em velocidades impostas, ou seja, trata-se de uma condio de

contorno do tipo Dirichlet e diversas tcnicas podem ser utilizadas para imp-la. Nesta tese de

doutoramento, prope-se uma formulao utilizando o Mtodo dos Elementos Finitos

Generalizados (GFEM) em conjunto com Multiplicadores de Lagrange para impor a condio

de contorno na interface do problema de IFE. Embora o GFEM tenha sido originalmente

concebido para aplicao em problemas de mecnica da fratura ( [65] e [66]), bons resultados

tambm tm sido obtidos em IFE. Essa formulao foi utilizada tambm em [67,68] e [69],

porm, nestas referncias, a interface foi discretizada definindo-se ns nas intersees desta

com os lados dos elementos da malha do fluido. O subespao de aproximao dos

Multiplicadores de Lagrange assim gerado desperta problemas de instabilidade numrica por

no satisfazer a condio LBB, conforme j havia sido mostrado por Mes et al. [70], que

classificam esta tcnica de discretizao da interface de desavisada. Mes et al. [70]

propem um algoritmo de discretizao da interface que reduz o subespao de aproximao

dos Multiplicadores de Lagrange atravs de uma seleo dos ns vencedores dentre aqueles

gerados de forma desavisada, garantindo assim estabilidade ao mtodo e taxa de

convergncia tima, porm pagando-se o preo da perda de simplicidade do mtodo. A

formulao proposta no Captulo 5 deste trabalho difere das referncias citadas h pouco por

construir a malha dos Multiplicadores de Lagrange na interface do problema de IFE de forma

independente da malha do fluido, podendo, na maioria das vezes, ser utilizada a prpria malha

da estrutura para este fim. Alm disso, prope-se tambm que as funes de interpolao dos

Multiplicadores de Lagrange sejam constantes dentro dos elementos e descontnuas entre

estes, de modo a melhorar ainda mais a estabilidade numrica do mtodo e simplificar o

clculo das integrais ao longo da interface. Uma ideia muito semelhante foi apresentada em

[71], praticamente simultnea a [72] e [73].

9

Diversas outras tcnicas vm sendo pesquisadas com o mesmo objetivo de impor condies

de contorno do tipo Dirichlet juntamente com mtodos de fronteiras imersas. Destacam-se

aqui o mtodo dos Multiplicadores de Lagrange Distribudos [74,75], no qual os

Multiplicadores de Lagrange so definidos nos ns dos elementos de fluido cortados pela

interface, variaes do mtodo de Nitsche [76,77], mtodos hbrido mistos nos quais os

elementos de fluido intersectados pela interface so enriquecidos com um campo de tenses

adicional para garantir a imposio da condio de Dirichlet na interface [78,79,80] e mtodos

de Galerkin descontnuos [81].

10

3. O Problema de Fluidos

Neste Captulo, ser apresentado o problema da Mecnica dos Fluidos descrito pelas

Equaes de Navier-Stokes para escoamentos incompressveis e sua soluo pelo Mtodo dos

Elementos Finitos. Apesar dos exemplos numricos estarem limitados a duas dimenses do

espao, toda a formulao ser apresentada para o caso geral em trs dimenses.

3.1. Equaes de Navier-Stokes

A conservao de momento linear do fluido dada por

divddt

u

T b , (3.1)

onde a densidade do fluido e b o vetor das foras de volume por unidade de massa. O

vetor das velocidades representado por u e /d dtu a derivada material da velocidade

com relao ao tempo e T o tensor das tenses de Cauchy. A equao de conservao da

massa, que no caso de escoamentos incompressveis se reduz equao de

incompressibilidade, pode ser escrita como

div 0u . (3.2)

Adotando um sistema de referncia Euleriano, a derivada material da velocidade na equao

de conservao do momento linear resultar em duas parcelas: derivada parcial com relao

ao tempo, ou acelerao local, e acelerao convectiva. Portanto, a eq. (3.1) pode ser reescrita

como

div u u u T b , (3.3)

onde u representa a derivada parcial da velocidade com relao ao tempo e o gradiente da

velocidade definido por

,i i u u e . (3.4)

11

Ao longo deste texto, utilizar-se- a conveno ,ia para denotar a derivada do vetor a com

relao coordenada i de um sistema cartesiano definido pelas bases ie . Tambm ser

adotada a conveno de somatria de ndices repetidos, ou mudos, de modo que o vetor

posio seja escrito como i ixx e . Para fluidos Newtonianos, o tensor das tenses de

Cauchy dado por

,2 ( )p T I u (3.5)

onde a viscosidade do fluido, p a presso e (u) o tensor taxa de deformao da

velocidade dado por

1( ) ( )2

T u u u . (3.6)

O divergente do tensor das tenses pode ser calculado a partir de (3.5), resultando em

div div2 ( )p T u u , (3.7)

onde a ltima parcela nula devido condio de incompressibilidade (3.2). Substituindo

ento a equao anterior em (3.3), chega-se finalmente a

em

em

div

2( )

0 ,

p

u u u u b

u

(3.8)

na qual a condio de incompressibilidade foi repetida aqui apenas por motivo de clareza da

exposio e o domnio do problema. Note que todos os termos da equao do momento

foram divididos por , originando a viscosidade e presso cinemticas definidas por

/ e /p p , respectivamente. A fim de que o sistema de equaes diferenciais

parciais governado por (3.8) seja bem posto, necessrio impor condies de contorno e

iniciais ao problema. Dessa forma, definem-se

12

em

em

em

0

( )

0,

u

tp

e t

u u

u I n t

u u

(3.9)

onde u a parte do contorno na qual so impostas as condies de contorno essenciais ou de

Dirichlet e t as condies de contorno naturais ou de Neumann, tais que u t e

u t . 0u a condio inicial do problema e precisa satisfazer a condio de

incompressibilidade, isto , div 0 0u . A barra sobre as variveis indica que seus valores

so prescritos.

O sistema de equaes (3.8) constitui as Equaes de Navier-Stokes para escoamento

incompressvel e sua deduo remonta ao sculo XIX. O leitor interessado em conhecer um

pouco do histrico da mecnica dos fluidos e seus principais personagens, desde Arquimedes

at Stokes, passando por grandiosas contribuies de Bernoulli, Euler e o curioso paradoxo de

dAlambert, encorajado a ler o captulo introdutrio de [82] e notas histricas de [83].

3.2. Forma fraca

Sejam e subespaos de Hilbert das funes de aproximao e peso do campo de velocidades, respectivamente, e cujas condies de contorno de Dirichlet so satisfeitas pelas

funes de aproximao1 enquanto que as funes peso so nulas em u . Esses subespaos

podem ser escritos como

em 1( ) u u u u (3.10)

e

em 01 1 ( )B u w w 0 . (3.11) 1 A rigor, esta definio s seria correta caso as condies de contorno de Dirichlet fossem homogneas, pois s neste caso a soma de duas funes seria ainda uma funo do mesmo subespao. Poder-se-ia definir a funo de aproximao como a soma de duas funes

0D u u u , a primeira satisfazendo as condies de contorno reais do problema e a segunda as condies de contorno homogneas. O

subespao seria ento formado pelas funes 0u e o rigor matemtico seria atendido, porm pagando-se o preo de uma notao menos concisa. Neste texto, o autor optou por cometer este abuso de linguagem uma vez que a imposio das condies de contorno essenciais em elementos finitos direta e trivial, tornando esta discusso um assunto de menor importncia.

13

Adicionalmente, seja 2( ) o espao das funes quadrado integrveis que contm as funes de aproximao e peso da presso. O sistema de equaes (3.8) pode ento ser

reescrito em sua forma integral valendo-se de funes teste arbitrrias ( , )q w da

seguinte forma

,

div

2( )

0.

d d d

pd d

q d

w u w u u w u

w w b

u

(3.12)

Integrando por partes os termos viscoso e da presso, obtm-se

,

div

2 : ( )

( ) .

d d d

pd pd p d

w u w u w u n

w w w n

(3.13)

Substituindo (3.13) em (3.12) e valendo-se da notao compacta definida logo a seguir, pode-

se reescrever o problema definido por (3.8) e (3.9) de forma equivalente conforme se segue.

Encontrar ( , ) 0,t T u x e ( , ) 0,p t T x , tal que, para qualquer par

( , )q w

div div, ; , , , , , ,t

c a p q w u u w u w u w u w t w b ,(3.14)

que a forma fraca da equao de Navier-Stokes para escoamentos incompressveis. Note que

o termo no contorno u nulo, pois w 0 nesta regio. Os termos convectivo e viscoso

so definidos atravs das formas trilinear e bilinear a seguir:

e ; , ( ) , :c d a d

u w u w u u u w w u . (3.15)

14

3.3. Tenses e pseudotenses no contorno

importante ressaltar que o vetor t aplicado no contorno t e apresentado na forma fraca

(3.14) dado por

( ) ,p t u I n (3.16)

estando em perfeita concordncia com (3.13) e (3.9). Ou seja, (3.16) de fato a condio de

contorno natural do sistema (3.8). Entretanto, por definio, o vetor tenso aplicado a uma

superfcie definida pelo vetor normal unitrio n vale

[2 ( ) - ] [ ( ) - ] ,Tp p t Tn u I n u I n (3.17)

que diferente de (3.16). Devido a esta inconsistncia, a tenso (3.16) muitas vezes

chamada de pseudotenso na literatura internacional [84]. Embora em grande parte dos

problemas, nos quais a ateno esteja voltada para os resultados de velocidade e presso do

escoamento, essa distino entre tenso e pseudotenso no seja relevante, nas aplicaes em

que as tenses no contorno do fluido so de suma importncia, a formulao precisa estar

consistente com as tenses verdadeiras dadas por (3.17). Problemas de IFE so um exemplo

de quando utilizar (3.17) e no (3.16).

Uma formulao do problema de escoamento incompressvel cuja condio de contorno

natural dada por (3.17) pode ser facilmente obtida combinando-se as equaes (3.3) e (3.5)

para escrever a equao de conservao do momento linear da seguinte forma

div(2 ) ,s p u u u u b (3.18)

onde novamente toda a equao foi dividida por . O termo s u o gradiente simtrico das

velocidades que, por definio, idntico ao tensor taxa de deformao das velocidades dado

por (3.6). Focando agora na integrao por partes dos termos viscoso e de presso, j na forma

fraca, de modo semelhante ao que foi feito em (3.13), obtm-se

15

div( ,

div

[ 2 )] 2 : 2 [ ( )]

( ) .

s s s sd d d

pd pd p d

w u w u w u n

w w w n

(3.19)

A forma fraca do problema ser ento dada pela mesma expresso de (3.14), ou seja,

div div, ; , , , , , , ,t

c a p q w u u w u w u w u w t w b

( , )q w , porm, com a seguinte definio da forma bilinear do termo viscoso

, 2 :s sa d

u w w u (3.20)

e a condio de contorno natural, ou de Neumann, dada por (3.17) em t .

3.4. Integrao no tempo

H vrios mtodos numricos que podem ser utilizados para discretizar no tempo a equao

(3.14). Ver, por exemplo, trabalho de Dettmer e Peric [85] no qual os autores promovem

extensa discusso e comparao de diversos algoritmos de integrao no tempo. Neste

trabalho, o mtodo de Newmark [86] ser apresentado e sua performance avaliada em funo

do parmetro atravs de vrios exemplos numricos a serem resolvidos mais adiante no

texto.

Mtodo de Newmark

Na discretizao temporal, o tempo discretizado em intervalos t de modo que

1 ,n nt t t

onde 1nt o instante de tempo subsequente a nt . No Mtodo de Newmark [86], a

acelerao do fluido no instante de tempo 1nt pode ser escrita como

1

1 1 ( ) (1 ) ,n n

n n

t

u uu u (3.21)

16

onde o parmetro de integrao de Newmark. interessante notar que, fazendo 1 ,

obtm-se

,1

1n n

n

t

u uu (3.22)

que o Mtodo de Euler Implcito. Este mtodo, apesar de ser de primeira ordem, tem

algumas vantagens como sua simplicidade, facilidade de implementao e convenincia

quando utilizado em problemas de IFE com fronteiras imersas. Por outro lado, ao fazer

1 / 2 , chega-se a

1

1 2( ) ,n n

n n

t

u uu u (3.23)

que, apesar de exigir o armazenamento da acelerao no instante anterior, nt , tem a vantagem

de possuir convergncia de segunda ordem. Utilizando, portanto, o Mtodo de Newmark de

discretizao do tempo, podemos reescrever as equaes de Navier-Stokes (3.8) e suas

condies de contorno e iniciais (3.9) de modo a obter

em

div em

11 1 2 1 1 1

1

1

(1 )( )

0

n nn n n n n n

n

n

pt t

u uu u u u b

u

u u

em

em 1 1 1( )

u

n n nt

n

p

u I n t

u em ,0 0

u

(3.24)

que o sistema de equaes semidiscreto do problema de escoamento viscoso incompressvel.

A forma fraca do problema descrito por (3.24) obtida de forma anloga a (3.14) e ser dada

por

div div

1

, ; , , , ,

, , , (1 ) , ,t

nn n

t c a p q

t

w u u w u w u w u

w b w t w u w u(3.25)

17

( , )q w .

3.5. Discretizao do espao por elementos finitos

Da forma como as equaes de Navier-Stokes foram apresentadas e, posteriormente, a forma

fraca foi deduzida, a discretizao espacial da eq. (3.25) acarretar numa formulao mista de

elementos finitos na qual a velocidade e a presso so as variveis primitivas do problema.

Adicionalmente, o campo de presses pode ser interpretado como Multiplicadores de

Lagrange da condio de incompressibilidade. Dessa forma, de se esperar que a escolha do

elemento finito misto para formulao do problema discreto no possa ser feita de forma

aleatria em relao aos pares de velocidade/presso. Insistindo na interpretao da presso

como um Multiplicador de Lagrange, a escolha de um elemento finito com muitos graus de

liberdade de presso e poucos de velocidade fatalmente causar travamento ou instabilidade

numrica devido ao excesso de restries ao campo de velocidades. A condio a qual o par

de presso/velocidade deve satisfazer a fim de evitar problemas de instabilidade numrica

conhecida na literatura por LBB ou inf-sup. Elementos finitos 2D com pares de

velocidade/presso populares na literatura so mostrados na Figura 3, onde os crculos

menores e cheios representam os ns com graus de liberdade de velocidade, enquanto os

crculos maiores e sem enchimento simbolizam ns de presso. Para outros elementos e

discusses mais aprofundadas da condio LBB, o leitor pode consultar [51], [87] e [84].

(a) Q2Q1 (Taylor-Hood)

(b) P2P1 (Taylor-Hood)

(c) Q2Q2

(d) P2P2

Figura 3. Exemplos de elementos finitos com diferentes pares de velocidade/presso.

Neste trabalho, o autor prope um novo elemento finito triangular que tem mostrado

resultados interessantes em diversas simulaes numricas [88,89,90]. Este elemento possui

seis ns com graus de liberdade de velocidade enquanto que a presso calculada apenas nos

ns dos pontos mdios dos lados, conforme ilustra a Figura 4. Portanto, as velocidades so

18

interpoladas por funes de forma quadrticas compatveis e a presso por funes lineares e

incompatveis, por isso o nome P2P1i atribudo a este elemento.

Figura 4. Elemento finito triangular P2P1i

A discretizao no espao da eq. (3.25) feita ao substituir as funes de aproximao e peso

definidas na Seo 3.2 por suas respectivas aproximaes. Por se tratar de uma aproximao

por elementos finitos, as funes escolhidas possuem dimenso finita e suporte compacto.

Essas funes podem ser escritas como

, ,

,

1 1

1 1

( ) ( )

( ) ( ) ,

Nel Nelh h

u e p ee e

Nel Nelh h

u e p ee e

p

q

u x x

w x x

N u N p

N w N q

(3.26)

onde Nel o nmero de elementos finitos da discretizao. O ndice h indica a aproximao

da funo em questo e Nu e Np so as matrizes que contm as funes de forma ou de

interpolao local do elemento finito e so dadas por

e ,1 2 1 2

u p

pu u u p pu pn n

N N N N N N I I I N N

onde un e pn so os nmeros de ns de velocidade e de presso por elemento,

respectivamente. Os vetores com ndice e referem-se aos valores nodais da varivel no

domnio local do elemento finito. Os vetores de velocidade e presso nodais, por exemplo, so

dados por

TT e 1 1 1 2 2 2 1 2 ,u u u pe n n n e nu v w u v w u v w p p p

u p

uu

u

,u p,u p

,u p

19

com ui, vi e wi sendo as componentes da velocidade do n i nas trs direes ortonormais do

sistema de referncia. Os subespaos que contm as funes de aproximao e peso agora

tero dimenso finita e a notao h , h e h ser adotada. O problema (3.25) pode ento ser enunciado, agora na sua forma discretizada, como: dados b , t , u e conhecida a soluo

do problema no instante de tempo n , encontrar ( )h hu x e ( )h hp x no instante 1n nt t t tal que, para qualquer par ( , )h h h hq w , tenha-se

div

div

1

1, ,

1

, ; , , ,

, , , (1 ) ,

, 0.

t

nh h h h h h h h h

nh h h h n h h n

nh h

t c a p

t t

q

w u u w u w u w

w b w t w u w u

u

(3.27)

Substituindo as expresses de (3.26) na equao anterior e aps algumas manipulaes

algbricas, obtm-se

T

1 1 1 1 1

1

(1 )( )

,

n n n n n n n

n

t t

M Mu C u K u Gp f u Mu

G u 0

(3.28)

onde u e p so os vetores de velocidade e presso nodais globais, que se relacionam com seus

respectivos vetores locais por

e .e e e e u A u p A p (3.29)

T, , , e M C K G G correspondem, respectivamente, matriz de massa, matriz convectiva,

matriz viscosa, operador gradiente e operador divergente. O vetor f contm as foras de

volume e tenses aplicadas. Esses vetores e matrizes globais so dados pelo espalhamento das

contribuies locais de cada elemento. Esta operao de espalhamento dada por

T T e 1 1

,Nel Nel

e e e e ee e

M A M A f A f (3.30)

20

no caso da matriz de massa e do vetor de foras aplicadas, respectivamente. As demais

matrizes so obtidas de forma anloga. Obviamente, as matrizes de espalhamento Ae no

precisam ser definidas no cdigo computacional, sua funo aqui apenas didtica na

exposio. As matrizes locais so definidas por

, ,

,

e

,

,( )

e e

e

e

e te

T Te u u e e u u e

Te u u i e i u e

Te u p e

T Te u e u te

d d

d

d

d d

e

b t

M N N K B B

C N N u N

G N N

f N N

(3.31)

com

,1,1 2,1 ,1

1,2 2,2 ,2

1,3 2,3 ,3

u

u

u

u u un

u u uu n

u u un

N N N

N N N

N N N

I I I

I I I

I I I

B

= e, 1, 2, ,u

u u uu i i i n i

N N N I I IN

.1,1 1,2 1,3 2,1 2,2 2,3 ,1 ,2 ,3u u uu n n nN N N N N N N N N N

Quando se utilizar a equao do momento linear na forma (3.18) e, consequentemente, a

forma fraca do termo viscoso for dada por (3.20), pode-se calcular a matriz viscosa eK

definindo antes o vetor taxa de deformao

,1 ,2 ,3 ,2 ,1 ,3 ,2 ,1 ,3( )T

u v w u v v w w u u (3.32)

e a matriz constitutiva dada por

21

2 0 0 0 0 0

0 2 0 0 0 0

0 0 2 0 0 0,

0 0 0 1 0 0

0 0 0 0 1 0

0 0 0 0 0 1

C (3.33)

de modo que

, ( ) ( ) .Ta d u w w u C (3.34)

A matriz viscosa local, neste caso, ser dada por

,e

Te u u ed

K B C B (3.35)

com

1,1 2,1 ,1

1,2 2,2 ,2

1,3 2,3 ,3

1,2 1,1 2,2 2,1 ,2 ,1

1,3 1,2 2,3 2,2 ,3 ,2

1,3 1,1 2,3 2,1 ,3 ,1

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0

0 0 0

0 0 0

u

u

u

u u

u u

u u

u u un

u u un

u u un

u u u u u u un n

u u u u u un n

u u u u u un n

N N N

N N N

N N N

N N N N N N

N N N N N N

N N N N N N

B .

A obteno das matrizes em (3.31) e (3.35) a partir da substituio das funes peso e

aproximao na forma fraca (3.27) detalhada no Apndice deste texto. Apresentaremos aqui

apenas a obteno da matriz de massa

T T T T

T T T .

,

e

e e

h h h h h he

e

u e u e e e u u e ee e

e e ee

d d

d d

w u w u w u

N A w N A u w A N N A u

w A M A u w Mu

(3.36)

22

Note que ao substituir o resultado da expresso anterior em (3.27), o vetor Tw ser

cancelado uma vez que este aparece em todas as parcelas e por ser arbitrrio. O sistema de

equaes algbricas (3.28) pode ser escrito tambm como

T

11 1

1

(1 )( ).

nn n n n

nt t

M MuC u K G f u Mu

p 0G 0 (3.37)

Note que o sistema de equaes algbricas acima no linear devido presena da matriz

convectiva. Nesta formulao, utilizado o mtodo de Newton-Raphson para soluo de

problemas no lineares e a matriz de rigidez tangente do problema ser obtida na Seo 3.8.

3.6. Estabilizao da condio LBB

Um caso particular de escoamento de fluidos ocorre quando se tem velocidades muito baixas.

Nesta situao possvel desprezar a conveco e a acelerao local do fluido, pois os efeitos

viscosos so predominantes no escoamento. Escoamentos com essa caracterstica so

conhecidos na literatura por escoamento de Stokes e podem ser descritos de forma

simplificada por

T .

K G u fp 0G 0

(3.38)

O escoamento de Stokes portanto um problema estacionrio descrito por um sistema linear

de equaes algbricas. Consequentemente, os possveis problemas de instabilidade numrica

esto restritos condio LBB. A demonstrao matemtica de que o elemento P2Pli satisfaz

ou no a condio LBB bastante complexa e foge ao escopo deste trabalho. Optou-se aqui

por fazer testes e observar possveis modos esprios de presso.

Problema 3.1: Escoamento de Stokes com soluo analtica

Este problema consta na referncia [84] e consiste em se determinar o campo de velocidades e

presso num domnio quadrado de lado unitrio [0,1] [0,1] com condies de contorno

do tipo Dirichlet homogneas em toda a fronteira. As foras de volume so dadas por

23

4 3 2 3 21 2 1 2 1 2 2 2 1

2 3 2 32 2 2 1 2 2 2

2 3 2 22 2 2 1 2 2 1

2 3 4 2 3 42 2 2 2 1 2 2 2

12 24 24 48 48 72 48 12

2 24 72 48 1 4 12 8

8 48 48 12 72 72

4 24 48 48 24 12 24 12 ,

b x x x x x x x x

x x x x x x x

b x x x x x x

x x x x x x x x

(3.39)

onde o vetor posio dado por i ixx e e a fora de volume i ibb e conforme orientao

do sistema de coordenadas mostrado na Figura 5.

Figura 5. Escoamento de Stokes com soluo analtica. Dados do Problema 3.1.

Este problema possui soluo analtica dada por

T

.

2 22 2 3 2 2 31 1 2 2 2 2 2 1 1 1

1 1

1 2 6 4 1 2 6 4

1

x x x x x x x x x x

p x x

u (3.40)

Note que neste problema o campo de presses s pode ser determinado a menos de uma

constante, pois a nica condio de contorno do tipo Dirichlet aplicada velocidade. De

modo a evitar que o sistema (3.38) seja singular, a presso imposta igual a zero no n do

canto superior esquerdo do domnio. Procedendo desta forma, os resultados obtidos com os

elementos Taylor-Hood Q2Q1 e P2P1 so mostrados na Figura 6. Na Figura 7, os resultados

de presso e velocidade vertical ao longo da linha mdia horizontal (x2=0,5) so comparados

com a soluo analtica, onde pode ser observada excelente concordncia dos mesmos. O

software Gid [91] foi utilizado para pr e ps-processamento dos dados do problema.

1e

2e

0,1 0,1

u 0 u 0

u 0

u 0

1

24

Figura 6. Malha, campo de velocidades e presso do escoamento de Stokes. Resultados obtidos usando os elementos Q2Q1 e P2P1.

Figura 7. Comparao dos resultados obtidos com os elementos Q2Q1 e P2P1 e a soluo analtica (exata) do Problema 3.1. Velocidade vertical (esquerda) e presso (direita) ao longo da linha mdia horizontal.

A boa concordncia dos resultados da Figura 7 j era esperada uma vez que os elementos de

Taylor-Hood satisfazem a condio LBB (ver [84] e [51]). Os elementos finitos com funes

de interpolao de mesma ordem na Figura 3, Q2Q2 e P2P2, no satisfazem a condio LBB

e a distribuio de presses obtida neste exemplo possui modos esprios. O elemento P2P1i,

por sua vez, apresenta modos esprios de presso em alguns tipos de malha. As figuras a

seguir mostram alguns resultados deste mesmo exemplo em funo da disposio da malha de

-0,015

-0,010

-0,005

0,000

0,005

0,010

0,015

0 0,2 0,4 0,6 0,8 1

v

x

exato Q2Q1 P2P1

0,00

0,05

0,10

0,15

0,20

0,25

0,30

0 0,2 0,4 0,6 0,8 1

p

x

exato Q2Q1 P2P1

25

elementos finitos. Note que no primeiro caso (Figura 8 e Figura 9), malha estruturada e

cruzada, o campo de presses no possui nenhum sinal de modos esprios ou instabilidade

numrica e pode-se observar excelente concordncia com a soluo analtica. No segundo

caso (Figura 10 e Figura 11), tem-se uma malha no estruturada e uma concordncia tambm

muito boa com a soluo analtica. Finalmente, no terceiro caso (Figura 12 e Figura 13), a

malha estruturada, porm no cruzada. Nesta disposio de malha, o resultado da

distribuio de presso obtido no possui significado fsico e, portanto, no tem qualquer

utilidade. Fica evidenciado assim que o elemento P2P1i no satisfaz a condio LBB e

alguma tcnica de estabilizao se faz necessria para que o mesmo possa ser empregado em

outras aplicaes. Neste trabalho, optou-se por estabilizar os elementos P2P1i, Q2Q2 e P2P2

usando o mtodo Galerkin/Least-squares [52] que ser apresentado na prxima Seo.

Figura 8. Malha, campo de velocidades e presso do Problema 3.1. Resultados obtidos usando o elemento P2P1i com malha estruturada e cruzada.

Figura 9. Comparao dos resultados obtidos com o elemento P2P1i e a soluo analtica (exata) do Problema 3.1 para o caso de malha estruturada e cruzada. Velocidade vertical (esquerda) e presso (direita) ao longo da linha mdia

horizontal.

-0,015

-0,010

-0,005

0,000

0,005

0,010

0,015

0 0,2 0,4 0,6 0,8 1

v

x

exato P2P1i

0,00

0,05

0,10

0,15

0,20

0,25

0,30

0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1

p

x

exato P2P1i

26

Figura 10. Malha, campo de velocidades e presso do Problema 3.1. Resultados obtidos usando o elemento P2P1i com malha no estruturada.

Figura 11. Comparao dos resultados obtidos com o elemento P2P1i e a soluo analtica (exata) do Problema 3.1 para o caso de malha no estruturada. Velocidade vertical (esquerda) e presso (direita) ao longo da linha mdia

horizontal.

Figura 12. Malha, campo de velocidades e presso do Problema 3.1. Resultados obtidos usando o elemento P2P1i com malha estruturada no cruzada.

-0,015

-0,010

-0,005

0,000

0,005

0,010

0,015

0 0,2 0,4 0,6 0,8 1

v

x

exato P2P1i

0,00

0,05

0,10

0,15

0,20

0,25

0,30

0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1

p

x

exato P2P1i

27

Figura 13. Comparao dos resultados obtidos com o elemento P2P1i e a soluo analtica (exata) do Problema 3.1 para o caso de malha estruturada no cruzada. Velocidade vertical (esquerda) e presso (direita) ao longo da linha

mdia horizontal.

Estabilizao da condio LBB usando Galerkin/Least-Squares

Esta tcnica de estabilizao da condio LBB permite a utilizao de elementos finitos com

qualquer par de velocidade/presso e foi desenvolvida originalmente por Hughes e Franca

[52]. A referncia [84] tambm aborda esta tcnica num contexto mais recente. A ideia central

est na minimizao do resduo quadrado da equao de conservao do momento em (3.8)

dado por

,2 2( , ) ( , )QR p p p u u b u b (3.41)

onde os termos transiente e convectivo foram desprezados. A minimizao deste funcional

pode ser feita exigindo que a derivada de Gteaux em u e p, com real, seja nula:

.0

( , )0 ( , )Q

dR p qq

d

u ww (3.42)

De (3.42) e (3.41), pode-se definir uma funo real ( ) dada por

,

.2 2 2 2

( ) ( , )

( ) ( , )

QR p q

p q p q

u w

u w b u w b(3.43)

-0,015

-0,010

-0,005

0,000

0,005

0,010

0,015

0 0,2 0,4 0,6 0,8 1

v

x

exato P2P1i

-0,20

-0,10

0,00

0,10

0,20

0,30

0,40

0,50

0,60

0 0,2 0,4 0,6 0,8 1

p

x

exato P2P1i

28

Como

,0

( , )( 0)Q

dR p q

d

u w (3.44)

e

,2 2( 0) 2( , )p q u b w (3.45)

a eq. (3.42) equivale a

.2 2( , ) 0 ( , )p q q u b w w (3.46)

Finalmente, a equao anterior deve ser adicionada forma fraca do problema. Procedendo

dessa forma, o problema de Stokes estabilizado e discretizado por elementos finitos pode ser

escrito como

div div

,2 21

( , ) ( , ) , ( , ) ( , )

( , )

t

e

h h h h h h h h

Nelh h h h h

ee

a p q

p q

w u w u w b w t

u b w (3.47)

( , )h h h hq w . A fim de evitar maiores exigncias de continuidade das funes peso

e de aproximao da velocidade, em razo da derivada de segunda ordem, o termo de

estabilizao admitido atuar apenas no interior dos elementos [84]. O parmetro e ,

chamado de parmetro de estabilizao, dado por [52]

,2

2e

eh

(3.48)

onde eh uma medida local do tamanho do elemento e a viscosidade cinemtica do fluido,

j definida anteriormente. A escolha do coeficiente depende das funes de forma da

velocidade e presso. 2,5 e eh definido como sendo a razo entre a rea e o permetro do

elemento tm fornecido bons resultados para o elemento P2P1i , como ser mostrado mais

adiante.

29

importante salientar aqui que, segundo Hughes e Franca [52], a eq. (3.47) vlida apenas

para elementos finitos com funes de aproximao da presso contnuas entre os elementos.

A utilizao de elementos com presso descontnua, caso do P2P1i, requereria a incluso de

um termo adicional nesta equao que necessita ser integrado ao longo das interfaces dos

elementos. Este termo adicional dado por

,

1

,2 e

Nelh he

e

hq p

onde hp

e hq

so os saltos da presso e da funo peso nas interfaces entre os elementos

adjacentes, respectivamente, e e a interface dos elementos. No caso particular do P2P1i,

apesar destes saltos serem no nulos ao longo de e , a contribuio deste termo parece no

ser significativa. Todos os resultados que sero apresentados nas prximas Sees foram

obtidos sem o termo adicional apresentado acima.

Outro aspecto interessante com relao estabilizao do elemento P2P1i dada pela eq. (3.47)

que, exceo das foras de volume, que dependem do problema a ser resolvido, os demais

termos da estabilizao so todos constantes, pois envolvem a segunda derivada da

velocidade, que quadrtica, e a primeira derivada da presso, que linear. Dessa forma, o

custo computacional ao se estabilizar o P2P1i muito baixo.

Utilizando agora a formulao estabilizada (3.47), o problema de Stokes descrito na Seo

anterior pode ser simulado novamente a fim de eliminar os modos esprios de presso

observados no terceiro caso da Figura 13. Os resultados obtidos so mostrados na Figura 14 e

Figura 15.

Figura 14. Malha, campo de velocidades e presso do escoamento de Stokes. Resultados obtidos utilizando o elemento P2P1i e formulao estabilizada com =2,5.

30

Figura 15. Presso ao longo da linha mdia horizontal obtida utilizando o elemento P2P1i e formulao estabilizada com =2,5.

A mesma formulao tambm foi utilizada para estabilizar os elementos Q2Q2 e P2P2. Os

resultados obtidos para esses elementos so mostrados da Figura 16 Figura 19.

Figura 16. Malha, campo de velocidades e presso do escoamento de Stokes. Resultados obtidos utilizando o elemento Q2Q2 e formulao estabilizada com =1.

Figura 17. Presso ao longo da linha mdia horizontal obtida utilizando o elemento Q2Q2 e formulao estabilizada com =1.

0,00

0,05

0,10

0,15

0,20

0,25

0,30

0 0,2 0,4 0,6 0,8 1

p

x

exato P2P1i

0,00

0,05

0,10

0,15

0,20

0,25

0,30

0 0,2 0,4 0,6 0,8 1

p

x

exato Q2Q2

31

Figura 18. Malha, campo de velocidades e presso do escoamento de Stokes. Resultados obtidos utilizando o elemento P2P2 e formulao estabilizada com =7.5.

Figura 19. Presso ao longo da linha mdia horizontal obtida utilizando o elemento P2P2 e formulao estabilizada com =7,5.

A Figura 20 mostra a comparao entre as taxas de convergncia da presso entre os

elementos finitos do tipo Taylor-Hood (Figura 3) e o elemento P2P1i (Figura 4) em funo do

nmero de divises do domnio do problema. O erro foi calculado atravs de

exato

exato

erro .2

2

( )hp p

p

0,00

0,05

0,10

0,15

0,20

0,25

0,30

0 0,2 0,4 0,6 0,8 1

p

x

exato P2P2

32

Figura 20. Comparao entre as taxas de convergncia da presso entre os elementos de Taylor-Hood e o P2P1i. Valores calculados para 4, 8, 16 e 32 divises do domnio.

Os elementos com funes quadrticas de aproximao da presso possuem taxa de

convergncia maior que os demais e foi omitido. Note que o elemento P2P1i converge mais

rapidamente que os elementos de Taylor-Hood, justificando seu uso apesar da necessidade de

estabilizao da condio LBB.

3.7. Estabilizao da conveco

Conforme mencionado na introduo deste texto, a soluo do problema de Navier-Stokes por

elementos finitos utilizando projees usuais de Galerkin, eq. (3.37), susceptvel a

instabilidades numricas em problemas com conveco dominante. Neste trabalho, utiliza-se

o conceito de SUPG para estabilizao da conveco com base nos trabalhos de Brooks e

Hughes [43] e Tezduyar e Osawa [50]. A formulao estabilizada do problema, discretizado

apenas no espao para simplificar a exposio, dada por

SUPG

div

+ , ,

div

1

, , ; , , , ,

( ) ( ) 0

, 0,

t

e

h h h h h h h h h h h

Nelh h h

e

h h

a c p

q

w u w u u w u w w b w t

w u u

u

(3.49)

1,00E-07

1,00E-06

1,00E-05

1,00E-04

1,00E-03

1,00E-021 10 100

Erro

na

pres

so

N divises

P2P1i

Q2Q1

P2P1

33

( , )h h h hq w . ( )hu o resduo da equao de conservao do momento linear

dado por

2( ) ( ) .h h h h h h hp u u u u u b (3.50)

De modo anlogo eq. (3.47), o termo de estabilizao da eq. (3.49) admitido atuar apenas

no interior dos elementos. O coeficiente SUPG calculado por (ver [50])

SUPG

1/

1 2 3

1 1 1r

r r rS S S

(3.51)

onde

, , 1 2 3 12e e

S S S Se e

tRe

C C

C M (3.52)

e

, e

.

2

( ) ,( ) ( ) ,e e

h h h h h h he e

he

e

Re

w u u u w u u

u

C M

C

C

(3.53)

Re o nmero de Reynolds definido localmente no elemento. Neste trabalho, foi utilizado

2r em (3.51) e as normas foram calculadas fazendo Ttr ( )( ) .

Note que a introduo do termo de estabilizao na eq. (3.49) acrescentar um novo termo ao

sistema de equaes algbricas. Este termo detalhado a seguir

SUPG SUPG

TSUPG

T TSUPG

,

.

,

,

,

( ) ( ) ( ) ( )

( ) ( )

( ) ( )

ee

e

e

h h h h hu i e i e

h hi u i e e

h he i u i e

d

d

d

w u u e u u

e u u

e u u

N w

N w

w N

(3.54)

34

O sistema de equaes (3.28) pode ser reescrito de modo a incluir o termo de estabilizao

anterior como

T ,

1 1 1 1 1 1

1

(1 )( ) )n n n n n n n n

n

t t

M Mu C u K u Gp r(u f u Mu

G u 0

(3.55)

onde

T TSUPG e .,1

( ) ( )

e

Nelh h

e e e i u i ee

d

e u ur A r r N (3.56)

3.8. O problema tangente

O sistema de equaes (3.55) no linear devido ao termo convectivo e de estabilizao. O

mtodo de Newton-Raphson ser utilizado para resolv-lo de forma iterativa e, portanto, se

faz necessria a determinao dos operadores tangentes. Reescrevendo o problema (3.55) na

forma

,( ) F d 0 (3.57)

onde d o vetor de deslocamentos generalizados dado por

.1

1

n

n

ud

p (3.58)

A soluo do problema (3.57) pode ser obtida de forma iterativa por

,1

( ) ( )k

k k

A d F dd d

(3.59)

onde k o nmero da iterao e

35

( ) ( )kk

A d F d

d (3.60)

a matriz de rigidez tangente do problema. Observe no sistema (3.55) que apenas os termos

convectivo e de estabilizao so no lineares. Os operadores tangentes desses termos podem

ser obtidos derivando-os em relao a u e p e, posteriormente, espalhando-os de modo

anlogo a (3.30) para obter ( )A d . Dessa forma, tem-se para o termo convectivo

T

T

,

,

,

( ) ( )

( )

( ) ( )

( ) ( ) ( )

e

e

e

h h h h h h

h h he

e

u e u i e i u e ee

h hu e u i e i u e e

e

d

d

d

d

w u u w u u

w u u

e

e u u

u u

u

N A w N A u N A uu

N A w N A I N A I

T T T T .,( ) ( )e

h he u i u i u u e e

e

d

e u uw A N N N N A

(3.61)

Observe que o termo de estabilizao em (3.49) definido localmente no domnio do

elemento. Portanto, calcularemos suas matrizes tangentes locais

SUPG

SUPG

,

,

( ) ( )

( ) ( ) ( ) ( )

e

e

h h h

e

h h h h h he

e e

d

w u u

w u u w u u

u

u u

(3.62)

e

SUPG SUPG, ,( ) ( ) ( ) ( )e

e

h h h h h he

e e

d

w u u w u up p (3.63)

36

onde as derivadas presentes no integrando das equaes anteriores so calculadas a seguir:

T T T ,

,

,

( ) ( ) ( ) ( )

( )

h h h hu i e i u

e

he u i i u

w u u e u

u e

N w N Iu

w N N

(3.64)

T T )

2

, , ,

,

( ) ( ) ( ) ( ) ( )

( ) ( ) ( )

( (

h h h h h h h h h

e e e e

h h huu i e i j u j u u jj

h ue u i i j

t

t

w u u w u u u u u

e u e u u

e u e

u u u u

NN w N N N

Nw N , ,) ( )

h hu j u u jj

u uN N N

(3.65)

e

T

T T

)

, ,

, ,

( ) ( ) ( )

( ) ( )

( ( )

h h h h h h

e e

hu i e i j p j

he u i i j p j

p

w u u w u

e u e

e u e

p p

N w N

w N N .

(3.66)

Substituindo as expresses (3.66) em (3.63) e (3.65) e (3.64) em (3.62) e, em seguida,

substituindo (3.62) e (3.61), aps espalhamento, em (3.60), pode-se escrever a matriz tangente

do problema (3.55) da seguinte forma

TT

,1( )n

t

MC u K G

A dG 0

(3.67)

onde TC e G so dadas pelo espalhamento de suas respectivas matrizes locais

37

T T T TSUPG

T )

, ,

, , ,

( ) ( ) ( )

( ( ) ( )

e

h h hTe u i u i u u u i i u

h h huu i i j u j u u jj edt

e u u u e

e u e u u

C N N N N N N

NN N N N

(3.68)

e

TSUPG ., ,( )( )e

he e u i i j p j ed

e u eG G N N (3.69)

3.9. Exemplos Numricos

Nesta seo sero apresentados alguns exemplos numricos com a finalidade de verificar a

formulao apresentada at o momento.

Problema 3.2: Escoamento numa cavidade

Este um problema clssico da literatura e consiste em se determinar o escoamento em

regime estacionrio em um domnio quadrado de lado unitrio [0,1] [0,1] e com

condies de contorno do tipo Dirichlet na velocidade em todo o contorno do problema,

conforme ilustra a Figura 21. O nmero de Reynolds dado por

,1Re

u l

onde l o lado do quadrado. A velocidade horizontal 1u imposta no contorno superior do

domnio corresponde, portanto, ao prprio nmero de Reynolds do problema. A viscosidade e

densidade do fluido so tambm mostradas na Figura 21. O domnio do problema foi

discretizado com 2339 elementos finitos triangulares P2P1i conforme Figura 22.

38

Figura 21. Dados do Problema 3.2: escoamento numa cavidade.

Figura 22. Malha de elementos triangulares P2P1i utilizada no problema de escoamento numa cavidade. 1104 elementos e 2339 ns.

As figuras seguintes mostram alguns resultados obtidos com este elemento para nmero de

Reynolds 1, 100, 400 e 1000. O campo de presses para este problema vale zero em quase

todo o domnio, possuindo duas singularidades nos cantos superiores. Por esta razo a malha

bastante refinada nestes cantos.

Figura 23. Vetores de velocidade, linhas de corrente e campo de presses para Re=1.

1e

2e

0,1 0,1

u 0 u 0

u 0

1 1uu e

1

39

Figura 24. Linhas de corrente superpostas ao campo de velocidades. Re=100, 400 e 1000.

Note que, medida que o nmero de Reynolds aumenta de 1 at 100, a posio do vrtice

principal se desloca para a direita e surge um vrtice secundrio no canto inferior direito.

Aumentando o Reynolds at 400 e 1000, o vrtice principal tende a retornar a linha mdia

vertical, porm numa posio mais abaixo comparada com Re=1. H tambm o aparecimento

de um terceiro vrtice no canto inferior esquerdo quando se atinge Re=1000. As posies dos

vrtices principais para Re=100, 400 e 1000 so mostradas na Tabela 1 e comparadas com

resultados da literatura. Na Figura 25 so mostrados os grficos da velocidade horizontal (u),

normalizada pelo nmero de Reynolds, ao longo da linha vertical mdia (x=0,5) para os

diferentes regimes de escoamento.

x1 x2

Re=100

Presente estudo 0,614 0,736

Donea e Huerta (2003) 0,620 0,740

Tuann e Olson (1978) 0,610 0,722

Re=400

Presente estudo 0,558 0,605

Donea e Huerta (2003) 0,568 0,606

Tuann e Olson (1978) 0,506 0,583

Re=1000

Presente estudo 0,533 0,565

Donea e Huerta (2003) 0,540 0,573

Ozawa (1975) 0,533 0,569

Tabela 1. Posio do vrtice principal para Re=100, 400 e 1000.

40

Figura 25. Velocidades horizontais ao longo da linha vertical mdia. esquerda, resultados obtidos com o elemento P2P1i; direita, resultados de Donea e Huerta (2003) [84].

A convergncia do mtodo de Newton-Raphson tambm foi investigada com o objetivo de

validar a formulao do problema tangente descrito na seo 3.8. Os grficos do erro da

velocidade e da presso em funo da iterao so mostrados na Figura 26 para Re=100. Os

dados utilizados nos grficos da Figura 26 so listados na Tabela 2.

Figura 26. Convergncia do mtodo de Newton-Raphson para Re=100.

Iterao Erro em u Erro em p

1 2.227E-01 4.601E-01

2 2.420E-02 4.630E-02

3 2.616E-04 1.200E-03

4 5.783E-08 1.697E-07

5 2.751E-15 5.213E-15

Tabela 2. Convergncia do mtodo de Newton-Raphson para Re=100.

0,0

0,1

0,2

0,3

0,4

0,5

0,6

0,7

0,8

0,9

1,0

-0,4 -0,2 0,0 0,2 0,4 0,6 0,8 1,0

StokesRe=100Re=400Re=1000

u

y

u

y

0,00

0,05

0,10

0,15

0,20

0,251 2 3 4 5

Erro

na

velo

cida

de

Iterao

0,00

0,10

0,20

0,30

0,401 2 3 4 5

Erro

na

pres

so

Iterao

41

Problema 3.3: Escoamento em torno de um cilindro

Este um exemplo clssico da literatura e pode ser encontrado em diversas referncias

bibliogrficas. Ver [50,43,92,93], por exemplo. Entretanto, para fins de verificao da

formulao apresentada at o momento, os parmetros da referncia [93] sero utilizados,

uma vez que esta referncia rene vrios resultados de simulaes do problema realizadas por

diversos grupos de pesquisa em atividade. O escoamento se d atravs de um canal estreito,

ver Figura 27, de comprimento m2,2L e altura m0,41H , com origem do sistema de

referncia no canto inferior esquerdo. Prximo entrada do canal, h um obstculo de seo

circular com raio m0,05r e centro posicionado nas coordenadas (0,20;0,20)C .

Figura 27. Escoamento em torno de um cilindro. Desenho esquemtico do Problema 3.3.

A velocidade na entrada do canal tem perfil parablico dado por

22 2 2(0, ) 4 ( ) / 0 ,T

mx U x H x H u

onde o vetor posio dado por i ixx e e mU a velocidade mxima na entrada do canal.

O nmero de Reynolds e de Strouhal so definidos neste exemplo, respectivamente, como

e ,UD fDRe StU

com a velocidade mdia U na entrada do canal definida por

2 (0, / 2)

,3H

U u

com D sendo o dimetro do cilindro e f a frequncia de desprendimento dos vrtices em s-1.

a viscosidade cinemtica do fluido, j definida anteriormente. A velocidade no topo e na

H

L(0,0)

C r+

42

base do canal vale zero. A sada do canal assumida como uma fronteira aberta, ou seja, com

tenso nula aplicada. A densidade e a viscosidade cinemtica do fluido valem 3kg/m1,0

e m /s3 210 , respectivamente. A simulao feita para duas situaes de escoamento:

estacionrio e transiente. Os parmetros do problema para essas condies so mostrados na

Tabela 3.

Estacionrio Transiente

3kg/m[ ] 1 1

m /s3 2[10 ] 1 1

m/s][mU 0,3 1,5 Re

20 100

Tabela 3. Parmetros do fluido para escoamento estacionrio e transiente do Problema 3.3.

Dois tipos de elementos foram utilizados para obteno dos resultados que sero mostrados

mais adiante: P2P1 e P2P1i. A malha utilizada, ver Figura 28, foi a mesma em todas as

simulaes.

Figura 28. Malha do fluido utilizada na simulao de escoamento em torno de um cilindro do Problema 3.3. 3848 elementos triangulares e 7952 ns.

O mtodo de Newton-Raphson foi forado a executar pelo menos duas iteraes para soluo

do sistema. Da Figura 29 Figura 31, so mostrados alguns resultados do escoamento

estacionrio (Re=20) para uma avaliao qualitativa dos resultados. Visualmente, no h

diferena entre os resultados obtidos com o elemento P2P1 e o P2P1i.

Figura 29. Campo de velocidades para Re=20. Resultados obtidos com o elemento P2P1.

43

Figura 30. Campo de presso para Re=20. Resultados obtidos com o elemento P2P1.

Figura 31. Vetores de velocidade prximos ao cilindro (esquerda) e tenses na superfcie do cilindro (direita). Re=20. Resultados obtidos com o elemento P2P1.

Para uma avaliao mais quantitativa, foram calculados os coeficientes de sustentao Lc e de

arrasto Dc definidos por

e2 2

,v hL DF F

c cUD UD

(3.70)

onde vF e hF so, respectivamente, a resultante vertical e horizontal da fora que atua no

cilindro. Foram calculados tambm a diferena de presso entre o ponto de estagnao

montante do cilindro (0,15 0,20)Tx e o bordo de fuga do mesmo em (0,25 0,20)Tx

e o comprimento da regio de recirculao ou de formao dos vrtices L . Esses resultados

esto listados na Tabela 4.

Lc Dc P L

Presente trabalho (P2P1) 0,0102 5,5921 0,1178 0,0845

Presente trabalho (P2P1i) 0,0110 5,6336 0,1183 0,0851

[R. Rannacher, S. Turek] 0,0106 5,5755 0,1173 0,0780

Tabela 4. Comparao dos resultados obtidos com os elementos P2P1 e P2P1i e resultados de [93] para o caso de escoamento estacionrio do Problema 3.3.

44

Percebe-se, dos resultados apresentados na Tabela 4, muito boa concordncia com os

resultados da literatura, exceo dos valores de L calculados com os elementos P2P1 e

P2P1i. O autor desconhece uma definio precisa do comprimento de recirculao, tendo o

mesmo sido estimado graficamente e, portanto, no se esperava boa concordncia deste

resultado com o resultado da literatura. No caso do escoamento transiente, o tempo foi

discretizado em diversos intervalos diferentes a fim de se analisar o desempenho do mtodo

de integrao apresentado em 3.4. O mtodo de Newton-Raphson novamente foi forado a

executar pelo menos duas iteraes por incremento de tempo. Da Figura 32 Figura 34, so

mostrados alguns resultados de velocidade, presso e tenso na superfcie do cilindro,

novamente para fins de anlise qualitativa. Os resultados mostram o escoamento aps o

desprendimento alternado de vrtices se estabilizar e assumir um comportamento peridico.

Mais uma vez, no h distino visual entre os resultados obtidos com os elementos P2P1 e

P2P1i.

Figura 32. Campo de velocidades para Re=100. Resultados obtidos com o elemento P2P1i.

Figura 33. Contornos de presso para Re=100. Resultados obtidos com o elemento P2P1i.

45

Figura 34. Vetores de velocidade prximos ao cilindro (esquerda) e tenses na superfcie do cilindro (direita) para Re=100. Resultados obtidos com o elemento P2P1i.

Os coeficientes de arrasto e de sustentao tambm foram computados ao longo do tempo e

seus grficos so mostrados na Figura 35. A Tabela 5, na sequncia, compara os coeficientes

de arrasto e sustentao mximos e o nmero de Strouhal aps estabilizar o regime peridico

de desprendimento de vrtices com resultados da literatura.

Figura 35. Coeficientes de arrasto (CD) e sustentao (CL) em funo do tempo para o caso de escoamento transiente do Problema 3.3 (Re=100).

CDmax CLmax St

Presente estudo (P2P1) 3,2450 1,0369 0,3008

Presente estudo (P2P1i) 3,3038 1,0746 0,3008

[R. Rannacher, S. Turek] 3,2300 1,0000 0,3000

Tabela 5. Comparao dos resultados obtidos com os elementos P2P1 e P2P1i e resultados de [93] para o caso de escoamento transiente do Problema 3.3.

2,95

3,00

3,05

3,10

3,15

3,20

3,25

3,30

3,35

1 3 5 7 9 11 13

CD

Tempo (s)

-1,50

-1,00

-0,50

0,00

0,50

1,00

1,50

1 3 5 7 9 11 13

CL

Tempo (s)

46

A convergncia dos resultados da Tabela 5, medida que o incremento de tempo decresce,

mostrada nos grficos das figuras a seguir.

Figura 36. Convergncia do coeficiente de arrasto (CD) em funo do incremento de tempo t e do parmetro de integrao do Mtodo de Newmark para o caso de escoamento transiente do Problema 3.3.

Figura 37. Convergncia do coeficiente de sustentao (CL) em funo do incremento de tempo t e do parmetro de integrao do Mtodo de Newmark para o caso de escoamento transiente do Problema 3.3.

Figura 38. Convergncia do nmero de Strouhal (St) em funo do incremento de tempo t e do parmetro de integrao do Mtodo de Newmark para o caso de escoamento transiente do Problema 3.3.

2,95

3,00

3,05

3,10

3,15

3,20

3,25

3,30

3,35

0,000 0,005 0,010 0,015 0,020

CD

T

P2P1i =1

P2P1i =1/2

P2P1 =1

P2P1 =1/2

0,00

0,20

0,40

0,60

0,80

1,00

1,20

0,000 0,005 0,010 0,015 0,020

CL

T

P2P1i =1

P2P1i =1/2

P2P1 =1

P2P1 =1/2

0,284

0,286

0,288

0,290

0,292

0,294

0,296

0,298

0,300

0,302

0,000 0,005 0,010 0,015 0,020

St

T

P2P1i =1

P2P1i =1/2

P2P1 =1

P2P1 =1/2

47

3.10. Concluso Parcial

Neste Captulo, a formulao do problema de escoamento incompressvel de fluidos

Newtonianos descrito pelas Equaes de Navier-Stokes e sua implementao computacional

usando o Mtodo dos Elementos Finitos foram apresentadas. Questes essenciais como a

estabilizao da condio LBB, estabilizao da conveco e integrao no tempo foram

discutidas e diversos exemplos numricos foram resolvidos com o objetivo de verificar o

contedo exposto. Os resultados obtidos nas simulaes numricas mostraram muito boa

concordncia com resultados da literatura e algumas concluses podem ser destacadas:

A etapa inicial deste trabalho obteno de um programa para simulao

computacional de problemas da Mecnica dos Fluidos foi concluda com sucesso;

Em problemas de IFE, deve-se utilizar o termo viscoso definido por (3.20) para

garantir consistncia no clculo das tenses provenientes do problema fluido na

interface com a estrutura;

O elemento triangular P2P1i, proposto neste trabalho, fornece resultados comparveis

aos elementos Taylor-Hood P2P1 e Q2Q1, porm a um custo de se estabilizar a

condio LBB. Este custo, entretanto, muito baixo devido ao grau dos polinmios de

interpolao deste elemento;

Ao discretizar o tempo, deve-se optar por mtodos de segunda ordem. Os resultados

de convergncia no te