UMA FORMULAÇÃO ESTÁVEL DE ELEMENTOS FINITOS …objdig.ufrj.br/60/teses/coppe_d/MiltonAlvesGoncalvesJunior.pdf · Milton Alves Gonçalves Junior Tese de Doutorado apresentada ao

UMA FORMULAÇÃO ESTÁVEL DE ELEMENTOS FINITOS COM PRECISÃO

TEMPORAL DE 2ª ORDEM PARA SIMULAÇÃO DE ESCOAMENTOS

INCOMPRESSÍVEIS

Milton Alves Gonçalves Junior

Tese de Doutorado apresentada ao Programa de

Pós-graduação em Engenharia Civil, COPPE, da

Universidade Federal do Rio de Janeiro, como

parte dos requisitos necessários à obtenção do

título de Doutor em Engenharia Civil.

Orientadores: Alvaro Luiz Gayoso de Azeredo

Coutinho

Paulo Augusto Berquó de Sampaio

Rio de Janeiro

Outubro de 2011



INCOMPRESSÍVEIS


TESE SUBMETIDA AO CORPO DOCENTE DO INSTITUTO ALBERTO LUIZ

COIMBRA DE PÓS-GRADUAÇÃO E PESQUISA DE ENGENHARIA (COPPE) DA

UNIVERSIDADE FEDERAL DO RIO DE JANEIRO COMO PARTE DOS

REQUISITOS NECESSÁRIOS PARA A OBTENÇÃO DO GRAU DE DOUTOR EM

CIÊNCIAS EM ENGENHARIA CIVIL.

Examinada por:

Prof. Alvaro Luiz Gayoso de Azeredo Coutinho, D.Sc.

Prof. Paulo Augusto Berquó de Sampaio, Ph.D.

Prof. José Luis Drummond Alves, D.Sc.

Prof. Paulo Roberto Maciel Lyra, Ph.D.

Prof. Eduardo Gomes Dutra do Carmo, D.Sc.

RIO DE JANEIRO, RJ – BRASIL

OUTUBRO DE 2011

ii

Gonçalves Junior, Milton Alves

Uma Formulação Estável de Elementos Finitos com

Precisão Temporal de 2ª Ordem para Simulação de

Escoamentos Incompressíveis/ Milton Alves Gonçalves

Junior. – Rio de Janeiro: UFRJ/COPPE, 2011.

XIII, 100, p.: il.; 29,7 cm.

Orientadores: Alvaro Luiz Gayoso de Azeredo

Coutinho


Tese (doutorado) – UFRJ/ COPPE/ Programa de

Engenharia Civil, 2011.

Referências Bibliográficas: p. 83–90

1. Formulações estabilizadas de elementos finitos. 2.

Fluidodinâmica computacional. 3. Estruturas de dados por

aresta. I. Coutinho, Alvaro Luiz Gayoso de Azeredo et al..

II. Universidade Federal do Rio de Janeiro, COPPE,

Programa de Engenharia Civil. III. Título.

iii

À minha esposa Luciana e à nossa amada filha Maria Luíza, a

grande obra de nossas vidas, que desde a sua chegada nos ilumina

e preenche nossos corações de muita felicidade.

iv

Agradecimentos

Eu chego ao fim do doutorado, após 4 anos de trabalho, com 2 trabalhos publicados

em periódicos, 4 trabalhos apresentados em congressos e uma apresentação inédita em

um congresso internacional realizado na Europa. Por isso, gostaria de fazer um breve

agradecimento a todos aqueles a quem pude recorrer e que tornaram esse feito possível.

Para alguns, essa é somente a ordem natural da vida. Para mim, uma emblemática e

incomensurável conquista social.

Aos meus orientadores, os pesquisadores Paulo Augusto Berquó de Sampaio e

Alvaro Luiz Gayoso de Azeredo Coutinho, pela motivação, apoio, compreensão e pela

garantia de acesso todas as vezes em que eu precisei. Pela oportunidade de realização de

pesquisas e conclusão desta importante etapa profissional em minha vida. Pela

realização de trabalhos e promoção de divulgação destes a comunidade científica

através de apresentações em congressos aqui e fora do Rio de Janeiro. E pelos grandes

exemplos de pessoas e de profissionais que são.

Aos meus pais, por tudo aquilo que eu não teria condições de retribuir por maior

que fosse o meu tempo de vida e esforço em fazê-lo.

Aos colegas de doutorado José Camata, Carlos Henao, Gabriel Guerra, José

Guilherme e Anderson Mendonça por toda ajuda desde o período de realização das

disciplinas até o final do curso, pelos momentos de descontração extra COPPE e pela

troca de idéias onde eu sempre fui o beneficiado. Ao colega de doutorado André Rossa

e à pesquisadora Paula Sesini pelo compartilhamento de trabalhos importantes para a

minha tese.

Ao professor e pesquisador Renato Nascimento Elias pela amizade, pela grande

contribuição na minha formação e pelo valioso suporte no entendimento do EdgeCFD®

e na utilização do cluster do NACAD. Também por todas as dicas que eu obtive com

esse profissional de tamanha competência.

Ao professor e pesquisador Paulo de Tarso Themistocles Esperança pela

oportunidade de realização de pesquisas de interesse da área naval, pelo acolhimento no

v

Hidrolab, pelos recursos disponibilizados para o trabalho e pelas conversas, às vezes de

cunho pessoal, durante o período em que trabalhamos juntos. Ao pesquisador Marcelo

Vitola pelo suporte na utilização do cluster do LabOceano, de onde obtive resultados

muito importantes para os meus estudos.

À Mara Prata, pelas prazerosas conversas todas as vezes em que me dirijo ao

NACAD, independente do motivo que me leva até o laboratório, e por toda a ajuda

nesses anos na solução de problemas institucionais. Sempre muito gentil, simpática,

solícita, prestativa e eficiente nos nossos atendimentos.

Ao engenheiro naval Cristiano H. P. Clemente, meu companheiro de Hidrolab,

responsável pela minha rápida adaptação ao laboratório e pela agradabilíssima parceria

e convivência.

À Comissão Nacional de Energia Nuclear pelo suporte financeiro, fundamental à

minha manutenção na UFRJ e para realização da pesquisa.

vi

Resumo da Tese apresentada a COPPE/UFRJ como parte dos requisitos necessários

para a obtenção do grau de Doutor em Ciências (D.Sc.)



INCOMPRESSÍVEIS


Outubro/2011

Orientadores: Alvaro Luiz Gayoso de Azeredo Coutinho


Programa: Engenharia Civil

Esta tese apresenta o desenvolvimento de uma formulação de elementos finitos

estabilizada com precisão temporal de 2ª ordem para as equações quase-incompressíveis

de Navier-Stokes. A implementação 3D do método foi realizada utilizando a estrutura

do software EdgeCFD®, onde características como reordenação de dados, estrutura de

dados por aresta e paralelismo foram mantidas.

Problemas típicos de mecânica de fluidos foram explorados para verificar a

qualidade dos resultados e o desempenho computacional do método.

vii

Abstract of Thesis presented to COPPE/UFRJ as a partial fulfillment of the

requirements for the degree of Doctor of Science (D.Sc.)

A SECOND-ORDER TIME ACCURATE STABLE FINITE ELEMENT

FORMULATION FOR SIMULATION OF INCOMPRESSIBLE FLOWS


October/2011

Advisors: Alvaro Luiz Gayoso de Azeredo Coutinho


Department: Civil Engineering

This thesis presents the development of a second-order time accurate stabilized

finite element formulation for quasi-incompressible Navier-Stokes equations. The 3D

implementation of the method was performed using the EgdeCFD® software frame,

where the features like data reordering, edge-based data structure and parallelism were

kept.

Representative problems of fluid dynamics were exploited to verify the quality of

results and the computational performance of the method.

viii

Sumário

Capítulo 1 ................................................................................................................................... 1

Introdução ................................................................................................................................. 1

1.1 Motivação ....................................................................................................................... 1

1.2 Objetivo .......................................................................................................................... 3

1.3 Contribuições da Tese ...................................................................................................... 3

1.4 Estrutura da Tese ............................................................................................................. 7

Capítulo 2 ................................................................................................................................... 9

Equações Governantes e Formulação de Elementos Finitos ..................................................... 9

2.1 Escoamento quase-incompressível .................................................................................... 9

2.2 Formulação com precisão temporal de 2ª ordem ............................................................... 10

2.3 A equação de atualização da pressão ............................................................................... 11

2.4 As equações para atualização da velocidade ..................................................................... 14

2.5 Passos de tempo locais e sincronização ............................................................................ 19

Capítulo 3 ................................................................................................................................. 24

Exemplos Numéricos Bidimensionais .................................................................................... 24

3.1 Escoamento induzido no interior de uma cavidade pelo movimento de uma tampa .............. 24

3.2 Escoamento com expansão abrupta ................................................................................. 29

3.3 Desprendimento periódico de vórtices em escoamento cruzado sobre um cilindro circular ... 33

3.4 Escoamento de Kim-Moin .............................................................................................. 41

Capítulo 4 ................................................................................................................................. 48

Implementação Tridimensional ............................................................................................... 48

4.1 Características básicas do EdgeCFD® .............................................................................. 48

4.2 Diferenças entre o método proposto e o método presente no EdgeCFD® ............................. 50

4.3 Estrutura de dados ......................................................................................................... 53

4.4 Matrizes de aresta .......................................................................................................... 54

Capítulo 5 ................................................................................................................................. 61

Exemplos Numéricos Tridimensionais ................................................................................... 61

5.1 Escoamento induzido no interior de uma cavidade pelo movimento de uma tampa – 3D ...... 61

Capítulo 6 ................................................................................................................................. 79

Conclusões .............................................................................................................................. 79

Trabalhos Futuros .................................................................................................................... 81

ix

Referências ............................................................................................................................... 83

Apêndices .................................................................................................................................... 91

Apêndice A .................................................................................................................................. 92

Apêndice B .................................................................................................................................. 95

Equações Governantes ............................................................................................................ 95

Formulação de Elementos Finitos ........................................................................................... 96

x

Lista de Figuras

Figura 1. Malha de base utilizada para 1 000L

Re . e 5 000L

Re . , contendo 2.601 nós e

5.000 elementos. .......................................................................................................................... 25

Figura 2. Escoamento induzido numa cavidade. Resultados obtidos com o método proposto:

Iso-linhas de pressão (a) e vorticidade (b) para 1 000L

Re . ; Iso-linhas de pressão (c) e

vorticidade (d) para 5 000L

Re . . ........................................................................................... 26

Figura 3. Velocidades na linha de centro da cavidade para 1 000L

Re . . Componente xu (a) e

componente yu (b). .................................................................................................................... 27


Re . . Componente xu (a)

e componente yu (b). .................................................................................................................. 28

Figura 5. Geometria e condições de contorno do problema de escoamento com expansão

abrupta. ........................................................................................................................................ 29

Figura 6. Malha de base estruturada utilizada no escoamento com expansão abrupta com 1.249

nós e 2.304 elementos triangulares. ............................................................................................ 30

Figura 7. Escoamento com expansão abrupta. Caso de referência: Malha estruturada com 74.497

nós e 147.456 elementos. Resultado do campo de velocidade obtido com o método proposto

para 100H

Re em 0

0 1 /t , H u . ........................................................................................ 31

Figura 8. Erro relativo da velocidade e tempo de CPU consumido para os métodos de 1ª e 2ª

ordem. Escoamento com expansão abrupta em diferentes níveis de refinamento. .................... 32

Figura 9. Escoamento com expansão abrupta. Refinamento adaptativo da malha intermediária (

min 0,03125h e máx 0,125h ) em 0

0 01 /t , H u . ........................................................... 33

Figura 10. Geometria e condições de contorno do problema de escoamento cruzado sobre um

cilindro circular. .......................................................................................................................... 34

Figura 11. Escoamento cruzado sobre um cilindro circular para 100d

Re . Campos de

velocidade (à esquerda) e pressão (à direita) nos instantes 0

4 /t d u ((a) e (b)), 0

12 /t d u

((c) e (d)), 0128 /t d u ((e) e (f)) e 0

188 /t d u ((g) e (h)). ............................................... 36

Figura 12. Evolução temporal dos coeficientes de arrasto (CD) e sustentação (CL) para o

escoamento cruzado sobre um cilindro circular para. ................................................................. 37

xi


escoamento cruzado sobre um cilindro circular para 125d

Re . ............................................ 38


escoamento cruzado sobre um cilindro circular para 150d

Re . ............................................ 38

Figura 15. Escoamento cruzado sobre cilindro. Malha refinada em 0

200 /t d u para

100d

Re . Tamanho característico mínimo de 0 02, d e máximo de 0 40, d . ........................ 39

Figura 16. Comparação entre os resultados do presente estudo com os dados obtidos em [35] e

[36]. ............................................................................................................................................. 40

Figura 17. Comparação entre as predições obtidas nas simulações para o número de Strouhal e

os resultados da correlação Strouhal-Reynolds [37]. .................................................................. 41

Figura 18. Malha estruturada triangular utilizada para o estudo do Escoamento de Kim-Moin.

(a) 20 20 , 441 nós e 800 elementos e (b) 68 68 , 4.761 nós e 9.248 elementos. ................. 43

Figura 19. Escoamento de Kim-Moin. Campo de velocidade (a) e pressão (b) no instante

1 0 t . s . .................................................................................................................................. 43

Figura 20. Evolução temporal do erro da componente de velocidade xu obtida com os métodos

de precisão temporal de 1ª e 2ª ordem para a 4 malhas menos refinadas. ................................... 46

Figura 21. Desmembramento da matriz de rigidez do elemento tetraédrico em suas matrizes de

aresta correspondentes................................................................................................................. 55

Figura 22. Desmembramento da matriz simétrica de rigidez do elemento tetraédrico em suas

matrizes de aresta e bloco diagonal correspondentes. X

representa os coeficientes “repetidos”

ou armazenados em outra sub-matriz. ......................................................................................... 59

Figura 23. Escoamento induzido no interior de uma cavidade 3D pelo movimento de uma

tampa. Malha não estruturada com 108.104 elementos tetraédricos e 20.589 nós...................... 62

Figura 24. Linhas de corrente no escoamento induzido no interior de uma cavidade 3D pelo

movimento de uma tampa para 1 000L

Re . . ......................................................................... 63

Figura 25. Velocidades na linha de centro da cavidade cúbica para 1 000L

Re . . Componente

xu (a) e componente

yu (b). Malha não estruturada de 108.104 elementos tetraédricos e 20.589

nós. .............................................................................................................................................. 64

Figura 26. Desempenho computacional. Tempo de processamento (a), speedup (b) e eficiência

computacional (c) para os métodos de 1ª e 2ª ordem. Escoamento induzido no interior de uma

cavidade 3D pelo movimento de uma tampa para 1 000L

Re . . Malha tetraédrica com

108.104 elementos e 20.589 nós. ................................................................................................ 66

xii


computacional (c) para os métodos de 1ª e 2ª ordem. Tolerância do solucionador linear de 10-6

para pressão e para velocidade. ................................................................................................... 68

Figura 28. Desempenho computacional. Tempo de processamento (a) e speedup (b) e eficiência


para pressão e de 10-3

para as componentes da velocidade. ........................................................ 69




para as componentes da velocidade. Malha tetraédrica com 1.613.455

elementos e 279.362 nós. ............................................................................................................ 70

Figura 30. Desempenho computacional. Tempo de processamento (a) e speedup (b) para os

métodos de 1ª e 2ª ordem e o EdgeCFD®. Malha não estruturada com 108.104 elementos

tetraédricos e 20.589 nós. ............................................................................................................ 76

xiii

Lista de Tabelas

Tabela 1. Resumo dos dados usados para comparação entre os métodos de precisão temporal de

1ª e 2ª ordem. ............................................................................................................................... 32

Tabela 2. Erro das componentes de velocidade e pressão em 1 0 t . s e 2 0 t . s para as

cinco malhas. ............................................................................................................................... 44

Tabela 3. Taxa de convergência do erro espacial em 1 0 t , s e 2 0 t , s considerando as

cinco malhas do exemplo. ........................................................................................................... 45

Tabela 4. Principais características entre os métodos implementados no EdgeCFD®. ............... 52

Tabela 5. Coeficientes armazenados para a matriz de rigidez para as estrutura de dados por

aresta e por elemento e para os métodos de 2ª ordem e o implementado originalmente no

EdgeCFD®. .................................................................................................................................. 56

Tabela 6. Tempo de processamento e número de iterações por passo de tempo para solução das

equações de pressão e velocidade em função da tolerância do solucionador linear para os

métodos de 1ª e 2ª ordem. Tolerância de 10-3

para a pressão e para a velocidade. ..................... 71




para a pressão e para a velocidade. ..................... 72




para a pressão e de 10-3

para a velocidade. ......... 72

Tabela 9. Características do cluster do LabOceano. ................................................................... 74

Tabela 10. Demanda de memória de cada método para o problema do escoamento induzido pelo

movimento de uma tampa em uma cavidade cúbica 3D. Malha não estruturada com 108.104

elementos tetraédricos e 20.589 nós. ........................................................................................... 75

Tabela 11. Tempo de processamento por passo de tempo para solução das equações de pressão e

velocidade em função do número de processadores para o EdgeCFD® original e o presente

método. Malha não estruturada com 108.104 elementos tetraédricos e 20.589 nós. .................. 77

1

Capítulo 1

Introdução

1.1 Motivação

A Fluidodinâmica Computacional – em inglês, Computational Fluid Dynamics

(CFD) – tem se tornado uma grande aliada da engenharia nos últimos anos. Tendo o seu

desenvolvimento paralelo ao de computadores cada vez mais potentes, algoritmos e

métodos numéricos mais eficientes, ela possibilita a análise de problemas extremamente

complexos, cuja reprodução é inviável devido a limitações que vão desde a

especificação de condições físicas adequadas até o alto custo financeiro. Técnicas

computacionais avançadas tem sido o alvo de pesquisas recentes que visam à redução

do custo computacional e produção de resultados confiáveis. Dentre elas podemos citar

solucionadores avançados lineares e não-lineares.

Outras áreas de pesquisa importantes para fluidodinâmica computacional são o

desenvolvimento de métodos de estabilização e modelagem da turbulência. Estes são

necessários à simulação de vários problemas típicos de engenharia. Métodos de

estabilização são fundamentais na simulação de escoamentos onde o fenômeno de

convecção é tão relevante quanto o de difusão. Em problemas de convecção dominante,

onde o número de Reynolds é elevado, além do método de estabilização, um modelo de

turbulência eficiente é indispensável para garantir a representação do escoamento em

todas as escalas de interesse.

2

O tratamento de escoamentos, que envolvem simplesmente mecânica de fluidos, é

regido por equações oriundas de balanço de massa e quantidade de movimento[1], [2].

Para fluidos com características específicas, essas equações são reduzidas às equações

de Navier-Stokes. Além de problemas puramente mecânicos, existem também

problemas que tratam do acoplamento entre o movimento do fluido e transporte de

energia, substâncias químicas ou propriedades físicas quaisquer. Tais problemas podem

ser tratados pela adição de equações características de transporte dessas propriedades às

equações de Navier-Stokes e estabelecendo uma relação entre elas. Projetos de

trocadores de calor, análise de desgaste e corrosão de tubulações devido à interação

fluido-estrutura ou estratificação térmica, análise de dispersão atmosférica de gases e

escoamento de correntes marítimas ao redor dos risers e amarras de plataformas são

