Metodo dos elementos finitos

  • View
    63

  • Download
    16

Embed Size (px)

DESCRIPTION

Tese mestrado

Text of Metodo dos elementos finitos

  • 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