exemplos de problemas reais presentes na indústria nuclear e offshore.

Métodos computacionais são utilizados há vários anos para a análise de segurança

de reatores nucleares. Apesar de confiáveis, vários códigos desenvolvidos para a análise

de sistemas primários de centrais nucleares e análise de acidentes severos foram

baseados em modelos físicos simplificados e no emprego de correlações experimentais

e empíricas. O desenvolvimento da fluidodinâmica computacional tem proporcionado o

estudo de problemas diversos de interesse da engenharia nuclear e pode ser tornar uma

ferramenta indispensável de projeto termo-hidráulico de reatores em um futuro próximo

[3]. Dentre os problemas para os quais a fluidodinâmica computacional pode trazer

benefícios, podemos citar corrosão de circuitos primários, secundários e núcleos de

reatores, diluição de boro em vasos de pressão, reações químicas, escoamento bifásico e

comportamento da interface gás/líquido e dispersão atmosférica de material radioativo.

Diante do grande desafio nacional de atender uma demanda energética cada vez

maior, movida principalmente pelo aumento populacional e pelo estágio de

desenvolvimento científico e econômico em que o país se encontra, algumas

ferramentas tecnológicas podem ser fundamentais nos estudos acerca da exploração de

petróleo e na geração de energia núcleo-elétrica, auxiliando na consolidação do Brasil

como uma futura potência mundial.

3

1.2 Objetivo

O objetivo do trabalho é o desenvolvimento de uma formulação estabilizada de

elementos finitos para simulação em 3D de escoamentos puramente mecânicos, mas que

possa ser estendida à solução de problemas com transferência de calor, em regimes de

convecção natural, forçada e mista, além do transporte de propriedades físicas.

Problemas dessa natureza podem ser encontrados tanto em análise de segurança de

sistemas nucleares atualmente em operação quanto em projetos de novos sistemas

avançados, em especial nos que prevêem remoção passiva de calor.

O EdgeCFD®

é um software desenvolvido pela COPPE / UFRJ para simulação de

problemas de mecânica de fluidos. As características do software são apresentadas no

Capítulo 4. Durante a pesquisa foram abordadas questões como a adequação da nova

formulação proposta às técnicas computacionais implementadas no EdgeCFD® e a

viabilidade de sua utilização alternativa no código computacional.

1.3 Contribuições da Tese

Um grande número de formulações estabilizadas de elementos finitos têm sido

propostas com o intuito de superar as deficiências que o método de Galerkin padrão

apresenta quando aplicado a dinâmica de fluidos. Uma dessas dificuldades está

associada à condição de incompressibilidade em escoamentos viscosos incompressíveis.

A condição de compatibilidade Babuška-Brezzi, (também conhecida como condição de

compatibilidade LBB – Ladyzhenskaya-Babuška-Brezzi) ou ainda, condição inf-sup,

diz que os espaços de interpolação para a pressão e para a velocidade não podem ser

escolhidos aleatoriamente, deve existir alguma relação entre eles. Em problemas

viscosos incompressíveis, a equação referente ao balanço de massa se reduz a um

divergente nulo do campo de velocidades [4] (daí, condição de incompressibilidade).

Como a pressão só aparece na equação de quantidade de movimento, o sistema

algébrico para a determinação dos valores nodais de pressão e velocidade, oriundo do

método clássico de Galerkin, resulta em uma matriz particionada com uma sub-matriz

nula em sua diagonal. A solução do sistema depende de uma escolha adequada dos

espaços de interpolação para a pressão e para a velocidade. Em [5], [6] e [7]

4

encontramos algumas das mais populares combinações de pressão-velocidade bem

sucedidas em atender a condição de compatibilidade de Babuška-Brezzi.

Existem também os métodos de penalidade – do inglês, penalty methods. O

conceito consiste de uma relaxação da condição de incompressibilidade, considerando

que o problema incompressível pode ser tratado como um problema levemente

compressível. Neste caso, a condição de incompressibilidade deixa de existir e pode ser

feita a escolha do mesmo espaço de interpolação para a pressão e velocidade. Há uma

vasta literatura sobre método de penalidade para escoamentos incompressíveis. Alguns

estudos podem ser encontrados em [8], [9], [10], [11], [12] e [13].

Em suma, algumas formulações de elementos finitos permitem contornar as

restrições de Babuška-Brezzi na escolha dos espaços de interpolação na aproximação

das equações incompressíveis de Navier-Stokes escritas em variáveis primitivas.

Formulações alternativas permitem tal feito, mesmo utilizando combinações de pressão-

velocidade que seriam instáveis na formulação clássica de Galekin [14], [15].

Outra dificuldade referente ao método de Galerkin aplicado a problemas de

dinâmica de fluido diz respeito às oscilações espúrias que aparecem normalmente na

simulação de problemas de convecção dominante. Essas oscilações são indesejadas e

incompatíveis com uma representação confiável do escoamento.

As formulações estabilizadas de elementos finitos correspondem basicamente ao

método padrão de Galerkin adicionado de termos extras responsáveis pela estabilização

do método. Ao mesmo tempo em que os termos extras contribuem na estabilização do

método, eles não afetam a sua consistência, tendendo a zero com o refinamento. As

formulações estabilizadas podem ser freqüentemente interpretadas como uma

aproximação de resíduos ponderados de Petrov-Galerkin, onde a função de ponderação

usual de Galerkin é modificada com a adição de uma perturbação. Os termos resultantes

da interação da perturbação com os resíduos produzem o efeito desejado de

estabilização, sem comprometer a consistência da aproximação. Usualmente, os termos

de estabilização aparecem multiplicados por parâmetros de estabilização, interpretados

como escalas de tempo intrínsecas que definem a quantidade de estabilização requerida,

5

dependendo do tamanho local da malha, velocidade e propriedades físicas. Uma revisão

da literatura sobre o tema pode ser vista em Donea e Huerta [15].

Certamente, as formulações estabilizadas são consideradas bem sucedidas em lidar

com problemas de convecção dominante possibilitando o controle das oscilações

espúrias, além de contornar as restrições de Babuška-Brezzi na escolha dos espaços de

interpolação para a pressão e para a velocidade [16], [17].

As formulações estabilizadas aplicadas a escoamentos transientes podem ser

caracterizadas quanto a seqüencia adotada para a discretização das equações

governantes com relação às variáveis independentes. Nas formulações convencionais,

primeiramente é realizada a discretização espacial das equações governantes e a

estabilização é introduzida de forma fraca (formulação variacional). Este tratamento

conduz a um conjunto de equações diferenciais ordinárias. As equações resultantes são

chamadas de equações semi-discretas e é necessária a utilização de uma técnica de

integração numérica para o avanço temporal [18]. Algumas definições para o parâmetro

de estabilização, presente nos termos extras propostos a priori, podem ser encontrados

em [19], [20], [21].

Algumas formulações estabilizadas seguem uma seqüencia diferente. A

discretização temporal é realizada e, somente depois, é realizada a discretização

espacial. Dentre elas podemos citar o método Taylor-Galerkin e a formulação de

Galerkin aplicada aos métodos da família θ – do inglês, θ family methods - [15], [22].

Em [23], De Sampaio apresentou uma formulação estabilizada na qual a estrutura

dos termos de estabilização surge naturalmente da aplicação de uma aproximação de

mínimos quadrados da discretização temporal do balanço da quantidade de movimento.

Ao contrário da prática usual, uma discretização temporal com base em diferenças

finitas precede a aproximação espacial de elementos finitos. Como resultado, o passo de

tempo t utilizado na discretização temporal do balanço da quantidade de movimento

desempenha o mesmo papel que o parâmetro de estabilização, chamado de escala de

tempo intrínseca, em outras formulações. O fato de que o parâmetro de estabilização é o

próprio passo de tempo implica que há um parâmetro a menos para definir em

aplicações transientes. Além disso, ele fornece uma indicação para a escolha do passo

6

de tempo, ou do parâmetro de estabilização, de acordo com a escala de tempo do

processo físico dominante. Maiores detalhes podem ser encontrados em [23] e [24].

Entretanto, com o objetivo de introduzir a quantidade correta de estabilização em

todo o domínio de análise, o passo de tempo pode ser definido localmente, conduzindo

a uma distribuição do passo tempo espacialmente variável. Nos trabalhos anteriores

[25], [26] e [27], algoritmos especiais foram empregados para considerar o uso de tal

distribuição e sincronizar a computação em problemas transientes. Neste trabalho, no

entanto, o procedimento adotado em [23] é seguido. Neste caso, a utilização do passo de

tempo local e o esquema de sincronização necessário estão inseridos no método. O

resultado é um método que se assemelha às bem conhecidas formulações estabilizadas

que empregam um único passo de tempo para todo o domínio e uma definição local dos

parâmetros de estabilização, mas que cujas origens são baseadas no uso do passo tempo

local combinado com o esquema de sincronização, como pode ser visto em [23], [24],

[28].

Nesta tese um método de elementos finitos para escoamentos quase-

incompressíveis é apresentado. Assim como nos trabalhos anteriores apresentados

acima, a discretização temporal precede a discretização espacial. Entretanto, aqui a

discretização temporal empregada é melhorada com uma precisão de 2ª ordem.

No método apresentado uma equação para a pressão é deduzida a partir de uma

aproximação de Taylor-Galerkin de 2ª ordem que combina as leis de conservação de

massa e quantidade de movimento. Em cada passo de tempo, uma vez que a pressão foi

determinada, o campo de velocidade é calculado a partir da solução das equações

discretizadas obtidas a partir de esquema específico de precisão de 2ª ordem e da

minimização do resíduo espacial da quantidade de movimento. Os termos que

estabilizam o método (controlando oscilações espúrias e contornando a condição de

Babuška-Brezzi) surgem naturalmente do processo, sem que sejam introduzidos a priori

na formulação variacional.

A nova formulação possibilita a construção de sistemas totalmente segregados ou

semi-segregados (ou parcialmente acoplados) dos graus de liberdade. No esquema de

solução semi-segregado, inicialmente são resolvidos os graus de liberdade de pressão e,

7

na seqüência, são resolvidos os graus de liberdade de velocidade, onde as suas

componentes formam um sistema acoplado. As matrizes oriundas são simétricas, o que

torna a nova formulação mais flexível em termos de solucionadores dos sistemas de

equações.

Outro fato importante é que o esquema de discretização utilizado na dedução da

nova formulação dispensa a utilização de um solucionador não linear. A não

linearidade, que está presente no termo convectivo da equação de quantidade de

movimento, é eliminada na nova formulação. O método de discretização lineariza esse

termo.

A linearização mencionada no parágrafo anterior torna dispensável a utilização de

um esquema especial para o avanço no tempo da computação. Em outras palavras, não é

realizada integração numérica temporal.

Além desses aspectos relacionados ao esquema de solução e que são dependentes

da formulação de elementos finitos, o estudo apresenta a avaliação do desempenho

computacional também a partir da estrutura de dados utilizada. Essa avaliação é

realizada através de comparações entre as estruturas de dados por elemento (element-by-

element - EBE) e por aresta (edge-based).

1.4 Estrutura da Tese

Nas próximas seções é realizada uma descrição acerca do estudo realizado para a

elaboração da tese. O Capítulo 2 apresenta o modelo físico referente às equações

governantes para o escoamento transiente quase-incompressível e as respectivas

condições de contorno. A formulação de elementos finitos proposta, de precisão

temporal de 2ª ordem, e o caminho percorrido na dedução das equações discretizadas

também são abordados no Capítulo 2. No Capítulo 3 são apresentados os resultados

obtidos com a nova formulação para simulações bidimensionais. Algumas

peculiaridades da adequação do método ao EdgeCFD® considerando as técnicas

presentes no software e a implementação 3D são tratadas no Capítulo 4. O Capítulo 5

apresenta os resultados de simulações de um problema típico tridimensional obtidos

8

com a nova formulação implementada no EdgeCFD®. O Capítulo 6 apresenta as

conclusões e sugestões para trabalhos futuros.

9

Capítulo 2

Equações Governantes e Formulação de

Elementos Finitos

2.1 Escoamento quase-incompressível

Consideremos um modelo contínuo para escoamentos viscosos quase-

incompressíveis. O problema é definido em um domínio limitado Ω, com contorno Γ,

contido no espaço Euclidiano n-dimensional. O escoamento é governado pelas equações

quase-incompressíveis de Navier-Stokes. Estas são apresentadas em coordenadas

cartesianas utilizando a convenção de soma a = 1, 2,..., n e b = 1, 2,..., n.

2

10a

a

p u

t xc (1)

0a a ab

bb b a

u u pu

t x x x (2)

onde / ( / )( / )t p t p e 2/ 1 /p c .

Neste modelo, as variáveis independentes são as componentes de velocidade e a

pressão, representadas por au e p , respectivamente. A tensão viscosa é dada por

10

ab a b b au x u x , onde é a viscosidade do fluido. A velocidade do

som no meio é denotada por c . A massa específica do fluido é indicada por .

O modelo é completado pelas condições iniciais e de contorno para o campo de

velocidade. As condições de contorno de velocidade e de forças por unidade de área

aplicadas ao contorno são prescritas nas partições das fronteiras não sobrepostas ua

e

ta, de modo que

ua ta:

onde ab

é o delta de Kronecker e bn denota as componentes cartesianas do vetor

normal ao contorno e que aponta para fora dele. As condições de contorno de pressão e

velocidade normal à fronteira são prescritas nas partições da fronteira não sobrepostas

p e G , de modo que p G :

2.2 Formulação com precisão temporal de 2ª ordem

A formulação de elementos finitos a seguir possui características semelhantes à

formulação apresentada anteriormente em [23], [24] e [28]. Assim, ela pertence a uma

classe de elementos finitos inerentemente estáveis sem que sejam estabilizados por

termos extras adicionados à formulação variacional. A grande diferença é que enquanto

as formulações apresentadas nos trabalhos anteriores mencionados apresentam precisão

temporal de 1ª ordem, aqui foi empregado um esquema de discretização para garantir a

precisão temporal de 2ª ordem, promovendo maior qualidade dos resultados em

problemas transientes.

, a a uau u tx x (3)

,ab ab b a ta

p n t tx x (4)

, p

p p tx x (5)

,b b Gu n G tx x (6)

11

2.3 A equação de atualização da pressão

Com o objetivo de deduzir a equação de atualização da pressão, consideramos a

expansão temporal da pressão em série de Taylor. Podemos escrever então:

120

2

nn n np p p t pt

t t t t (7)

onde 0 1 . Os subscritos n e 1n indicam o nível do tempo e t denota o

passo de tempo. A variação da pressão durante o passo de tempo é representado por

1n np p p .

Substituindo o balanço de massa dado pela Equação (1) na Equação (7), obtemos:

2

2

10

2

nn

a a

a a

u up tt

t x t xc (8)

Somando e dividindo o segundo termo do lado direito pela massa específica do

fluido e alternando ax e t , a Equação (8) é reescrita como:

2

2

10

2

nn

a a

a a

u up tt

t x x tc (9)

Nesta etapa estamos prontos para introduzir o balanço de momento dado pela

Equação (2) na Equação (9). Considerando que o termo do gradiente de pressão na

Equação (2) é tomado no nível de tempo 1n , o que é equivalente a adotar 1 para

o termo do gradiente de pressão e 0 para os demais termos, obtemos a seguinte

equação:

2

2

10

2 aa a

p t pF t

t x xc (10)

12

Como 1 para o gradiente de pressão, temos que

1n n

a a a

p p p

x x x. E,

portanto:

2

n

a a ab

a ba b b b a

u u ptF u

x x x x x (11)

Note que o procedimento descrito acima começa com a expansão temporal da

pressão em série de Taylor e depois incorpora os balanços de massa e de quantidade de

movimento.

Empregando o método clássico de Galerkin na Equação (10), usando a identidade

de Green e introduzindo as condições de contorno dadas pela Equação (6), obtemos um

sistema simétrico de equações para calcular os valores nodais de atualização da pressão:

2

1

ˆ ˆ1

2

ˆ ˆ ˆˆ

2

ˆ

2G

ii

a a

n n n

i a abn

ba b b a

n

a n nii

a

Np t pN d d

t x xc

t N u pu d

x x x x

u NN d G G d

x

(12)

Para as variáveis do problema temos ˆn na j aju N u , ˆn n

j jp N p e ˆ

j jp N p .

Note que iN e

jN representam as funções de forma do elemento finito e as variáveis

com os subscritos i e j são valores nodais.

De forma mais compacta, a Equação (12) pode ser reescrita utilizando notação

matricial:

pM fp (13)

Em termos da contribuição de cada elemento, a matriz M é definida como:

13

M M

nel

e

eA (14)

M Me e

ij (15)

M = 2

1

2e e

je e iij i j

a ae

Nt NN N

x xc t (16)

O índice e e nel indicam um elemento específico e o número de elementos de uma

dada malha, respectivamente. A matriz Me é a matriz referente a cada elemento e e

A é o operador de assemblagem dos elementos – do inglês, assembly operator. Meij

são os coeficientes da matriz referente ao elemento e .

O vetor pf é definido como:

p pf f

nel

e

eA (17)

p pf fe e

i (18)

pf =

1

ˆ ˆ ˆˆ

2

ˆ

2

e

eG

n n ni a abe n

i ba b b a

na n ni

ia

t N u pu

x x x x

u NN G G

x

(19)

No presente método foi usada a mesma ordem de interpolação para velocidade e

pressão. Entretanto, é importante ressaltar que quando a formulação mista padrão de

Galerkin é usada para aproximar escoamentos incompressíveis, deve-se escolher

interpolações compatíveis para velocidade e pressão. Mais precisamente, os espaços de

aproximação para velocidade e pressão devem satisfazer a condição de Babuška-Brezzi

[14]. Em particular, a escolha da mesma ordem de interpolação para velocidade e

pressão não funciona. Essa limitação é bem conhecida e surge do fato de que no método

padrão de Galerkin o campo de pressão não é diretamente relacionado à conservação de

massa. A pressão simplesmente não aparece na equação do balanço de massa que, para

14

escoamentos incompressíveis, se reduz a 0a au x , configurando o campo de

velocidade com divergência nula [4].

Porém esse não é o caso com o presente método. Note que a Equação (13)

considera as equações de conservação de massa e quantidade de movimento, e envolve

a variação do campo de pressão durante o passo de tempo. Certamente, ela é usada para

calcular a variação do campo de pressão, sem levar em conta a escolha dos espaços de

interpolação empregados para aproximar velocidade e pressão.

No trabalho de Hughes et al. para Stokes [29] e nos trabalhos de De Sampaio para

Navier-Stokes [23] e [30], uma equação análoga a Equação (13) foi obtida adicionando

a equação de conservação de massa alguns termos extras de estabilização que surgem

naturalmente do método dos mínimos quadrados aplicado ao balanço de quantidade de

movimento. Aqui, no entanto, o mesmo termo de estabilização surge naturalmente a

partir do método de Taylor-Galerkin empregado para deduzir a Equação (13).

2.4 As equações para atualização da velocidade

Uma vez que o campo de pressão foi determinado, empregamos uma discretização

temporal de 2ª ordem ao balanço de quantidade de movimento para obter as equações de

atualização da velocidade.

As variáveis da Equação (2) são tomadas no nível de tempo 1/ 2n :

1/2 1/2 1/2 1/2

0

n n n n

a a abb

b b a

u u pu

t x x x (20)

onde

1/21

20

nn n

a a au u u

tt t

(21)

e

15

1/21

2110

2

nn n

n na a ab b b

b b b

u u uu u u tx x x

(22)

Logo, a discretização temporal da Equação (2) é dada por:

21 10

2 2

nna a a abb b a

b b b

u u uu u Q t

t x x x (23)

onde

1/2nn nn a ab

a bb b a

u pQ u

x x x (24)

No Apêndice A é mostrado que a Equação (23) possui precisão temporal de 2ª

ordem. Uma importante característica da discretização temporal de 2ª ordem acima é o

fato de que ela envolve uma linearização do termo convectivo. Portanto, a variação da

solução durante o passo de tempo pode ser calculada diretamente, isto é, sem recorrer a

um processo iterativo dentro do passo de tempo.

Como anteriormente mencionado, os subscritos n e 1n indicam o nível de

tempo e t denota o passo de tempo. A variação de velocidade durante o passo de

tempo t é representada por 1n n

a a au u u . O campo de pressão no nível de tempo

1/ 2n (o que significa que, neste caso, 1/ 2 ) é escrito como

1/2 1

2n np p p .

Considere a seguinte discretização espacial das variáveis do problema: ˆn na j aju N u

e a j aju N u . Novamente,

jN denota as funções de forma do elemento finito e as

variáveis com subscrito j são valores nodais. Usando o campo das variáveis

discretizadas, podemos escrever a seguinte expressão para o resíduo quadrado gerado

pela discretização da Equação (23):

16

ˆ â aS S

(25)

onde é uma parâmetro de escala para ser definido posteriormente e âS é o resíduo da

discretização da Equação (23) e é escrito como:

ˆ ˆ ˆ ˆ1 1ˆ ˆˆ ˆ2 2

nna a a ab

a b b ab b b

u u uS u u Q

t x x x (26)

Minimizando , dado pela Equação (25), em relação aos valores nodais livres

aiu , temos:

livreˆˆ ˆˆ 0

2 2

nn i b

i b a i c cib c

N ut tN u S N S u

t x x (27)

Os índices a , b e c representam coordenadas cartesianas.

Aqui, definimos t com o intuito de normalizar (e adimensionalizar) as

funções de peso da Equação (27). Perceba que a função de peso, presente no primeiro

termo da Equação (27), possui a mesma estrutura da função de peso encontrada no

método SUPG concebido por [29]. A função de peso remanescente, que atua sobre o

segundo termo, é específica do método apresentado aqui.

Tanto no presente método quanto em [23], [24] e [30], as estruturas das funções de

ponderação obtidas são diretamente dependentes do tempo de discretização adotado

originalmente. Assim, se ao invés da Equação (23), tivéssemos adotado uma

discretização temporal explícita da velocidade (com o campo de velocidade tomado no

tempo n ) obteríamos a função de ponderação do método clássico de Galerkin e nenhum

efeito de ponderação a montante – em inglês, upwind – por menor que fosse.

Note que a Equação (27) se refere somente aos resíduos da quantidade de

movimento no interior do domínio de análise. Neste ponto precisamos introduzir as

condições de contorno dadas pela Equação (4). Elas são adicionadas à Equação (27), na

forma de resíduos ponderados, como mostrado a seguir:

17

2 ˆˆ ˆˆ2 2

, 0

ta

nn i b

i b a i c

b c

i ab ab b a ci

N ut tN u S N S

x x

N p n t t u x livre (28)

Usando a identidade de Green obtemos um sistema simétrico de equações para

calcular os valores nodais para atualização da velocidade. Para problemas 2D, por

exemplo, temos:

uu uv uT

vuv vv

A A f

f A A x

y

u

u (29)

Considerando as contribuições de todos os elementos, a matriz uu

A é definida

como:

uu uuA A

nel

e

eA (30)

uu uuA A

ij

e e (31)

uuA =

2

ˆ ˆ1 1

2 2 2 2

ˆ 1

4 2

ije

e

je e x e i e x ei b j b

e b b

y je i ii j

Nt u t N t u tN u N u

t x x x x

u Nt N NN N

x x x

e

jN

y y

(32)

Da mesma forma, a matriz uv

A :

uv uvA A

nel

e

eA (33)

uv uvA A

ij

e e (34)

18

uvA =

ˆ ˆ1

2 2 2

ˆ ˆ1

2 2 2 2

ije

e e

e e x e i xi b j

b

y j y je e ij b i

b

t u t N uN u N

x x y

u N u Nt t NN u N

y x x y x

(35)

E, finalmente, o vetor u

f :

u uf f

nel

e

eA (36)

u uf fe e

i (37)

uf

ˆ ˆ1ˆ ˆ 2

2

ˆˆ

2

ie e e

e e

e x i i xi b

b

yi x i

eb

u N N uN u p p

x x x x

uN u N

y y y x

t Nu

ˆ ˆ ˆ 1ˆ

2

ˆ ˆ ˆ 1 1ˆ

2 2 2

e

ex

i x xi b

b b

y yei b i x x

b t

u u pN u p

x x x x

u ut pN u p N t tx x x

(38)

As matrizes vv

A e vu

A e o vetor v

f são definidas de forma análoga às matrizes

uuA e uv

A e ao vetor u

f . A única diferença é que os índices, variáveis e derivadas em

x e y são alternados.

É importante mencionar que a discretização temporal é realizada antes da

discretização espacial. E esta, utilizando elementos finitos de classe padrão 0C . Os

termos multiplicados por t na Equação (28) são responsáveis por controlar oscilações

espúrias em escoamentos com convecção dominante e por estabilizar o cálculo,

independentemente das restrições de Babuška-Brezzi na escolha dos espaços de

interpolação de pressão e velocidade. Em particular, o uso da mesma ordem de

interpolação para a pressão e para a velocidade adotada se torna possível pela adequada

escolha de t . É importante registrar que os termos de estabilização não são propostos

a priori, e sim, emergem naturalmente da minimização do resíduo quadrado da

19

discretização temporal do balanço da quantidade de movimento em relação aos graus de

liberdade de pressão e velocidade (valores nodais livres).

Note que a presente abordagem conduz a um sistema de equações parcialmente

acoplado, onde a solução dos graus de liberdade de pressão é obtida primeiro, e em

seguida a solução dos graus de liberdade das componentes de velocidade é realizada.

Quando comparado ao trabalho apresentado em [23], observamos que a precisão

temporal de 2ª ordem do presente método é alcançada ao custo de resolver as

componentes de velocidade acopladas, enquanto que em [23] as componentes de

velocidade são calculadas de maneira totalmente segregada.

Por outro lado, com um procedimento de solução eficiente, como o método dos

gradientes conjugados pré-condicionado, e estruturas de dados adequadas, um código

baseado no presente método é provavelmente tão bem escalável em aplicações paralelas

quanto um código baseado em um procedimento totalmente segregado, como

apresentado em [23]. Isso porque uma das grandes similaridades entre os métodos é o

cálculo do campo de pressão que é realizada de forma segregada. Como o cálculo da

pressão necessita de muito mais iterações, o sistema de equações de velocidade é

calculado de forma muito rápida em ambos os métodos. Os resultados referentes a esse

estudo serão mostrados mais a frente.

2.5 Passos de tempo locais e sincronização

Para elementos lineares, uma quantidade apropriada do efeito upwind é introduzida

no balanço de quantidade de movimento, Equação (28), se o passo de tempo for

definido como:

ee n

ht

u (39)

onde

20

Re 2coth

2 Reh

h

(40)

Na Equação (39), n n na au uu é o módulo da velocidade local e

eh o

tamanho característico do elemento. Em malhas 2D, eh pode ser obtido calculando a

raiz quadrada da área do elemento. Em malhas 3D ele pode ser estimado pela raiz

cúbica do volume do elemento. O número de Reynolds do elemento é

Reh e

hnu . A Equação (40) define , chamado de parâmetro ótimo de upwind.

Definido dessa forma conduz a soluções nodalmente exatas para problemas de

convecção e difusão, estacionários e unidimensionais [29].

O passo de tempo dado pela Equação (39) é adequado para seguir a evolução

temporal dos processos de convecção e difusão resolvíveis numa malha com elementos

de tamanho eh [30]. Para convecção forte (Re 1

h), obtemos

n

e et h u ,

enquanto que para difusão pura (Re 0h

), temos 2 6

e et h . Entretanto,

adotamos o modo alternativo de calcular o passo de tempo usado em [17], que é

equivalente a aproximar o parâmetro dado pela Equação (40). Assim, em vez de

utilizar as Equações (39) e (40), o passo de tempo é escolhido de acordo com o menor

valor das escalas de tempo características de convecção e difusão. Dessa forma o passo

de tempo é definido por mine

t t , onde min ,e e et tc td ,

n

e etc h u é a

escala de tempo característico de convecção e 2 6

e etd h é a escala de tempo

característico de difusão. O que é equivalente a aproximar o parâmetro na Equação

(40) por ' dado por:

Re 6 if Re 6 Re

1 if Re 6h h

hh

(41)

E então obter o passo de tempo de acordo com n

e et h u . O parâmetro ' ,

dado pela Equação (41), é precisamente uma aproximação assintótica para

21

introduzido em [29] para simplificar o cálculo da função de ponderação no método

SUPG.

Diante do que foi exposto acima está claro que a resolução espacial no presente

método, representada pelo tamanho do elemento eh , e a resolução temporal, associada

ao passo de tempo dado por n

e et h u , não são independentes. Em outras

palavras, o tamanho do passo de tempo diminui à medida que a malha é refinada. Como

resultado, não encontramos as instabilidades relacionadas a pequenos passos de tempo

investigadas por Bochev et al. [31] no caso do escoamento de Stokes. Na verdade, como

o passo de tempo local et é o próprio parâmetro de estabilização, a razão entre o

passo de tempo e o parâmetro de estabilização é mantida / 1t independente de

eh e

et , como apontado em [31].

Note que o passo de tempo calculado a partir de n

e et h u varia

espacialmente no domínio de acordo com a velocidade local, propriedades físicas e

tamanho da malha. Por conta do ótimo passo de tempo que varia com a posição,

precisamos recorrer a um esquema especial para sincronizar o avanço temporal da

computação. Para este fim, foi adotado o procedimento introduzido em [23]. Ele se

baseia na seleção do passo de tempo de sincronização *t , que deverá ser o mesmo

para todas as variáveis do escoamento e não variar no espaço, que é o conceito

convencional para o passo de tempo. O passo de tempo de sincronização é escolhido

para ser o menor passo de tempo do problema e é atualizado a cada passo calculando

* 0,999t t .

Considere que au e p sejam as variações das variáveis obtidas quando

utilizamos o passo de tempo local para resolver as Equações (13) e (29). Por outro lado,

permita que as variações das variáveis entre os instantes de tempo nt e

*tt n (tempo

de sincronização) sejam denotadas como ˆ *au e ˆ *p . Assim, mantendo a mesma

taxa de variação, obtemos as seguintes relações:

22

ˆ ˆ *

*

p p

t t (42)

ˆ ˆ *

*a au u

t t (43)

onde ˆ* *j j

p N p e ˆ * *a j aju N u . Claramente, ˆ *p e ˆ *

au são os valores

nodais das variações das variáveis do problema do instante de tempo nt ao

*tt n.

Na prática, o cálculo baseado no passo de tempo local e a fase de sincronização não

devem ser realizados separadamente. Certamente, a fase de sincronização, representada

nas Equações (42) e (43), podem ser inseridas nas Equações (13) e (29). Assim, a

solução sincronizada pode ser obtida diretamente resolvendo os seguintes sistemas

simétricos de equações:

p

M* * fp (44)

uu uv uT

vuv vv

A * A * * f

* f A * A * x

y

u

u (45)

O esquema de sincronização, introduzido por De Sampaio em [23] e [24], pode ser

estendido a problemas com transferência de energia e transporte de radionuclídeos [28],

por exemplo, onde são calculadas outras escalas de tempo de acordo com os fenômenos

físicos presentes.

O procedimento de solução é semi-segregado uma vez que a pressão é calculada

primeiro, e de forma independentemente em cada passo, mas a solução da componentes

de velocidade é acoplada, como mostra a Equação (29). Depois de cada passo as

variáveis do problema são atualizadas e o procedimento de cálculo continua até o tempo

de análise final especificado ser alcançado.

As Equações (44) e (45) geram matrizes simétricas positivas definidas. O código

computacional onde a formulação foi inicialmente implementada utiliza um método de

gradientes conjugados pré-condicionado como solucionador linear e estrutura de dados

por elemento para resolver os sistemas de equações. Tal código foi desenvolvido para o

estudo de problemas 2D, onde triângulos lineares são usados para interpolar todas as

23

variáveis independentes. O código utiliza ainda refinamento adaptativo baseado em

estimativa de erro [32], e programação paralela otimizada para alto desempenho em

sistemas de computação paralela com memória distribuída.

É importante considerar que o presente método foi desenvolvido para soluções

muito precisas em escoamentos transientes, para os quais soluções não estacionárias são

esperadas. É por isso que o passo de tempo de sincronização é o mínimo passo de tempo

local na malha. Por outro lado, no caso de resolver problemas estacionários, um método

mais eficiente pode ser desenvolvido adotando o procedimento descrito em [26]. Lá, o

passo de tempo de interpolação é escolhido para ser muito próximo ao maior passo de

tempo na malha e as variáveis que atingem o instante de tempo de interpolação são

temporariamente removidos da análise, acelerando o cálculo. Mais detalhes podem ser

encontrados em [26].

Para a simulação de problemas turbulentos, a presente formulação de elementos

finitos introduz naturalmente uma modelagem para as escalas não resolvíveis. Elas

podem ser comparadas ao método de Simulação de Grandes Escalas – em inglês, Large

Egde Simulation (LES) – e são equivalentes ao método de discretização de Galerkin das

equações filtradas no espaço, onde é aplicado um modelo sub-malha particular,

proporcional ao resíduo de discretização. Desta forma, o modelo implícito é seletivo,

pois atua preponderantemente na escala sub-malha e possui efeito desprezível nas

escalas resolvíveis [28]. Estes tópicos não serão abordados nesta tese onde o foco é

aumento da precisão temporal. A relação entre modelos sub-malha e formulações

estabilizadas de elementos finitos foi elucidada por Hughes, quando propôs o

Variational MultiScale Method em [33].

No Apêndice B é detalhada uma formulação de elementos finitos com precisão

temporal de 1ª ordem, utilizada para comparação com a formulação presente em alguns

problemas numéricos apresentados adiante [23].

24

Capítulo 3

Exemplos Numéricos Bidimensionais

O método proposto foi usado na análise de quatro problemas bem conhecidos em

CFD: o escoamento induzido no interior de uma cavidade pelo movimento de uma

tampa, a fase inicial transiente de um escoamento com uma expansão abrupta, a

simulação do problema transiente do desprendimento periódico de vórtices em

escoamento cruzado sobre um cilindro circular e o problema do escoamento de Kim-

Moin. Em todos os exemplos bidimensionais apresentados a tolerância para o critério de

parada do solucionador linear e de 10-3

. O critério de parada baseia-se no resíduo

relativo.

3.1 Escoamento induzido no interior de uma cavidade pelo

movimento de uma tampa

Inicialmente parece estranho mostrar um problema estacionário para avaliar a

desempenho de um método concebido para soluções precisas de escoamentos

transientes. No entanto, como mostrado no Capítulo 2, a estrutura da função de

ponderação no presente método é deduzida da discretização temporal de 2ª ordem dada

pela Equação (23). Assim, as funções de ponderação são exclusivas do presente método

e nunca foram testadas antes. É por isso que o método é utilizado em um exemplo

estacionário, onde podemos comparar as suas predições com resultados publicados

anteriormente. Em particular, este exemplo demonstra que os termos de estabilização

25

que são inerentes à formulação possibilitam calcular corretamente os campos de

velocidade e pressão, apesar da escolha do mesmo espaço de interpolação.

São mostrados os resultados da simulação do problema clássico do escoamento

induzido no interior de uma cavidade pelo movimento de uma tampa para 1 000L

Re .

e 5 000L

Re . . O número de Reynolds é definido como 0L

Re u L / , onde L é

a medida do lado da cavidade quadrada e 0u é a velocidade de referência do

escoamento.

A Figura 1 mostra a malha de base usada neste exemplo.

Figura 1. Malha de base utilizada para 1 000L

Re . e 5 000L

Re . , contendo

2.601 nós e 5.000 elementos.

São prescritas as velocidades 1xu e 0 0

yu , no topo do domínio para simular o

movimento da tampa da cavidade para a direita, sentido positivo do eixo x . As laterais

e a base da cavidade são consideradas paredes sólidas, ou seja, nelas são prescritas as

velocidades 0 0xu , e 0 0

yu , .

Os resultados para 1 000L

Re . e 5 000L

Re . são mostrados na Figura 2.

x

y

26

(a) (b)

(c) (d)

Figura 2. Escoamento induzido numa cavidade. Resultados obtidos com o método

proposto: Iso-linhas de pressão (a) e vorticidade (b) para 1 000L

Re . ; Iso-linhas de

pressão (c) e vorticidade (d) para 5 000L

Re . .

, uma vez que o

método introduz uma quantidade adequada de estabilização.

As Figuras 3 e 4 apresentam as componentes de velocidades na linha de centro da

cavidade.

27

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

-0,6

-0,4

-0,2

0,0

0,2

0,4

0,6

0,8

1,0

Ghia et. al.

Presente

Co

mp

on

en

te u

x

Coordenada y

Componente ux de velocidade na linha de centro vertical para Re = 1.000

(a)

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

-0,6

-0,4

-0,2

0,0

0,2

0,4

0,6

0,8

1,0

Ghia et. al.

Presente

Co

mp

on

en

te u

y

Coordenada x

Componente uy de velocidade na linha de centro horizontal para Re = 1.000

(b)


Re . .

Componente xu (a) e componente

yu (b).

28

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

-0,6

-0,4

-0,2

0,0

0,2

0,4

0,6

0,8

1,0


Ghia et. al.

Presente

Co

mp

on

en

te u

x

Coordenada y

(a)

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

-0,6

-0,4

-0,2

0,0

0,2

0,4

0,6

0,8

1,0

Ghia et. al.

Presente

Co

mp

on

en

te u

y

Coordenada x

Componente uy de velocidade na linha de centro horizontal para Re = 5.000

(b)


Re . .


yu (b).

29

Os resultados obtidos, apresentados nas Figuras 2, 3 e 4 para 1 000L

Re . e

5 000L

Re . , são tão boas quanto aos resultados da solução de Guia et al. apresentado

em [34], utilizado como referência para esta análise. Significa, portanto, que a presente

formulação se mostra apropriada para a simulação de problemas de mecânica de fluidos.

3.2 Escoamento com expansão abrupta

Aqui foi considerada a simulação dos estágios iniciais de um escoamento transiente

com expansão abrupta. O desempenho do presente método de precisão temporal de 2ª

ordem foi comparado ao método de precisão temporal de 1ª ordem apresentado

anteriormente por De Sampaio [23] na análise do transiente rápido.

Figura 5. Geometria e condições de contorno do problema de escoamento com

expansão abrupta.

A Figura 5 mostra a geometria e as condições de contorno para o problema do

escoamento com expansão abrupta. A velocidade na entrada é uniforme com 1,0xu e

0,0yu . O fluido no interior do domínio é considerado parado no instante 0t .

Para comparar as soluções obtidas com os métodos de precisão temporal de 1ª e 2ª

ordem, foi adotado o ponto de monitoração da velocidade no ponto 0,125; 0x y ,

no instante 00 1 /t , H u , onde

0u é a velocidade de referência do escoamento. Assim,

as simulações transientes realizadas abrangem os estágios iniciais do transiente desde o

início, quando o fluido estava em repouso.

Ponto Monitorado

ux = uy = 0

ux = uy = 0

ux = uy = 0 p = 0

X

Y

H/2

ux = 1

uy = 0 (0,125 ; 0)

H 4H

H

30

A malha de base de discretização espacial mais grosseira utilizada neste exemplo é

apresentada na Figura 6. Ela possui 1.249 nós e 2.304 elementos triangulares. O

tamanho característico dos elementos dessa malha é 0,125eh , onde 2H .

Figura 6. Malha de base estruturada utilizada no escoamento com expansão

abrupta com 1.249 nós e 2.304 elementos triangulares.

Uma solução discretizada, obtida com o esquema de precisão de 2ª ordem e uma

malha estruturada com 74.497 nós e 147.456 elementos, foi utilizada com solução de

referência (benchmark). A avaliação de desempenho foi feita sobre ambos os métodos,

em discretizações espaciais mais grosseiras, comparando suas predições com a solução

de referência. É importante lembrar que em ambos os métodos a discretização especial e

temporal não são independentes, mas relacionadas de acordo com une et h .

Assim, quando uma malha é refinada, o passo de tempo de interpolação é também

refinado. O tamanho característico do elemento na malha de referência é 0,015625h

e o passos de tempo de interpolação na simulação são da ordem de 2E-5.

O erro relativo da velocidade obtido nas diversas rodadas é calculado de acordo

com be u , onde r be u u . Note que ru é a velocidade no ponto monitorado

obtido numa rodada particular e bu é a velocidade correspondente a solução de

referência. As velocidades no ponto e instante monitorados são 0,756340xbu e

0,995090ybu . Os resultados são mostrados na Figura 7.

31

Figura 7. Escoamento com expansão abrupta. Caso de referência: Malha

estruturada com 74.497 nós e 147.456 elementos. Resultado do campo de velocidade

obtido com o método proposto para 100H

Re em 0

0 1 /t , H u .

As simulações foram realizadas utilizando versões seqüenciais dos códigos com as

formulações de 1ª e 2ª ordem no tempo, executadas em um computador pessoal com

processador AMD Athlon 64 X2 Dual Core 4400+. Três níveis de refinamento foram

considerados no estudo. Os níveis de refinamento são especificados pelos tamanhos

característicos mínimo e máximo escolhidos para o refinamento adaptativo de cada

malha. Quando os tamanhos mínimo e máximo são iguais, o refinamento adaptativo é

automaticamente desativado e a malha de base não sofre alterações.

O método de 2ª ordem, embora mais preciso, necessita de mais tempo de

processamento que o método de 1ª ordem para um mesmo caso. Para uma comparação

justa entre os dois métodos, foram confrontados os seus desempenhos mostrando o erro

relativo da velocidade contra o tempo gasto de CPU. A Tabela 1 apresenta os tamanhos

mínimo e máximo dos elementos empregados nas simulações, a velocidade no ponto

monitorado no instante 0

0 1 /t , H u , o erro relativo da velocidade e o tempo de CPU

correspondente.

32

Tabela 1. Resumo dos dados usados para comparação entre os métodos de precisão

temporal de 1ª e 2ª ordem.

1ª Ordem

hmin hmax ux uy Erro Tempo de CPU (s)

0,125000 0,125000 0,245690 -0,207000 0,751314 3,28

0,031250 0,125000 0,381740 -0,560610 0,458973 132,36

0,015625 0,062500 0,640660 -0,856500 0,144431 2.966,78

2ª Ordem

hmin hmax ux uy Erro Tempo de CPU (s)

0,125000 0,125000 0,420740 -0,426830 0,528009 12,79

0,031250 0,125000 0,740120 -0,958580 0,031963 287,55

0,015625 0,062500 0,760330 -1,016110 0,017118 5.817,79

A Figura 8 representa o erro relativo da velocidade e o tempo de CPU,

demonstrando o desempenho superior do método de 2ª ordem no exemplo de uma

transiente rápido.

Figura 8. Erro relativo da velocidade e tempo de CPU consumido para os métodos

de 1ª e 2ª ordem. Escoamento com expansão abrupta em diferentes níveis de

refinamento.

0,00

0,10

0,20

0,30

0,40

0,50

0,60

0,70

0,80

1 10 100 1000 10000

Erro

Re

lati

vo d

a V

elo

cid

ade

Tempo de CPU (s)

1ª Ordem

2ª Ordem

33

Através dos dados da Tabela 1 e do gráfico apresentado na Figura 8 podemos

verificar que a medida com que a malha utilizada é mais refinada, o erro relativo da

velocidade diminui e o tempo de CPU aumenta. Esse comportamento é esperado para

ambos os métodos.

Verificamos ainda que com o método de 2ª ordem é possível atingir uma solução

mais precisa, com uma malha mais grosseira e utilizando menos tempo de CPU. Mesmo

com a malha mais refinada, o método de 1ª ordem não fornece uma solução tão precisa

quanto à obtida pelo o método de 2ª ordem com a malha intermediária. A malha

intermediária é a malha de base refinada adaptativamente. A Figura 9 mostra a malha

intermediária após a simulação.

Figura 9. Escoamento com expansão abrupta. Refinamento adaptativo da malha

intermediária ( min 0,03125h e máx 0,125h ) em 0

0 01 /t , H u .

3.3 Desprendimento periódico de vórtices em escoamento cruzado

sobre um cilindro circular

Nesta seção são apresentados os resultados da simulação do problema de

desprendimento periódico de vórtices. O fenômeno é observado quando o escoamento

de um fluido incompressível é perturbado pela presença de um cilindro circular em seu

caminho. A geometria e as condições de contorno do problema são apresentadas na

Figura 10.

34

Figura 10. Geometria e condições de contorno do problema de escoamento cruzado

sobre um cilindro circular.

O número de Reynolds é dado por 0d

Re u d / , onde d é o diâmetro do

cilindro. Novamente, 0u é a velocidade de referência do escoamento.

Simulações para 100d

Re , 125d

Re e 150d

Re foram realizadas. As forças

de arrasto (drag) DF e sustentação (lift)

LF foram calculadas ao longo do

processamento. A freqüência f do desprendimento de vórtices é a freqüência de

oscilação da força de sustentação. Essa grandeza é analisada pelo número adimensional

de Strouhal dado por 0

/St fd u .

Os coeficientes de arrasto e sustentação são dados por 2

02 /

D DC F u d e

2

02 /

L LC F u d , respectivamente. As forças de arrasto e sustentação são calculadas a

partir dos dados do escoamento de acordo com as seguintes equações:

2

c

yx xD y x x

uu uF n n pn d

y x x (46)

e

6 d 12 d

ux = uy = 0

ux = uy = 0

p = 0

ux = 1

uy = 0

ux = uy = 0

d x

y

35

2

c

yx xL x y y

uu uF n n pn d

y x x (47)

onde an é o vetor normal à superfície do cilindro

c .

Para realização do estudo foi utilizada uma malha de base não estruturada com

4.083 nós e 7.652 elementos triangulares, bastante refinada em torno do cilindro.

A Figura 11 mostra os campos de pressão e velocidade do problema para

100d

Re .

36

t=4d/u0

(a)

(b)

(c)

t=12d/u0

(d)

(e)

t=128d/u0

(f)

(g)

t=188d/u0

(h)

Figura 11. Escoamento cruzado sobre um cilindro circular para 100d

Re .

Campos de velocidade (à esquerda) e pressão (à direita) nos instantes 0

4 /t d u ((a) e

(b)), 012 /t d u ((c) e (d)), 0

128 /t d u ((e) e (f)) e 0188 /t d u ((g) e (h)).

37

A evolução temporal dos coeficientes de arrasto e sustentação para os casos

simulados são apresentados nas Figuras 12, 13 e 14. O tempo de simulação nos 3 casos

é de 0

200 /t d u .

0 50 100 150 200-0,4

0,0

0,4

0,8

1,2

1,6 Re = 100

Co

efi

cie

nte

s

Tempo

CD

CL

Figura 12. Evolução temporal dos coeficientes de arrasto (CD) e sustentação (CL)

para o escoamento cruzado sobre um cilindro circular para.

O comportamento observado para os coeficientes de arrasto de sustentação é, de

fato, o comportamento esperado. A comparação entre os coeficientes calculados na

simulação e os dados presentes na literatura é feita mais adiante.

38

0 50 100 150 200

-0,4

0,0

0,4

0,8

1,2

1,6

CD

CL

Re = 125

Co

efic

ien

tes

Tempo


para o escoamento cruzado sobre um cilindro circular para 125d

Re .

0 50 100 150 200

-0,4

0,0

0,4

0,8

1,2

1,6

CD

CL

Re = 150

Co

efic

ien

tes

Tempo


para o escoamento cruzado sobre um cilindro circular para 150d

Re .

39

Para uma melhor representação do escoamento, foi utilizado o refinamento

adaptativo do código. A Figura 15 mostra a malha refinada em 0

200 /t d u para

100d

Re . O tamanho característico mínimo e máximo do elemento escolhidos para o

refinamento são 0 02, d e 0 40, d , respectivamente. Esses valores correspondem aos

tamanhos mínimo e máximo dos elementos encontrados na malha de base.

Figura 15. Escoamento cruzado sobre cilindro. Malha refinada em 0

200 /t d u

para 100d

Re . Tamanho característico mínimo de 0 02, d e máximo de 0 40, d .

O refinamento adaptativo é guiado pelos gradientes de velocidades mais intensos,

baseado em uma estimativa de erro [32]. Daí o refinamento adaptativo também mais

intensivo na região da esteira, onde há o desprendimento de vórtices.

As Equações (48) e (49) representam uma correlação experimental para o cálculo

do valor r.m.s. do coeficiente de sustentação [35]. Esta referência é utilizada como

benchmark para este exemplo:

47

47d

Re

(48)

2

30 90LC (r.m.s)

(49)

A comparação entre os resultados obtidos na simulação para o valor r.m.s do

coeficiente de sustentação e o calculado de acordo com as Equações (48) e (49) e para o

40

coeficiente médio de arrasto e os dados experimentais encontrados em [36] são

mostrados na Figura 16.

Figura 16. Comparação entre os resultados do presente estudo com os dados

obtidos em [35] e [36].

O número de Strouhal é comparado com os dados obtidos através da correlação

experimental de Strouhal-Reynolds [37]:

3 32650 1816 0 00016

dd

.St . . Re

Re

(50)

A Figura 17 mostra os resultados obtidos.

0,0

0,2

0,4

0,6

0,8

1,0

1,2

1,4

1,6

100 110 120 130 140 150

Co

efi

cie

nte

s

Número de Reynolds

CD mean exp

CD mean simul

CL rms exp

CL rms simul

41

Figura 17. Comparação entre as predições obtidas nas simulações para o número de

Strouhal e os resultados da correlação Strouhal-Reynolds [37].

Como mostrado nas Figuras 16 e 17, os resultados conseguidos no presente estudo

estão em ótima concordância com os dados encontrados na literatura.

3.4 Escoamento de Kim-Moin

O escoamento de Kim-Moin [39] é um problema bidimensional e possui solução

analítica para escoamentos viscosos incompressíveis numa cavidade quadrada de

domínio [0,1]x[0,1]. O escoamento de Kim-Moin é um problema fortemente

indicado para verificação da taxa de convergência de um determinado método. Além

disso, ele pode ser utilizado para analisar a evolução temporal do erro de uma grandeza

calculada no escoamento. Algumas análises utilizando o escoamento de Kim-Moin

podem ser encontradas em [40], [41] e [42].

O erro para o escoamento de Kim-Moin é baseado na norma da diferença entre a

solução aproximada, a solução obtida na simulação, e a solução exata.

A solução exata, utilizada para a comparação com os resultados numéricos, é dada

por:

0,15

0,16

0,17

0,18

0,19

0,20

100 110 120 130 140 150

Nú

me

ro d

e S

tro

uh

al

Número de Reynolds

Strouhal simul

Strouhal exp

42

2 2( , , ) cos( )sen( )exp( 2 )xu x y t a x a y a t (51)

2 2( , , ) sen( )cos( )exp( 2 )yu x y t a x a y a t (52)

2 21( , , ) (cos(2 ) cos(2 ))exp( 4 )

4p x y t a x a y a t (53)

A solução apresentada nas Equações (51), (52) e (53) é estacionária no espaço e

diminui monotônicamente no tempo. As componentes de velocidade e a pressão são

representadas por xu , yu e p , respectivamente. x , y e t são as coordenadas espaciais

e o tempo. Para o estudo apresentado, temos o parâmetro 2a e a viscosidade

cinemática do fluido 0.01.

A solução exata é um produto entre as soluções temporal e espacial [39]. Isso

possibilita o redimensionamento do erro a fim de remover o decaimento temporal deste,

característico da solução. Dessa forma, o erro absoluto espacial pode ser definido como:

2 2

2: exp(2 )h

uerr u u a t (54)

2 2

2: exp(2 )h

verr v v a t (55)

2 2

2: exp(4 )h

perr p p a t (56)

onde o sobrescrito h indica a solução numérica e 2

. representa a norma usual 2L .

Foram realizados testes em malhas estruturadas de 20 20 , 32 32 , 44 44 ,

56 56 e 68 68 células quadradas, onde cada uma delas foi dividida ao meio pela sua

diagonal, dando origem a elementos triangulares. A seqüência das divisões para as

malhas não seguem nenhum esquema especial, tendo sido adotada de forma aleatória.

A Figura 18 apresenta as malhas de 20 20 elementos e de 68 68 células. Para

este exemplo elas representam a malha mais grosseira e a mais refinada,

respectivamente. Não foi utilizado o recurso de refinamento adaptativo em nenhum dos

casos apresentados.

43

(a) (b)

Figura 18. Malha estruturada triangular utilizada para o estudo do Escoamento de

Kim-Moin. (a) 20 20 , 441 nós e 800 elementos e (b) 68 68 , 4.761 nós e 9.248

elementos.

A Figura 19 apresenta os campos de velocidade e pressão em 1.0 t s obtidos com

o método de precisão temporal de 2ª ordem na malha 68 68 numa nova demonstração

da adequação do método.

As condições iniciais do problema são prescritas a partir da solução analítica do

problema, representada pelas Equações (51), (52) e (53).

(a) (b)

Figura 19. Escoamento de Kim-Moin. Campo de velocidade (a) e pressão (b) no

instante 1 0 t . s .

44

Perceba que os campo de velocidade e pressão são devidamente representados, sem

a evidência de oscilações espúrias, destacando mais uma vez a estabilização da solução

pela presente formulação. Note que é utilizada a mesma ordem de interpolação para

velocidade e pressão, como mencionado anteriormente.

A Tabela 2 mostra os erros definidos pelas Equações (54), (55) e (56) para as cinco

malhas e obtidos a partir das simulações realizadas com os métodos de precisão

temporal de 1ª e 2ª ordem em 1 0 t . s e 2 0 t . s .

Tabela 2. Erro das componentes de velocidade e pressão em 1 0 t . s e

2 0 t . s para as cinco malhas.

Malha t = 1,0 s 1ª Ordem 2ª Ordem

Tamanho helem dt passos errou errov errop errou errov errop

20 x 20 0,035355 2,08E-04 4808 4,08E-04 3,58E-04 6,30E-04 4,06E-04 3,67E-04 3,50E-04

32 x 32 0,022097 8,13E-05 12300 9,34E-05 8,89E-05 2,42E-04 9,30E-05 8,92E-05 9,39E-05

44 x 44 0,016071 4,30E-05 23256 3,50E-05 3,44E-05 1,27E-04 3,50E-05 3,42E-05 3,81E-05

56 x 56 0,012627 2,65E-05 37679 1,67E-05 1,68E-05 7,86E-05 1,68E-05 1,66E-05 1,92E-05

68 x 68 0,010399 1,80E-05 55556 9,30E-06 9,40E-06 5,33E-05 9,30E-06 9,20E-06 1,10E-05

Malha t = 2,0 s 1ª Ordem 2ª Ordem

Tamanho helem dt passos errou errov errop errou errov errop

20 x 20 0,035355 2,08E-04 9615 5,50E-04 4,86E-04 7,25E-04 4,78E-04 4,49E-04 5,12E-04

32 x 32 0,022097 8,13E-05 24600 1,21E-04 1,18E-04 2,54E-04 1,12E-04 1,09E-04 1,34E-04

44 x 44 0,016071 4,30E-05 46512 4,47E-05 4,50E-05 1,31E-04 4,27E-05 4,19E-05 5,47E-05

56 x 56 0,012627 2,65E-05 75358 2,13E-05 2,18E-05 7,97E-05 2,07E-05 2,03E-05 2,78E-05

68 x 68 0,010399 1,80E-05 111111 1,17E-05 1,21E-05 5,38E-05 1,15E-05 1,13E-05 1,62E-05

Como mencionado anteriormente, em ambos os métodos o tamanho característico

do elemento e o passo de tempo de sincronização não são independentes. Pode-se

verificar que à medida que a malha é refinada, o erro, num determinado instante, vai

ficando cada vez menor. Por outro lado, o número de passos de tempo e,

conseqüentemente o número de iterações, aumenta consideravelmente até o instante de

tempo escolhido para análise. Apesar do caráter altamente não linear do problema de

Kim-Moin, o aumento do número de interações não interfere de forma relevante na

convergência espacial do método.

45

O erro da discretização espacial pode ser sintetizado através de uma efetiva taxa de

convergência [41] obtida pela Equação (57):

: iP

i ierr h (57)

onde iP é a taxa de convergência do erro espacial para cada variável i .

A Equação (57) pode ser linearizada da seguinte forma:

log( ) log( ) log( )i i ierr P h (58)

Dessa forma, iP é facilmente identificado como o coeficiente angular no gráfico

log( ) x log( )ierr h .

A Tabela 3 apresenta as taxas de convergência para cada variável nos instantes

1 0 t , s e 2 0 t , s .

Tabela 3. Taxa de convergência do erro espacial em 1 0 t , s e 2 0 t , s

considerando as cinco malhas do exemplo.

Taxa de Convergência

Tempo Método Pu Pv Pp

t = 1,0 s 1ª ordem 2,475 2,483 2,171

2ª ordem 2,476 2,484 2,461

t = 2,0 s 1ª ordem 2,414 2,417 2,161

2ª ordem 2,428 2,438 2,370

Ambos os métodos apresentam excelente taxa de convergência espacial. Vale

lembrar que os dois apresentam a mesma ordem de discretização no espaço se

diferenciado essencialmente na precisão temporal.

À medida que o tempo vai passando a tendência é que a taxa de convergência vá

diminuindo. Como podemos verificar na Tabela 3, a taxa de convergência é

ligeiramente reduzida de 1 0 t , s para 2 0 t , s . Percebe-se ainda que o esquema de

46

precisão temporal de 2ª ordem apresenta taxa de convergência em torno de 10% maior

para a pressão, o que não acontece com as componentes de velocidade.

A fim de investigar a evolução temporal do erro, foram realizadas novas simulações

com as malhas originais até o instante 8 0 t , s . A Figura 20 mostra os resultados da

evolução temporal dos erros da velocidade obtidos com os dois métodos.

Figura 20. Evolução temporal do erro da componente de velocidade xu obtida com

os métodos de precisão temporal de 1ª e 2ª ordem para a 4 malhas menos refinadas.

Com ambos os métodos o erro da velocidade se torna cada vez menos acentuado de

acordo com o refinamento da malha utilizada na simulação. O erro da velocidade obtido

com o método de 1ª ordem é sempre crescente, enquanto que o método de 2ª ordem

0,0E+00

5,0E-04

1,0E-03

1,5E-03

2,0E-03

2,5E-03

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

Erro

r u

Time (s)

20x20 - 2ª ordem

20x20 - 1ª ordem

32x32 - 2ª ordem

32x32 - 1ª ordem

44x44 - 2ª ordem

44x44 - 1ª ordem

56x56 - 2ª ordem

56x56 - 1ª ordem

47

impede o crescimento do erro da velocidade a partir de um determinado instante,

tornando o erro praticamente constante ao longo do tempo.

Ainda assim, os gráficos para uma mesma malha apresentam resultados muito

próximos no intervalo que se estende de 0 0 t , s até o tempo em torno de 1 8 t , s .

Este resultado é observado para a malha de 20x20x2 elementos. O intervalo de tempo

inicial em que os gráficos são semelhantes para ambos os métodos aumenta à medida

que a malha é mais refinada.

O estudo do escoamento de Kim-Moin apresentado na tese segue a referências

citadas no início desta seção [40], [41] e [42]. Na determinação do erro espacial, a

norma da diferença entre as soluções exata e aproximada é multiplicada por uma função

exponencial crescente, que tende ao infinito com o tempo. Por outro lado, à medida que

o problema avança no tempo, é esperado que a solução exata e a solução aproximada se

tornem mais próximas, fazendo com que a diferença entre elas também diminua. E

como a solução exata decai monotônicamente, tendendo à zero com o tempo, a

diferença entre as soluções exata e aproximada também tenderia à zero.

A questão é que a precisão da solução aproximada é limitada pela tolerância do

solucionador linear. Com o avanço do tempo, a diferença entre as soluções tende à

tolerância do solucionador linear e o erro calculado passa a ser a tolerância do

solucionador linear vezes uma exponencial crescente. Em outras palavras, os gráficos de

todos os erros apresentados na Figura 20 devem crescer. Esse comportamento só não é

observado devido ao tempo de simulação ( 8 0 t , s ) que, no presente estudo, é

insuficiente para tal. O crescimento monotônico do erro, conforme esperado, pode ser

visto em [42].

A principal suspeita para o amortecimento apresentado no método de 2ª ordem fica

a cargo do procedimento de solução adotado, onde as velocidades são resolvidas de

forma acoplada. Nesse caso, parece que há uma maior exigência do solucionador linear,

o que retarda o crescimento do gráfico. Esse assunto merece uma investigação mais

detalhada

48

Capítulo 4

Implementação Tridimensional

Nos capítulos anteriores foram apresentados o modelo físico para escoamentos

quase-incompressíveis, a formulação de elementos finitos com precisão temporal de 2ª

ordem com estabilização inerente ao método de discretização utilizado e os resultados

obtidos para problemas de mecânica de fluidos bidimensionais.

O presente método foi implementado em um código 3D, de onde herdou algumas

das suas características originais, com objetivo de ser disponibilizado como uma

alternativa na solução de problemas transientes com boa precisão atrelada a um bom

desempenho computacional.

Neste capítulo é apresentada uma breve descrição do código 3D utilizado como

estrutura para o presente método, apontando suas principais características e diferenças

entre a formulação original e a formulação proposta.

4.1 Características básicas do EdgeCFD®

O EdgeCFD®

é um software desenvolvido pela COPPE / UFRJ para simulação de

problemas em 3D de mecânica de fluidos com base em elementos finitos estabilizados e

estrutura de dados por aresta. Além de estrutura de dados por aresta, o software conta

com paralelismo híbrido com MPI e OpenMP, ordenação de dados, acoplamento

escoamento-transporte, dentre outras facilidades [47], [48], [49], [50], [51], [52], [53] e

[54].

49

O método original do EdgeCFD® utiliza estabilização PSPG – do inglês, Pressure-

Stabilizing/Petrov-Galerkin – para contornar o problema da escolha do mesmo espaço

de interpolação para pressão e velocidade, restrição de Babuška-Brezzi, [14], [19]. Para

outra fonte de problema, relacionada a oscilações espúrias verificadas em problemas de

convecção dominante e discretização com o método clássico de Galerkin, o EdgeCFD®

utiliza estabilização SUPG – do inglês, Streamline-Upwind/Petrov-Galerkin [29], [19].

Para problemas de escoamentos com número de Reynolds muito alto, o software dispõe

ainda da estabilização LSIC – do inglês, Least Squares Incompressibility Constraint

[43], para atenuar possíveis flutuações da solução.

A formulação de elementos finitos que este software utiliza leva as equações

governantes do movimento do fluido à um sistema acoplado de 4 graus de liberdade (3

de velocidade e 1 de pressão). As matrizes resultantes são não simétricas e o método

utilizado na solução é o GMRES – do inglês, Generalized Minimum RESidual [55],

[56]. O método de Newton-Inexato é utilizado para a solução do problema não-linear a

cada passo de tempo.

É adotado um esquema de marcha no tempo implícito realizado através do

algoritmo preditor-multicorretor de 2ª ordem. Alternativamente ao passo de tempo fixo,

o EdgeCFD® dispõe de um passo de tempo adaptativo onde o controle do passo utiliza

um controlador do tipo PID (Proporcional-Integral-Derivativo) [49].

A técnica de DD (Desativação Dinâmica) é opcional em problemas de

conhecimento prévio das regiões de variação e “congelamento” da solução, bem como o

algoritmo de conservação de massa, em problemas onde há “perda” de fluido [57].

Além de problemas de transferência de calor, cujo acoplamento com as equações de

Navier–Stokes é feito através de aproximação Bousinesq, o EdgeCFD® permite a

simulação de problemas de superfície livre através do método VOF – em inglês, Volume

of Fluid. A formulação de elementos finitos, aplicada às equações de transporte, utiliza

estabilização SUPG e operador de captura de descontinuidades YZβ [44]. Em condições

específicas, o operador de captura de descontinuidades YZβ reproduz o CAU – do

inglês, Consistent Approximated Upwind [45]. O método de estabilização CAU é um

operador de captura de descontinuidades e atua em regiões de intensos gradientes em

direções diferentes das linhas de corrente. A solução das equações de transporte é

50

segregada, ou seja, os graus de liberdade de temperatura e da função marcadora do VOF

são resolvidos separadamente dos graus de liberdade de pressão e velocidade.

O modelo de turbulência implementado no software é baseado na simulação das

grandes escalas através do modelo de turbulência clássico de Smagorinsky [46],

Smagorinsky dinâmico e do método RB-VMS, do inglês Residual–Based Variational

Multi-Scale [58].

O EdgeCFD®

reproduz uma gama extensa de problemas de mecânica de fluidos,

puramente mecânicos ou acoplados com transferência de calor e superfície livre,

escoamentos internos e externos, escoamento de fluidos newtonianos e não-

newtonianos, correntes gravitacionais, dispersão de gases, escoamento bifásico, etc.

com qualidade e desempenho estupendos registradas em vários trabalhos científicos.

4.2 Diferenças entre o método proposto e o método presente no

EdgeCFD®

A formulação de elementos finitos implementada no EdgeCFD® é realizada da

seguinte forma: Primeiramente é realizada a discretização espacial das equações

governantes. A discretização espacial utiliza elementos tetraédricos lineares. No caso

das equações de Navier-Stokes, o modelo de turbulência é incluído ainda no modelo

físico e os termos de estabilização SUPG e PSPG são propostos a priori. Esses termos

permitem a simulação de problemas de convecção dominante e a escolha dos mesmos

espaços de interpolação para pressão e velocidade, respectivamente.

Esse método dá origem a um sistema de equações diferenciais ordinárias não

lineares, chamadas de equações semi-discretas. Em função da semi-discretização das

equações governantes, a formulação depende de um esquema de integração temporal

que completa a discretização. Para tal, o EdgeCFD® utiliza o método preditor-

multicorretor, que pertence a família dos métodos trapezoidais generalizados. Após a

discretização das equações de Navier-Stokes, a formulação apresenta um sistema

acoplado, com 4 graus de liberdade, de equações algébricas não lineares, que precisam

ser resolvidas a cada passo de tempo. Para a solução do sistema não linear, é a realizada

a linearização através do método de Newton-Inexato, que é um solucionador avançado

51

implementado. A cada problema não linear, são resolvidos vários problemas lineares.

Como as matrizes não são simétricas, o EdgeCFD® utiliza o método GMRES com pré-

condicionamento bloco diagonal, que é um método baseado em subespaços de Krylov e

que necessita de avaliação do produto matriz vetor, sendo esta a parte mais custosa do

algoritmo. Porém, o elevado custo do GMRES é compensado pela estrutura de dados

que o software utiliza, que produz ganhos significativos em utilização de memória,

endereçamento indireto e operações de ponto flutuante.

A estratégia de discretização da nova formulação de elementos finitos proposta foi

apresentada anteriormente. Entretanto, as diferenças básicas entre a formulação

implementada no EdgeCFD® e a nova formulação proposta começam na seqüência

adotada para da discretização. Na presente formulação, é realizada primeiramente a

discretização temporal e, em seguida, a discretização espacial das equações

governantes. A discretização espacial é realizada através da minimização do resíduo

quadrado das equações discretizadas no tempo em relação aos graus de liberdade livres.

No caso das equações de Navier-Stokes, o acoplamento entre as equações de massa e

quantidade de movimento é feito de forma consistente, de modo que o mesmo espaço de

interpolação possa ser usado para pressão e velocidade.

Os termos de estabilização para problemas de convecção dominante, de estrutura

similar a estrutura SUPG, surgem naturalmente da aplicação do método. A mesma coisa

acontece com o modelo de turbulência. Como subproduto da nova formulação, temos

um modelo de turbulência baseado em LES, onde o modelo sub-malha – em inglês, sub-

grid – não precisa ser incluído, a priori, no modelo físico. Portanto, ele é considerado

um modelo implícito [28].

Esse tratamento resulta em um sistema de equações algébricas lineares (o problema

da não linearidade do termo convectivo da equação de quantidade de movimento é

contornado na discretização temporal através da linearização do termo convectivo),

simétricas, onde, em cada passo de tempo, é resolvido um problema “estacionário”.

Como as matrizes são simétricas, além do próprio GMRES, a nova formulação pode

utilizar um solucionador linear baseado em Gradientes Conjugados.

Outra característica importante do ponto de vista do desempenho computacional

está relacionada ao sistema de equações que precisa ser resolvido. Esse resulta em

52

sistemas semi-segregados, ou totalmente segregado dependendo da precisão temporal

do método, utilizando matrizes menores que as matrizes presentes no EdgeCFD®, já que

este resolve um único sistema totalmente acoplado dos graus de liberdade de pressão e

velocidades, para problemas puramente mecânicos.

Com base nas características do presente método, técnicas e ferramentas presentes

no EdgeCFD® foi realizada uma nova implementação da formulação de elementos

finitos descrita no Capítulo 3. O novo método possui precisão temporal de 2ª ordem,

discretização espacial com tetraedros lineares, ordenação de dados, estrutura de dados

por aresta, sistema de equações algébricas lineares resultando em matrizes simétricas

positivas definidas, solucionador linear baseado em gradientes conjugados pré-

condicionado e paralelismo com MPI e OpenMP.

A Tabela 4 sintetiza as principais diferenças entre as formulações e técnicas

computacionais comuns às duas implementações.

Tabela 4. Principais características entre os métodos implementados no EdgeCFD®.

Características Método original - EdgeCFD® Método presente

Discretização Espacial Temporal Temporal Espacial

Termos de estabilização SUPG / PSPG / LSIC (a priori) Surgem naturalmente

Modelo de turbulência LES – Smagorinsky / VMS LES – Implícito

Sistemas de equações EDOs não lineares totalmente acopladas u-p (4 graus de liberdade por nó). Não simétricas.

Algébricas lineares parcialmente acopladas p u. Simétricas positivas definidas.

Integração numérica temporal Preditor / Multi-corretor Não possui

Solucionador não linear Newton Inexato Não possui

Solucionador linear GMRES PCG

Estrutura de dados Aresta por aresta Aresta por aresta

Implementação paralela MPI / OpenMP / Híbrida MPI / OpenMP / Híbrida

As técnicas computacionais presentes no EdgeCFD®

e julgadas mais importantes

para a realização do presente estudo são a estrutura de dados por aresta e a programação

paralela em MPI/OpenMP. Como será visto adiante, a utilização de matrizes simétricas,

bem como o procedimento de solução parcialmente segregado, possibilitam uma

considerável economia de memória, redução de busca e gravação de informações na

memória (endereçamentos indiretos) e redução de operações de pontos flutuantes.

53

4.3 Estrutura de dados

A formulação de elementos finitos apresentada, aplicada às equações de Navier-

Stokes quase-incompressíveis, dá origem a um sistema simétrico de equações algébricas

lineares. Na prática, é como se tivéssemos um problema estacionário para ser resolvido

a cada passo de tempo. Portanto, não é necessária a utilização um esquema especial de

integração numérica para o avanço no tempo. Como o sistema é simétrico positivo-

definido, é utilizado o método dos gradiente conjugados, com pré-condicionamento

diagonal (para um grau de liberdade por nó) ou bloco diagonal, onde o grande esforço

computacional é atribuído ao produto matriz-vetor de cada iteração.

Desta forma, o foco para melhorar o desempenho computacional, tendo em vista o

consumo de memória, endereçamento indireto e operações de ponto flutuante, é a

redução do custo da operação de multiplicação matriz-vetor e passa necessariamente

pela forma com que a matriz é armazenada.

Na implementação 2D da presente formulação, foi utilizada uma estrutura de dados

elemento por elemento para armazenamento dos coeficientes da matriz produzida pelo

método. Lançar mão da utilização de uma estrutura de dados adequada é fundamental

quando se trata de elementos finitos devido ao elevado grau de esparsidade das matrizes

geradas por esse método de discretização. É importante minimizar o desperdício de

memória com armazenamento de uma quantidade grande de “zeros”. E o procedimento

adotado neste trabalho, em síntese, consiste na redução da matriz global em matrizes

menores como uma alternativa eficiente da redução do custo computacional.

Utilizando o EdgeCFD®

como base para a implementação 3D, foi herdada a

estrutura de dados por aresta. A opção pela implementação aresta por aresta é justificada

pelos bons resultados performances obtidos com o EdgeCFD® [48].

Uma descrição detalhada sobre a teoria de grafos e a estrutura de dados por aresta

podem ser encontradas em [59], [60].

54

4.4 Matrizes de aresta

Na prática, o EdgeCFD® constrói o lado esquerdo (matriz de “rigidez”) de cada

elemento e depois então a desmembra em matrizes de aresta [60]. Os coeficientes da

matriz de “rigidez” de um determinado elemento serão utilizados na construção das

matrizes das arestas pertencentes ao elemento.

Numa implementação tridimensional com elementos tetraédricos são armazenados

16 coeficientes da matriz de rigidez para cada elemento, para um grau de liberdade por

nó, já que a matriz de elemento é quadrada de ordem 4 (ngl noel , onde ngl é número

de graus de liberdade por nó e noel é número de nós por elemento). Os coeficientes

associados a um determinado nó guardam as contribuições de todos os elementos que

contêm esse nó. Quando é feito o desmembramento em matrizes de aresta, cada matriz

de elemento dá origem a 6 novas matrizes de aresta, quadradas de ordem 2. A forma

genérica de uma matriz de aresta é dada por:

e eij ije

ij e eij ij

k kT

k k (59)

onde eijk é o coeficiente correspondente à aresta ij do elemento e . Para mais de um

grau de liberdade, eijk corresponde a blocos de matrizes quadradas de ordem ngl . O

sinal negativo é uma convenção utilizada pelos arranjos de incidência do programa.

Como uma aresta pertence, em geral, a mais de um elemento, a matriz de aresta é

construída a partir das contribuições dos seus elementos concorrentes. A matriz de

aresta com as suas devidas contribuições pode ser escrita como:

1 1

1 1

n ne eij ij

a e eij n n

e eij ij

e e

k k

T

k k

(60)

onde n é o número total de elementos que compartilham a aresta ij e eijk é a sub-

matriz do elemento e correspondente aos graus de liberdade dos nós i e j da aresta a .

55

O desmembramento da matriz de rigidez de um elemento tetraédrico em matrizes

de aresta é mostrado na Figura 21:

11

22

3 3

13

31

23

32

24

42

12

2

4

43

3

44

4

1

1

41

e

k

k

k

k

k

k

k

k

k kk

kK

kk

k

k

Matriz de rigidez do elemento – um grau de

liberdade por nó, quatro nós por elemento.

11

22

12

21

eA

kT

k

kk

Matriz da aresta A

22

33

23

32

eB

kT

k

kk

Matriz da aresta B

11

33

13

31

eC

kT

k

kk

Matriz da aresta C

14

41

11

44

eD

kT

k

k

k

Matriz da aresta D

24

42

22

44

eE

kT

k

k

k

Matriz da aresta E

34

43

33

44

eF

kT

k

k

k

Matriz da aresta F

Figura 21. Desmembramento da matriz de rigidez do elemento tetraédrico em suas

matrizes de aresta correspondentes.

As matrizes eNT são sub-matrizes da matriz de rigidez do elemento e .

A dedução apresentada no Capítulo 2 é genérica para qualquer número de

dimensões espaciais. As alterações nas equações de atualização da pressão e da

velocidade ficam por conta da terceira dimensão agora presente. O sistema de equações

para atualização da pressão tem a sua forma preservada enquanto que o sistema de

equações para atualização da velocidade sofre uma ligeira modificação. Em 3D, as

novas matrizes, levando em conta o esquema de sincronização, são apresentadas como:

F

A

C

B

D E

1

2

3

4

56

p

M* * fp (61)

uu uv uw uTuv vv vw vT T T

wuw vw ww

A * A * A * * f

A * A * A * * f

* f A * A * A *

x

y

z

u

u

u

(62)

Ou simplesmente:

pM* * fp (63)

UA* U* f (64)

A Tabela 5 apresenta a quantidade de coeficientes armazenados para a matriz de

rigidez levando em conta as estruturas de dados por aresta e por elemento e os métodos

de 2ª ordem e o implementado originalmente no EdgeCFD®.

Tabela 5. Coeficientes armazenados para a matriz de rigidez para as estrutura de

dados por aresta e por elemento e para os métodos de 2ª ordem e o implementado

originalmente no EdgeCFD®.

EdgeCFD®

Método Proposto

Graus de

liberdade por nó 4 1 3

Elemento

192 coeficientes x número

de elementos

(+ 16 x número de nós)


de elementos



de elementos


Aresta


de arestas


1 coeficiente x número

de arestas



de arestas


Se considerarmos a montagem das matrizes na Equação (63) para cada elemento, a

matriz M* de pressão, de ordem noel , totalizaria 16 coeficientes, como já foi visto.

Aqui, no entanto, as matrizes são simétricas. Portanto, para cada elemento, são

armazenados 10 coeficientes. Para matrizes de aresta, são armazenados 3 coeficientes

ao invés de 4.

Para a montagem das matrizes na Equação (64) para cada elemento, a matriz

acoplada de velocidade A* , de ordem 12 (ngl noel ), teríamos 144 coeficientes

57

armazenados. Se as matrizes são simétricas, que é o que temos, ficamos com 90

coeficientes por elemento e 3 coeficientes por matriz de aresta. Sendo que nesse caso,

são 18 matrizes de aresta para cada matriz de elemento desmembrada.

Vale ressaltar que cada sub-matriz A * em (64), onde , u, v, w , não é

simétrica se considerada de forma independente.

No EdgeCFD® encontramos ainda um esquema especial para o armazenamento dos

coeficientes da diagonal de cada matriz de elemento. A utilização da matriz bloco

diagonal é importante no armazenamento dos coeficientes mencionados e na

determinação do pré-condicionador empregado no solucionador linear. A matriz bloco

diagonal possui também um papel fundamental na realização do produto matriz-vetor

implementado no EdgeCFD®, como pode ser visto em [60]. A matriz bloco diagonal

tem a forma:

( )blockB A (65)

A matriz B é a matriz bloco diagonal nodal obtida através do operador block

aplicado a matriz de rigidez. A matriz bloco diagonal é quadrada e possui ordem ngl .

Portanto, para a matriz de pressão, a matriz B possui apenas um coeficiente por nó.

Enquanto que, para a matriz de velocidade, a matriz B possui 9 coeficientes por nó.

Na prática, B

é declarado um vetor de ordem 3 3 nnodes , onde nnodes é o

número total de nós da malha, e que é “subaproveitado” durante o procedimento de

solução da pressão. Para a pressão, na verdade, bastaria um vetor de ordem nnodes

para armazenamento dos coeficientes da matriz diagonal.

Para um grau de liberdade por nó (como é o caso da pressão no método

apresentado), cada matriz de elemento possui 16 coeficientes. Como a matriz é

simétrica, são 10 coeficientes, como mencionado anteriormente. Dessa forma, 4

coeficientes da matriz do elemento são armazenados no vetor referente à matriz

diagonal (1 coeficiente para cada nó do elemento).

Para 3 graus de liberdade por nó (matriz de velocidade), cada matriz de elemento

possui 144 coeficientes. Destes, 36 são armazenados no vetor referente à matriz bloco

58

diagonal (9 coeficientes para cada nó do elemento). Cada matriz de elemento possui 4

matrizes bloco diagonais de ordem ngl , ao invés de 4 coeficientes como ocorre para um

grau de liberdade.

Utilizando a matriz bloco diagonal como recurso auxiliar para armazenamento dos

coeficientes da matriz de rigidez e considerando matrizes simétricas, temos então que,

para um grau de liberdade por nó, 4 coeficientes são armazenados na matriz diagonal e

os outros 6 na própria matriz de rigidez de elemento. Para 3 graus de liberdade por nó,

36 coeficientes são armazenados na matriz bloco diagonal e os outros 54 na própria

matriz de rigidez do elemento.

Por aresta, temos 4 coeficientes armazenados na sub-matriz diagonal do elemento e

apenas um coeficiente para cada sub-matriz de rigidez de aresta. A Figura 22 remonta a

estrutura apresentada na Figura 21 para o caso de um grau de liberdade por nó de

acordo com a formulação proposta:

59

11

22

3 3

141

43

44

3

2 13

42 2e

X

X XX

X

X

k

kk

k

k

k

Kk

k

kk

Matriz de rigidez do elemento – um grau de

liberdade por nó, quatro nós por elemento.

12eA

X

XX

kT

Matriz da aresta A

23eB

X

XX

kT

Matriz da aresta B

13eC

X

XX

kT

Matriz da aresta C

14eD

X

XX

kT

Matriz da aresta D

24eE

X

XX

kT

Matriz da aresta E

34eF

X

XX

kT

Matriz da aresta F

11

22

33

44

e X

X

X

X

X

X

XX

X

XX

X

k

kB

k

k

Matriz bloco diagonal do elemento e

Coeficientes do elemento e armazenados no vetor

global de aresta: 13 24 31 22 414 3,, ,,, kk k kk k

Coeficientes do elemento e armazenados no vetor

global bloco diagonal: 11 22 33 44, , ,k k k k

Figura 22. Desmembramento da matriz simétrica de rigidez do elemento tetraédrico

em suas matrizes de aresta e bloco diagonal correspondentes. X

representa os

coeficientes “repetidos” ou armazenados em outra sub-matriz.

F

A

C B

D E

1

2

3

4

60

Para a solução do sistema linear Ax = b , o pré-condicionamento diagonal (para

um grau de liberdade por nó) ou bloco diagonal nodal utilizado, tanto na implementação

original do EdgeCFD®

quanto na nova formulação, consiste da pré multiplicação

-1 -1B Ax = B b , onde a matriz -1B é a inversa da matriz B . A matriz B definida pela

Equação (65) é positivo-definida e de fácil inversão.

Para a solução dos graus de liberdade de pressão, ou seja, para a solução de um

grau de liberdade por nó, o pré-condicionamento é diagonal, como já foi mencionado.

Para a matriz ii nnodesB k , a sua inversa é definida como

1

ii nnodes

-1Bk

.

O sinal negativo do coeficiente iik , atribuído à convenção utilizada nos arranjos

de incidência do programa, foi suprimido por conveniência.

61

Capítulo 5

Exemplos Numéricos Tridimensionais

Como exemplo do que foi feito nos estudos bidimensionais, a implementação 3D,

utilizando a estrutura do EdgeCFD®, foi testada no escoamento induzido em uma

cavidade pelo movimento de uma tampa.

Após a verificação do método com exemplos bidimensionais, neste capítulo são

apresentados os resultados de verificação do método implementado no EdgeCFD®, com

ênfase na avaliação do desempenho computacional do método.

5.1 Escoamento induzido no interior de uma cavidade pelo movimento

de uma tampa – 3D

Nesta seção são apresentados os resultados para o escoamento induzido no interior

de uma cavidade 3D pelo movimento de uma tampa para 1 000L

Re . . O número de

Reynolds é definido como 0L

Re u L / , onde L é comprimento da aresta da

cavidade cúbica e 0u é a velocidade de referência do escoamento.

A Figura 23 mostra o modelo geométrico da cavidade cúbica e a malha utilizada no

exemplo.

62

Figura 23. Escoamento induzido no interior de uma cavidade 3D pelo movimento

de uma tampa. Malha não estruturada com 108.104 elementos tetraédricos e 20.589

nós.

A tampa, plano superior perpendicular ao eixo z , se move na direção e no sentido

positivo do eixo x . O movimento da tampa é simulado pela prescrição das velocidades

1 0xu , , 0 0

yu , e 0 0

zu , .

Os dois planos perpendiculares ao eixo y definem a parte da frente e a parte de trás

da cavidade. Nestes planos é prescrita a condição de impenetrabilidade ( 0 0yu , ).

Conseqüentemente, xu e z

u são deixadas livres.

Os planos laterais da cavidade, planos perpendiculares ao eixo x , e sua base são

consideradas paredes sólidas. A condição de não deslizamento é necessária para

verificação do efeito da viscosidade do fluido.

A Figura 24 apresenta as linhas de corrente da solução para 1 000L

Re . .

63

Figura 24. Linhas de corrente no escoamento induzido no interior de uma cavidade

3D pelo movimento de uma tampa para 1 000L

Re . .

A Figura 25 mostra a comparação entre as componentes de velocidades na linha

central da cavidade cúbica, obtidas com o método original do EdgeCFD®

, previamente

validado [53], e o presente método.

64

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

-0,6

-0,4

-0,2

0,0

0,2

0,4

0,6

0,8

1,0

Presente

EdgeCFD

Co

mp

on

en

te u

x

Coordenada y


(a)

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

-0,6

-0,4

-0,2

0,0

0,2

0,4

0,6

0,8

1,0

Componente uy de velocidade na linha de centro vertical para Re = 1.000

Presente

EdgeCFD

Co

mp

on

en

te u

y

Coordenada x

(b)

Figura 25. Velocidades na linha de centro da cavidade cúbica para 1 000L

Re . .


yu (b). Malha não estruturada de 108.104 elementos

tetraédricos e 20.589 nós.

65

A simulação do problema do escoamento induzido por uma tampa numa cavidade

3D apresentado foi realizada pelo o EdgeCFD®

com o seu método original e com o

método proposto. Bons resultados para esse problema e a conseqüente validação do

EdgeCFD®

podem ser encontrados em [60].

O presente método é adequado para a simulação de problemas transientes. Para

problemas estacionários, o EdgeCFD® é mais indicado. A afirmação se deve à sua

formulação original, que requer uma estratégia de solução totalmente implícita. Além

do recurso de utilização de passos de tempos fixos, maiores que os definidos pela escala

de tempo local apresentados na seção 3.3. O tempo de CPU gasto pelo EdgeCFD®

na

solução do problema apresentado é bem inferior ao conseguido com a presente

formulação.

Para a comparação do desempenho computacional do presente método, portanto,

decidiu-se pela utilização de um método de características semelhantes às suas. A

formulação apresentada no Apêndice B foi implementada no EdgeCFD® mantendo

todas as características que tinham sido mantidas quando da implementação do método

de precisão temporal de 2ª ordem. Sendo assim, as principais diferenças entre as duas

novas versões do EdgeCFD®

são a ordem de precisão temporal e o procedimento de

solução. Enquanto que no método presente temos precisão temporal de 2ª ordem e o

procedimento de solução parcialmente acoplado, na versão obtida a partir da

implementação da formulação descrita no Apêndice B, temos precisão temporal de 1ª

ordem e procedimento de solução totalmente segregado.

A Figura 26 apresenta os gráficos do tempo de processamento, speedup e eficiência

computacional em função do número de processadores para os métodos de elementos

finitos de 1ª e 2ª ordem temporal. Vale lembrar que ambos os métodos, uma vez

implementados no EdgeCFD®, utilizam reordenação de dados e estrutura de dados por

aresta.

As características do exemplo anterior foram mantidas e o tempo de duração da

simulação para este estudo foi fixado em 05t L / u .

66

(a)

(b)

(c)

Figura 26. Desempenho computacional. Tempo de processamento (a), speedup (b)

e eficiência computacional (c) para os métodos de 1ª e 2ª ordem. Escoamento induzido

no interior de uma cavidade 3D pelo movimento de uma tampa para 1 000L

Re . .

Malha tetraédrica com 108.104 elementos e 20.589 nós.

0,00

0,20

0,40

0,60

0,80

1,00

1,20

1 2 4 8 16 32

Tem

po

(h

ora

s)

Número de CPUs

1ª Ordem

2ª Ordem

0,00

2,00

4,00

6,00

8,00

10,00

12,00

14,00

0 10 20 30 40

Spe

ed

up

Número de CPUs

1ª Ordem

2ª Ordem

0,00

0,20

0,40

0,60

0,80

1,00

1,20

1 2 4 8 16 32

Efic

iên

cia

Número de CPUs

1ª Ordem

2ª Ordem

67

Até aqui, as simulações 3D apresentadas foram feitas com a tolerância de 10-3

do

solucionador linear (gradientes conjugados pré-condicionado) tanto para a pressão

quanto para as velocidades. No caso do método de 1ª ordem, a mesma tolerância de 10-3

foi adotada para a pressão e para as componentes de velocidade, que são resolvidas

separadamente. Para o método de 2ª ordem, o mesmo foi feito. Sendo que, nesse caso,

as componentes de velocidade são resolvidas juntas.

A fim de verificar o desempenho computacional a partir da alteração da tolerância

do solucionador linear, algumas combinações de tolerância para pressão e velocidade

foram realizadas.

Nas Figuras 27 e 28 são apresentados os gráficos de tempo de processamento e de

speedup em função do número de processadores para duas dessas combinações.

Também foi realizado um teste com uma malha mais refinada, com 1.613.455

elementos tetraédricos e 279.362 nós. A tolerância do solucionado linear foi ajustada

em 10-3

para a pressão e para as componentes de velocidade para ambos os métodos.

Foi mantido 1 000L

Re . e o tempo de simulação foi de 0

t L / u . Esses resultados

são mostrados na Figura 29.

68

(a)

(b)

(c)


e eficiência computacional (c) para os métodos de 1ª e 2ª ordem. Tolerância do

solucionador linear de 10-6

para pressão e para velocidade.

0,00

0,20

0,40

0,60

0,80

1,00

1,20

1 2 4 8 16 32

Tem

po

(h

ora

s)

Número de CPUs

1ª Ordem

2ª Ordem

0,00

2,00

4,00

6,00

8,00

10,00

12,00

0 10 20 30 40

Spe

ed

up

Número de CPUs

1ª Ordem

2ª Ordem

0,00

0,20

0,40

0,60

0,80

1,00

1,20

1 2 4 8 16 32

Efic

iên

cia

Número de CPUs

1ª Ordem

2ª Ordem

69

(a)

(b)

(c)

Figura 28. Desempenho computacional. Tempo de processamento (a) e speedup (b)




para as componentes da velocidade.

0,00

0,20

0,40

0,60

0,80

1,00

1,20

1 2 4 8 16 32

Tem

po

(h

ora

s)

Número de CPUs

1ª Ordem

2ª Ordem

0,00

2,00

4,00

6,00

8,00

10,00

12,00

14,00

0 10 20 30 40

Spe

ed

up

Número de CPUs

1ª Ordem

2ª Ordem

0,00

0,20

0,40

0,60

0,80

1,00

1,20

1 2 4 8 16 32

Efic

iên

cia

Número de CPUs

1ª Ordem

2ª Ordem

70

(a)

(b)

(c)





para as componentes da velocidade.

Malha tetraédrica com 1.613.455 elementos e 279.362 nós.

0,00

0,10

0,20

0,30

0,40

0,50

0,60

1 2 4 8 16 32

Tem

po

(h

ora

s)

Número de CPUs

1ª Ordem

2ª Ordem

0,00

2,00

4,00

6,00

8,00

10,00

12,00

14,00

0 10 20 30 40

Spe

ed

up

Número de CPUs

1ª Ordem

2ª Ordem

0,00

0,20

0,40

0,60

0,80

1,00

1,20

1 2 4 8 16 32

Efic

iên

cia

Número de CPUs

1ª Ordem

2ª Ordem

71

Perceba que o tempo de processamento é quase sempre menor para o método de 2ª

ordem em rodadas paralelas, independente da malha utilizada ou da tolerância do

solucionador linear.

O speedup e a eficiência computacional são melhores para o método de 2ª ordem.

Na verdade, para a malha com mais elementos o método de 1ª ordem se torna melhor.

Com 2 e 8 processadores o speedup e a eficiência são ligeiramente melhores para o

método de 1ª ordem em todos os casos apresentados.

Diante dos resultados apresentados nesse exemplo podemos concluir que os dois

métodos apresentam boa escalabilidade e bom desempenho computacional.

As Tabelas 6, 7 e 8 mostram o tempo de CPU e o número de iterações para solução

dos sistemas de pressão e velocidade para os métodos de 1ª e 2ª ordem de acordo com o

número de processadores e a tolerância do solucionador linear. A malha utilizada possui

108.104 elementos e 20.589 nós.

Tabela 6. Tempo de processamento e número de iterações por passo de tempo para

solução das equações de pressão e velocidade em função da tolerância do solucionador

linear para os métodos de 1ª e 2ª ordem. Tolerância de 10-3

para a pressão e para a

velocidade.

toler p = 10-3; toler U = 10-3 1ª Ordem 2ª Ordem

Número de CPUs

Dados de Desempenho

p ux uy uz Tempo de CPU por passo de tempo (s)

p U Tempo de CPU por passo de tempo (s)

1 Iterações 108 3 3 3

3,45 113 3

3,00 Tempo de CPU (s) 0,156 0,008 0,008 0,008 0,668 0,044


1,77 80 5

1,59 Tempo de CPU (s) 0,052 0,004 0,004 0,004 0,356 0,044


0,94 70 7

0,78 Tempo de CPU (s) 0,024 0,004 0,004 0,004 0,132 0,036


0,51 66 7

0,48 Tempo de CPU (s) 0,020 0,008 0,008 0,008 0,124 0,044


0,38 53 7

0,26 Tempo de CPU (s) 0,048 0,008 0,008 0,008 0,064 0,028


0,30 43 7

0,18 Tempo de CPU (s) 0,076 0,012 0,016 0,012 0,076 0,024

72




para a pressão e para a

velocidade.


Número de CPUs

Dados de Desempenho




3,55 181 7

3,44 Tempo de CPU (s) 0,252 0,012 0,008 0,012 1,068 0,084


1,83 159 9

1,95 Tempo de CPU (s) 0,092 0,004 0,004 0,008 0,696 0,076

4 Iterações 140 11 11 11

0,98 141 13

0,92 Tempo de CPU (s) 0,044 0,008 0,004 0,004 0,256 0,064

8 Iterações 121 12 12 12

0,52 121 14

0,60 Tempo de CPU (s) 0,036 0,008 0,008 0,008 0,216 0,072

16 Iterações 98 12 12 12

0,39 99 13

0,34 Tempo de CPU (s) 0,084 0,016 0,012 0,012 0,120 0,044

32 Iterações 82 13 13 14

0,33 83 15

0,27 Tempo de CPU (s) 0,128 0,024 0,020 0,024 0,140 0,040




para a pressão e de 10-3

para a

velocidade.


Número de CPUs

Dados de Desempenho




3,54 181 3

3,40 Tempo de CPU (s) 0,248 0,008 0,008 0,008 1,056 0,044


1,82 158 5

1,92 Tempo de CPU (s) 0,100 0,008 0,008 0,004 0,684 0,044


0,96 141 7

0,90 Tempo de CPU (s) 0,048 0,004 0,004 0,004 0,260 0,036


0,52 121 7

0,56 Tempo de CPU (s) 0,032 0,008 0,008 0,008 0,216 0,044


0,43 99 7

0,38 Tempo de CPU (s) 0,088 0,008 0,008 0,008 0,140 0,080


0,30 83 8

0,26 Tempo de CPU (s) 0,124 0,020 0,020 0,020 0,140 0,026

73

É importante considerar que o tempo de processamento e o número de iterações na

solução de cada sistema não são os mesmos ao longo de toda a simulação. Foram

selecionados os menores tempos de processamento de cada caso.

Em ambos os métodos, a solução da pressão é quem precisa de mais iterações até a

convergência. A velocidade converge muito mais rápido. Em alguns casos, verificamos

que a solução da pressão consome 15 vezes mais tempo que a solução da velocidade.

Isso pode ser constado na rodada serial para o método de 2ª ordem apresentada na

Tabela 6. O restante do tempo é utilizado, principalmente, na construção das equações.

Note a razão entre o tempo de CPU gasto na solução das equações em um passo de

tempo, e o tempo total de CPU por passo tempo, muda de acordo com o número de

processadores e a tolerância do solucionador linear para cada sistema.

O programa utilizado para a implementação bidimensional da presente formulação

possui refinamento adaptativo da malha de elementos finitos. Para tal, é necessária

atualização dos gradientes das funções de interpolação considerando os novos dados dos

nós gerados e da conectividade entre eles, sempre que é feito o refinamento. O

EdgeCFD® apresenta, alternativamente, uma formulação estabilizada de elementos

finitos integrada em um domínio Langrangiano-Euleriano Arbitrário [63]. Como nessa

formulação há a alteração das coordenadas dos nós da malha, este recurso também

requer o cálculo dos gradientes das funções de interpolação a cada atualização dos

dados da malha. Na implementação tridimensional dos métodos de precisão temporal de

1ª e 2ª ordem, é realizado o cálculo dos gradientes das funções de interpolação toda a

vez que os coeficientes dos sistemas de equações são determinados. Dessa forma, temos

um elevado tempo total de CPU por passo de tempo, comparado com o tempo de CPU

gasto somente na solução dos sistemas de equações a cada passo de tempo. O cálculo do

passo de tempo local, inerente às formulações de 1ª e 2ª ordem, mesmo sendo

paralelizado, representa um custo de 10% do tempo de CPU por passo de tempo.

Ainda assim, a razão entre o tempo de CPU gasto em cada passo de tempo e o

tempo gasto para a solução dos sistemas de equações é muito grande em ambos os

métodos, incompatível com o que é normalmente observado nos testes realizados com

EdgeCFD® e no código utilizado para os exemplos bidimensionais apresentado nesta

tese. Esse gasto excessivo de tempo precisa ser mais bem estudado.

74

O tempo de processamento para a solução da pressão no método de 2ª ordem é, em

geral, maior do que no método de 1ª ordem. Por outro lado, como somente um sistema é

construído para a solução da velocidade e todas as componentes são resolvidas de uma

única vez, o tempo global de processamento por passo de tempo é menor para o método

de 2ª ordem. Isso ocorre mesmo sendo o tempo de solução do sistema acoplado de

velocidade maior que o tempo gasto na solução segregada.

Os dados de desempenho computacional apresentados na presente seção, foram

obtidos no cluster do Laboratório de Tecnologia Oceânica (LabOceano), que faz parte

do Programa de Engenharia Naval e Oceânica da COPPE/UFRJ. A Tabela 9 mostra

uma breve descrição das características do cluster.

Tabela 9. Características do cluster do LabOceano.

SGI / Altix XE 1200

Processadores 11 CPUs Quad Core Intel® Xeon

® CPU X5355 @ 2,66 GHz: 88 cores

Memória 16 GBytes RAM por nó

Compiladores Intel Fortran-90 e C/C++ com suporte OpenMP e MPI

Rede Infiniband e gigabit

Swicht Voltaire

Armazenamento 47 Gbytes + 840 Gbytes

Sistema Operacional Linux Suse Enterprise V10 SP1

Sistema de Filas Torque

Gerenciador do Cluster Scalimanager

A Tabela 10 mostra o uso de memória com os métodos utilizados no estudo para

uma malha com 8 partições do problema anterior: o método originalmente

implementado no EdgeCFD®, o método de precisão temporal de 2ª ordem e o método

de precisão temporal de 1ª ordem.

75

Tabela 10. Demanda de memória de cada método para o problema do escoamento

induzido pelo movimento de uma tampa em uma cavidade cúbica 3D. Malha não

estruturada com 108.104 elementos tetraédricos e 20.589 nós.

Método Memória (MB) / processador

1ª Ordem 2,6

2ª Ordem 3,8

EdgeCFD®

20,7

Além do ganho na implementação por aresta quando comparado a estrutura de

dados elemento por elemento, verifica-se também que os métodos com discretização

temporal de 1ª e 2ª ordem, que utilizam métodos de solução totalmente ou parcialmente

segregado, possuem economia de memória muito maior que o método totalmente

acoplado do EdgeCFD®

. O método de 1ª chega a utilizar cerca de 12% apenas da

memória utilizada pelo EdgeCFD®

para o problema em questão.

Por outro lado, para a solução de problemas estacionários, o método totalmente

acoplado do EdgeCFD®

é muito mais rápido. A Figura 30 apresenta os gráficos de

tempo de processamento e de speedup em função do número de processadores para os

métodos de 1ª e 2ª ordem implementados na estrutura do EdgeCFD®

e para o

EdgeCFD® e o seu método original.

76

(a)

(b)

Figura 30. Desempenho computacional. Tempo de processamento (a) e speedup

(b) para os métodos de 1ª e 2ª ordem e o EdgeCFD®. Malha não estruturada com

108.104 elementos tetraédricos e 20.589 nós.

Como podemos verificar na Figura 30, os tempos de CPU obtidos são

incomparáveis. Como já era previsto, o EdgeCFD® original é muito mais veloz. Com o

programa configurado com um passo de tempo fixo apropriado para o problema

resolvido, maior que o passo local calculado a partir das escalas de tempo de convecção

e difusão, o EdgeCFD®

atinge o tempo final de simulação em relativamente poucos

passos de tempo. No entanto, o método com precisão temporal de 2ª ordem aproveita

0,00

0,20

0,40

0,60

0,80

1,00

1,20

1 2 4 8 16 32

Tem

po

(h

ora

s)

Número de CPUs

1ª Ordem

2ª Ordem

EdgeCFD

0,00

2,00

4,00

6,00

8,00

10,00

12,00

14,00

0 10 20 30 40

Spe

ed

up

Número de CPUs

1ª Ordem

2ª Ordem

EdgeCFD

77

melhor a implementação paralela, o que é refletido em sua melhor escalabilidade. Isso

acontece por causa da demanda de memória por processador. Como o uso de memória

por processador nos métodos de discretização temporal de 1ª e 2ª ordem é menor que o

método totalmente acoplado do EdgeCFD®, a presente formulação é beneficiada pelo

recurso da memória cache dos processadores, o que reduz a busca e gravação de

informações na memória do computador.

Apesar de o tempo total de CPU do EdgeCFD®

para problemas estacionários ser

muito menor que o tempo de CPU do método de 2ª ordem, o seu tempo de CPU por

passo de tempo é bem superior. A Tabela 11 apresenta os dados de tempo de

processamento por passo de tempo em função do número de processadores para o

EdgeCFD® original e o presente método.

Tabela 11. Tempo de processamento por passo de tempo para solução das

equações de pressão e velocidade em função do número de processadores para o

EdgeCFD® original e o presente método. Malha não estruturada com 108.104

elementos tetraédricos e 20.589 nós.

Número

de CPUs

EdgeCFD®

2ª Ordem Toler_linear = 10-3; Toler_n-linear = 10-3. Toler_linear = 10-3

Tempo de CPU por passo de tempo (s) Tempo de CPU por passo de tempo (s)

1 5,33 3,00

2 3,37 1,59

4 1,89 0,78

8 1,75 0,48

16 1,04 0,26

32 0,90 0,18

As tolerâncias do solucionador linear e do solucionador não linear do EdgeCFD®

foram ajustadas em 10-3

. Outros parâmetros inerentes ao método original do EdgeCFD®

também foram configurados de forma adequada para o exemplo em questão. A

tolerância do solucionador linear para o método de 2ª ordem foi ajustado também em

10-3

.

Os dados apresentados Figura 30 e na Tabela 11 mostram que, na simulação de um

problema transiente, onde a qualidade da solução está diretamente ligada à precisão da

78

escala de tempo adotada, o método estável de elementos finitos com precisão temporal

de 2ª ordem pode se mostrar altamente competitivo.

A expressão “tempo de CPU” é usada sempre que é feita referência ao tempo gasto

no processamento. No entanto, no processamento paralelo o tempo é monitorado e

impresso por um único processador, o mais ocioso deles, que espera todos os outros

finalizarem as suas tarefas e compartilharem os dados processados. Portanto, “tempo de

parede” seria uma expressão mais adequada para o tempo de processamento monitorado

em execuções paralelas.

79

Capítulo 6

Conclusões

Uma formulação de elementos finitos com precisão temporal de 2ª ordem foi

apresentada. Os balanços de massa e de quantidade de movimento são combinados em

uma expansão em Série de Taylor para a pressão. Esta é discretizada no espaço com o

método de Galerkin e resulta numa equação adequada para a atualização da pressão.

A discretização temporal do balanço da quantidade de movimento é realizada com

diferenças finitas. A minimização do resíduo espacial é feita, utilizando o método dos

mínimos quadrados, para obter as equações para atualização das componentes de

velocidade.

O método proposto introduz automaticamente os termos de estabilização

necessários para controlar as oscilações espúrias em problemas com convecção

dominante e responsáveis por contornar as restrições de Babuška-Brezzi na escolha dos

espaços de interpolação para pressão e velocidade. Em outras palavras, a formulação é

estabilizada de forma natural, ou seja, os termos de estabilização surgem na dedução das

equações ao invés de serem propostos a priori como nas formulações tradicionais. O

método conduz a um sistema parcialmente acoplado, onde os graus de liberdade de

pressão são resolvidos primeiro e, na seqüência, os graus de liberdade de velocidade são

calculados.

O método proposto utiliza ainda um esquema especial de sincronização,

considerando que cada elemento possui o seu próprio passo de tempo. O passo de tempo

local é necessário para que o método funcione de forma apropriada e ele é definido a

80

partir das escalas de tempo dos fenômenos físicos em cada região do domínio

garantindo, dessa forma, que uma quantidade adequada de estabilização seja introduzida

de acordo com as características do escoamento em cada elemento. A sincronização é

realizada adotando o menor passo de tempo calculado entre as escalas de convecção e

difusão de cada elemento do domínio computacional.

A precisão transiente superior do método presente é mostrada na solução dos

estágios iniciais do escoamento com expansão abrupta, onde são comparadas as

predições do presente método com as predições realizadas com um método similar,

porém, de 1ª ordem.

A eficiência do método na estabilização da solução é mostrada na simulação de

escoamentos incompressíveis (condição limite do quase-incompressível), como é feito

claramente nos exemplos que envolvem o escoamento induzido pelo movimento de uma

tampa no interior de uma cavidade e no escoamento de Kim-Moin.

A solução do desprendimento periódico de vórtices em escoamento cruzado sobre

um cilindro circular também mostra o bom desempenho do método na simulação de

escoamentos viscosos transientes incompressíveis.

A alta taxa de convergência espacial do método para as variáveis de pressão e

componentes de velocidade é mostrada no escoamento de Kim-Moin, bem como a

evolução temporal do erro. O crescimento do erro temporal nesse método para as

componentes de velocidade é limitado, ficando relativamente estável ao longo do tempo

de análise.

O método de elementos finitos proposto foi implementado no software EdgeCFD®.

O programa utiliza elementos tetraédricos lineares. Dentre as principais características

do programa, foram mantidas reordenação de dados, a estrutura de dados por aresta e o

paralelismo. Foi feita a comparação dos resultados obtidos com o EgdeCFD® original e

a formulação proposta para o problema do escoamento induzido pelo movimento de

uma tampa numa cavidade.

Também foi feita a implementação 3D no EdgeCFD® da formulação estável de 1ª

ordem. O escoamento induzido pelo movimento de uma tampa no interior de uma

cavidade foi utilizado para a análise do desempenho computacional com as formulações

81

estabilizadas com discretização temporal de 1ª e de 2ª ordem. Verificamos que a

economia de memória proporcionada pelas formulações de 1ª e 2ª ordem, cujos métodos

de solução são segregados (total ou parcial), associadas à estrutura de dados por aresta,

é muito maior que o método totalmente acoplado implementado no EdgeCFD®.

Trabalhos Futuros

Como sugestão de trabalhos futuros, podemos elencar as seguintes possibilidades:

1) Investigação mais minuciosa da qualidade dos resultados obtidos em problemas

típicos tridimensionais: A verificação da qualidade dos resultados pode ser mais bem

explorada, utilizando uma faixa maior do número de Reynolds para o problema aqui

apresentado e respectivos resultados de referência. Bem como para problemas

transientes 3D, já que a ênfase dada aqui foi ao desempenho computacional do método

de precisão temporal de 2ª ordem numa implementação aresta por aresta paralela.

2) Identificação do modelo implícito de turbulência: O presente método, assim como o

método descrito no Apêndice B, se destina a simulação de escoamentos turbulentos. A

validade dessa afirmação fica por conta da realização de testes de escoamentos com

elevados números de Reynolds e a identificação do modelo de turbulência na

formulação, como foi feito em [28].

3) Obtenção de mais dados comparativos para a análise do desempenho

computacional do método: Uma vez que o método foi implementado no EdgeCFD®, a

análise do desempenho do método pode justificar a sua viabilidade como uma

ferramenta de maior precisão temporal dentro do programa, devendo ser utilizada de

acordo com as características do problema a ser resolvido.

4) Acoplamento térmico-mecânico: Inclusão da equação de energia nas equações

governantes e aplicação da formulação de precisão temporal de 2ª ordem à temperatura.

Tal integração deverá dar origem a um sistema parcialmente acoplado dos graus de

liberdade de velocidade e temperatura, onde a pressão deverá ser resolvida primeiro e

separadamente. Um programa 2D com acoplamento térmico-mecânico e formulação

estável de precisão temporal de 2ª ordem tem sido desenvolvido em um estudo paralelo

82

realizado por estudantes do Programa de Pós Graduação do Instituto de Engenharia

Nuclear e o desafio nesse momento seria a sua implementação 3D, utilizando técnicas

de reordenação de dados e estrutura de dados por aresta.

5) Aplicação do método aos problemas de engenharia nuclear: A pesquisa realizada ao

longo dos últimos 4 anos tinha como objetivo inicial o desenvolvimento de um código

computacional 3D para simulação de problemas de mecânica de fluidos com

acoplamento térmico, típicos de projeto de reatores nucleares. Apesar de conclusão da

implementação 3D do método para simulação de problemas puramente mecânico,

alguns testes ainda devem ser feitos. E como sugestão fica a aplicação efetiva do

método às pesquisas referentes ao desenvolvimento de reatores inovadores.

6) Aplicação do método a problemas de escoamento bifásico, de superfície livre ou

escoamento em meios porosos: Na formulação estável de 1ª ordem temporal

apresentada no apêndice B, a capacidade de simulação desses problemas se resume

simplesmente a inclusão de uma nova equação às equações governantes do escoamento.

Quando se pretende preservar a precisão temporal de 2ª ordem, o procedimento de

dedução das equações resultantes da formulação de elementos finitos é um pouco mais

complexo, impactando também em outros aspectos. Por exemplo, se a “habilitação” do

software para simulação de problemas envolvendo superfície livre que dependem, por

exemplo, da implementação do método Volume of Fluids seguir a mesma linha no

acoplamento térmico-mecânico, isso implica em mais uma variável acoplada. E aí se

deve verificar até que ponto o investimento vale a pena em detrimento do desempenho

computacional.

83

Referências

[1] Bird, R. B., Stewart, W. E., Lightfoot, E. N., 1960, Transport Phenomena, New

York, John Wiley & Sons ltd.

[2] Fox, R. W., McDonald, A. T., 2001, Introdução à Mecânica dos Fluidos, 5ª edição,

Rio de Janeiro, LTC.

[3] INTERNATIONAL ATOMIC ENERGY AGENCY (IAEA), 2003, Use of

Computational Fluid Dynamics Codes for Safety Analysis of Nuclear Reactor

Systems, IAEA-TECDOC-1379.

[4] Zienkiewicz, O.C., Qu, S., Taylor, R.L., Nakazawa, S., 1986, “The patch test for

mixed formulation”, International Journal for Numerical Methods in

Engineering, v. 23, pp. 1873–1883.

[5] Taylor, C. Hood, P., 1973, “A Numerical solution of the Navier-Stokes equations

using the finite element technique”, Computer & Fluids, v. 1 (1), pp. 73–100.

[6] Coruzeix, M., Raviart, P.-A., 1973, “Conforming and nonconforming finite element

methods for solving the stationary Stokes equations I”, RAIRO Analyse

Numerique, v. 7 (R-3), pp. 33–75.

[7] Arnold, D. N., Brezzi., Fortin, M., 1984, “A stable finite element for Stokes

equations”, Calcolo, v. 21(4), pp. 337–344.

[8] Bercovier, M., 1978, “Perturbation of a mixed variational problems. Application to

mixed finite element methods”, RAIRO Analyse Numerique, v. 12 (3), pp. 211–

236.

[9] Malkus D. S., Hughes, T. J. R., 1078, “Mixed finite element methods – reduced and

selective integration techniques: a unification of concepts”, Communications in

Numerical Methods in Engineering, v. 15 (1), pp. 63–81.

84

[10] Zienkiewics, O. C., Taylor, R. L., 2000, The Finite Element Method, Volume 1: The

Basis, 5th

edition, Oxford, Butterworth-Heinemann.

[11] Bercovier, M., Engelman, M., 1979, “A finite element for the numerical solution of

viscous incompressible flows”, Journal of Computational Physics, v. 30 (2),

pp. 181 – 201.

[12] Girault, V., Raviart, P.-A, 1986, Finite Element Methods for Navier-Stokes

Equations. Theory and Algorithms, Berlin, Springer-Verlag.

[13] Terman, R., 2001, Navier-Stokes Equations. Theory and Numerical Analysis, 1984,

Providence, RI, AMS Chelsea Publishing. Correct reprint of the 1984 edition

[North-Holland, Amsterdam, 1984].

[14] Brezzi, F., Fortin, M., 1991, Mixed and Hybrid Finite Element Methods, New York,

Springer.

[15] Donea, J., Huerta, A., 2003, Finite Elements Methods for Flow Problems,

Chichester, England, John Wiley & Sons ltd.

[16] Tezduyar, T. E., 1992, “Stabilized finite element formulations for incompressible

flow computations”, Advances in Applied Mechanics, v. 28.

[17] Tezduyar, T.E., Mittal, S., Ray, S.E., Shih, R., 1992, “Incompressible flow

computations with stabilized bilinear and linear equal-order-interpolation

velocity-pressure elements”, Computer Methods in Applied Mechanics and


[18] Franca, L. P., Frey, S. F., 1992, “Stabilized Finite Element Methods: II. The

Incompressible Navier-Stokes Equations”, Computer Methods in Applied

Mechanics and Engineering, v. 99, pp. 209–233.

[19] Tezduyar, T. E., Osawa, Y., 2000, “Finite element stabilization parameters

computed from element matrices and vectors”, Computer Methods in Applied

Mechanics and Engineering, v. 190, pp. 411 – 430.

[20] Tezduyar, T. E., 2001, Adaptative determination of the finite element stabilization

parameters, Rice University – Mechanical Engineering and Materials Science.

85

[21] Codina, R., 2000, “On stabilized finite element method for linear systems of

diffusion-convecção-reaction equation”, Computer Methods in Applied

Mechanics and Engineering, v. 156 (1–4), pp. 185–210.

[22] Zienkiewics, O. C., Taylor, R. L., 2000, The Finite Element Method, Volume 3:

Fluid Dynamics, 5th

edition, Oxford, Butterworth-Heinemann.

[23] De Sampaio, P.A.B., 2005, “A finite element formulation for transient

incompressible viscous flows stabilized by local time-steps”, Computer

Methods in Applied Mechanics and Engineering, v. 194, pp. 2095–2108.

[24] De Sampaio, P.A.B, 2006, “A stabilized finite element method for incompressible

flow and heat transfer: a natural derivation based on the use of local time-

steps”, Computer Methods in Applied Mechanics and Engineering, v. 195 (44–

47), pp.6177–6190.

[25] De Sampaio, P.A.B., 1993, “Transient solutions of the incompressible Navier–

Stokes equations in primitive variables employing optimal local time

stepping”, In: Proceedings of the 8th International Conference on Numerical

Methods for Laminar and Turbulent Flow, pp. 1493–1504, Swansea, Jul.

[26] De Sampaio, P.A.B., Coutinho, A.L.G.A., 1999, “Simulation of free and forced

convection incompressible flows using an adaptive parallel/vector finite

element procedure”, International Journal for Numerical Methods in Fluids, v.

29, pp. 289–309.

[27] De Sampaio, P.A.B., Hallak, P.H., Coutinho, A.L.G.A., Pfeil, M.S., 2004, “A

stabilized finite element procedure for turbulent fluid–structure interaction

using adaptive time–space refinement”, International Journal for Numerical

Methods in Fluids, v. 44 (6), pp. 673–693,Chichester, UK.

[28] De Sampaio, P.A.B., Gonçalves Jr., M.A., Lapa, C.M.F., 2008, “A CFD approach

to the atmospheric dispersion of radionuclides in the vicinity of NPPs”.

Nuclear Engineering and Design, v. 238, pp.250–273.

[29] Brooks, A., Hughes, T.J.R., 1982, “Streamline–upwind/Petrov–Galerkin

formulations for convection dominated flows with particular emphasis on the

86

incompressible Navier–Stokes equations”. Computer Methods in Applied

Mechanics and Engineering, v. 32, pp. 199–259.

[30] De Sampaio, P.A.B., 1991, “A Petrov–Galerkin formulation for the incompressible

Navier–Stokes equations using equal order interpolation for velocity and

pressure”, International Journal for Numerical Methods in Engineering, v. 31,

pp. 1135–1149.

[31] Bochev, P.B., Gunzburger, M.D., Lehoucq, R.B., 2007, “On stabilized finite

element methods for the Stokes problem in the small time step limit”,

International Journal for Numerical Methods in Fluids, v. 53, pp. 573–597.

[32] Zienkiewicz, O.C., Zhu, J.Z., 1987, “A simple error estimator and adaptive

procedure for practical engineering analysis”, Journal for Numerical Methods

in Engineering, v. 24, pp. 337–357.

[33] Hughes, T.J.R., 1995, “Multiscale Phenomena - Greens-Functions, The Dirichlet-

To-Neumann Formulation, Subgrid Scale Models, Bubbles And The Origins

Of Stabilized Methods”, Computer Methods In Applied Mechanics And

Engineering, v. 127 (1-4), pp. 387–401.

[34] Ghia, U., Ghia, K.N., Shin, C.T., 1982, “High-Re solutions for incompressible-flow

using the Navier-Stokes equations and a multigrig method”, Journal of

Computational Physics, v. 48 (3), pp. 387–411.

[35] Norberg, C., 2001, “Flow around a circular cylinder: Aspects of fluctuating lift”,

Journal of Fluids and Structures, v. 15, pp. 459–469.

[36] Baranyi, L., 2003, “Computation of unsteady momentum and heat transfer from a

fixed circular cylinder in laminar flow”, Journal of Computational and Applied

Mechanics, v. 4, pp. 13-25.

[37] Williamson, C.H.K., 1988, “Defining a universal and continuous Strouhal–

Reynolds number relationship for the laminar vortex shedding of a circular

cylinder”, Physics of Fluids, v. 31, pp. 2742–2744.

87

[38] Erturk, E., Corke, T.C., Gökçöl, C., 2005, “Numerical solutions of 2-D steady

incompressible driven cavity flow at high Reynolds numbers”, International

Journal for Numerical Methods in Fluids, v. 48, pp. 747–774.

[39] Föster, C., Wall, W. and Ramm, E., 2008, “Stabilized finite element formulation for

incompressible flow on distorted meshes”, International Journal for Numerical

Methods in Fluids, v. 60, pp. 1103–1126.

[40] Rossa, A.L., Camata J.J., Coutinho, A.L.G.A., 2009, “Stabilized

SUPG/PSPG/LSIC FEM formulation using libMesh”, In: Proceedings of 30th

Iberian-Latin-American Congress on Computational Methods in Engineering –

CILAMCE, Armação dos Búzios.

[41] Rossa, A.L., Camata J.J., Coutinho A.L.G.A., 2010, “Hourglass Control for the

Stabilized Finite Element Solution of Coupled Incompressible Viscous Flows

and Heat Transfer”, In: Proceedings of 9th

Argentinian Congress on

Computational Mechanics – MECOM, and 31st Iberian-Latin-American

Congress on Computational Methods in Engineering – CILAMCE, Buenos

Aires.

[42] Sesini, P.A., Coutinho, A.L.G.A., Cruz A.G.B., Rochinha, F.A., 2009,

“Verification of a Stabilized and Residual-Based Multiscale Finite Elment

Code for Incompressible Flows”, In: Proceedings of 30th

CILAMCE, Armação

dos Búzios.

[43] Tezduyar, T.E., 2001, “Finite Element Methods for Flow Problems with Moving

Boundaries and Interfaces”, Archives of Computational Methods in


[44] Bazilevs, Y., Calo, V.M., Tezduyar, T.E., Hughes, T.J.R., 2007, “YZ

Discontinuity Capturing for Advection-Dominated Processes with Application

to Arterial Drug Delivery”, International Journal for Numerical Methods in

Fluids, v. 54, pp. 593–608.

88

[45] Galeão, A.C. and do Carmo, E.G.D., 1998, “A Consistent Approximate Upwind

Petrov-Galerkin Method for Convection-Dominated Problems”, Computer

Methods in Applied Mechanics and Engineering, v. 68, pp. 83–95.

[46] Smagorinsky, J., 1963, “General Circulation Experiments with Primitive Equations,

part i: the Basic Experiment”, Monthly Weather Review, v. 91, pp. 99–152.

[47] Elias, R.N., Coutinho, A.L.G.A., 2007, “Stabilized Edge-Based Finite Element

Simulation of Free-Surface Flows”, International Journal for Numerical

Methods in Fluids, v. 54, pp. 965–993.

[48] Elias, R.N., Martins, M.A.D., Coutinho, A.L.G.A., 2006, “Parallel Edge-Based

Solution of Viscoplastic Flows with the SUPG/PSPG Formulation”,

Computational Mechanics, v. 38, pp. 365–381.

[49] Valli, A.M.P., Carey, G.F., Coutinho, A. L. G. A., 2005, “Control Strategies for

Timestep Selection in FE Simulation of Incompressible Flows and Coupled

Reaction-Convection-Diffusion Processes”, International Journal for

Numerical Methods in Fluids, v. 47, pp. 201–231.

[50] Elias, R. N., Coutinho, A. L. G. A., Martins, M. A. D., 2006, “Inexact Newton-

Type Methods for the Solution of Steady Incompressible Viscoplastic Flows

with the SUPG/PSPG Finite Element Formulation”, Computer Methods in

Applied Mechanics in Engineering, v. 195, pp. 3145–3167.

[51] Coutinho, A.L.G.A., Martins, M.A.D., Sydenstricker, R.M., Elias, R.N., 2006,

“Performance Comparison of Data-Reordering Algorithms for Sparse Matrix-

Vector Multiplication in Edge-Based Unstructured Grid Computations”.

International Journal for Numerical Methods in Engineering, v. 66, pp. 431–

460.

[52] Martins, M.A.D., Elias, R.N., Coutinho, A.L.G.A., 2007 “Edgepack: A Parallel

Vertex and Node Reordering Package for Optimizing Edge-Based

Computations in Unstructured Grids”, Lecture Notes in Computer Science, v.

4395, pp. 292–304.

89

[53] Elias, R.N., Martins, M.A.D., Coutinho, A.L.G.A., 2005, “Parallel Edge-Based

Inexact Newton Solution of Steady Incompressible 3D Navier-Stokes

Equations”, Lecture Notes in Computer Science, v. 3648, pp. 1237–1245.

[54] Elias, R.N., Gonçalves Jr., M.A., Coutinho, A.L.G.A., Esperança, P.T.T., Martins,

M. A. D. and Ferreira, M. D. A. S., 2009, “Computational Techniques for

Stabilized Edge-Based Finite Element Simulation of Nonlinear Free-Surface

Flows”, Journal of Offshore Mechanics and Arctic Engineering, v. 131, pp.

041103-1-041103-7.

[55] Saad, Y., 1996, Iterative Methods for Sparse Linear Systems, Boston, PWS

Publishing Company.

[56] Shakib, F., Hughes, T. J. R., Johan, Z., 1989, “A Multi-element Group

Preconditioned GMRES Algorithm for Nonsymmetric Systems Arising in

Finite Element Analysis”, Computer Methods in Applied Mechanics in


[57] Lohner, R., Camelli F., 2004, “Dynamic Deactivation for Advection-Dominated

Contaminant Transport”, Communications in Numerical Methods in


[58] Lins, E.F., Elias, R.N., Guerra, G.M., Rochinha, F.A., Coutinho, A.L.G.A., 2009,

“Edge-based Finite Element Implementation of the Residual-based Variational

Multiscale Method”, International Journal for Numerical Methods in Fluids, v.

61, pp. 1–20.

[59] Catabriga, L., 2000, Soluções Implícitas das Equações de Euler Empregando

Estrutura de Dados por Arestas. Tese de D.Sc., COPPE/UFRJ, Rio de Janeiro,

RJ, Brasil.

[60] Elias, R.N., 2007, Estruturas de Dados por Arestas para a Simulação Paralela de

Escoamentos Incompressíveis pelo Método Estabilizado de Elementos Finitos.

Tese de D.Sc., COPPE/UFRJ, Rio de Janeiro, RJ, Brasil.

90

[61] De Sampaio, P.A.B., Gonçalves Jr., M.A, 2010, “A second-order time accurate

finite element method for quasi-incompressible viscous flows”, International

Journal for Numerical Methods in Fluids, DOI: 10.1002/fld.2450.

[62] Ribeiro, F.L.B., Coutinho A.L.G.A., 2005, “Comparison between Element, Edge

and Compressed Storage Schemes for Iterative Solutions in Finite Element

Analysis”, International Journal for Numerical Methods in Engineering, v.

36(4), pp. 659–685.

[63] Elias, R.N., Gonçalves Jr., M.A., Coutinho, A.L.G.A., Esperança, P.T.T., Martins,

M. A. D. and Ferreira, M. D. A. S., 2008, “Wave and Free-Surface Flow

Simulation Using Stabilized Edge-Based Finite Element Methods”. In: 29th

Iberian-Latin-American Congress on Computational Methods in Engineering –

CILAMCE, Búzios.

[64] van Gijzen, M.B., 1995, “Large Scale Finite Element Computations with

GMRESlike Methods on a CRAY Y-MP”, Applied Numerical Mathematics, v.

19, pp. 51–62.

http://onlinelibrary.wiley.com/doi/10.1002/fld.2450/abstract

http://onlinelibrary.wiley.com/doi/10.1002/fld.2450/abstract

91

Apêndices

92

Apêndice A

Neste apêndice mostramos que a discretização temporal do balanço de massa

apresentado na Equação (23) possui precisão de 2ª ordem. A demonstração pode ser

encontrada também em [61].

Usando uma Série de Taylor para au , obtemos:

21 3( )

2

nnn n a aa a

u utu u t O t

t t t (66)

Como temos que:

1/2 2

2

nn n

a a au u ut

t tt t t t

(67)

Podemos reescrever a Equação (66) como:

1/21 3( )

nn n aa a

uu u t O t

t (68)

Multiplicando e dividindo a Equação (68) por e t , respectivamente:

1 1/220 ( )

n n n

a a au u u

O tt t

(69)

Usando a Equação (2) no instante de tempo 1/2n :

1/21

20 ( )

nn n

aba a ab

b b a

pu u uu O t

t x x x (70)

Entretanto, o termo

1/2n

ab

b

uux

é não–linear.

93

Considere que:

1/2 11 21

( )2

n n nn na a a

b b bb b b

u u uu u u O tx x x

(71)

Fazendo b

A u e a

b

uB

x, perceba que:

1/2 1/2

1 1 1/2 1/2

1/2 1/2

1/2 1/2

1 1

2 2 2 2

2 2

n n

n n n n n n

n n

n n

t A t BA B A B A B

t t

t A t BA B

t t

(72)

Multiplicando os termos no lado direito da Equação (72):

21 1

21/2 2

1 1

2 2 2 2 4

( ) ( )2 2 4

n n n n

n

t A t B t A BA B A B AB B A

t t t t

t A t B t A BAB B A AB O t

t t t t

(73)

Trocando as variáveis A e B por suas definições originais, demonstramos a

relação apresentada na Equação (71). Usando a Equação (71) na Equação (70),

obtemos:

1/2 1/21 11 21

0 ( )2

n nn n n nabn na a a a

b bb b b a

pu u u uu u O t

t x x x x (74)

A introdução de 1n n

a a au u u e 1/2 1

2n n

ab ab ab nos conduz a:

1/2

2

1 1

2 2

( )

nabna a a

b bb b b

n nnabn a

bb b a

u u uu u

t x x x

puu O tx x x

(75)

94

Que é exatamente a equação (23). Perceba que a discretização do termo convectivo

é linear.

95

Apêndice B

O Apêndice B apresenta uma formulação de elementos finitos para as equações

incompressíveis de Navier-Stokes. Esta formulação, de precisão temporal de 1ª ordem,

foi apresentada em [17], [24] e [28] e é utilizada para comparação da qualidade dos

resultados e de dados de desempenho com a formulação proposta em implementações

2D e 3D.

Quando implementada no EdgeCFD® ela se apresenta como uma boa opção em

função do procedimento de solução totalmente segregado inerente ao método e,

principalmente, à estrutura de dados por aresta.

Equações Governantes

Consideramos um modelo contínuo para escoamentos viscosos incompressíveis. O

problema é definido sobre um domínio limitado aberto , com contorno , contido no

espaço Euclidiano n-dimensional.

O escoamento é modelado pelas equações de Navier-Stokes incompressíveis. Estas

são escritas em coordenadas Cartesianas usando a convenção de soma, com a =1, 2,..., n

e b=1, 2,...,n:

0a a abb

b b a

u u pu

t x x x (76)

0a

a

x

u

(77)

96

As variáveis dependentes incluem os campos de velocidade e de pressão

representadas por au e p . Note que a tensão viscosa é dada por

ab a b b au x u x , onde é a viscosidade do fluido.

Completamos o modelo introduzindo as condições iniciais e de contorno dos

campos de velocidade e pressão. As condições de contorno de velocidade e tensão são

prescritas através de valores fornecidos sobre partições da fronteira não sobrepostas ua

e ta , de tal forma que ua ta :

uaaa tuu xx,

(78)

taababab ttnp xx,

(79)

onde ab é o delta de Kronecker e bn denota as componentes Cartesianas do vetor

normal ao contorno e que apontam para fora.

As condições de contorno de pressão e velocidade normal à fronteira são associadas

ao balanço de massa. Elas são prescritas por valores fornecidos sobre partições da

fronteira não sobrepostas p e G , de tal forma que p G :

ptpp xx,

(80)

Gbb tGnu xx,

(81)

As Equações (76) e (77) envolvem o gradiente de pressão, mas não a pressão.

Assim, ao menos um valor de referência de pressão deve ser prescrito para que se tenha

um campo de pressão único.

Formulação de Elementos Finitos

As Equações (76) e (77) são discretizadas em relação ao tempo como:

97

a

ab

an

ba F

x

p

x

uu

t

u

(82)

a

n

a

a

a

x

u

x

u

(83)

onde

n

a aba b

b b a

u pF u

x x x (84)

Os subscritos n e 1n indicam o nível de tempo e t é o passo de tempo. As

mudanças das variáveis durante o passo de tempo t são representadas por

1n np p p e 1n n

a a au u u . Os campos de velocidade e pressão no tempo n

são representados por 1 1n n n

a a au u u e 1 1n n np p p ,

respectivamente, onde 0 1 .

A discretização temporal apresentada possui precisão de 1ª ordem.

As equações acima são discretizadas no espaço usando elementos finitos

convencionais da classe 0C . A mesma ordem de interpolação é usada para aproximar

velocidade e pressão.

Consideremos a seguinte discretização espacial das variáveis do problema:

ˆn n

a j aju N u , ˆ n n

j jp N p , â j aju N u , ˆ

j jp N p . Note que jN representa as

funções de forma do elemento finito e as variáveis com o subscrito j são valores

nodais. Usando os campos das variáveis discretizadas, podemos escrever as seguintes

expressões para os resíduos quadrados gerados pela discretização:

ˆ â aR R d (85)

onde o parâmetro são parâmetros de escala a serem definidos posteriormente e âR é

o resíduo da discretização, dado por:

98

a

ab

an

ba

a Fx

p

x

uu

t

uR ˆˆˆ

ˆˆˆ

(86)

Minimizando , dado pela Equação (85), em relação aos ip livres e aos valores

nodais aiu livres, obtemos:

ˆˆ 0 livren ii b a ai

b

NN t u R d u

t x (87)

ˆ 0 livreia i

a

NR d p

x (88)

O parâmetro é escolhido como o objetivo de normalizar (e adimensionalisar) as

funções de peso na Equação (87).

Combinando a Equação (87) e a condição de contorno de tensão, dada pela

Equação (79), obtemos o seguinte balanço discretizado de quantidade de movimento:

livre

0ˆˆˆˆ

ai

ab

n

abab

n

ia

b

in

bi

u

dtnpNdRx

NutN

ta

(89)

Uma equação para atualização da pressão é obtida combinando a Equação (88) e o

balanço discretizado de massa, Equação (83):

livre 0î

a

n

a

a

aia

a

i pdx

u

x

uNdR

x

Nt

(90)

Introduzindo as equações de contorno dadas pela Equação (81) e usando a

identidade de Green, a Equação (90) acima produz:

1

ˆ

ˆ ˆ ˆˆˆ

G

i

a a

n n nnn n ni a ab ab i i

a b b a a

Nt pd

x x

N u ut pu d N d N G G d

x x x x x

(91)

99

Na equação acima ˆn

ab é calculado no interior dos elementos de acordo com

ˆ ˆ ˆn n n

ab a b b au x u x .

Note que a Equação (91) combina o método dos mínimos quadrados da quantidade

de movimento em relação à pressão, Equação (88), com o balanço discretizado de

massa, Equação (83). Por essa razão a chamamos de equação da pressão-continuidade.

Da Equação (89) obtemos a equação governante de variação da velocidade durante

o passo de tempo. Ela é dada por:

ˆˆ ˆ ˆ

ˆ ˆˆ ˆ

ˆ ˆˆ ˆ ˆ

ta

n ni ai c a b

c b

nn ni a

i c b

c b a

n nn n n ni ab i ic ab i a

c b a b a

N uN t u u t u d

x t x

N u pN t u u d

x x x

N N Npt u d d p d N t d

x x x x x

(92)

Da mesma forma como no método apresentado nesta tese, os termos multiplicados

por t nas Equações (91) e (92) são responsáveis pelo controle das oscilações espaciais

em escoamentos com convecção dominante, e por estabilizar o cálculo,

independentemente das restrições de Babuška-Brezzi na escolha dos espaços de

interpolação para velocidade e pressão. Em particular, a mesma ordem de interpolação

para todas as variáveis adotadas aqui se torna possível através de uma escolha adequada

de t . É importante observar que ao contrário de serem propostos a priori, os termos de

estabilização no presente método surgem naturalmente do método dos mínimos

quadrados da discretização temporal com respeito aos graus de liberdade (valores

nodais livres).

Utilizando o esquema de sincronização adotado na presente tese, podemos

reescrever as Equações (91) e (92) para 3 dimensões espaciais como:

' ' '

pM* p* f (93)

' ' '

uu uA * * f (94)

100

' ' '

vv vA * * f (95)

' ' '

ww wA * * f (96)

As matrizes resultantes são simétricas positivas definidas. O procedimento de

solução adotado para o método descrito no Apêndice B é totalmente segregado.

Primeiro é resolvida a pressão e depois as componentes de velocidade, uma a uma.

Documents

UMA FORMULAÇÃO ESTÁVEL DE ELEMENTOS FINITOS …objdig.ufrj.br/60/teses/coppe_d/MiltonAlvesGoncalvesJunior.pdf · Milton Alves Gonçalves Junior Tese de Doutorado apresentada ao