UNIVERSIDADE FEDERAL DE PERNAMBUCO CENTRO DE TECNOLOGIA E GEOCIÊNCIAS

PROGRAMA DE RECURSOS HUMANOS DA ANP – PRH 26

TESE DE DOUTORAMENTO

Um Sistema Computacional Utilizando uma Formulação de Passo Fracionado e o Método

dos Elementos Finitos por Arestas para a Análise de Escoamentos Incompressíveis Tridimen-

sionais usando Computação Paralela

Bolsista DSc. AlessandroRomario Echevarria Antunes

Orientador: Prof. Paulo Roberto Maciel Lyra (PhD)

Co-Orientador: Prof. Ramiro Brito Willmersdorf (PhD)

RECIFE – 2008

ii

Universidade Federal de Pernambuco Centro de Tecnologia e Geociências Departamento de Engenharia Civil

UM SISTEMA COMPUTACIONAL UTILIZANDO UMA FORMULAÇÃO DE PASSO FRACIONADO E O MÉTODO DOS ELEMENTOS FINITOS POR A-RESTAS PARA A ANÁLISE DE ESCOAMENTOS INCOMPRESSÍVEIS TRI-

DIMENSIONAIS USANDO COMPUTAÇÃO PARALELA

Alessandro Romário Echevarria Antunes

TESE SUBMETIDA AO CORPO DOCENTE DO CURSO DE PÓS-GRADUAÇÃO EM ENGENHARIA CIVIL DA UNIVERSIDADE FEDERAL DE PERNAMBUCO COMO PARTE DOS REQUISITOS NECESSÁRIOS À OB-TENÇÃO DO GRAU DE DOUTOR EM CIÊNCIAS DE ENGENHARIA CIVIL.

ÁREA DE CONCENTRAÇÃO: ESTRUTURAS

Paulo Roberto Maciel Lyra, PhD

(Orientador)

Recife, Pernambuco, Brasil. ©Alessandro Romário Echevarria Antunes, Dezembro de 2008.

iii

A636u Antunes, Alessandro Romário Echevarria.

Um sistema computacional utilizando uma formulação de passo fracionado e o método dos elementos finitos por arestas para a análise de escoamentos incompressíveis tridimensionais usando computação paralela / Alessandro Romário Echevarria Antunes. - Recife: O Autor, 2008.

x, 166 folhas, il : tabs. grafs., figs Tese (Doutorado) – Universidade Federal de Pernambuco. CTG.

Programa de Pós-Graduação em Engenharia Civil. Inclui Referências Bibliográficas e Apêndices. 1. Engenharia Civil. 2. Equações de Navier-Stokes Incompressí-

veis. 3. Método dos Elementos Finitos. 4. Método do Passo Fracionado. 5. Computação Paralela I. Título

UFPE 624 BCTG/ 2009-061

iv

v

DEDICATÓRIA

Aos meus pais, João Batista Melo Antunes e Maria das Graças Eche-varria Antunes, pela educação, apoio e exemplos de força e perseverança im-prescindíveis na conquista desta realização, sem os quais este caminho não te-ria sido percorrido.

Aos meus irmãos, que com carinho e generosidade foram essenciais para

a conclusão desta jornada. À Fernanda Maria Bezerra de Mello, minha amada esposa, que acom-

panhou esta jornada com paciência e confiança, tendo sempre a certeza de que os nossos objetivos seriam alcançados.

vi

AGRADECIMENTOS

A todos aqueles que de alguma forma, direta ou indiretamente, contri-buíram para a realização deste trabalho.

A todos os amigos que fizeram parte desta jornada, dentre eles, Erwin,

Demétrio, Jacek, Jorge, Yuri, Saul, Sílvio Alan, Sílvio Leandro, Seniz, Flávio e tantos outros.

Ao amigo Darlan Karlo Elisiário de Carvalho, pelas valiosas e inúme-

ras conversas científicas/acadêmicas que tivemos, e principalmente pelos as-suntos aleatórios tratados durante nossos “café com cultura” diários.

Ao amigo Rogério Soares da Silva, pelas valiosas discussões sobre

Computação Paralela que muito auxiliaram no desenvolvimento deste traba-lho e pela disposição em auxiliar e contribuir para o fechamento deste traba-lho.

Ao amigo Daniel Ventura, que acompanhou grande parte deste traba-

lho, e contribuiu com o desenvolvimento da estrutura do programa desenvol-vido.

À amiga Eliane Alves da Silva, que desde o início da minha chegada à

Recife, sempre foi muito prestativa nas questões burocráticas da pós gradua-ção.

À amiga Roesane Teixeira de Lima que sempre esteve disposta a aju-

dar nos assuntos referentes aos auxílios financeiros.

Aos professores Paulo R. Maciel Lyra e Ramiro Brito Willmersdorf pelos valiosos ensinamentos imprescindíveis para a elaboração desta tese.

Aos professores, funcionários e amigos do Departamento de Enge-

nharia Mecânica e do Programa de Pós-Graduação em Engenharia Civil da Universidade Federal de Pernambuco.

Ao Núcleo de Atendimento em Computação de Alto Desempenho

(NACAD) da COPPE/UFRJ pelo suporte computacional durante o desenvol-vimento deste trabalho.

Ao Programa de Recursos Humanos da Agência Nacional do Petróleo

(PRH-26) e ao Conselho Nacional de Pesquisa (CNPq) pelo suporte financeiro durante o período de realização desta tese.

Aos meus sogros Ayrton Bezerra de Mello e Nair Mello por terem me

acolhido como um filho, tendo sido um porto seguro nos momentos difíceis.

vii

RESUMO O objetivo do presente trabalho é desenvolver um sistema computacional capaz de simu-

lar numericamente escoamentos laminares de fluidos governados pelas Equações de Navier-Stokes escritas na sua forma incompressível, com variáveis primitivas em domínios tridimen-sionais. Para tal, o fluido é considerado viscoso e incompressível, e no domínio fluido empre-ga-se uma formulação com estabilização da convecção e da pressão do tipo SU (Streamline Upwind), o Método dos Elementos Finitos e o Método do Passo Fracional para realizar o fra-cionamento dos sistemas de equações algébricas no tempo, resultando em uma formulação monolítica, e com estabilidade de segunda ordem no tempo. Foi utilizada uma estrutura de dados baseada nas arestas dos elementos tetraédricos lineares, que apresentam como principal vantagem o ganho de tempo computacional, visto que cada termo pode ser obtido via pré-processamento na forma de coeficientes de contribuição de cada aresta. O programa foi escri-to em linguagem FORTRAN 90. Foram aplicados os conceitos de computação paralela em computadores com memória distribuída para permitir a análise de problemas de grande porte, comuns na Dinâmica dos Fluidos Computacional. Utilizou-se a decomposição de domínio via PARMETIS (Parallel Graph Partitioning and Sparse Matrix Ordering) para particionar a ma-lha computacional e MPI (Message Passing Interface) para viabilizar a comunicação paralela. Também foi utilizado o pacote PETSc (Portable Extensible Toolkit for Scientific Computati-ons), que possui uma vasta gama de estruturas paralelas úteis no desenvolvimento de algorit-mos paralelos para computação distribuída. Diversas aplicações foram analisadas de forma a verificar e validar o sistema computacional desenvolvido. Os resultados obtidos nas análises estão em conformidade com dados experimentais, teóricos e numéricos disponíveis na litera-tura. Palavras Chave: Equações de Navier-Stokes Incompressíveis, Método dos Elementos Finitos, Método do Passo Fracionado, Computação Paralela

viii

ABSTRACT

The main objective of the present work is to develop a computational system capable of numerically simulate three-dimensional, laminar incompressible Navier-Sokes fluid flow problems written in primitive variable form, using parallel computing. For this purpose, the fluid is considered viscous and incompressible. At the fluid domain we are using a Finite Element Method formulation with stabilization for pressure and advective terms and the Fractional Step Method in ordering to obtain a second ordering stable monolithic formulation. The edge-based data structure is used. It’s main advantage refer to the reduction of the computational time requested as all coefficients are calculated in a pre-processing step. A FORTRAN 90 program is developed and parallel computing concepts are applied to allow the computation of large scale tridimensional applications frequent in Computational Fluid Dynamics. Domain decomposition is obtained via PARMETIS (Parallel Graph Partitioning and Sparse Matrix Ordering) and MPI (Message Passing Interface) is the standard communicator. We are using the PETSc (Portable Extensible Toolkit for Scientific Computations) toolkit, which is an efficient set of solvers, data structures (sequential and parallel), pre-conditioners, etc, need in parallel computing developments. Several applications were analyzed in ordering to verify and validate the computational system developed. The obtained results are in agreement with the experimental, theoretical and numerical data available in the literature. Keywords: Incompressible Navier-Stokes Equations, Finite Element Method, Fractional-Step Method, Parallel Computing

ix

SUMÁRIO

1. INTRODUÇÃO..................................................................................................... 1

1.1 Motivação e Considerações Gerais: Escoamentos Incompressíveis................. 1

1.2 Motivação e Considerações Gerais: Método dos Elementos Finitos.................3

1.3 Motivação e Considerações Gerais: Computação Paralela................................5

1.4 Objetivos e Organização Geral da Tese.............................................................7

2. FORMULAÇÃO MATEMÁTICO-NUMÉRICA...............................................9

2.1 Introdução...........................................................................................................9

2.2 Equações de Navier-Stokes Incompressíveis.....................................................9

2.2.1 Condições de Contorno para as Equações de Navier-Stokes...................12

2.2.2 Condições de Contorno de Neumann para a Velocidade nas Equações

de Navier-Stokes.......................................................................................17

2.2.3 Dificuldades Associadas aos Fluxos Incompressíveis..............................20

2.2.4 Formulação de Elementos Finitos Estabilizada........................................25

2.3 Método do Passo Fracionado............................................................................. 25

2.4 Formulação Segregada Final..............................................................................45

2.5 Termos de SUPG para Estabilização Espacial...................................................47

2.5.1 Estabilização Espacial da Equação de Poisson para a Pressão..................47

2.5.2 Estabilização Espacial das Equações de Quantidade de Movimento........48

2.6 Parâmetros Importantes......................................................................................48

2.6.1 Termo de Preservação da Monotonicidade................................................48

2.6.2 Passo de Tempo para Estabilidade Temporal............................................49

2.6.2.1 Estabilização Temporal das Equações de Quantidade de Movimen-

to.......................................................................................................50

2.6.2.2 Estabilização Temporal da Equação de Poisson para a Pressão.......50

3. MODELAGEM NUMÉRICO-COMPUTACIONAL.........................................52

3.1 Introdução e Considerações Gerais....................................................................52

3.2 Discretização Espacial via Método dos Elementos Finitos............................... 54

3.3 Implementação Computacional..........................................................................55

x

4. COMPUTAÇÃO PARALELA..............................................................................73

4.1 Conceitos Fundamentais.....................................................................................73

4.2 Aspectos Computacionais...................................................................................74

4.2.1 Particionamento da Malha.........................................................................75

4.2.2 Gerenciamento da Comunicação Paralela via MPI...................................81

4.2.3 Utilização do Tollkit PETSC.....................................................................82

4.3 Programa NS_INCOMP3D................................................................................85

4.3.1 Considerações Gerais.................................................................................85

5. APLICAÇÕES NUMÉRICAS...............................................................................92

5.1 Escoamento em Tubo........................................................................................ .93

5.2 Backfacing Step........................................……………..................................... .96

5.2.1 Número de Reynolds 100...........................................................................97

5.2.2 Número de Reynolds 1000.........................................................................99

5.3 Escoamento sobre Cilindro Circular.................................................................100

5.3.1 Número de Reynolds 40...........................................................................101

5.3.2 Número de Reynolds 100.........................................................................102

5.4 Escoamento em Canal com Cavidade...............................................................105

5.5 Escoamento em Canal com Expansão Brusca...................................................110

5.6 Estudos de Desempenho do Programa Paralelo................................................117

6. CONCLUSÕES E SUGESTÕES PARA TRABALHOS FUTUROS..............130

6.1 Conclusões........................................................................................................130

6.2 Trabalhos Futuros.............................................................................................132

- Interação Fluido-Estrutura..............................................................................132

- Modelos de Turbulência.................................................................................133

- Fluidos Não-Newtonianos..............................................................................134

- Adaptação de Malha.......................................................................................134

- Otimização Estrutural com Interação Fluido-Estrutura..................................135

REFERÊNCIAS BIBLIOGRÁFICAS......................................................................136

Apêndice A: Matriz de Adjacências..........................................................................144 Apêndice B: Mapeamento dos Nós............................................................................153

xi

Apêndice C: Comparação entre Estruturas de Dados por Elementos, Aresta Pura e Arestas CSR.......................................................................................................157

Apêndice D: Resultados Preliminares Obtidos para Difusão Pura e Convecção Difu-são ransiente.......................................................................................................162

1

Capítulo 1

Introdução

Neste capítulo serão expostos alguns dos principais conceitos que definem os conteúdos

que são objetos de estudo neste trabalho, bem como as suas características intrínsecas que tor-

nam o tema desafiador e ao final apresentamos os principais objetivos, assim como a organi-

zação desta tese. Dentre os principais temas envolvidos no trabalho, primeiramente será

apresentada uma breve explanação sobre os escoamentos que serão abordados, evidenciando

suas características e as hipóteses que foram adotadas. A seguir é apresentado o método nu-

mérico utilizado para discretizar as equações governantes dos fenômenos físicos, evidencian-

do as principais dificuldades associadas à aplicação deste método. Por fim, alguns conceitos

sobre computação paralela são abordados, identificando aqueles que foram utilizados neste

trabalho.

1.1 Motivação e Considerações Gerais: Escoamentos Incompressíveis

Escoamentos fluidos onde as alterações na massa específica são desprezíveis podem ser

considerados escoamentos incompressíveis. Em geral estes tipos de escoamentos são caracte-

rísticos dos líquidos, porém mesmo os gases quando em escoamentos a baixas velocidades, ou

com o número de Mach menor ou igual a 0.3 e sem a presença de grandes gradientes de tem-

peratura, podem ser modelados como escoamentos incompressíveis. Muitos exemplos de es-

2

coamentos incompressíveis podem ser citados, dentre eles, fluxo de fluidos ao redor de estru-

turas submarinas, escoamento de ar sobre estruturas como edifícios, pontes e torres, escoa-

mento sanguíneo nas artérias humanas, escoamentos líquidos envolvendo dispersão de

poluentes, escoamentos internos em dutos e canais, etc. Estes escoamentos são governados

pelas equações de Navier-Stokes com a restrição de incompressibilidade, que possuem solu-

ção analítica apenas quando várias hipóteses simplificadoras são adicionadas à modelagem do

problema, por exemplo, escoamentos permanentes sobre geometrias simples. Dado isso,

quando escoamentos sobre geometrias complexas, com fluxos transientes são considerados,

estas equações não possuem solução analítica, sendo necessária utilização de técnicas numéri-

cas para a obtenção de soluções, ou o uso de experimentação em laboratório, sendo esta últi-

ma opção, em geral, menos flexível e mais cara em termos de custo e tempo.

Ao longo dos anos, muitos esquemas numéricos foram desenvolvidos para tratar as equa-

ções de Navier-Stokes, dentre estes esquemas citamos: Esquemas Monolíticos, Esquemas de

Pré-Condicionamento das Equações de Navier-Stokes Compressíveis, Formulações de Pena-

lidade e Esquemas Baseados no Fracionamento das Equações [Codina, 2001; Chorin, 1967;

Costa, 2004; Gresho e Sani, 1998; Temam, 1969]. Todos os esquemas citados apresentam

vantagens e desvantagens, tais como: os Esquemas Monolíticos satisfazem a restrição de in-

compressibilidade com o campo de velocidades obtido pelas equações de Quantidade de Mo-

vimento, mas forçam a utilização do mesmo método de solução para o sistema de equações,

desconsiderando o caráter diferente de cada equação; os Esquemas de Pseudo-

Compressibilidade adicionam um termo de variação no tempo da massa específica e no limite

da incompressibilidade fazem este termo tender a zero, satisfazendo assim a restrição de in-

compressibilidade, porém estes esquemas são fortemente dependentes de matrizes de condi-

cionamento, muitas vezes de obtenção complexa, baixando, assim, a eficiência computacional

e numérica destes esquemas; esquemas que envolvem formulações de penalidade são geral-

mente complicados, por envolverem parâmetros livres na formulação; esquemas de fraciona-

mento das equações resultam em formulações acuradas no tempo e com características de

estabilidade para a pressão, porém adicionam novas variáveis ao sistema, introduzindo erros

devido ao fracionamento, e com isso introduzem novas dificuldades a serem tratadas. Como

visto, todos esses esquemas são capazes de fornecer formulações estáveis para as equações de

Navier-Stokes [Henriksen e Holmem, 2002], ficando sob critério do pesquisador, definir o

método que vai utilizar. Neste trabalho optou-se por utilizar um esquema de fracionamento

das equações de Navier-Stokes, conhecido como “Fracitonal Step Method” [Chang et al,

2002; Codina, 2000, 2001 e 2002; Codina e Soto, 2003; Donea et al., 1982; Guermond et al.,

3

2006; Henriksen e Holmem, 2002; Perot, 1993; Löhner, 2004; Löhner et al., 2006], que pode

ser visto de duas maneiras diferentes, utilizando conceitos físicos baseados na decomposição

de Helmholtz [Gresho e Sani, 1998, vol 2], e utilizando conceitos algébricos baseados na fato-

ração LU do sistema algébrico obtido a partir das equações de Navier-Stokes. Esse esquema

fornece uma formulação implícita, monolítica (a cada passo de tempo), de segunda ordem, e

com estabilidade para a pressão. Para a discretização espacial foi aplicado o Método dos Ele-

mentos Finitos.

1.2 Motivação e Considerações Gerais: Método dos Elementos Finitos

O Método dos Elementos Finitos foi introduzido na Dinâmica dos Fluidos Computacional

no início da década 70, tendo sido utilizado no final da década de 50 e início da década de 60

pela indústria aeronáutica, que estava em crescente expansão no período logo após a Segunda

Guerra Mundial, na análise de tensões sobre a asa de aeronaves. Muitas qualidades do método

chamaram a atenção dos pesquisadores, dentre estas qualidades podemos citar a sólida base

matemática, a possibilidade de utilização através de discretização com malhas não estrutura-

das, que podem ser aplicadas a domínios com modelamentos geométricos complexos, o tra-

tamento consistente das condições de contorno envolvendo diferenciais, e a possibilidade de

ser programado de forma geral e flexível. Aproximações de Elementos Finitos padrão são ba-

seados no Método dos Resíduos Ponderados [Zienkiewicz e Morgan, 1983, Zienliewicz e Ta-

ylor, 1988] (ou na formulação variacional correspondente) através da formulação de Galerkin.

Essa formulação já tinha mostrado grande sucesso em aplicações na Mecânica dos Sólidos e

também em outras áreas, como fenômenos envolvendo condução de Calor, cujas equações

governantes são equações diferenciais parciais que apresentam termos de segunda derivada,

caracterizando equações elípticas. A razão para o sucesso dessa formulação quando aplicada a

problemas governados por equações com caráter difusivo é que estes problemas são governa-

dos por equações diferenciais elípticas e/ou parabólicas que possuem termos auto-adjuntos, e

desta forma o Método de Galerkin produz matrizes simétricas. Nesse caso, a diferença entre a

solução de elementos finitos e a solução exata é minimizada com respeito à norma de energia,

e essa é a denominada “propriedade da melhor aproximação” do Método de Galerkin. Nesses

casos, existe um funcional quadrático que quando minimizado corresponde a satisfazer as e-

quações diferenciais governantes destes problemas. Por exemplo, em elasticidade linear, a po-

sição de equilíbrio de uma estrutura corresponde ao mínimo de um funcional quadrático que

4

expressa a energia potencial total do sistema. Similarmente, em problemas de condução de ca-

lor em regime permanente, o equilíbrio térmico resultante da solução das Equações de Lapla-

ce ou Poisson corresponde ao mínimo de um funcional quadrático expresso em termos do

fluxo de calor, que representa fisicamente a energia térmica total do sistema.

O sucesso do Método dos Elementos Finitos de Galerkin para problemas de Mecânica Es-

trutural, Condução de Calor e outros problemas do tipo potencial, motivou a aplicação do mé-

todo na simulação de problemas na Dinâmica dos Fluidos Computacional [Anderson, 1995;

Donea, 2003]. Esperava-se obter os mesmos resultados obtidos nos fenômenos para os quais o

método vinha sendo aplicado. Atualmente sabemos que este foi um ponto de vista excessiva-

mente otimista, visto que fenômenos envolvendo escoamentos fluidos apresentam caracterís-

ticas diferentes, e talvez a principal delas seja o fato de que existem termos convectivos nas

equações governantes, que obviamente são termos não-simétricos, e que fazem com que o

Método dos Elementos Finitos de Galerkin perca a propriedade da melhor aproximação, pois

estes termos não são auto-adjuntos. Na prática, soluções para problemas que apresentam ca-

racterísticas de transporte convectivo dominante utilizando o Método dos Elementos Finitos

de Galerkin apresentam soluções com oscilações espúrias e/ou instabilidades [Brooks e Hu-

ghes, 1982; Catabriga, 2000; Catabriga e Coutinho, 2002; De Sampaio, 1991, Zienkiewicz e

Taylor, 1988]. Esses infortúnios podem ser resolvidos aplicando-se malhas com alto grau de

refinamento e com passos de tempo muito pequenos, quando considerando-se esquemas de

avanço explícito no tempo. Porém essas medidas tornam o problema computacionalmente

impraticável. Estas dificuldades motivaram o desenvolvimento de formulações alternativas

que possibilitassem a aplicação do Método dos Elementos Finitos sem a necessidade deste

grau de refinamento da malha e aplicação de passos de tempo computacionalmente inviáveis.

Essas formulações alternativas ficaram conhecidas como formulações estabilizadas e repre-

sentaram um grande avanço para a utilização do Método dos Elementos Finitos na Dinâmica

dos Fluidos Computacional. Neste trabalho foi utilizada uma formulação de Elementos Fini-

tos Mistos e Estabilizada para tratar as equações de Navier-Stokes Incompressíveis. A estabi-

lização é do tipo SUPG (Streamline Upwind Petrov Galerkin) [Almeida e Galeão, 1996,

Brooks e Hughes, 1982; Catabriga, 2000; Catabriga e Coutinho, 2002] que tem por finalidade

adicionar difusão suficiente para evitar o surgimento de oscilações espúrias em escoamentos

que apresentam transporte convectivo dominante.

Em termos de implementação computacional do Método dos Elementos Finitos na Dinâ-

mica dos Fluidos Computacional, classicamente utiliza-se uma estrutura de dados baseada nos

elementos que compõem a malha computacional. Este tipo de estrutura de dados é muito efi-

5

caz, porém necessita que sejam efetuados “loops” em todos os elementos da malha a cada

passo de tempo e as contribuições sejam acumuladas considerando-se os nós dos elementos.

Desta forma tem-se grande quantidade de endereçamento indireto, visto que, por exemplo, pa-

ra elementos tetraédricos lineares, cada elemento precisa identificar quatro nós pertencentes a

este elemento. Uma estrutura de dados alternativa, baseada nas arestas da malha [Peraire et

al., 1993, Löhner, 2001, Löhner e Galle, 2002, Lyra et al., 1993], começou a ser utilizada para

soluções explícitas de escoamentos compressíveis utilizando-se o Método dos Elementos Fi-

nitos em malhas não estruturadas, que teve por inspiração o uso semelhante no contexto do

Método dos Volumes Finitos [Barth, T. J., 1992]. A vantagem de uma estrutura de dados ba-

seada nas arestas dos elementos reside no fato de que propicia uma grande redução nas opera-

ções de endereçamento indireto, na quantidade de memória requerida, e no número de

operações de ponto flutuante necessárias. As operações de endereçamento indireto ainda po-

dem ser reduzidas se estruturas de arestas mais elaboradas são utilizadas, tais como “stars, su-

per-edges, chains” [Löhner, 2001], porém a complexidade algoritmica associada a essas

alternativas reduzem a atratividade das mesmas. Neste contexto, a estrutura de dados baseada

nas arestas pode ser interpretada como a representação de um grafo nodal de malhas tetraédri-

cas em domínios tridimensionais. Outra vantagem diz respeito à possibilidade de lidar com

qualquer tipo de elementos na etapa de processamento, uma vez os dados estando pré-

processados, além de possibilitar de forma aproximada estender técnicas eminentemente uni-

dimensionais num contexto multidimensional.

1.3 Motivação e Considerações Gerais: Computação Paralela

Até a metade do século XX a pesquisa e resolução de problemas científicos e tecnológicos

baseavam-se em duas metodologias: a teórica e a experimental. Com o advento do computa-

dor surge uma terceira: a metodologia computacional, que abre perspectivas próprias de de-

senvolvimento, complementando e suplementando as metodologias de pesquisa tradicionais.

Não é surpresa que o desenvolvimento do computador, que desde seu surgimento até os dias

de hoje vem apresentando impressionantes ganhos nas capacidades de velocidade e de memó-

ria, resulte em grande impacto nos processos de desenvolvimento científico e tecnológico. Na

área de mecânica computacional, por exemplo, o advento da supercomputação, atingida atra-

vés de processamento paralelo, permite que projetos sejam quase que completamente desen-

6

volvidos no computador: protótipos virtuais antes dos protótipos reais. É desnecessário frisar

a importância deste ganho computacional para atender à velocidade de inovações tecnológicas

que sociedade e mercado demandam. Atualmente diversos grupos de pesquisa, espalhados pe-

lo mundo, utilizam “supercomputadores” no auxílio ao desenvolvimento de pesquisas e proje-

tos de alto impacto econômico e tecnológico, porém, o alto custo de instalação e manutenção

destes “supercomputadores” vem dificultando a ampliação de novas instalações e também

diminuindo a vida útil dos atuais. Para sobrepor esta dificuldade, surgiu uma alternativa com

custos reduzidos, capaz de proporcionar ganhos de desempenho utilizando computadores pes-

soais. Esta alternativa ficou conhecida como “cluster de PC’s”, ou agrupamento de computa-

dores pessoais (PC’s) conectados em rede de altíssima velocidade, com um conjunto de

processadores podendo ser utilizados para resolver problemas de grande porte. O conceito ba-

seia-se no fato de que muitas das tarefas executadas em computador podem ser divididas em

partes e cada parte pode ser resolvida por um processador diferente. Com este novo conceito,

houve um espalhamento global de grupos instalando seu próprio “cluster”, e desta forma di-

fundindo os conhecimentos e desenvolvimentos nesta área. Essa alternativa foi tão bem aceita

que os maiores computadores paralelos, em termos de velocidade e capacidade de processa-

mento atualmente são “clusters” de computadores pessoais, tendo o “RoadRunner” da IBM

instalado no “Departament of Energy’s Los Alamos National Laboratory” que possui 129.600

núcleos, alcançou a performance de 1.026 petaflop/s, sendo inclusive o primeiro a alcançar a

marca dos petaflop/s. 1 petaflop/s (peta é um jargão da computação que significa o número 1

seguido de 15 zeros) equivale a 1 quatrilhão de operações de ponto flutuante por segundo. Há

apenas um ano antes, o primeiro lugar pertencia ao hoje quarto colocado “BlueGene/L” da

IBM, que possuía 131.072 processadores alcançando em média 280.6 teraflops com picos de

340 teraflops, e que hoje possui 212.992 núcleos. O “BlueGene” está instalado no LLNL

(Lawrence Livermore National Laboratory) (www.top500.org, acessado em 15/11/2008).

Neste contexto, o objetivo principal deste trabalho foi desenvolver um programa computacio-

nal baseado na utilização de “clusters” de computadores pessoais de forma a analisar proble-

mas da Dinâmica dos Fluidos Computacional de grande porte envolvendo escoamentos

incompressíveis sobre domínios tridimensionais. Ainda, para determinados tipos de aplica-

ções na atualidade, investiga-se o uso do chamado “Grid Computing”, onde diversos “clus-

ters” são usados para resolver uma única tarefa.

7

1.4 Objetivos e Organização Geral da Tese

O objetivo central desse trabalho foi o estudo de problemas de escoamentos fluidos in-

compressíveis, e para tal foi desenvolvido um sistema computacional para tratar as Equações

de Navier-Stokes Incompressíveis, em variáveis primitivas, em domínios tridimensionais, a-

través do Método dos Elementos Finitos, utilizando uma malha computacional não estrutura-

da de elementos tetraédricos lineares e uma estrutura de dados baseada nas arestas dos

elementos e implementado em computadores com memória distribuída. A seguir apresenta-se

uma breve descrição dos passos seguidos neste trabalho:

• Estudo analítico das equações governantes do problema fluido, descritas em uma

formulação Euleriana.

• Estudo das técnicas numéricas e computacionais necessárias para a viabilização das

soluções, por exemplo: Método dos Elementos Finitos, Estrutura de Dados por Ares-

tas, Algoritmos de Solução, Pré-Condicionadores.

• Obtenção de modelos numéricos simplificados visando o aprofundamento dos con-

ceitos de estrutura de dados e Fortran 90, por exemplo: Convecção-Difusão 2D.

• Obtenção do modelo discreto das Equações de Navier-Stokes Incompressíveis via

Método dos Elementos Finitos por Arestas: obtenção dos coeficientes das arestas,

montagem das matrizes e dos sistemas de equações algébricas resultantes.

• Estudo das técnicas de programação paralela, tais como: particionamento de malha,

comunicação paralela, vetores e matrizes definidas em ambientes paralelos.

• Desenvolvimento do programa principal de Navier-Stokes Incompressível utilizando

Fortran 90 e definido para computadores com memória distribuída.

• Obtenção de resultados necessários para verificação e validação da ferramenta com-

putacional desenvolvida.

Desta forma, ao final deste trabalho obteve-se um programa computacional capaz de tra-

tar escoamentos fluidos incompressíveis em domínios tridimensionais, utilizando computado-

res com memória distribuída através de uma formulação numérica que permite uma solução

monolítica a cada passo de tempo, com precisão temporal de segunda ordem.

8

O restante desta tese está organizado da seguinte maneira: no capítulo 2 são apresenta-

das as equações governantes e a formulação matemática utilizadas juntamente com os princi-

pais parâmetros de estabilização necessários à formulação. No capítulo 3, essa formulação

matemática é discretizada espacialmente utilizando-se o Método dos Elementos Finitos dando

origem à formulação numérica, e também é apresentada a estrutura de dados baseada em ares-

tas que foi empregada, e os coeficientes associados às arestas para cada termo discreto das

equações contínuas. No capítulo 4, alguns tópicos sobre computação paralela são definidos,

dentre eles, a partição de domínio e suas conseqüências, a comunicação entre processadores

em ambientes paralelos gerenciada via MPI, as ferramentas paralelas utilizadas via PETSc e o

programa principal com detalhes sobre as rotinas desenvolvidas. No capítulo 5, alguns exem-

plos numéricos são mostrados para verificação e validação do código desenvolvido com res-

peito à eficiência na captura dos fenômenos envolvidos e à precisão temporal e/ou espacial da

formulação utilizada. No capítulo 6, as principais conclusões são apresentadas, e finalizando

esta etapa de trabalho, citamos novas investigações e desenvolvimentos que podem ser efetu-

ados como continuidade da presente tese. Por último, apresentamos uma lista de referências

bibliográficas utilizadas e alguns apêndices descrevendo alguns tópicos complementares de

interesse e que surgiram ao longo do desenvolvimento desta tese.

9

Capítulo 2

Formulação Matemático-Numérica

2.1 Introdução

Neste trabalho estamos interessados em obter soluções para problemas de escoamentos

fluidos governados pelas equações de Navier-Stokes, e neste capítulo será apresentada a for-

mulação matemático-numérica utilizada. Primeiramente o sistema de equações de Navier-

Stokes é apresentado, a seguir procede-se um fracionamento deste sistema para obter estabili-

zação temporal para a pressão e é obtida uma formulação fracionada com características de

estabilidade para a pressão e a velocidade. Também são apresentados nesta seção os termos de

estabilização espacial para o termo convectivo, e os parâmetros importantes na obtenção desta

estabilização.

2.2 Equações de Navier-Stokes Incompressíveis

Se for necessário um grande incremento de pressão para produzir um pequeno decrés-

cimo de volume em um fluido, é razoável assumir total incompressibilidade do fluido. Esta

condição é normalmente adotada para a simulação de líquidos, porém, escoamentos de gases

10

com velocidades baixas quando comparadas com a velocidade local de propagação do som e

com pequenas variações de temperatura também podem ser assumidos como incompressíveis

(Curie, 1974 e White, 1991).

A hipótese de incompressibilidade requer que a massa específica ρ seja constante, e a

equação da continuidade, pode ser escrita, usando notação indicial, como,

0=∂

∂

j

j

x

u (2.1)

que representa uma restrição no campo de velocidades, que obrigatoriamente deve ser livre de

divergente. Nessa equação, ju representa a componente da velocidade na direção j , e jx re-

presenta a componente espacial na direção j .

Fazendo uso da Eq.(2.1), as componentes do tensor de tensões cisalhantes ijτ , conside-

rando-se um fluido Newtoniano, podem ser escritas como,

∂

∂+

∂

∂=

i

j

j

i

ijx

u

x

uµτ (2.2)

onde µ é a viscosidade dinâmica do fluido.

Nas situações de interesse no presente estudo, onde desconsideramos gradientes de tem-

peratura, a viscosidade dinâmica pode ser considerada constante, e isso juntamente com a e-

quação da continuidade permite uma simplificação da equação da conservação da quantidade

de movimento. Desse modo uma grande simplificação pode ser feita no grupo de equações de

Navier-Stokes sempre que os efeitos de transferência de calor não são relevantes para a análi-

se do escoamento em questão, visto que a equação da conservação da energia fica desacopla-

da das equações da continuidade e da conservação da quantidade de movimento, isto é, os

campos de velocidade e pressão podem ser obtidos sem a necessidade de fazer referência à

equação da conservação da energia. O campo de temperatura pode ser determinado posteri-

ormente, se necessário, a partir do campo de velocidade já obtido. Considerando-se ainda um

escoamento monofásico, com viscosidade constante, essa equação pode então ser escrita co-

mo,

11

( ) ( ) ( )i

jj

i

j

ij

j

jii fxx

u

x

p

x

uu

t

uρµ

δρρ+

∂

∂=

∂

∂+

∂

∂+

∂

∂ 2

(2.3)

onde ρ é a massa específica do fluido, p é a pressão, t é o tempo, ijδ é o delta de Kronecker

e if é a componente do vetor de carregamentos externos.

As Eqs. (2.1) a (2.3) podem ser utilizadas para a obtenção dos campos de velocidade e

pressão conhecendo-se a massa específica e viscosidade do fluido e condições de contorno e

iniciais. Estas equações ainda podem ser escritas em uma forma não conservativa, como uma

função das variáveis primitivas ρ , ju e p , como segue,

0=∂

∂

j

j

x

u (2.4)

i

jj

i

jj

i

j

i fxx

u

x

p

x

uu

t

uρµρρ +

∂

∂=

∂

∂+

∂

∂+

∂

∂ 2

(2.5)

onde a equação da Continuidade, Eq. (2.4), foi utilizada para simplificar a Eq. (2.3) e chegar a

Eq. (2.5).

As Eqs. (2.4) e (2.5) podem ser escritas usando-se notação vetorial e o operador de Euler

∇ , que representa o divergente de um campo vetorial ou o gradiente de um campo escalar re-

sultando

fτuuu

u

=⋅∇−∇+

∇⋅+

∂

∂

=⋅∇

pt

ρ

0

(2.6)

Estas equações descrevem o movimento de um fluido Newtoniano livre de divergente

com viscosidade constante. A hipótese de divergente nulo tem muitas implicações, não apenas

nos conceitos físicos envolvidos e no modelo matemático que descreve o fluxo, mas também,

na natureza dos algoritmos numéricos desenvolvidos para este tipo de simulação.

12

2.2.1 Condições de Contorno para as Equações de Navier-Stokes

Há duas formas de interpretar as Equações de Navier-Stokes com relação aos balanços de

quantidade de movimento: a formulação convencional velocidade-pressão e a formulação ten-

são-divergente [Gresho e Sani, 1998, vol 2]. Essas diferentes formulações ocasionam a neces-

sidade de diferentes condições de contorno de Neumann, ou seja, do fluxo de quantidade de

movimento através das fronteiras do domínio. Esta discussão é importante devido ao fato de

que imposições de condições de contorno de Newman em domínios tridimensionais necessi-

tam que as informações sejam obtidas para as faces do contorno, que representam triângulos

na malha computacional de tetraedros, e como a formulação numérica neste trabalho é obtida

por arestas, não é desejável impor condições sobre as faces, visto que isto demandaria a ne-

cessidade de manter-se a estrutura de elementos das faces ao longo da simulação.

A seguir apresentamos as Equações de Navier-Stokes, escritas na sua forma incompressí-

vel, ou seja, assumindo que o campo de velocidade é livre de divergente, e o que isso ocasio-

na em temos de condições de contorno. Seja então a equação de conservação de Quantidade

de Movimento escritas na forma convencional

( )pt

υ∂

+ ⋅∇ + ∇ = ∇ ⋅ ∇∂

uu u u (2.7)

onde υ é a viscosidade cinemática.

Por simplicidade, e sem perda de generalidade, analisamos apenas a equação no plano

x y− , sendo desta forma

2u u u pu v u

t x y xυ

∂ ∂ ∂ ∂+ + + = ∇

∂ ∂ ∂ ∂ (2.8)

2v v v pu v v

t x y yυ

∂ ∂ ∂ ∂+ + + = ∇

∂ ∂ ∂ ∂ (2.9)

onde jviu ˆˆ +=u , é o vetor velocidade do fluido, associado ao sistema de coordenadas cartesi-

anas.

13

Completando o sistema de equações, tem-se a equação de Conservação da Massa, dada

por

0=∂

∂+

∂

∂

y

v

x

u (2.10)

Existem diversas formas de representar a Eq. (2.7), em relação aos termos viscoso e con-

vectivo, e estas diferentes representações conduzem a diferentes formas espacialmente discre-

tas, que geralmente não são equivalentes, e algumas vezes apresentam vantagens ou

desvantagens sobre a forma da Eq. (2.7), denominada convencional. A seguir, escreveremos a

Eq. (2.7) na forma conhecida como Tensão-Divergente e analisaremos as condições contorno

de Newman para a velocidade neste caso.

Seja então a Eq. (2.7) escrita na forma Tensão-Divergente

( )pt

ρ∂

+ ⋅∇ = ∇ ⋅ ≡ ∇ ⋅ − ∂

uu u σ τ I (2.11)

onde σ é o tensor de tensões total, e µ≡ ∇τ u , ou equivalentemente,

j iij

i j

u u

x xτ µ

∂ ∂= + ∂ ∂

(2.12)

e I é o tensor identidade, e τ é o tensor das tensões cisalhantes.

Note que a Eq. (2.11) é realmente uma equação de balanço da quantidade de movimento,

e desta forma σ é o tensor de tensões total. A razão para muitos autores preferirem esta forma

de escrever as equações de Navier-Stokes, Eq. (2.11), esta relacionada à forma fraca dessa

equação, quando utilizando o Método dos Elementos Finitos, pois esta representação leva à

condições de contorno naturais que representam realmente as forças físicas sobre porções do

contorno que necessitam condições de Newman.

Para ilustrar este fato, consideremos as condições de contorno necessárias para as equa-

ções de Navier-Stokes Incompressíveis para um caso genérico. Se o fluido for considerado

viscoso, então as condições de parede sólida devem ser prescritas como não penetração e não

deslizamento, ou seja, as componentes normal e tangencial da velocidade do fluido devem ser

14

idênticas às componentes normal e tangencial da velocidade da parede sólida. Em regiões de

entrada do escoamento, é possível e razoável prescrever as componentes da velocidade. Am-

bas as condições de contorno citadas anteriormente são denominadas condições de contorno

de Dirichlet. Em regiões de saída do escoamento, pode-se utilizar condições de contorno aber-

to, ou seja, deixar a velocidade livre. Em geral, nessa região do contorno prescreve-se pelo

menos um ponto de pressão para garantir a unicidade da solução. As condições de contorno

relatadas são bem conhecidas e geralmente não levantam dúvidas sobre a sua aplicação. Po-

rém o mesmo não ocorre em regiões do domínio onde deve ser avaliado um balanço de forças,

denominado “traction”.

Seja a componente normal do tensor de tensões total, dado no lado direito da Eq. (2.11),

i.e.,

pµ⋅ = ∇ ⋅ − =σ n u n n F (2.13)

onde n é o versor normal apontando para fora do domínio, F é a força por unidade de área

aplicada no contorno, que é válida para superfícies de qualquer forma, porém, analisamos a-

penas superfícies planas. Desta forma, a Eq. (2.13) pode ser escrita como

pn

µ∂

− =∂

un F (2.14)

Seja então o vetor das forças aplicadas escrito como

( ) ( )( ) ( )

1

2

/ /

/ /x x y

y x y

F n u x p n u y

F n v x n v y p

µ µ

µ µ

∂ ∂ − + ∂ ∂⋅ ⋅ = = = ∂ ∂ + ∂ ∂ −⋅ ⋅

e σ nF

e σ n (2.15)

onde 1e e 2e são os versores das direções x e y , respectivamente, e 0⋅ =n t , onde n e t são

as direções normal e tangencial à superfície de contorno, respectivamente. Desta forma, usan-

do as definições anteriores, e n x x y yF n F n F= + , t x x y yF t F t F= + , na Eq. (2.15), temos

n x x y y x x y y x y

u u v vF n F n F n n p n n n n p

x y x yµ µ µ µ

∂ ∂ ∂ ∂ = + = − + + + − = ∂ ∂ ∂ ∂

15

2 2x x y x y y

u u v vn p n n n n n p

x y x yµ µ µ µ

∂ ∂ ∂ ∂ − + + + − =

∂ ∂ ∂ ∂

( )2 2 2 2x x y y x y

u u v vn n n n p n n

x y x yµ µ µ

∂ ∂ ∂ ∂+ + + − +

∂ ∂ ∂ ∂

como ( ) 22 2 1x yn n+ = =n , resulta

2 2n x x y y

u u v vF n n n n

x y x yµ µ µ

∂ ∂ ∂ ∂= + + +

∂ ∂ ∂ ∂

Definindo agora t x x y yF t F t F= + , temos

x x y y x x y y x y

u u v vF F F t n p n t n n p

x y x yτ τ τ µ µ µ µ

∂ ∂ ∂ ∂ = + = − + + + − = ∂ ∂ ∂ ∂

x x x x x y y x y y y y

u u v vt n t n p t n t n t n t n p

x y x yµ µ µ µ

∂ ∂ ∂ ∂− + + + − =

∂ ∂ ∂ ∂

x x x y y x y y x x y y

u u v vt n t n t n t n t n p t n p

x y x yµ µ µ µ

∂ ∂ ∂ ∂+ + + − − =

∂ ∂ ∂ ∂

x x x y y x y y

u u v vt n t n t n t n

x y x yµ µ µ µ

∂ ∂ ∂ ∂+ + +

∂ ∂ ∂ ∂

pois ( ) ( ) 0x x y y x x y yt n p t n p t n t n p p− − = − + = − ⋅ =t n , pois t e n são os versores das direções

tangente e normal à superfície, respectivamente, logo são ortogonais. Desta forma o vetor das

forças aplicadas pode ser escrito nas direções tangencial e normal à superfície como dado na

Eq. (2.16) a seguir:

( )2 2/ / / /

/ / / /

x x y yn

t x x y y x y y

n u x n n u y v x n v y pF

F t n u x tn u y t n v x t n v y

µ

µ

∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ −⋅ ⋅ = = = ⋅ ⋅ ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂

n σ nF

t σ n (2.16)

Retornando para a Eq. (2.11), e olhando apenas para a equação de Conservação da Quan-

tidade de Movimento na direção x, podemos escrevê-la como:

16

x

Du

Dt= ∇ ⋅σ (2.17)

onde Dt

D é a derivada total, ou material da variável u , e com

x x y

u up

x yυ υ

∂ ∂ = − +

∂ ∂ σ e e (2.18)

Aplicando-se o Método dos Resíduos Ponderados na Eq. (2.17), com função de pondera-

ção φ , e integrando por partes o segundo termo da Eq. (2.17) obtém-se a forma fraca

x x x

Dud d d

Dtφ φ φ φ

Ω Ω Γ Ω

Ω = ∇ ⋅ Ω = ⋅ Γ − ∇ ⋅∫ ∫ ∫ ∫σ n σ σ (2.19)

ou

x x

Dud d

Dtφ φ φ

Ω Γ

+ ∇ ⋅ Ω = ⋅ Γ

∫ ∫σ n σ (2.20)

Expandindo a equação anterior, fazendo o uso da Eq. (2.16), temos

Du u up d

Dt x x y y

φ φφ υ υ

Ω

∂ ∂ ∂ ∂ + − + Ω =

∂ ∂ ∂ ∂ ∫ x y

u un p n d

x yφ υ υ

Γ

∂ ∂ − + Γ

∂ ∂ ∫ (2.21)

ou ainda,

Duu p d

Dt x

φφ υ φ

Ω

∂ + ∇ ⋅∇ − Ω =

∂ ∫ x y

u un p n d

x yφ υ υ

Γ

∂ ∂ − + Γ

∂ ∂ ∫ (2.22)

Note que o termo no contorno é uma condição de contorno natural para esta forma de es-

crever as equações de Conservação da Quantidade de Movimento:

17

x x y x

u un p n F

x yυ υ

∂ ∂ ⋅ = − + =

∂ ∂ n σ (2.23)

onde xF é a força aplicada, por unidade de área, na direção x, que é prescrita. Logo a integral

no contorno resulta

∫Γ

ΓdFxφ (2.24)

Então xF é a verdadeira componente na direção x, do vetor “traction”, que representa uma

força física aplicada no contorno.

Se 0=xF em Γ , então nada precisa ser feito no contorno, ou seja, a variável fica livre.

Diante do que foi visto anteriormente, é preciso cuidado na prescrição das condições de

Neuman nas Equações de Navier-Stokes. Vamos agora obter as condições que devem ser a-

plicadas para as diferentes formulações.

2.2.2 Condições de contorno de Newman para a velocidade nas Equações de Navier-

Stokes:

As condições de contorno têm que ser aplicadas de acordo com a forma que as equações

estão escritas, por exemplo:

1) Forma Tradicional das Equações de Navier-Stokes

uuuu 2∇+−∇=

∇⋅+

∂

∂µρ p

t (2.25)

0=⋅∇ u (2.26)

Aplicando-se o Método dos Elementos Finitos para obter a forma fraca e integrando por

partes o termo viscoso, surge a seguinte integral no contorno, considerando a direção x.

18

∫Γ

Γ⋅ dxτnφ (2.27)

com

x x y

u u

x yµ µ

∂ ∂= +

∂ ∂τ e e (2.28)

onde xτ é o tensor de tensões cisalhantes na direção x.

As condições de contorno de Newman para este termo precisam ser funções das derivadas

das componentes do vetor velocidade, na direção normal ao contorno, ou seja:

n

u

∂

∂ e

n

v

∂

∂

2) Forma “tensão-divergente” das Equações de Navier-Stokes

( )Iτσuuu

pt

−⋅∇=⋅∇=

∇⋅+

∂

∂ρ (2.29)

com

( )µ= ∇τ u (2.30)

Neste caso a integração por partes aplica-se também ao gradiente de pressão, logo a inte-

gral no contorno resulta

( )∫ ∫Γ Γ

Γ⋅−=Γ⋅ dpd nIτnσ (2.31)

onde σ é conhecido como “total stress tensor”, pois considera as tensões cisalhantes e nor-

mais. Esta forma é preferida por muitos autores porque representa realmente um balanço das

forças aplicadas dado pelas equações de Conservação da Quantidade de Movimento, e as

19

condições de contorno, por sua vez são as forças físicas aplicadas a estes contornos. Então, as

condições de contorno apropriadas resultam:

Para u:

x y x

u un p n F

x yµ µ

∂ ∂ − + =

∂ ∂ em N

uΓ (2.32)

uu = em D

uΓ (2.33)

Para v:

y x y

v vn p n F

y yµ µ ∂ ∂

− + = ∂ ∂

em N

vΓ (2.34)

vv = em D

vΓ (2.35)

onde N

uΓ é a parte do contorno onde são prescritas as condições de contorno de Newman para

a componente da velocidade na direção x, D

uΓ é a parte do contorno onde são prescritas as

condições de contorno de Dirichlet para a componente da velocidade na direção x, N

vΓ é a par-

te do contorno onde são prescritas as condições de contorno de Newman para a componente

da velocidade na direção y, e D

vΓ é a parte do contorno onde são prescritas as condições de

contorno de Dirichlet para a componente da velocidade na direção y. Sendo D N

u u uΓ = Γ Γ∪ e

D N

v v vΓ = Γ Γ∪ o contorno do domínio de interesse. u é o valor da velocidade na direção x

prescrito em D

uΓ e v é o valor da velocidade na direção y prescrito em D

vΓ .

Concluindo, note que quando o p∇ é integrado por partes, basta prescrever xF e yF , e

quando o p∇ não é integrado por partes, é preciso prescrever n

u

∂

∂ e

n

v

∂

∂.

Neste trabalho utilizou-se a forma tradicional, portanto o gradiente de pressão não foi in-

tegrado por partes, e dessa forma as condições no contorno de Newman para a velocidade são

funções apenas das tensões cisalhantes, da seguinte maneira:

20

dΓ

⋅ Γ∫ τ n (2.36)

com

jiij

j i

uu

x xτ µ

∂∂= + ∂ ∂

(2.37)

Considerando-se escoamentos incompressíveis, as condições de fluxo prescrito podem ser

impostas diretamente no campo de velocidades, como condições de Dirichlet, e as condições

de escoamento livre em regiões do contorno longe de obstáculos, podem ser prescritas utili-

zando o conceito de “traction free”, da seguinte maneira:

0=∂

∂

n

un (2.38)

E desta forma, o sistema de equações de Navier-Stokes fica bem posto com as condições

de contorno de Newman (prescrição de fluxo) destacadas acima, e as condições de contorno

de Dirichlet (prescrição do valor da variável), sendo D NΓ = Γ Γ∪ , e D NΓ Γ = ∅∩ .

2.2.3 Dificuldades Associadas aos Fluxos Incompressíveis

Algumas dificuldades associadas à simulação numérica de escoamentos de fluidos incom-

pressíveis são discutidas nesta seção.

A primeira dificuldade diz respeito à presença do termo convectivo nas equações de con-

servação da quantidade de movimento, visto que este termo é não simétrico e não linear, e o

método de Galerkin padrão é instável quando aplicado a tal termo. Estas instabilidades au-

mentam com o aumento do número de Reynolds, pois os fluxos com convecção dominante

determinam maior importância do termo convectivo. Estas instabilidades podem ser evitadas

usando técnicas de estabilização, tais como SUPG (Streamline Upwind Petrov-Galerkin)

[Brooks e Hughes, 1982], GLS (Galerkin Least Squares) [Hughes et al., 1989, Franca et al.,

1988], SGS (Subgrid Scale) [Codina et al., 2006]. Estes métodos são essenciais na obtenção

21

de soluções numéricas de escoamentos incompressíveis na presença de fenômenos convecti-

vos dominantes.

Outra dificuldade numérica é imposta pela própria condição de incompressibilidade, que

para fluxos incompressíveis resulta na imposição de ocorrência de um campo de velocidade

solenoidal, ou seja, livre de divergente. Neste caso a pressão deve ser considerada como uma

variável sem nenhuma relação com equações constitutivas. Sendo assim, a presença de termos

de pressão na equação de conservação da quantidade de movimento introduz um grau de li-

berdade adicional de forma a satisfazer a restrição de incompressibilidade. Esta formulação

matemática pode ser entendida se considerarmos a pressão como um multiplicador de La-

grange, caracterizando o acoplamento entre as variáveis de pressão e velocidade.

Neste trabalho estamos utilizando uma formulação baseada nas variáveis primitivas, velo-

cidade e pressão, das equações de Navier-Stokes, caracterizando um método dos Elementos

Finitos Mistos [Zienkiewicz e Taylor, 1988]. Estes métodos apresentam dificuldades numéri-

cas devido ao fato da ocorrência de ponto de sela que resulta do problema variacional associ-

ado. Esta ocorrência faz com que o sistema algébrico para os valores nodais das variáveis de

velocidade e pressão, quando utilizando-se uma formulação de Galerkin, seja governado por

uma matriz particionada que contém uma submatriz nula sobre uma diagonal. Desta forma a

solução do sistema algébrico depende da escolha dos espaços de elementos finitos para velo-

cidade e pressão. Estes espaços devem satisfazer uma condição de compatibilidade, conhecida

como “LBB Condition” (Ladyzhenskaya-Babuska-Brezzi) [De Sampaio, 1991, Gresho e Sani,

1998, Hughes et al., 1996, Zienkiewicz e Taylor, 1988], que identifica as condições que de-

vem ser satisfeitas pelos espaços de elementos finitos de forma a tornar a formulação numéri-

ca estável.

Condição de Ladyzhenskaya, Babuska e Brezzi (LBB Condition)

Dado que o sistema algébrico resultante possui uma submatriz nula sobre a diagonal, é

preciso definir sob quais condições este sistema pode ser resolvido. Obviamente estas condi-

ções precisam estar associadas às variáveis de velocidade e pressão, pois, dado o sistema

=

h

f

p

u

0G

GKT

(2.39)

22

vemos que se o núcleo da matriz G é nulo, 0N(G) = , então a matriz global é não singular e

desta forma as variáveis de velocidade e pressão são unicamente definidas. Entretanto se o

núcleo, ou espaço nulo da matriz G não for nulo, então as soluções obtidas apresentam con-

vergência para as variáveis de velocidade, mas resultam em oscilações no campo de pressão.

Neste contexto, Ladyzhenskaya em 1969, Babuska em 1970/1071 e Brezzi em 1974 determi-

naram uma condição de compatibilidade, conhecida como “LBB condition”, que os espaços

discretos devem satisfazer para garantir a estabilidade dos métodos dos elementos mistos. A

não satisfação desta condição implica no aparecimento de oscilações espaciais no campo de

pressão, usualmente denominadas modos de oscilações espúrias. Estas oscilações não físicas

surgem devido aos resíduos oriundos do fato de que o campo de velocidade obtido solução

das equações de conservação da quantidade de movimento não garantem um campo solenoi-

dal.

A “LBB condition” diz que os espaços discretos utilizados para as variáveis de velocidade

e pressão não podem ser escolhidos de forma arbitrária, é preciso que estejam relacionados.

Para ilustrar este fato, consideremos o sistema matricial resultante das equações de Navier-

Stokes dado pelas Eq. (2.39), onde eqeq nn ×K é uma matriz quadrada de ordem eqn ,

eqeq nn ˆ×G é

uma matriz retangular, e u é o vetor das variáveis de velocidade, p é o vetor das variáveis de

pressão, f e h são os vetores das contribuições dos carregamentos e valores prescritos no

contorno, com as respectivas dimensões relativas às matrizes definidas, sendo eqn o número

de variáveis de velocidade, e eqn o número de variáveis de pressão. Note que para garantir

que o sistema seja compatível e determinado, ou seja, possua solução única, as matrizes do

sistema precisam possuir posto completo, ou seja, a dimensão do espaço das linhas de cada

matriz deve ser igual à dimensão da matriz quadrada associada a cada matriz. Neste caso, co-

mo K é uma matriz regular, então eqnrank =)(K , implicando que os vetores linha de TG

devem ser linearmente independentes para garantir o posto completo do sistema algébrico.

Então eq

T nrank ˆ)( =G . Como a matriz TG possui eqn linhas, uma condição necessária é que

eqeq nn ≤ˆ . Isso significa que uma condição necessária, mas não suficiente para a obtenção da

solução do sistema algébrico é que o espaço discreto que define as variáveis de pressão deve

possuir dimensão menor ou igual ao espaço discreto que define as variáveis de velocidade.

A condição suficiente que relaciona os espaços discretos definidos é dada pela “LBB

condition” de compatibilidade, que diz que: “a existência de uma solução aproximada de ele-

23

mentos finitos estável ( )hh pu , depende da escolha dos espaços discretos hV e hQ , tal que a

seguinte condição seja satisfeita:

( ),

inf sup 0h h hh h h

h h

h

q Q Q V

q

qα

∈ ∈

∇ ⋅≥ >

w V

w

w (2.40)

onde inf, é o maior limitante inferior, sup é o menor limitante superior, α é um parâmetro in-

dependente do tamanho da malha, i é a norma 2L , ( ),a b a bdΩ

= ⋅ Ω∫ é uma forma bilinear

contínua de Q× →V . Se a “LBB condition” de compatibilidade é satisfeita, então existe

um único hh Vu ∈ e um hh Qp ∈ que resolve o sistema algébrico. Desta forma, a “LBB con-

dition” é uma condição necessária e suficiente para a resolução satisfatória do sistema.

Perceba que verificar a “LBB condition” não é tarefa simples, porém notamos que a satis-

fação desta condição diz respeito aos espaços de interpolação das variáveis de pressão e velo-

cidade que são definidos de forma algébrica. Sendo assim, podemos encontrar as restrições na

definição destes espaços de forma que a condição seja satisfeita [De Sampaio, 1991]. Seja o

sistema da Eq. (2.39), onde u e p são as variáveis livres associados aos campos de velocidade

e pressão, respectivamente. Assumindo que existem nu variáveis de velocidade e np variá-

veis de pressão. Desta forma as matrizes K , G e TG possuem dimensão n n×u u , n n×u p e

n n×p u , respectivamente. Note que a matriz K corresponde a discretização dos termos de

massa, convecção e difusão, e que se considerarmos o fluxo de Stokes, então esta matriz é si-

métrica e positiva definida, portanto possuindo todos os autovalores e pivots positivos e assim

garantindo a não singularidade. Então é possível escrever a partir do sistema (2.39)

+ =Ku Gp f (2.41)

então, isolando o termo Gp na Eq. (2.41), resulta

= −Gp f Ku (2.42)

Multiplicando agora a Eq. (2.42) por 1−K , visto que K é inversível resulta

24

1 1 1− − −= −K Gp K f K Ku (2.43)

e como 1− =K K I , e ainda multiplicando a Eq. (2.43) por TG temos

1 1T T T− −= −G K Gp G K f G u (2.44)

Agora, substituindo T =G u h na Eq. (2.44) resulta

1 1T T− −= −G K Gp G K f h (2.45)

e definindo

1T −=H G K G (2.46)

tem-se

1T −= −Hp G K f h (2.47)

Note que a matriz H não pode ter posto maior do que nu , que é o posto da matriz K , e

portanto, da matriz 1−K . Por outro lado, para se obter uma solução única para a pressão, o

posto da matriz H deve ser np , logo a seguinte desigualdade é necessária para garantir a es-

tabilidade do sistema.

n n≤p u (2.48)

Porém, a condição dada em (2.48) não é suficiente para garantir a não singularidade da

matriz H . Para tal, o espaço nulo da matriz G deve ser zero, 0N(G) = , ou seja

=Gp 0 se e somente se =p 0 . (2.49)

Assim, a “LBB condition” fica satisfeita se as Eqs. (2.48) e (2.49) forem satisfeitas.

25

2.2.5 Formulação de Elementos Finitos Estabilizada

As formulações numéricas baseadas no método de Galerkin resultam em aproximações

instáveis para as variáveis de velocidade e pressão quando os escoamentos apresentam con-

vecção dominante, pois aproximações de Galerkin são aproximações do tipo centradas. Entre-

tanto é possível obter formulações de elementos finitos, denominadas estabilizadas, que

permitem o uso de espaços discretos de mesma dimensão para as variáveis de velocidade e

pressão, satisfazendo a “LBB condition”, e ainda possibilitando a solução de problemas de es-

coamentos com elevado número de Reynolds.

Segundo Tezduyar (1992), formulações estabilizadas de elementos finitos para as equa-

ções de Navier-Stokes são obtidas com a adição à formulação padrão de Galerkin de termos

de estabilização, baseados no fato de que em escoamentos com fenômenos convectivos domi-

nantes, as aproximações não podem exceder primeira ordem, pelo menos não em regiões do

escoamento que apresentam gradientes de pressão acentuados, desta forma evitando oscila-

ções espúrias principalmente no campo de pressão. Neste caso, então a redução da ordem da

aproximação, de segunda ordem, referente ao método de Galerkin, para primeira ordem, pode

ser obtida adicionando-se termos de difusão adicional, numérica, controladas por algum pa-

râmetro que identifique as regiões onde esta redução da ordem de aproximação faz-se neces-

sária, e apenas nestas regiões reduza a ordem da aproximação, e mantendo segunda ordem no

restante do domínio. Os termos associados com os parâmetros de estabilização permitem o

uso de elementos finitos mistos com mesma ordem de interpolação para as variáveis de velo-

cidade e pressão, sendo ainda estes termos, resíduos ponderados que garantem a consistência

da formulação.

2.3 Método do Passo Fracionado

Um método bastante popular utilizado para discretizar as equações de Navier-Stokes no

tempo é o Método do Passo Fracional (Fractional Step Method) (Donea e Huerta, 2003) , que

consiste em obter o avanço no tempo através da decomposição do passo no tempo em uma

seqüência de dois ou mais passos. Estes métodos foram desenvolvidos para as equações de

Navier-Stokes independentemente por Chorin, 1968 e Teman, 1969. Esta decomposição per-

mite superar as dificuldades numéricas associadas ao problema de ponto de sela que origina-

se da formulação variacional das Equações de Navier-Stokes. A idéia básica é dividir o trata-

26

mento numérico dos vários operadores nas equações, decompondo o sistema original, de

complicada obtenção da solução, em problemas mais simples, facilitando a obtenção da solu-

ção. Há muitas maneiras de realizar este fracionamento das equações, e desta forma, há tam-

bém diversos Métodos de Passo Fracional para as equações de Navier-Stokes. Em geral, os

Métodos de Passo Fracional realizam a discretização no tempo antes da discretização espacial,

e quando isso é feito, surge uma controvérsia entre os autores a respeito das condições de con-

torno em cada passo de tempo, visto que o fracionamento implica em variáveis adicionadas ao

sistema, e estas novas variáveis demandam algum tipo de condição para que o problema con-

tinue bem posto matematicamente. Outra característica relacionada aos Métodos de Passo

Fracional é a precisão com relação à discretização temporal, pode ser de primeira ou segunda

ordem, dependendo do fracionamento.

O fracionamento das equações de Navier-Stokes implica em definir o tratamento dos ter-

mos em relação à discretização temporal, ou seja, tratar todos os termos de forma explícita,

implícita, ou semi-implícita. A escolha desta discretização implica em conseqüências na ob-

tenção da solução. Muitos autores preferem o tratamento semi-implícito, devido ao fato de

que o termo convectivo pode facilmente ser linearizado, e também a formulação resultar in-

condicionalmente estável, por outro lado uma formulação totalmente implícita também é in-

condicionalmente estável, porém necessita grande esforço computacional, e um tratamento

explícito fica restrito em relação ao tamanho do passo de tempo, que neste caso deve obedecer

os limites de estabilidade.

Em geral, para um esquema de passo fracional, tem-se resultados acurados para pressão

integrando os termos convectivos de forma explícita, e consequentemente passo de tempo im-

posto pelo limite de estabilidade.

Há duas formas de se obter o fracionamento das Equações de Navier-Stokes, dando ori-

gem a duas visões do mesmo fato:

Chorin (1967 e 1968) e Teman (1969) utilizam os conceitos de Métodos de Projeção, onde

é possível computar os campos de velocidade e pressão separadamente através da obtenção de

uma velocidade intermediária, que é então projetada no subespaço das funções vetoriais que

são solenoidais. A idéia básica desta visão dos Métodos de Passo Fracional é o fato que todo

vetor pode ser decomposto em sua parte livre de divergente e a sua parte que é gradiente de

uma função escalar, conhecida como Princípio da Decomposição de Helmholtz, que baseia-se

no teorema da Decomposição Ortogonal de Ladyzhenskaya. Este teorema diz que qualquer

campo vetorial w admite uma decomposição ortogonal única da forma

27

φ∇+= uw (2.50)

isto é, em um campo solenoidal, u , livre de divergente e o gradiente de alguma função escalar

φ .

Os Métodos de Projeção de Chorin e Teman podem ser vistos como uma decomposição

do sistema de equações algébricas originado das Equações de Navier-Stokes (Chorin, 1967,

Donea e Huerta, 2003, Teman, 1969), dando origem a outra visão dos Métodos de Passo Fra-

cional. Nesta visão, uma decomposição LU é realizada no sistema original, resultando em um

fracionamento puramente algébrico das equações. Muitos autores preferem este tipo de de-

composição porque como nenhum conceito físico é utilizado, as variáveis adicionais que sur-

gem do fatoramento das equações não têm nenhum significado físico e desta forma não

demandam condições de contorno.

Neste trabalho foi utilizado um método de projeção, ou “Fractional Step Method” (Codi-

na, 2002), que baseia-se no fracionamento algébrico do sistema de equações, dividindo o sis-

tema resultante em sistemas menores, o que neste caso foi obtido a partir de uma fatoração

LU do sistema matricial resultante (Codina, 2002, Henriksen et al., 2002, Perot, 1993). Abai-

xo mostramos a idéia central do método.

Sejam as equações de Navier-Stokes na sua forma Incompressível

0

2

=⋅∇

∇+−∇=∇⋅+∂

∂

u

uuuu

υpt (2.51)

Discretizando de forma implícita os termos de convecção, difusão e de pressão no tempo,

tipo “Euler Backward” a Eq. (2.51), temos

11 1 1 2 1

n nn n n n

pt

υ+

+ + + +−+ ⋅∇ = −∇ + ∇

∆

u uu u u (2.52)

ou

11 1 2 1 1

n nn n n n

pt t

υ+

+ + + ++ ⋅∇ − ∇ + ∇ =∆ ∆

u uu u u (2.53)

28

e com a condição de incompressibilidade

01 =⋅∇ +nu (2.54)

Adicionado à Eq. (2.53) o termo ( )n np p∇ − ∇ , que será utilizado para determinar o in-

cremento de pressão, e escrevendo as Eqs. (2.53) e (2.54) resultam matricialmente em

+

=

+

2

1

bc

bc

0

r

q

u

0D

GA nn 1

(2.55)

onde

→A Operador Implícito para os termos de Massa, Conveção e Difusão

→G Operador Gradiente

→D Operador Divergente

→nr Vetor das forças aplicadas modificado pela contribuição da pressão

→q Contém as variáveis de pressão, sendo q = pn+1 - pn (Chang et al., 2002)

→1bc Condições de contorno para as equações de Conservação da Quantidade de Mo-

vimento

→2bc Condições de contorno para a restrição de incompressibilidade

A estrutura da matriz A depende do esquema de avanço no tempo, da malha e do elemen-

to finito adotado. A formulação mais comum é manter a convecção explícita e a difusão im-

plícita. Desta forma, ao invés das Eqs. (2.53) e (2.54) temos

12 1 1

n nn n n n

pt t

υ+

+ +− ∇ + ∇ = − ⋅∇∆ ∆

u uu u u (2.56)

01 =⋅∇ +nu (2.57)

Se a convecção for explícita na Eq. (2.55) υLI

A −=→∆t

(2.58)

29

Se a convecção for implícita na Eq. (2.55) NLI

A +−=→ υ∆t

(2.59)

Onde

→L Operador Laplaciano

→N Operador Convectivo Linearizado

Com convecção implícita a matriz deixa de ser uma matriz simétrica, bloco diagonal, de-

vido ao operador N . Porém mesmo com convecção implícita, o termo N pode ser especifi-

cado de forma aproximada, tal que a matriz A seja ainda bloco diagonal, mas obviamente

isso reduz a precisão temporal para 1ª ordem.

Note que a dificuldade em resolver o sistema acima vem do fato da existência de uma

submatriz zero na matriz global, então, como para que o sistema tenha solução é necessário

que a matriz tenha posto completo, cuidados devem ser tomados com respeito às ordens de in-

terpolação das variáveis de velocidade e pressão (Zienkiewicz and Taylor, 1988). Entretanto

isto pode ser contornado dividindo o sistema inicial em sistemas menores, que é a idéia cen-

tral dos métodos de projeção. Um exemplo destes métodos é apresentado a seguir.

Efetuando uma fatoração LU do sistema dado na Eq. (2.55) resulta

− →

−

− −

GDAD

0A

0D

GA GA 1

1

CC 12 (2.60)

ou seja,

−=

−

− I0

GAI

GDAD

0A

0D

GA 1

1 (2.61)

onde

→2C Segunda coluna da matriz A

→1C Primeira coluna da matriz A

→I Matriz Identidade

30

O sistema de equações (2.55), usando a decomposição LU dada em (2.61) fica

+

=

−

+−

−2

1

bc

bc

0

r

q

u

I0

GAI

GDAD

0A n1n1

1 (2.62)

Fazendo

=

+=

−++−

Iq

w

Iq

GqAIu

q

u

I0

GAI 11n1n1

(2.63)

O sistema (2.62) pode ser reescrito como

+

=

− −2

1

bc

bc

0

r

Iq

w

GDAD

0A n

1 (2.64)

Que resulta

1bcrAw += n (2.65)

e

2bcGIqDADw =− −1 (2.66)

ou

2bcDwGqDA −=−1 (2.67)

Usando a primeira linha do sistema (2.63) podemos reescrever a Eq. (2.65) em termos de

1+nu , i.e.

AwbcrGqAu 1 =+=++ n1n (2.68)

31

Multiplicando a Eq. (2.68) por 1−A leva a

GqAAwAu 111n −−+ −= (2.69)

resultando

n 1 1+ −= −u w A Gq (2.70)

Fazendo 1−A , em (2.62), na matriz triangular inferior igual a H1 e na matriz triangular

superior igual a H2 resulta

( )

= = −−

22

2 11

A AH GA 0A G I H G

D D H H GD DH GD 0 0 I (2.71)

e diferentes aproximações da matriz 1−A , ou seja, para H1 e H2, levam a diferentes tipos de

métodos de projeção (Chang et al., 2002), como pode ser observado abaixo:

1−== AHH 21 - Sistema Inicial (2.72)

1−≠= AHH 21 - Apenas massa é conservada (2.73)

1−=≠ AH;HH 221 - Apenas quantidade de movimento é conservada (2.74)

1−≠≠ AHH 21 - Nenhuma grandeza física é conservada (2.75)

Neste trabalho optou-se por conservar massa ao preço de perturbar as equações de quanti-

dade de movimento, desta forma a mesma aproximação da inversa da matriz A é empregada

em todas as equações do sistema.

Fazendo 1t

− ≠ = = ∆1 2A H H I tem-se o chamado “Classic Fractional Step Method”

(Chang et al., 2002), e nenhuma matriz inversa precisa ser calculada.

Reescrevendo o sistema dado por (2.65), (2.67) e (2.70), após as aproximações, fazendo

1t

− ≠ = ∆ =1 2Α H I H , resulta

32

GqHwu

bcDwGqDH

bcHrHwbcrAw

2

21

1111

−=

−=

−=⇒+=

+1n

nn

(2.76)

Nesta formulação, estamos fazendo

wu =~ (2.77)

nn fr = (2.78)

1bcruA += n~ (2.79)

21 bcuDGqDH −= ~ (2.80)

1n 1+ = −u u H Gq (2.81)

Suprimindo por enquanto as condições de contorno, isto é, 1 2= =bc bc 0 , o sistema de

equações dado por (2.36), com a nomenclatura definida de (2.76) a (2.81), pode ser expresso,

de forma analítica, como

1n θ n θ n n θ n θp υ

t

+ + + + +∂+ ⋅∇ + ∇ − ∇ =

∂2u

u u u f (2.82)

( ) 112 ~ ++ ⋅∇=−∇∆ nnn ppt u (2.83)

( )1 1

1 0n n

n np p

t

+ ++−

+ ∇ − =∆

u u (2.84)

para uma discretização temporal em dois passos, dada pelo método Trapezoidal, onde θ con-

trola o tipo da aproximação, onde se 1θ = , tem-se uma discretização do tipo “Euler Back-

ward” padrão que possui erro na discretização temporal de ordem ( )O t∆ , e se 1/ 2θ = , tem-

se uma discretização do tipo “Crank-Nicholson”, que possui erro na discretização temporal de

ordem 2( )O t∆ , sendo portanto, mais acurado no tempo.

Comparando as Eqs. (2.82-2.84) com as Eqs. (2.65,2.67,2.70), tem-se

( )

( )

1 2

2 1 1

1 11 10

n n n n n n

n n n 1

n nn n n 1

pt

t p p

p pt

θ θ θ θυ+ + + + +

+ + −

+ ++ + −

∂+ ⋅∇ + ∇ − ∇ = ⇒ =

∂∆ ∇ − = ∇ ⋅ ⇒ =

− + ∇ − = ⇒ = − ∆

uu u u f Aw r

u DA Gq Dw

u uu w A Gq

(2.85)

33

Note que estamos fazendo

IA t∆=−1 (2.86)

tυ= − +

∆

IA L N (2.87)

Implicando que

qADGqAGqDA 2∇== −−− 111 (2.88)

Lembrando que

→N Operador Convectivo

→L Operador Laplaciano

→D Operador Divergente

→G Operador Gradiente

Então

nn pp −= +1q (2.89)

uw ~= (2.90)

1~ +⋅∇= nuDw (2.91)

( )1 1

1 1 1 1 1 0n n

n n n n nt p p

t

+ ++ − + + +−

− = − ⇒ − = −∆ ∇ ⇒ + ∇ − =∆

u uu w A Gq u u I q

(2.92)

Obtendo assim uma formulação seguindo o Método do Passo Fracionado. Sendo que os

super–índices n e θ , referem-se ao passo de tempo e à discretização temporal dada pela regra

trapezoidal, respectivamente. u~ é a velocidade intermediária, que é introduzida na formula-

ção para permitir o fracionamento das equações de conservação da quantidade de movimento.

34

Erro da Aproximação

É possível estimar o erro nas aproximações dos termos convectivo e difusivo, devido ao

fato destes serem tratados implicitamente, conforme apresentamos a seguir. Sendo

( )nnnn ppt −∇∆+= +++ 111~ uu (2.93)

e

( ) nnn uuu ~1~~ 1 θθθ −+= ++ (2.94)

Implicando que

( )[ ] ( ) ( )[ ]111 1~ −+++ −∇∆+−+−∇∆+= nnnnnnn pptppt uuu θθθ (2.95)

Utilizando-se as Eqs. (2.93) e (2.94) na primeira expressão da Eq. (2.85) temos

( )[ ] ( )[ ] ( ) ( )[ ]

( )[ ] ( ) ( )[ ] ( )[ ] ( ) ( )[ ] θθθυ

θθ

θθ

+−++

−++

−++++

=−∇∆+−+−∇∆+∆

−∇+−∇∆+−+−∇∆+∇

⋅−∇∆+−+−∇∆++∆

−−∇∆+∴

nnnnnnn

nnnnnnn

nnnnnnnnnn

pptppt

ppptppt

pptpptt

ppt

fuu

uu

uuuu

111

111

11111

1

1

1

(2.96)

Arrumando os termos acima,

( ) ( ) ( ) ( ) ( ) 1

1 1 1 11 1n n

n n n n n n n np p t p p t p p

tθ θ θ θ

++ + + −−

+ ∇ − + + − + ∆ ∇ − + − ∆ ∇ − ⋅ ∆

u uu u

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )

1 1 1

1 1 1

1 1

1 1

n n n n n n

n n n n n n n

t p p t p p

t p p t p pθ

θ θ θ θ

υ θ θ θ θ

+ + −

+ + − +

∇ + − + ∆ ∇ − + − ∆ ∇ − −

∆ + ∆ ∇ − + − + − ∆ ∇ − =

u u

u u f (2.97)

Como ( )baba +∇=∇+∇ , implica

( ) ( ) ( )

( ) ( )( )( )

1 1

1 1

1

1

n n n

n

n n n n n

t p p t p p

t p p p p t pθ

θ θ

θ θ δ

+ −

+ − +

∆ ∇ − + − ∆ ∇ − =

∆ ∇ − + − − = ∆ ∇ (2.98)

35

Então, a expressão final resulta

( ) ( )1

1n n

n n n n np t p t p

t

θ θ θ θδ δ+

+ + + + +−+ ∇ + + ∆ ∇ ⋅∇ + ∆ ∇

∆

u uu u

( )n n nt pθ θ θυ δ+ + +− ∆ + ∆ ∇ =u f (2.99)

Ou ainda,

( ) ( ) ( )

11

n nn n n n

n n n n n n n

pt

t p t p t p t p

θ θ θ

θ θ θ θ θ θ θ

υ

δ δ δ υ δ

++ + + +

+ + + + + + +

−+ ∇ + ⋅∇ − ∆ =

∆

− ∆ ∇ ⋅∇ + ∆ ∇ − ⋅∇ ∆ ∇ + ∆ ∆ ∇

u uu u u

f u u (2.100)

E, finalmente

( ) ( )( ) ( )θθθθθθθ

θθθ

δδυδδ

υ

+++++++

+++++

∇∆∇⋅−∇∇⋅∇∆+∇∆+∇⋅∇∆−

=∆−∇⋅+∇+∆

−

nnnnnnn

nnnnnn

ptptptpt

pt

uuf

uuuuu 1

1

(2.101)

Os primeiros cinco termos da equação anterior representam as equações de conservação

da quantidade de movimento discretizadas no tempo usando o método θ . Os termos restantes

representam o erro do esquema implícito. Note que o erro está associado ao fato de que os

termos convectivo e difusivo foram fracionado, e este erro é função do gradiente de pressão,

logo, utilizando oincremento de pressão é possível anular o erro, mantendo a ordem da apro-

ximação temporal inicial.

Formulação Variacional

Uma formulação variacional discreta pode ser obtida para as equações deduzidas anteri-

ormente e dadas em (2.85): dados ( )n

h

n

h p,u , encontrar ( )11, ++ n

h

n

h pu em hh Q×V tal que

( ) ( ) ( ), , ,n n n

h h h h h h ht

θ θ θυ+ + +∂+ ⋅∇ + ∇ ∇

∂u v u u v u v

36

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

h h h h hpθ θυ + + +

Γ− ∇ + ∇ =u v v f v (2.102)

( )( ) ( )11 ~,, ++ ⋅∇−=∇−∇∆ n

hh

nnqqppt u (2.103)

( ) ( )1 11, ( ),n n n n

h h h h h hp pt

+ +− = ∇ −∆

u u v v (2.104)

onde a Eq. (2.102) foi integrada por partes no termo difusivo, e note que termo de contorno

trata-se das “pseudo-tractions” definidas nas seções 2.2.1 e 2.2.2 deste trabalho, e em geral

são nulas.

As equações são válidas ( ) hhhh Qq ×∈∀ Vv , , onde hV e hQ são os espaços de elementos

finitos padrão, onde as componentes da velocidade e a pressão estão sendo aproximados, res-

pectivamente, e hv e hq são as respectivas funções de ponderação. Tem-se que hh uu ,~ e hp

são as aproximações de elementos finitos para a velocidade intermediária, a velocidade e a

pressão, respectivamente, e o subescrito h refere-se ao problema discreto. Finalmente

( ) ∫Ω

Ω⋅= dbaba, (2.105)

e

( ) ∫Γ

Γ Γ⋅= dbaba, (2.106)

para a e b vetores quaisquer.

Voltando ao sistema de Eqs. (2.102-2.104), introduzindo uma estrutura matricial, para fa-

cilidade de manipulação, e suprimindo os termos de contorno

( )1 11)n n n n n n

t

θ θ θ+ + + + +− + + =∆

M U U K(U U GP F (2.107)

( ) 11 ~ ++ −=−∆ nnnt UDPPL (2.108)

( ) ( ) 0PPGUUM =−+−∆

+++ nnnn

t

111 ~1 (2.109)

37

onde

→M Matriz de Massa

→K Matriz das contribuições Difusivas e Convectivas (termo de convecção linearizado)

→G Matriz Gradiente

→F Vetor das forças aplicadas

→L Matriz Laplaciana

→D Matriz Divergente, sendo TGD =

→U Vetor contendo as variáveis nodais hu

→U~

Vetor contendo as variáveis nodais hu~

→P Vetor contendo as variáveis nodais hp

Isolando 1~ +nU na Eq. (2.109), resulta

( )nnnn t PPGMUU −∆+= +−++ 1111~ (2.110)

Substituindo a Eq. (2.110), na Eq. (2.107), temos

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

tt

θ θ θ+ − + + + + ++ ∆ − − + + =∆

M U M G P P U K(U U GP F (2.111)

ou

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

t

θ θ θ+ + + + + +− + + = − −∆

MU U K(U U GP F G P P (2.112)

Substituindo, agora a Eq. (2.110) na Eq. (2.108),

( ) ( )( )1 1 1 1n n n n nt t

+ + − +∆ − = − + ∆ −L P P D U M G P P (2.113)

ou ainda

38

( )( ) 0PPGDMLDU =−−∆+ +−+ nnn t 111 (2.114)

Obs.: Como ( )nnnn t PPGMUU −∆+= +−++ 1111~, então ( )nnnn t PPGMUU −∆+= +−++ θθθ 1~

, e

como ( ) nnn PPP θθθ −+= ++ 11 , implica que:

( )[ ]nnnnn t PPPGMUU −−+∆+= +−++ θθθθ 1~ 11 (2.115)

( ) θθθθ θ +++−++ +=−∆+= nnnnnn t SUPPGMUU 11~ (2.116)

Sendo assim, como K é linear, pois o termo difusivo é linear e o termo convectivo foi li-

nearizado, então

θθθθθθθθ

θθθθθθ

++++++++

++++++

+++=

++=nnnnnnnn

nnnnnn

SK(SUK(SSK(UUK(U

SUSK(UUUK(

))))

))(~

)~

(2.117)

onde

( )nnn t PPGMS −∆= +−+ 11 θθ (2.118)

Logo,

))~

)~ θθθθθ +++++ += nnnnn E(UUK(UUUK( (2.119)

Onde

θθθθθθθ +++++++ ++= nnnnnnn SK(SUK(SSK(UE(U )))) (2.120)

Sendo assim, o sistema formado pelas Eq. (2.112) e Eq. (2.114), monolítico resulta

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

t

θ θ θ θ+ + + + + +− + + = +∆

M U U K(U U GP F E(U (2.121)

39

( )( ) 0PPGDMLDU =−−∆+ +−+ nnn t 111 (2.122)

Neste sistema de equações, a estabilização da pressão é obtida do termo

( ) 11 +−−∆ nt PGDML que é uma matriz positiva semi-definida. O termo ( )θ+nUE representa o

erro resultante do tratamento implícito dos termos de convecção e difusão, que é da ordem do

quadrado do passo de tempo, ( )2tO ∆ . Note que a estabilização da pressão é dependente do

passo de tempo, logo esta formulação não apresenta-se estável nos termos de pressão, para

passos de tempo muito pequenos, contudo, a formulação pode ser reinterpretada como um sis-

tema monolítico estabilizado, definindo-se um parâmetro de estabilização τ , que tem como

função ponderar a quantidade de estabilização necessária em cada região do domínio.

Esquema Monolítico Estabilizado

Uma formulação monolítica estabilizada do sistema dado pelas Eq. (2.112), (2.113) e

(2.105) pode ser escrito da seguinte maneira: dado n

hu , encontrar ( )111 ,, +++ n

h

n

h

n

h p ξu em

hhh Q VV~

×× , tal que

( ) ( ) ( ) ( ) ( )h

n

h

n

hh

n

hh

n

h

n

hhn pt

vfvvuvuuvu ,,,,,1 1 θθθθ υδ +++++ =∇−∇∇+∇⋅+∆

(2.123)

( ) ( ) 0,, 11 =∇−∇+⋅∇ ++h

n

h

n

hh qpq ξu τ (2.124)

( ) ( )h

n

hh

n

h p vvξ ~,~, 11 ++ ∇= (2.125)

Para ( ) hhhhhh Qq VVvv~~,, ××∈∀ e n

h

n

hh uuu −= +1δ . O espaço funcional hV~

é o mesmo hV ,

porém sem as condições de contorno, e onde τ é o passo de tempo crítico, necessário para e-

vitar oscilações espúrias da pressão. Obviamente este passo de tempo crítico deve ser o mes-

mo utilizado para integrar no tempo os termos difusivo e convectivo explicitamente.

Como,

40

( ) ( )n

h

n

hhhtt

UUMvu −∆

=∆

+11,

1δ (2.126)

É possível escrever as Eqs. (2.125-2.127) em forma matricial, resultando:

( )1 11)n n n n n n

t

θ θ θ+ + + + +− + + =∆

M U U K(U U GP F (2.127)

( )1 1 1n n nτ+ + ++ − =DU LP DΕ 0 (2.128)

1 1n n+ +=MΕ GP (2.129)

onde D~

, M~

e G~

são as mesmas matrizes D , M e G , sem as condições de contorno.

Escrevendo 1n+Ε em função de 1+nP , e usando a Eq. (2.129),

1 1 1n n+ − +=Ε M GP (2.130)

E inserindo na Eq. (2.128), resulta no sistema

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

t

θ θ θ+ + + + +− + + =∆

M U U K U U GP F (2.131)

( ) 0PGMDLDU =−+ +−+ 111 ~~~ nn τ (2.132)

Que é muito similar ao sistema fracionado obtido anteriormente Eqs. (2.121) e (2.122). A

principal diferença é que agora, a estabilidade é independente de t∆ , e que o erro devido ao

tratamento implícito dos termos viscosos e convectivos ( )θ+nUΕ é automaticamente elimina-

do.

O sistema de Eqs. (2.131-2.132) pode ser resolvido de forma bastante satisfatória usando

um algoritmo monolítico para obter as variáveis de velocidade e pressão, porém neste traba-

lho, queremos resolver este sistema de forma segregada, aproveitando-se do melhor resolve-

dor linear para cada tipo de sistema linear resultante de cada equação. Infelizmente, em

41

muitos casos, o sistema de Eqs. (2.131-2.132) apresenta-se muito mal condicionado, impedin-

do a convergência do resolvedor. Desta forma o sistema precisa ser pré-condicionado, o que

pode ser resolvido, fracionando-se novamente o sistema de Eqs. (2.131-2.132), obtendo uma

nova formulação estabilizada.

Seja o sistema, dado nas Eqs. (2.127-129)

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

t

θ θ θ+ + + + +− + + =∆

M U U K U U GP F (2.133)

( )1 1 1n n nτ+ + ++ − =DU LP DΕ 0 (2.134)

1 1n n+ +− =MΕ GP 0 (2.135)

Este deve ser fracionado novamente. Então, escrevendo as matrizes

1

1

1

n n n

n

n

τ τ

θ+ +

+

+

+ − = −

A G 0 U r F

D L D P 0

0 G M Ε 0

(2.136)

onde A é a matriz contendo as contribuições da massa, convecção e difusão, D é a matriz

divergente, G é a matriz a gradiente, L é a matriz Laplacian, M , D e G são as matrizes de

massa, divergente e gradiente sem considerar condições de contorno, respectivamente, nr é o

vetor das condições de contorno, n θ+F é o vetor dos carregamentos externos, e 1n+U , 1n+P e

1n+E são as variáveis de velocidade, pressão e projeção do gradiente de pressão no espaço de

elementos finitos, respectivamente.

Fazendo a fatoração da matriz do sistema de Eqs. (2.136), resulta

( )

12 1

13 2

C C 1

C C 1

1

τ τ τ τ

τ Lower

τ

τ

−

−

− −

− − −

−

− → − − − −

→ − = − −

A G

Β D

A G 0 A 0 0

D L D D L DA G D

0 G M 0 G M

A 0 0

D L DA G 0

0 G M GΒ D

(2.137)

42

e,

Upperτ1

1

=

−

−−

−

000

DΒI0

0GAI~

(2.138)

onde ( )1τ −= −Β L DA G

Agora escrevendo novamente o sistema, desta vez com a matriz fatorada

1

1

1

1 n n n

1 1 n

1 n

τ τ

τ

θ− + +

− − +

− +

+ − − = − −

A 0 0 I A G 0 U r F

D L DA G 0 0 I Β D P 0

0 G M GΒ D 0 0 0 Ε 0

(2.139)

Como

1

1

1

1 n

n

n n 1

− +

− +

+ +

− =

11

2

I A G 0 U w

0 I Β τD P w

0 0 0 Ε Ε

(2.140)

resulta

θ++= nn FrAw1 (2.141)

( )1τ −+ − =1 2Dw L DA G w 0 (2.142)

( )1 n 1τ − +− + − =2Gw M GΒ D Ε 0 (2.143)

usando a definição de 2w dada na Eq. (2.140) na Eq. (2.143), temos

( ) ( )n 1 1 n 1 1 n 1τ τ+ − + − +− − + − =G IP Β DΕ M GΒ D Ε 0 (2.144)

ou

43

1 1 1 1 1 1n n n nτ τ+ − + + − +− + + − =GP GΒ DΕ MΕ GΒ DΕ 0 (2.145)

e como GG =~

sem as condições de contorno, a Eq. (2.145) fica

1 1n n+ +=GP MΕ (2.146)

Após alguma manipulação, usando as Eqs. (2.133-2.135), o novo sistema fracionado re-

sulta:

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

t

θ θ θ+ + + + +− + = − +∆

M U U K U U GP F (2.147)

( ) ( )1 1 1 1n n n n nt τ+ + + +∆ − + − = −L P P LP DΕ DU (2.148)

( ) ( )1 1 11 n n n n

t

+ + +− − − =∆

M U U G P P 0 (2.149)

1 1n n+ +=MΕ GP (2.150)

Para escrever o sistema monolítico, 1~ +nU é obtido da Eq. (2.150) e inserido na Eq. (2.148)

( )nnnn t PPGMUU −∆−= +−++ 1111~ (2.151)

e, inserindo na Eq. (2.147), temos

( ) ( )1 1 1 11 n n n n n n n nt

t

θ θ θ+ − + + + + + − ∆ − − + + = ∆M U M G P P U K U U GP F (2.152)

Que após o rearranjo dos termos resulta

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

t

θ θ θ+ + + + +− + + =∆

M U U K U U GP F (2.153)

44

Que juntamente com as Eq. (2.148), e Eq. (2.150) formam a formulação monolítica, escri-

ta novamente abaixo:

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

t

θ θ θ+ + + + +− + + =∆

M U U K U U GP F (2.154)

( ) ( )1 1 1 1n n n n nt τ+ + + +∆ − + − = −L P P LP DΕ DU (2.155)

1 1n n+ +=MΕ GP (2.156)

Note que o sistema formado pelas Eqs. (2.154-2.156) é o mesmo sistema formado pelas

Eqs. (2.133-2.135), a menos do termo ( )nnt PPL −∆ +1 na Eq. (2.155). Este termo é positivo

definido (possui todos os autovalores e todos os pivots positivos), pois L é a matriz Laplaci-

ana, portanto, pode ser vista como um pré-condicionador do sistema formado pelas Eqs.

(2.154-2.156), permitindo que uma solução desacoplada, utilizando um método iterativo, do

tipo Gauss-Seidel seja aplicado, e o sistema segregado converge para a solução do sistema

monolítico com a convergência do método iterativo Gauss-Seidel. Para soluções transientes, o

termo adicional, na Eq. (2.155) leva à soluções diferentes das obtidas com a Eq. (2.134), pois

note que para soluções transientes, não está garantido o fato de que a restrição de incompres-

sibilidade seja satisfeita a cada passo de tempo, obviamente em casos permanentes este pro-

blema não ocorre, pois em regime permanente 1n n+ =P P . Estas diferenças podem ser

minimizadas, se o termo for re-escrito da seguinte maneira: ( )1,1,1 −++ −∆ inint PPL , onde o su-

per-índice i refere-se à iteração Gauss-Seidel no intervalo de tempo atual. Desta forma, o sis-

tema monolítico desacoplado resultante é:

( ) ( )1, , 1 , 1, 11 n i n n i n i n i n

t

θ θ θ+ + − + + − +− + + =∆

M U U K U U GP F (2.157)

( ) ( )1, 1, 1 1, 1, 1 1,n i n i n i n i n it τ+ + − + + − +∆ − + − = −L P P LP DΕ DU (2.158)

1, 1,n i n i

l

+ +=M Ε GP (2.159)

onde lM~

é a matriz de massa “lumped”, e note que o termo convectivo foi linearizado.

45

Note que neste ponto, a equação de Poisson para a pressão está estabilizada, mas ainda

falta estabilizar as equações de conservação da quantidade de movimento, o que pode ser fei-

to, seguindo o mesmo raciocínio utilizado para a pressão, ou seja, utilizar uma estabilização

de alta ordem do tipo SUPG (Streamline Upwind Petrov Galerkin), baseada na estabilização

de alta ordem utilizada para a pressão.

2.4 Formulação Segregada Final

O problema variacional discreto final, incluindo estabilização nas diversas equações [Soto

et al., 2001], pode ser escrito como:

Dado n

hu , encontrar ( )1111 ,,, ++++ n

h

n

h

n

h

n

h p ξπu em hhhh Q VVV~~

××× tal que

( ) ( )h

in

h

in

hh

n

h

in

ht

vuuvuu ,,1 ,1,1,1 θ

δ+−++ ∇⋅+− + ( ) ( ), 1, 1, ,n i n i

h h h hpθυ + + −∇ ∇ + ∇u v v

+ ( )( )h

in

h

in

h

inin

h

in

h vuπuu ∇⋅−∇⋅ −+−+−++−+ 1,1,1,1,1, , θθθθ βτ = ( ) ( )nu

h

in

h

n

Γ

−++ ⋅+ vnσvf ,, 1,θθ (2.160)

( ) ( )( )h

in

h

inin

hh

in

h

in

h qpqppt ∇−∇+∇∇−∇ −+−++−++ ,, 1,11,1,11,1,1 ξβτδ = - ( )h

in

h q,,1+⋅∇ u (2.161)

( ) ( )h

in

h

in

hh

in

h vuuvπ ~,~, ,,, θθθ +++ ∇⋅= (2.162)

( ) ( )h

in

hh

in

h p vvξ ~,~, ,1,1 ++ ∇= (2.163)

( ) hhhhhhhh Qq VVVvvv~~~,~,, ×××∈∀ (2.164)

onde 1 1 1 1, , ,n n n n

h h h hp

+ + + +u π ξ são o vetor velocidade, a pressão, a projeção do termo convectivo e a

projeção do gradiente de pressão no espaço de elementos finitos, respectivamente, e

hhhh Q VVV~~

××× são os espaços vetoriais correspondentes às variáveis acima citadas, respecti-

vamente. Os sobrescritos n e i referem-se ao passo de tempo e ao contador de iterações em

cada passo de tempo, respectivamente, o subescrito h refere-se à variáveis discretas, δt é o ta-

manho do passo de tempo, β é o parâmetro de estabilização na presença de gradientes acentu-

46

ados, controlado pelo gradiente de pressão, cuja variação restringe-se ao intervalo I=[0 1],

tomando valores próximos de 1 em regiões onde o campo de pressão é suave, e próximo de 0

em regiões com fortes gradientes de pressão (Löhner, 2001, Soto et. al., 2004), τ é o parâme-

tro que controla a estabilidade e a convergência do sistema, θ é o controlador do método de

integração no tempo (por exemplo, θ = 1.0 “Backward Euler” e θ = 0.5 “Crank-Nicholson”),

σ é o tensor das tensões viscosas, f é o vetor das forças externas aplicadas, Гnu é a porção do

contorno com condições de contorno de Newman, υ é a viscosidade cinemática do fluido, e

hv , são as funções de ponderação para a velocidade, hv~ são as funções de ponderação para as

projeções da convecção e do gradiente de pressão e, hq são as funções de ponderação para a

pressão, respectivamente nos espaços de elementos finitos.

Note que variáveis π e ξ são auxiliares e não necessitam condições de contorno e inici-

ais.

Sendo assim, o sistema de equações final dado pelas Eqs. (2.160-2.163), é resolvido da

seguinte maneira:

Enquanto o número máximo de iterações não for alcançado, faça

Obtém passo de tempo

Enquanto o resíduo for maior do que a tolerância especificada, faça

Passo 1: As equações de balanço da quantidade de movimento são resolvidas de

forma segregada para as componentes do vetor velocidade, Eq. (2.160);

Passo 2: A equação de Poisson para a pressão é resolvida com as variáveis de ve-

locidade obtidas no Passo 1, Eq. (2.161);

Passo 3: As equações da projeção da convecção são resolvidas de forma segregada

para as componentes da convecção, utilizando as variáveis de velocidade obtidas no Passo

1, Eq. (2.162);

Passo 4: As equações da projeção do gradiente de pressão são resolvidas de forma

segregada para as componentes do gradiente de pressão, utilizando a pressão obtida no

Passo 2, Eq. (2.163);

Passo 5: A convergência do sistema Gauss-Seidel é testada com o resíduo da pres-

são e da velocidade;

47

2.5 Termos de Estabilização Espacial

Considerando as dificuldades associadas à obtenção da solução das Equações de Navier-

Stokes discutidas anteriormente, foi necessário incorporar ao sistema de equações os termos

de estabilização das equações de Quantidade de Movimento, baseado na mesma estabilização

da pressão obtida através da formulação fracionada, como pode ser verificado a seguir.

2.5.1 Estabilização Espacial da Equação de Poisson para a Pressão

O termo que caracteriza a estabilização da equação de Poisson para a Pressão, Eq. (2.161),

é dado por:

( )( )h

in

h

inin

h qp ∇−∇ −+−++ ,1,11,1,1 ξβτ (2.165)

Este termo foi obtido de forma algébrica através do fracionamento das equações. Note que

este termo, em sua forma discreta resulta em uma matriz positiva-definida simétrica, portanto

servindo aqui de pré-condicionamento da equação da Pressão. A quantidade de estabilização

está sendo controlada pelo parâmetro artificial β que determina as regiões de gradientes de

pressão elevados e desta forma controla o condicionamento do sistema. Obviamente este pa-

râmetro deverá estar contido em [ ]1,0=I .

2.5.2 Estabilização Espacial das Equações de Quantidade de Movimento

O termo que caracteriza a estabilização da equação de Quantidade de Movimento, Eq.

(2.160) é dado por:

( )( )h

in

h

in

h

inin

h

in

h vuπuu ∇⋅−∇⋅ −+−+−++−+ 1,1,1,1,1, , θθθθ βτ (2.166)

Note que este termo é o mesmo termo de SUPG utilizado por muitos autores [Brooks e

Hughes, 1982, Cebral et al., 2001], que pode ser obtido de diversas formas. Este termo consis-

te em adicionar difusão total, diminuindo para primeira ordem a aproximação, e retirar nas re-

48

giões onde a solução mostra-se suave o suficiente para permitir soluções de segunda-ordem.

A quantidade de difusão necessária é controlada pelo mesmo parâmetro de captura de gradi-

entes elevados β , utilizado na equação da Pressão, que será definido mais adiante. Neste pon-

to é importante frisar que uma formulação do tipo SUPG, de ordem 1/ 20( )ph + , envolve a

ponderação de todo o resíduo, incluindo o termo transiente e o termo viscoso (uma formula-

ção do tipo Petrov-Galerkin requer uma função de ponderação modificada aplicada em toda a

equação). Portanto neste caso, como estamos tratando apenas o termo convectivo, esta formu-

lação caracteriza-se como uma formulação Streamline Upwind, ou seja, que apresenta estabi-

lização ao longo das linhas de corrente.

2.6 Parâmetros Importantes

Nesta seção serão desenvolvidos os parâmetros controladores da estabilização nas equa-

ções discretizadas. Estes parâmetros são obtidos de forma automática, sem a interferência do

usuário, e são totalmente baseados na solução em cada passo de tempo e nos parâmetros físi-

cos e numéricos do problema analisado.

2.6.1 Termo de Preservação da Monotonicidade

A formulação final apresentada nas Eqs. (2.162-2.165) apresenta alta ordem, portanto não

é capaz de evitar oscilações espúrias que ocorrem em fluxos com convecção dominante. Estas

oscilações afetam de forma destrutiva a convergência da formulação, sendo necessária a in-

clusão de termos que detectem regiões no escoamento que apresentam acentuados gradientes

de pressão, e nestas regiões baixar a ordem da aproximação evitando desta forma a ocorrência

destas oscilações. Para localizar as regiões com altos gradientes foi utilizado um parâmetro

controlador, baseado na diferença entre o gradiente de pressão e a sua projeção no espaço de

elementos finitos [Soto et al., 2004]. Este parâmetro é obtido para cada direção, para uma es-

tutura de dados baseada nas arestas dos elementos tetraédricos da seguinte maneira:

, , ,

,

, , ,

0.5 ( )1

0.5 ( )

i j k ij k i k j

k ij

i j k ij k i k j

p p l

p p l

ξ ξβ

ξ ξ

− + += −

− + + (2.167)

49

onde ip é a pressão no ponto nodal i , ijkl , é a dimensão da aresta ij na direção k , e ik ,ξ -e a

k-ésima componente da projeção do gradiente de pressão no ponto nodal i . É possível verifi-

car que 10 ≤≤ β , assumindo valores próximos de 1 onde o campo de pressão é suave, e pró-

ximo de 0 em regiões de elevados gradientes de pressão. Como verificado nas Eqs. (2.154), o

parâmetro funciona como um controlador de anti-difusão, ou seja, em regiões de elevados

gradientes de pressão, o termo SUPG é incluído nas equações de conservação da quantidade

de movimento e a estabilização da pressão é mantida na equação de Poisson da pressão. É im-

portante ressaltar que sensores de pressão são utilizados por muitos autores para evitar oscila-

ções numéricas em fluxos que apresentam convecção dominante.

2.6.2 Passo de Tempo para Estabilidade Temporal

O passo de tempo, obtido por aresta, é calculado seguindo a expressão apresentada em Zi-

enkiewicz et al. (1999), dada a seguir

t tt

t t

σ ν

σ ν

δ δδ

δ δ=

+ (2.168)

com

E

E

ht

uσδ = (2.169)

2

2Eh

tνδµ

= (2.170)

onde as Eqs. (2.169) e (2.170) são os limites de estabilidade dos fenômenos advectivos e difu-

sivos, respectivamente. Neste trabalho estamos utilizando uma estrutura de dados baseada nas

arestas dos elementos, então, Eh é a dimensão característica do elemento, neste caso o com-

primento da aresta, Eu é a velocidade absoluta na aresta, e µ é viscosidade dinâmica. Este

intervalo de tempo é usado como referência, pois como a formulação é implícita, em diversas

situações o avanço no tempo pode ser mais rápido sem prejudicar a qualidade da solução.

50

2.6.2.1 Estabilização Temporal das Equações de Quantidade de Movimento

O parâmetro τ que controla o termo de estabilização nas equações de Quantidade de Mo-

vimento é obtido de forma a garantir a consistência da aproximação, sendo assim, este parâ-

metro deve ser obtido como dependente apenas do nó i de cada aresta. Se isto não

acontecesse, ou seja, se este parâmetro fosse obtido por aresta, com dependência também do

no j da aresta, então valores nodais exatos da velocidade não iriam satisfazer a equação dis-

creta. Isto pode ser evidenciado pelo fato de que a soma considerando-se o nó j não reproduz

o termo convectivo se as derivadas parciais das funções de forma jN forem multiplicadas por

um fator que depende também de j [Codina, 2000]. Neste trabalho este parâmetro foi obtido

da seguinte maneira:

2

4 2i

i

i i i

h

hτ

υ=

+ a (2.171)

onde ih é o comprimento mínimo considerando todas as arestas que têm conectividade com o

nó i, iυ é a viscosidade cinemática no ponto nodal i, e ia é a velocidade de convecção no

ponto nodal i.

2.6.2.2 Estabilização Temporal das Equações de Poisson para a Pressão

Por razões de consistência, o parâmetro que controla a estabilidade e a convergência do

sistema, τ , deveria ser obtido por nó, e não por aresta, na Eq. (2.161), como foi obtido para as

equações de Quantidade de Movimento. Entretanto, isto destruiria a simetria do sistema. Des-

ta forma como os termos de estabilização da pressão são no mínimo de 4ª ordem, o que signi-

fica que os seus efeitos na precisão do sistema final são muito pequenos, optou-se por garantir

simetria e conservação ao preço de perder alguma consistência. Sendo assim, neste trabalho

este parâmetro foi obtido da seguinte maneira:

( ) ( )

2 2

4 2

i j

ij

i j i i j j

h h

h hτ

υ υ

+=

+ + +a a (2.172)

51

onde as variáveis na expressão acima são as mesmas definidas anteriormente para a

estabilização temporal das equações de Quantidade de Movimento. Note que neste caso, o

termo na aresta está sendo obtido como uma média dos termos de estabilização nos nós.

52

Capítulo 3

Aspectos Numérico-Computacionais

3.1 Introdução e Considerações Gerais

Neste capítulo apresentaremos alguns aspectos computacionais relacionados à solução de

problemas de escoamentos fluidos incompressíveis, bem como considerações sobre soluções

numéricas de problemas de escoamentos fluidos e suas principais características computacio-

nais.

Antes de apresentar uma descrição detalhada das técnicas da Dinâmica dos Fluidos Com-

putacional, apresentaremos alguns aspectos importantes sobre técnicas de simulação. Técnicas

de simulação computacional são utilizadas por engenheiros e cientistas das mais diversas á-

reas do conhecimento para obter informações sobre manufaturas de produtos, ou sobre o

comportamento físico dos mais variados fenômenos quando submetidos a determinadas con-

dições de contorno medidas ou assumidas, a estados iniciais, a carregamentos, etc. Uma gran-

de variedade de razões pode ser citada para justificar a importância que as técnicas de

simulação têm alcançado nas recentes décadas, dentre outras, temos:

(a) Necessidade de Melhoramento de Desempenho – a falta de habilidade em prever precisa-

mente o desempenho de um novo produto pode ser crucial para causar um efeito devastador

53

em empresas. Um fabricante de avião ou carro, por exemplo, precisa construir protótipos e

submetê-los a exaustivos e custosos testes de desempenho antes de lançá-lo no mercado. Sem

as simulações computacionais tais testes tornar-se-iam impraticáveis, o que aumentaria os ris-

cos de rejeição dos usuários. A forma mais eficaz de minimizar os riscos de uma performance

inesperada e obter o máximo de informações sobre o produto a ser lançado, consiste na utili-

zação de técnicas de simulação computacional (“prototipagem virtual”).

(b) Custo dos Experimentos – experimentos, fórmulas empíricas e analíticas constituem ou-

tras alternativas além de simulações computacionais, porém estas podem tornar-se custosas

e/ou imprecisas e até mesmo impossíveis.

(c) Impossibilidade de Experimentos – alguns experimentos são impossíveis de serem realiza-

dos, por exemplo, eventos galáticos ou solares, explosões atmosféricas nucleares, situações

biomédicas que colocam em risco a vida de pacientes, e muitas outras situações.

(d) Avanço Computacional – Nos últimos anos estamos presenciando um avanço tecnológico

em termos de velocidade de processamento e capacidade de memória dos computadores. Es-

tes avanços associados aos avanços dos métodos numéricos permitem que as simulações se-

jam cada vez mais realísticas, garantindo a precisão dos resultados.

Desta forma fica evidenciada a necessidade e importância da realização de simulações

computacionais no estudo de fenômenos físicos e no desenvolvimento de novas tecnologias.

Aliado a isso, quando procuram por uma descrição quantitativa de fenômenos físicos, os en-

genheiros e cientistas estabelecem modelos matemáticos, em geral, dados por sistemas de e-

quações diferenciais ordinárias ou parciais válidas em uma certa região, ou domínio,

submetidos à certas condições de contorno e iniciais. Neste estágio o modelo matemático está

completo, e para as aplicações práticas é necessário encontrar um conjunto de dados que satis-

façam este sistema juntamente com suas condições de contorno e iniciais para que se obtenha

a solução para a aplicação estudada. Entretanto, somente para poucas situações simples e mui-

to simplificadas é possível resolver analiticamente os sistemas de equações resultantes que

descrevem fenômenos físicos com alto grau de realismo. Para sobrepor esta dificuldade, e

possibilitar o estudo de problemas complexos, utiliza-se a ferramenta mais poderosa desen-

volvida no século XX, o computador digital, que resolve problemas de forma puramente algé-

brica, envolvendo apenas as operações aritméticas básicas. Para alcançar isso, várias formas

de discretização dos problemas descritos como contínuos e definidos pelas equações diferen-

ciais podem ser usadas. Em tais discretizações um conjunto infinito de números representando

54

as variáveis do problema é trocado por um grupo finito de variáveis desconhecidas, e esse

processo envolve diversas formas de aproximação.

Neste trabalho utilizou-se o método de discretização por elementos finitos, onde o contí-

nuo é dividido em elementos e as equações governantes do problema em questão são integra-

das espacialmente em cada elemento. O procedimento de obtenção da solução consiste em

acumular as contribuições para os nós de cada elemento em uma matriz, resultando em um

sistema algébrico de equações, o qual é resolvido utilizando-se algum método de solução de

sistemas lineares ou não lineares.

A seguir descreveremos a discretização espacial utilizando-se o Método dos Elementos

Finitos, evidenciando o fato de que neste trabalho estamos utilizando uma estrutura de dados

baseada nas arestas dos elementos tetraédricos computacionais e ressaltando as características

dessa estrutura de dados em temos de implementação computacional.

3.2 Discretização Espacial via Método dos Elementos Finitos

É importante notar que nas equações de balanço da quantidade de movimento, Eq. (2.6) há

simultaneamente termos convectivos (primeira ordem) e termos difusivos (segunda ordem). A

presença do termo convectivo quebra a simetria do operador diferencial, tornando-o não auto-

adjunto, e neste caso a propriedade da melhor aproximação (Hughes, 1987) é perdida quando

o método de Galerkin padrão é utilizado para realizar a discretização espacial, e oscilações

espaciais podem surgir quando o método de Galerkin é aplicado a problemas com convecção

dominante. Na Eq. (2.6) a contribuição convectiva é não-linear, envolvendo o produto da ve-

locidade com o seu gradiente, porém na prática, em formulações explícitas, esta equação é li-

nearizada tomando-se, por exemplo, o campo de velocidade no passo anterior e adotando-se

um processo iterativo.

Uma segunda dificuldade refere-se à necessidade de escolher interpolações compatíveis

para velocidade e pressão, a fim de satisfazer a condição de Babuska-Brezzi (De Sampaio,

1991) e obter soluções convergentes com formulações mistas para escoamentos viscosos in-

compressíveis. Contudo, as restrições para a escolha das interpolações podem ser evitadas u-

tillizando-se métodos estabilizados como os propostos por Hughes et al (1986), De Sampaio

(1991), e Zienkiewicz e Wu (1991). As dificuldades mencionadas acima levam a adoção de

formulações mistas tipo Petrov-Galerkin do método dos elementos finitos, que diferem da

formulação padrão de Galerkin, e tentam através da mudança das funções de ponderação, com

55

a introdução de funções que causem efeito “upwinding”, estabilizar o método. Esta nova fun-

ção de ponderação é composta pela função clássica de Galerkin mais uma perturbação respon-

sável pelo efeito “upwinding”. Neste capítulo será apresentada brevemente uma formulação

de Elementos Finitos Estabilizada, proposta por Soto et al., (2001), que introduz explicita-

mente termos de estabilização para a convecção, dispensando a necessidade de atender a con-

dição de Babuska-Brezzi permitindo a mesma ordem de interpolação para velocidade e

pressão, o que se torna atrativo do ponto de vista de implementação computacional.

3.3 Implementação Computacional

As matrizes resultantes da discretização espacial por elementos finitos são obtidas medi-

ante uma estrutura de dados baseada nas arestas dos elementos tetraédricos (Löhner, 2001),

que tem como principal vantagem sobre uma estrutura baseada nos elementos, o ganho de

tempo computacional, pois muitos termos discretos são computados no início, como pré-

processamento, considerando-se malhas fixas, e apenas algumas poucas operações necessitam

serem implementadas para que sejam re-computados os diferentes termos não-lineares do sis-

tema. Também, através de uma implementação por arestas, é possível garantir a conservação

local e a simetria em nível discreto (Soto et al., 2004). Duas formas de estrutura de dados por

arestas básicas foram implementadas e comparadas no caso bidimensional (veja Apêndice C),

sendo uma delas a estrutura tradicional, onde o laço é feito na aresta e nos nós que compõem

tal aresta, e outra, uma estrutura alternativa, tipo estrela, onde o laço é realizado nos nós e, a

seguir, nos nós que têm conectividade com aquele nó (Löhner, 2001, Soto et al., 2004). Em

termos de eficiência computacional, nenhuma das duas estruturas apresentou vantagem sobre

a outra, e, desta forma neste trabalho optou-se pela tradicional, ou seja, o loop é realizado nas

arestas. Outras estruturas de dados por arestas são ainda possíveis, tais como “chains”, “super

edges” (Löhner, 2001), e que podem ter vantagens em termos de redução de perda de cache,

porém são de manipulação complexa. A seguir evidenciamos o fato de que a mudança de uma

estrutura de dados por elementos para uma estrutura de dados por arestas é simplesmente uma

forma de visualizar as contribuições de cada elemento como as contribuições de cada aresta

deste mesmo elemento.

Seja a malha computacional dada na Fig.(1)

56

Figura 1. Representação de uma malha computacional composta de tetraédros.

onde

1IJKL

EΩ = - elemento formado pelos nós IJKL

2IJLM

EΩ = - elemento formado pelos nós IJLM

3IJMP

EΩ = - elemento formado pelos nós IJMP

4IJPK

EΩ = - elemento formado pelos nós IJPK

Assim, por exemplo, considerando a aresta IJ na Fig.(1), notamos que esta aresta está

contida nos elementos IJKLE , IJLME , IJMPE e IJPKE . A partir da formulação padrão de elemen-

tos, define-se um desmembramento dos coeficientes gerando as contribuições da aresta IP ,

que podem ser generalizada para todas as arestas da malha computacional.

Consideremos um único elemento tetraédrico da malha conforme Fig.(2).

57

Figura 2. Elemento tetraédrico padrão.

O elemento da Fig.(2) dá origem a uma matriz que considera as contribuições das arestas

deste elemento para cada nó I, J, K, L, dada conforme:

0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0

0 0

0 0

0 0 0 0

IJKL IJ IK IL

JK

E A A A

A

• • • • × × × × × × • • • • × × = + + • • • • × × • • • • × ×

× × + + × ×

0 0 0 0 0

0 0 0 0

0 0 0 0 0 0 0 0

0 0 0 0JL KLA A

× × × × +

× × × ×

onde • , × e 0 representam as submatrizes de ordem 6. Além disso, sobre a aresta ip incidem

informações dos elementos adjacentes a ela, como pode ser observado na Fig.(2). Assim, to-

das as contribuições referentes à aresta IJ , estão presentes nos elementos IJKLE , IJLME , IJMPE

e IJPKE .

Seja então o termo convectivo das equações de Quantidade de Movimento dada pela Eq.

(2.6), considerando a contribuição de todos os elementos discretos da malha computacional,

sendo N a função de forma de elementos finitos geral, e , , ,I J K LN N N N as funções de forma

referente aos nós I, J, K e L, respectivamente:

58

( ) ( ) ˆ, ( )E

J J I E J

E J

N Nd N N d⊃Ω Ω

⋅∇ = ⋅∇ Ω = ⋅∇ Ω =∑∫ ∫u u u u u u

( ) ( ) ( ) ( )( ) ( ) ( ) ( )( ) ( ) ( ) ( )( ) ( ) ( ) ( )

ˆ

ˆ

ˆ

ˆE

I I I J I I K I I L I I I

I J J J J J K J J L J J J

E

I K K J K K K K K L K K K

I L L J L L R L L L L L L

N N N N N N N N

N N N N N N N Nd

N N N N N N N N

N N N N N N N N

Ω

⋅∇ ⋅∇ ⋅∇ ⋅∇ ⋅∇ ⋅∇ ⋅∇ ⋅∇ Ω ⋅∇ ⋅∇ ⋅∇ ⋅∇

⋅∇ ⋅∇ ⋅∇ ⋅∇

∫

u u u u u

u u u u u

u u u u u

u u u u u

(3.1)

Esta matriz pode agora ser vista como uma soma de todas as matrizes que contém as con-

tribuições das arestas da seguinte maneira:

( ) ( )

( ) ( )

/ 3 0 0

0 0 0 0

0 0 0 0

0 0 / 3

I I I L I I

I L L L L L

N N N N

N N N N

⋅∇ ⋅∇ +

⋅∇ ⋅∇

u u

u u

( ) ( )

( ) ( )

/ 3 0 0

0 0 0 0

0 / 3 0

0 0 0 0

I I I K I I

I K K K K K

N N N N

N N N N

⋅∇ ⋅∇ + ⋅∇ ⋅∇

u u

u u

( ) ( )( ) ( )

+

∇⋅∇⋅

∇⋅∇⋅

0000

0000

003/

003/

JJJJJI

IIJIII

NNNN

NNNN

uu

uu

( ) ( )( ) ( )

0 0 0 0

0 0 0 0

0 0 / 3

0 0 / 3K K K L K K

K L L L L L

N N N N

N N N N

+ ⋅∇ ⋅∇

⋅∇ ⋅∇

u u

u u

( ) ( )( ) ( )

0 0 0 0

0 / 3 0

0 / 3 0

0 0 0 0

J J J K J J

J K K K K K

N N N N

N N N N

⋅∇ ⋅∇ + ⋅∇ ⋅∇

u u

u u

( ) ( )

( ) ( )

0 0 0 0

0 / 3 0

0 0 0 0

0 0 / 3

J J J L J J

J L L L L L

N N N N

N N N N

⋅∇ ⋅∇

⋅∇ ⋅∇

u u

u u

( ) ( ) ( ) ( )( ) ( ) ( ) ( )( ) ( ) ( ) ( )( ) ( ) ( ) ( )

I I I J I I K I I L I I

I J J J J J K J J L J J

I K K J K K K K K L K K

I L L J L L K L L L L L

N N N N N N N N

N N N N N N N N

N N N N N N N N

N N N N N N N N

⋅∇ ⋅∇ ⋅∇ ⋅∇

⋅∇ ⋅∇ ⋅∇ ⋅∇ = ⋅∇ ⋅∇ ⋅∇ ⋅∇

⋅∇ ⋅∇ ⋅∇ ⋅∇

u u u u

u u u u

u u u u

u u u u

(3.2)

Desta forma é possível verificar que a estrutura de dados por aresta requer um “loop” nas

arestas dos elementos tetraédricos, acumulando-se a contribuição desta matriz no sistema al-

gébrico final. Isto requer também que todos os termos das equações, que são classicamente

escritos em função dos elementos, devem ser obtidos como contribuição das arestas.

59

A seguir é apresentado como é realizada a mudança de estrutura de dados, partindo do

termo contínuo para o termo discreto com loop nos elementos e modificando para loop nas a-

restas.

Seja então o termo da Eq. (2.160) que possui caráter de transporte difusivo, eliminando os

sub/super-índices por simplicidade, e assumindo que E, S, I, J, IJ, k, ndim, Ω , EΩ , kx , m,

E I⊃∑ ,

E IJ⊃

∑ , 1

m

S =

∑ , referem-se a: elemento, aresta, nó inicial da aresta, nó final da aresta, ares-

ta do nó I a J, direção coordenada, dimensão espacial, domínio total, domínio do elemento,

coordenada espacial na direção k, número de arestas, somatório sobre todos os elementos que

contém o nó I, somatório sobre todos os elementos que contém a aresta IJ, e somatório sobre

todas as arestas, respectivamente. Seja então o termo:

ˆ ˆ( , )E E

E E E E E

J I E I J E J

E J E J

N Nd N d N N dυ υ υ υ⊃ ⊃Ω Ω Ω

∇ ∇ = ∇ ∇ Ω = ∇ ⋅∇ Ω = ∇ ⋅∇ Ω∑ ∑∫ ∫ ∫u u u u

ˆE

E E EE E EEJ J JI I I

E J

E J

N N NN N Nd

x x y y z zυ

⊃ Ω

∂ ∂ ∂∂ ∂ ∂= + + Ω

∂ ∂ ∂ ∂ ∂ ∂ ∑ ∫ u (3.3)

onde υ foi considerado constante, e ˆ E

Ju é o valor nodal aproximado da variável de velocida-

de, também considerado constante no elemento.

Para um elemento tetraédrico linear como o da Fig.(2), a Eq. (3.3) resulta:

ˆ

ˆ

ˆ

ˆE

E E E E E E E E E

I I I J I K I L I

E E E E E E E E E

J I J J J K J L J

EE E E E E E E E EE K I K J K K K L K

E E E E E E E E E

L I L J L K L L L

N N N N N N N N

N N N N N N N Nd

N N N N N N N N

N N N N N N N N

υΩ

∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ Ω ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇

∑ ∫

u

u

u

u

(3.4)

que é a matriz viscosa para o elemento tetraédrico linear.

Então para a malha apresentada na Fig(1), temos

1 2

1 1 1 2 2 21 2

E

E l m l l m l

E I

Nd N N d N N dυ υ υ⊃ Ω Ω Ω

∇ ⋅∇ Ω = ∇ ⋅∇ Ω + ∇ ⋅∇ Ω +∑ ∫ ∫ ∫u u u

3 4

3 3 3 4 4 43 4l m l l m lN N d N N dυ υ

Ω Ω

∇ ⋅∇ Ω + ∇ ⋅∇ Ω∫ ∫u u (3.5)

60

e escrevendo as matrizes de cada elemento,

1

1 1 1 1 1 1 1 1 1

1 1 1 1 1 1 1 1 1

11 1 1 1 1 1 1 1 1

1 1 1 1 1 1 1 1 1

ˆ

ˆ

ˆ

ˆ

I I I J I K I L I

J I J J J K J L J

K I K J K K K L K

L I L J L K L L L

N N N N N N N N

N N N N N N N Nd

N N N N N N N N

N N N N N N N N

υΩ

∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ = Ω + ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇

∫

u

u

u

u

2

2 2 2 2 2 2 2 2 2

2 2 2 2 2 2 2 2 2

22 2 2 2 2 2 2 2 2

2 2 2 2 2 2 2 2 2

ˆ

ˆ

ˆ

ˆ

I I I J I L I M I

J I J J J L J M J

L I L J L L L M L

M I M J M L M M M

N N N N N N N N

N N N N N N N Nd

N N N N N N N N

N N N N N N N N

υΩ

∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇

∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ Ω + ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇

∫

u

u

u

u

3

3 3 3 3 3 3 3 3 3

3 3 3 3 3 3 3 3 3

13 3 3 3 3 3 3 3 3

3 3 3 3 3 3 3 3 3

ˆ

ˆ

ˆ

ˆ

I I I J I M I P I

J I J J J M J P J

M I M J M M M P M

P I P J P M P P P

N N N N N N N N

N N N N N N N Nd

N N N N N N N N

N N N N N N N N

υΩ

∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇

∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ Ω + ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇

∫

u

u

u

u

4

4 4 4 4 4 4 4 4 4

4 4 4 4 4 4 4 4 4

14 4 4 4 4 4 4 4 4

4 4 4 4 4 4 4 4 4

ˆ

ˆ

ˆ

ˆ

I I I J I P I K I

J I J J J P J K J

P I P J P P P K P

K I K J K P K K K

N N N N N N N N

N N N N N N N Nd

N N N N N N N N

N N N N N N N N

υΩ

∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ Ω = ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇ ∇ ⋅∇

∫

u

uG

u

u

(3.6)

e montando a matriz global para a malha dada, fazendo:

1 1 2 2 3 3 4 4; ; ;Ω = Χ Ω = Χ Ω = Χ Ω = Χ

Logo, a matriz G resulta

61

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

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

1 4 1 4 1 4 1 4

1 2 1 2 1 1 2 2

2 3 2 3 2 2 3 3

3 4 3 4 4 3 3 4

ˆ

ˆ

ˆ

ˆ

ˆ

I

J

K

L

Χ + Χ + Χ + Χ Χ + Χ + Χ + Χ Χ + Χ Χ + Χ Χ + Χ Χ + Χ Χ + Χ + Χ + Χ Χ + Χ + Χ + Χ Χ + Χ Χ + Χ Χ + Χ Χ + Χ Χ + Χ Χ + Χ Χ + Χ Χ − Χ

= Χ + Χ Χ + Χ Χ Χ + Χ Χ −

Χ + Χ Χ + Χ − Χ Χ + Χ Χ

Χ + Χ Χ + Χ Χ − Χ Χ + Χ

u

u

uG

u

u

ˆM

P

u

(3.7)

Note que a matriz G é simétrica, o que já era esperado, devido à simetria do próprio termo

viscoso. Escrevendo a equação para o nó I da malha tem-se:

( ) ( ) ( ) ( )

( ) ( )1 2 3 4 1 2 3 4 1 4 1 2

2 3 3 4

ˆ ˆ ˆ ˆ

ˆ ˆ

I J K L

M P

⇒ Χ + Χ + Χ + Χ + Χ + Χ + Χ + Χ + Χ + Χ + Χ + Χ +

Χ + Χ + Χ + Χ

u u u u

u u (3.8)

e expandindo os termos da expressão acima

1 2 3 4

1 1 2 2 3 3 4 41 2 3 4 ˆ

I I I I I I I I IN N d N N d N N d N N dυ υ υ υΩ Ω Ω Ω

⇒ ∇ ⋅∇ Ω + ∇ ⋅∇ Ω + ∇ ⋅∇ Ω + ∇ ⋅∇ Ω +

∫ ∫ ∫ ∫ u

1 2 3 4

1 1 2 2 3 3 4 41 2 3 4 ˆ

I J I J I J I J JN N d N N d N N d N N dυ υ υ υ

Ω Ω Ω Ω

∇ ⋅∇ Ω + ∇ ⋅∇ Ω + ∇ ⋅∇ Ω + ∇ ⋅∇ Ω +

∫ ∫ ∫ ∫ u

1 4

1 1 4 41 4 ˆ

I K I K KN N d N N dυ υ

Ω Ω

∇ ⋅∇ Ω + ∇ ⋅∇ Ω +

∫ ∫ u

1 2

1 1 2 21 2 ˆ

I L I L LN N d N N dυ υ

Ω Ω

∇ ⋅∇ Ω + ∇ ⋅∇ Ω +

∫ ∫ u

1 3

2 2 3 32 3 ˆ

I M I M MN N d N N dυ υ

Ω Ω

∇ ⋅∇ Ω + ∇ ⋅∇ Ω +

∫ ∫ u

62

3 4

3 3 4 43 4 ˆ

I P I P PN N d N N dυ υ

Ω Ω

∇ ⋅∇ Ω + ∇ ⋅∇ Ω

∫ ∫ u (3.9)

Agora, escrevendo os termos da Eq. (3.9) de forma compacta, considerando a expressão

para o nó I, tem-se

ˆ ˆE E

E E E E

I I j E j I I E I

j I E Ij E I

r N N d N N dυ≠ ⊃ ⊃Ω Ω

= ∇ ⋅∇ Ω + ∇ ⋅∇ Ω

∑ ∑ ∑∫ ∫u u (3.10)

onde, para a malha dada, e considerando o nó I,

1 2 3 4

1 1 2 2 3 3 4 41 2 3 4

ˆ

ˆ

E

E E

I I E I

E I

I I I I I I I I I

N N d

N N d N N d N N d N N d

⊃ Ω

Ω Ω Ω Ω

∇ ⋅∇ Ω =

∇ ⋅∇ Ω + ∇ ⋅∇ Ω + ∇ ⋅∇ Ω + ∇ ⋅∇ Ω

∑ ∫

∫ ∫ ∫ ∫

u

u (3.11)

e,

+

Ω∇⋅∇+Ω∇⋅∇

+

Ω∇⋅∇+Ω∇⋅∇

+

Ω∇⋅∇+Ω∇⋅∇

+

Ω∇⋅∇+Ω∇⋅∇

+

Ω∇⋅∇+Ω∇⋅∇+Ω∇⋅∇+Ω∇⋅∇

=

Ω∇⋅∇

∫∫

∫∫

∫∫

∫∫

∫∫∫∫

∑ ∑ ∫

ΩΩ

ΩΩ

ΩΩ

ΩΩ

ΩΩΩΩ

≠ ⊃ Ω

PPIPI

MMIMI

LLILI

KKIKI

JJIJIJIJI

j

Ij IjE

E

E

j

E

I

udNNdNN

udNNdNN

udNNdNN

udNNdNN

udNNdNNdNNdNN

udNN

E

ˆ

ˆ

ˆ

ˆ

ˆ

ˆ

13

32

21

11

1321

444

333

333

222

222

111

444

111

444

333

222

111

(3.12)

Notando que a partir da propriedade de conservação para as funções de forma, tem-se

63

∑=

=4

1

1i

E

iN (3.13)

e, portanto, derivando em relação a kx

4 4

1 1

0E

E ii

i ik k

NN

x x= =

∂∂ = = ∂ ∂

∑ ∑ (3.14)

podemos ainda escrever

∑≠ ∂

∂−=

∂

∂

ij k

E

j

k

E

i

x

N

x

N (3.15)

e, desta forma, o segundo termo de Ir , Eq. (3.10), pode ser escrito como

I

IE

E

Ij

E

j

E

II

IE

E

E

I

E

I udNNudNN

EE

ˆˆ ∑ ∫ ∑∑ ∫⊃ Ω ≠⊃ Ω

Ω

∇−⋅∇=Ω∇⋅∇ (3.16)

que para o nó I, resulta

( )

( ) ( )

( )

1

2 3

4

1 1 1 11

2 2 2 2 3 3 3 32 3

4 4 4 44

ˆ

ˆ

E

E E

I j E I I J K L

E I j I

I J L M I J M P

I J P K I

N N d u N N N N d

N N N N d N N N N d

N N N N d u

⊃ ≠Ω Ω

Ω Ω

Ω

∇ ⋅ − ∇ Ω = − ∇ ⋅ ∇ + ∇ + ∇ Ω +

− ∇ ⋅ ∇ + ∇ + ∇ Ω + − ∇ ⋅ ∇ + ∇ + ∇ Ω +

− ∇ ⋅ ∇ + ∇ + ∇ Ω

∑ ∑∫ ∫

∫ ∫

∫

64

1 1 1

2 2 2

3 3 3

4 4 4

1 1 1 1 1 11 1 1

2 2 2 2 2 22 2 2

3 3 3 3 3 33 3 3

4 4 4 4 4 44 4 4

I J I K I L

I J I L I M

I J I M I P

I J I P I K

N N d N N d N N d

N N d N N d N N d

N N d N N d N N d

N N d N N d N N d

Ω Ω Ω

Ω Ω Ω

Ω Ω Ω

Ω Ω Ω

= − ∇ ⋅∇ Ω + ∇ ⋅∇ Ω + ∇ ⋅∇ Ω +

∇ ⋅∇ Ω + ∇ ⋅∇ Ω + ∇ ⋅∇ Ω +

∇ ⋅∇ Ω + ∇ ⋅∇ Ω + ∇ ⋅∇ Ω +

∇ ⋅∇ Ω + ∇ ⋅∇ Ω + ∇ ⋅∇ Ω

∫ ∫ ∫

∫ ∫ ∫

∫ ∫ ∫

∫ ∫ ∫ ˆIu

(3.17)

substituindo as expressões (3.12) e (3.17) em Ir , dado pela Eq. (3.10), temos

1 2 3

4 1 4

1 2 2

3

1 1 2 2 3 31 2 3

4 4 1 1 4 44 1 4

1 1 2 2 2 21 2 2

3 33

ˆ ˆ ˆ

ˆ ˆ ˆ

ˆ ˆ ˆ

ˆ

I I J J I J J I J J

I J J I K K I K K

I L L I L L I M M

I M M

r N N d N N d N N d

N N d N N d N N d

N N d N N d N N d

N N d

υΩ Ω Ω

Ω Ω Ω

Ω Ω Ω

Ω

= ∇ ⋅∇ Ω + ∇ ⋅∇ Ω + ∇ ⋅∇ Ω +

∇ ⋅∇ Ω + ∇ ⋅∇ Ω + ∇ ⋅∇ Ω +

∇ ⋅∇ Ω + ∇ ⋅∇ Ω + ∇ ⋅∇ Ω +

∇ ⋅∇ Ω + ∇

∫ ∫ ∫

∫ ∫ ∫

∫ ∫ ∫

∫

u u u

u u u

u u u

u3 4

1 1 1

2 2 2

3 3

3 3 4 43 4

1 1 1 1 1 11 1 1

2 2 2 2 2 22 2 2

3 3 3 33 3

ˆ ˆ

ˆ ˆ ˆ

ˆ ˆ ˆ

ˆ ˆ

I P P I P P

I J I I K I I L I

I J I I L I I M I

I J I I M I

N N d N N d

N N d N N d N N d

N N d N N d N N d

N N d N N d

υ

Ω Ω

Ω Ω Ω

Ω Ω Ω

Ω Ω

⋅∇ Ω + ∇ ⋅∇ Ω −

∇ ⋅∇ Ω + ∇ ⋅∇ Ω + ∇ ⋅∇ Ω +

∇ ⋅∇ Ω + ∇ ⋅∇ Ω + ∇ ⋅∇ Ω +

∇ ⋅∇ Ω + ∇ ⋅∇ Ω +

∫ ∫

∫ ∫ ∫

∫ ∫ ∫

∫ ∫

u u

u u u

u u u

u u3

4 4 4

3 33

4 4 4 4 4 44 4 4

ˆ

ˆ ˆ ˆ

I P I

I J I I P I I K I

N N d

N N d N N d N N d

Ω

Ω Ω Ω

∇ ⋅∇ Ω +

∇ ⋅∇ Ω + ∇ ⋅∇ Ω + ∇ ⋅∇ Ω

∫

∫ ∫ ∫

u

u u u (3.18)

e podemos agrupar os termos acima da seguinte maneira:

65

( )

( )

( )

1 2 3 4

1 4

1 2

2

1 1 2 2 3 3 4 41 2 3 4

1 1 4 41 4

1 1 2 21 2

2 22

ˆ ˆ

ˆ ˆ

ˆ ˆ

I I J I J I J I J I J

I K I K I K

I L I L I L

I M

r N N d N N d N N d N N d

N N d N N d

N N d N N d

N N d

υ

υ

υ

υ

Ω Ω Ω Ω

Ω Ω

Ω Ω

Ω

= ∇ ⋅∇ Ω + ∇ ⋅∇ Ω + ∇ ⋅∇ Ω + ∇ ⋅∇ Ω − + +

∇ ⋅∇ Ω + ∇ ⋅∇ Ω − + +

∇ ⋅∇ Ω + ∇ ⋅∇ Ω − + +

∇ ⋅∇ Ω + ∇

∫ ∫ ∫ ∫

∫ ∫

∫ ∫

∫

u u

u u

u u

( )

( )

3

3 4

3 33

3 3 4 43 4

ˆ ˆ

ˆ ˆ

I M I M

I P I P I P

N N d

N N d N N dυ

Ω

Ω Ω

⋅∇ Ω − + +

∇ ⋅∇ Ω + ∇ ⋅∇ Ω − +

∫

∫ ∫

u u

u u

(3.19)

ou ainda em uma forma compacta, ver (3.12), temos

( )ˆ ˆE

E E

I I J E I J

J I E IJ

r N N dυ≠ ⊃ Ω

= ∇ ⋅∇ Ω − +

∑ ∑ ∫ u u jI ≠∀ (3.20)

e considerando uma malha computacional com estrutura de dados baseada nas arestas dos e-

lementos tetraédricos, podemos escrever, como “loop” nas arestas, o acúmulo das contribui-

ções dos elementos para a formação do sistema resultante

( )1

ˆ ˆE

mE E

I J E J I

S E IJ

N N dυ= ⊃ Ω

∇ ⋅∇ Ω −

∑ ∑ ∫ u u (3.21)

ou

( )dim

1 1

ˆ ˆE

EEm nJI

IJ E J I

S E IJ k k k

NNd

x xυ

= ⊃ =Ω

∂∂ Ω −

∂ ∂ ∑ ∑ ∑∫ u u (3.22)

e, chamando de IJ

C o coeficiente associado à aresta IJ dado por

66

E

IJE

n

k k

E

J

k

E

IIJ d

x

N

x

NC

E

Ω

∂

∂

∂

∂= ∑ ∫ ∑

⊃ Ω =

dim

1

(3.23)

resulta

( )1

ˆ ˆm

IJ IJ J I

S

Cυ=

−∑ u u (3.24)

onde IJυ é a viscosidade cinemática do fluido, que foi considerada constante no elemento te-

traédrico, e é agora computada por aresta. Neste trabalho, 2

JI

IJ

υυυ

+= . Note que a simetria

é mantida em nível discreto. E ainda u é a velocidade do fluido, N é a função de forma do

elemento tetraédrico considerado, neste caso esta função é linear e idêntica à função de pon-

deração, e ˆIu é o valor da velocidade do fluido no nó considerado.

Note que o termo final da Eq. (2.152) é naturalmente conservativo, mesmo em nível dis-

creto, o que do ponto de vista físico significa que o transporte difusivo do ponto I para o ponto

J, é igual ao transporte difusivo do ponto J para o ponto I. Essa propriedade, na presente im-

plementação, é garantida para todos os termos de cada equação do sistema resultante, fazendo

com que os termos que não são naturalmente conservativos, tornem-se conservativos, através

da obtenção os termos da diagonal principal de cada matriz resultante, de cada termo do sis-

tema, como a negativa da soma dos elementos que estão fora da diagonal principal, porém na

mesma linha (Soto et al., 2004).

Note, assim, que em uma estrutura de dados baseada nas arestas, cada aresta está associa-

da a todos os nós dos elementos que contém esta aresta, e a mudança de estrutura de dados é

realizada considerando-se os dados referentes à malha de elementos já existente. Desta forma

então, para a aresta IJ , tem-se as contribuições

3 31 1 2 2 4 4

3 31 1 2 2 4 4

IJ IJKL IJLM IJPKIJMPA E E EE

× ×× × × × × ×• • = + + + × ×× × × × × ×• •

onde • e i× identificam as contribuições de cada elemento para a formação da matiz da ares-

ta.

67

As informações topológicas e geométricas da malha são manipuladas para gerar a nova es-

trutura baseada nas arestas dos elementos. Isto é realizado em uma etapa de pré-

processamento.

Todos os termos das equações do sistema (2.151-2.154) foram escritos em função das

contribuições das arestas dos elementos de forma semelhante ao termo difusivo deduzido an-

teriormente.

Considerando a Equação de Conservação da Quantidade de Movimento:

a) Termo Transiente:

1

,20

E

mJE I

E

S E IJ

VN Nd

t t t t= ⊃Ω

∂∂∂ ∂ = Ω = + ∂ ∂ ∂ ∂

∑ ∑∫uuu u

(3.25)

onde EV é o volume do elemento tetraédrico, definido como E

EV dΩ

= Ω∫ .

b) Termo Convectivo:

( ) ( )dim

1 1

,k

E E

m nJ

E IJ I E J I

S k E IJ k

NN Nd a N d

x= = ⊃Ω Ω

∂ ⋅∇ = ⋅∇ Ω = Ω −

∂ ∑∑ ∑∫ ∫a u a u u u (3.26)

onde a é o vetor da velocidade convectiva, que deve ser definido por arestas.

c) Termo do Gradiente de Pressão:

( ) ( )1

,E E

mJ

E I E J I

S E IJ k

Np N pNd N d p p

x= ⊃Ω Ω

∂ ∇ = ∇ Ω = Ω −

∂ ∑ ∑∫ ∫ (3.27)

d) Termo de Estabilização da Convecção:

( ) ( )( ),E

EN N dτ τΩ

⋅∇ ⋅∇ = ⋅∇ ⋅∇ Ω =∫a u a a u a

68

( )dim dim

1 1 1k n

E

m n nJI

IJ IJ IJ E J I

S k n E IJ k n

NNa a d

x xτ

= = = ⊃ Ω

∂∂ Ω − ∂ ∂

∑ ∑ ∑ ∑ ∫ u u (3.28)

e) Termo de Projeção de Alta Ordem:

( ) ( )dim

1 1

,k k k

E E

m nI

E IJ IJ IJ J E J I

S E IJ k k

NN Nd a N d

xτβ τβ τ β π π

= ⊃ =Ω Ω

∂ ⋅∇ = ⋅∇ Ω = Ω −

∂ ∑ ∑ ∑∫ ∫π a πa (3.29)

Considerando a Equação de Poisson para a Pressão:

a) Termo Laplaciano:

( )( ) ( )

( ) ( )

,E

E

E

E

I J E I J

E IJ

t p q t p Nd

t N N d p p

τ τ

τ

Ω

⊃ Ω

∆ + ∇ ∇ = ∆ + ∇ ⋅∇ Ω =

∆ + ∇ ⋅∇ Ω − +

∑ ∫

∑ ∫ (3.30)

b) Termo de Estabilização de Alta Ordem:

( ) ( )dim

1

,k k

E E

nI

E J E J I

E E IJ k k

Nq Nd N d

xτ τ τ ξ ξ

⊃ =Ω Ω

∂ ∇ = ⋅∇ Ω = Ω −

∂ ∑ ∑ ∑∫ ∫ξ ξ (3.31)

c) Termo Divergente

( ) ( )dim

1

,k k

E E

nI

E J E J I

E E IJ k k

Nq Nd N d u u

x⊃ =Ω Ω

∂ ∇ ⋅ = ∇ ⋅ Ω = Ω −

∂ ∑ ∑ ∑∫ ∫u u (3.32)

Para as equações da Projeção do Gradiente de Pressão e da Projeção do Termo Convecti-

vo no espaço de Elementos Finitos utilizou-se a matriz de massa consistente de Elementos Fi-

nitos e os termos do lado direito são exatamente idênticos ao termo do gradiente de pressão e

do termo convectivo, respectivamente, já discutidos anteriormente.

Nos termos definidos acima, as contribuições de cada elemento são avaliadas consideran-

do-se as arestas dos elementos, e desta forma é possível obter o sistema algébrico final acu-

69

mulando-se as contribuições de cada aresta, portanto as contribuições de todas as variáveis

que pertencem a estes termos devem também ser obtidas em função das arestas dos elemen-

tos. Então, de acordo com a formulação utilizada neste trabalho, a velocidade convectiva

1, −+= in

h

θua deve ser aproximada de forma a manter consistência com a formulação por aresta.

A primeira idéia é utilizar o mesmo raciocínio que foi utilizado para o termo viscoso υ , onde

a viscosidade cinemática foi obtida nas arestas simplesmente como uma média das viscosida-

des dos nós associados à aresta, da seguinte maneira:

2JI

IJ

υυυ

+= (3.33)

isto é, simplesmente utilizar uma média dos valores nodais da aresta e garantir conservação

obtendo a diagonal da matriz do termo convectivo como a negativa do somatório dos termos

de fora da diagonal, na mesma linha. Neste passo, os valores nodais já foram obtidos como

média dos valores nos elementos que contém o respectivo nó. Entretanto, o uso de procedi-

mento semelhante no termo convectivo destrói a principal característica da aproximação do

termo convectivo, que é manter-se de segunda ordem (Galerkin), pois representa uma apro-

ximação por diferenças centradas, quando elementos lineares são utilizados. Esta característi-

ca é refletida em nível discreto pelo fato do termo da diagonal principal da matriz ser zero

para pontos do domínio, quando utiliza-se o Método dos Elementos Finitos padrão, com ele-

mentos lineares ou Diferenças Finitas Centradas. Esse fato é demonstrado a seguir.

Seja o termo convectivo, escrito considerando as contribuições por arestas dos elementos

tetraédricos, como dado abaixo:

( ) ( )∫Ω

=Ω∇⋅=∇⋅ NdN uaua , ( )∑ ∑ ∑ ∫= ⊃ = Ω

−

Ω

∂

∂m

S

IJ

IJE

ndof

k

E

k

JIIJk

E

dx

NNa

1 1, ˆˆ uu (3.34)

onde

→N função de forma de elementos finitos linear;

→m número de arestas da malha;

→ndof número de graus de liberdade em cada nó.

70

E considere uma malha unidimensional representada na Fig.(3).

Figura 3. Representação de malha unidimensional.

onde

22N → Função de forma de elementos finitos referente ao nó 2 e elemento 2;

23N → Função de forma de elementos finitos referente ao nó 3 e elemento 2;

33N → Função de forma de elementos finitos referente ao nó 3 e elemento 3;

34N → Função de forma de elementos finitos referente ao nó 4 e elemento 3.

Montando o termo da diagonal principal da equação do nó 3, considerando a velocidade

convectiva avaliada na aresta IJ como a média dos valores nodais, isto é,

223

32

aaa

+= (3.35)

e

243

34

aaa

+= (3.36)

resulta na matriz dada abaixo:

( ) ( )

Ω∇⋅−Ω∇⋅− ∫∫ΩΩ ...

ˆ

...

.........

......

.........

33434322323

32

uaaEE

EE dNNdNN (3.37)

71

Note que na matriz acima, o termo da diagonal principal não resulta zero devido ao fato

das velocidades 32a e 34a serem diferentes. Desta forma o uso da média aritmética para o cál-

culo da velocidade convectiva na aresta leva a uma aproximação que não é de segunda ordem,

justificando assim, a avaliação da velocidade convectiva de forma alternativa.

Sendo assim, a única forma de garantir segunda ordem na aproximação deste termo, é a-

nular o somatório na diagonal da matriz da Eq. (3.37), e isto pode ser realizado fazendo IJa

ser função apenas do ponto nodal I , ou seja IIJ aa = [Soto et al., 2004], implicando segunda

ordem e garantindo conservação em nível discreto. Neste trabalho, a velocidade convectiva

foi obtida da seguinte maneira:

Pontos de Domínio

∑≠

=ptsn

IJ

IJ

kJIJ

Ik

m

uma , J∀ (3.38)

Pontos de Contorno

∑=

=ptsn

J

IJ

kJIJ

Ik

m

uma

1

, J∀ (3.39)

onde

→Ika , componente da velocidade convectiva do nó I, na direção k;

→IJm termo na posição IJ da matriz de massa consistente de elementos finitos;

→Jku , componente da velocidade do nó I na direção k;

→ptsn número de pontos da malha computacional.

Obs.: O somatório no denominador das equações acima apenas deve ser realizado sobre os

pontos nodais J que têm conectividade não nula com o ponto nodal I, pois 0=IJm se I não

está conectado a J.

72

Computando a velocidade convectiva como dado acima, resulta

( ) ( )

Ω∇⋅−Ω∇⋅− ∫∫ΩΩ ...

ˆ

...

.........

......

.........

334332233

32

uaaEE

EE dNNdNN (3.40)

E como, estamos utilizando elementos lineares, para elementos de domínio resulta

( ) ( )∫∫ΩΩ

Ω∇⋅−=Ω∇⋅32

34332233

EE

EE dNNdNN aa (3.41)

Garantindo assim segunda ordem na aproximação do termo convectivo.

O mesmo procedimento é efetuado para arestas de domínio e de contorno.

73

Capítulo 4

Computação Paralela

4.1 Conceitos Fundamentais

O interesse por problemas cada vez mais complexos tem levado à necessidade de compu-

tadores cada vez mais potentes, em termos de memória e velocidade de processamento, para

resolvê-los. Entretanto limitações físicas têm restringido o aumento da velocidade dos compu-

tadores seqüenciais, ou seja, computadores que executam instruções em série. Por outro lado,

problemas computacionais usualmente podem ter algumas de suas partes divididas, de tal

forma que essas partes poderiam ser solucionadas ao mesmo tempo, ou processadas em para-

lelo. Processamento paralelo é então uma forma pela qual a demanda computacional é suprida

através do uso simultâneo de recursos computacionais como processadores para a solução do

problema. Desta forma, o uso de computadores paralelos oferece a possibilidade de resolver

problemas de tamanhos praticamente impossíveis de serem tratados por máquinas com apenas

uma CPU. Para máquinas de memória distribuída, o tamanho dos problemas que podem ser

tratados pelo computador (cluster de PC’s por exemplo) pode ser aumentado, uma vez que

cada processador passe a tratar uma parte do problema. O tempo de tratamento de um deter-

minado problema de grande porte, pode também, se os devidos cuidados forem efetuados, ser

consideravelmente reduzido.

74

Os problemas envolvendo escoamentos de fluidos possuem natureza extremamente com-

plexa e, portanto, requerem um enorme esforço computacional para realizar simulações de

larga escala. As características tridimensionais, não lineares e transientes das equações gover-

nantes dos vários campos envolvidos, implicam em um alto custo computacional relacionado

ao armazenamento de informações e tempo de processamento a cada instante de tempo.

Neste contexto, o uso da computação paralela é uma importante ferramenta a ser utilizada

na simulação de fenômenos cada vez mais complexos e de forma completa, levando em con-

sideração todos os fenômenos envolvidos.

4.2 Aspectos Computacionais

Computação paralela é caracterizada pelo uso de várias unidades de processamento ou

processadores para executar uma computação de forma mais rápida. A computação paralela

baseia-se no fato de que o processo de resolução de um problema pode ser dividido em tarefas

menores, que podem ser realizadas simultaneamente através de algum tipo de coordenação.

Uma maneira clássica de se classificar arquiteturas de computadores é pela forma como o

fluxo de instruções e os dados se apresentam. Essa classificação é conhecida como taxonomia

de Flynn (Pacheco, 1997). Em 1966, Michael Flynn classificou os sistemas computacionais de

acordo com o número de instruções efetuadas e o número de dados analisados como: SISD

(Single Instruction-Single Data), SIMD (Single Instruction Multiple Data), MISD (Multiple

Instruction-Single Data) e MIMD (Multiple Instruction-Multiple Data).

Características importantes de programas paralelos incluem:

- Eficiência: tempo de execução de uma tarefa em um único processador dividido pelo

tempo de execução da mesma tarefa utilizando-se múltiplos processadores;

- Sobrecarga de Paralelismo: trabalho extra, associado à versão paralela de uma solução

comparado a sua versão seqüencial, envolvendo o tempo de CPU e a quantidade de memória

para a sincronização, comunicação e criação do ambiente paralelo.

- Sincronização: coordenação de tarefas simultâneas;

- Tarefa: pedaço de trabalho computacional discreto e independente.

Uma alternativa aos Supercomputadores vetoriais, muito utilizados até a alguns anos atrás

no desenvolvimento de ferramentas na computação paralela são os “clusters” de estações de

Trabalho ou de microcomputadores (PC’s), que atualmente vem recebendo grande aceitação

75

devido ao custo relativamente baixo de instalação e manutenção, à facilidade de atualizações,

à utilização de programas com códigos abertos e à independência de fornecedores.

4.2.1 Particionamento de Malha

Uma das etapas importantes no desenvolvimento de ferramentas que utilizem computação

de alto desempenho é o particionamento de malhas. Para isso foi utilizado o PARMETIS (Pa-

rallel Graph Partitioning and Sparse Matrix Orderinging) [Karypis et al., 2003], que é uma bi-

blioteca paralela baseada em MPI (Message Passing Interface) [Pacheco, P. S., 1995 e 1997]

que possui uma grande variedade de algoritmos para particionamento e reparticionamento de

malhas não estruturadas, de forma a reduzir a esparsidade de matrizes. Em se tratando de pro-

blemas de grande escala, ou seja, envolvendo malhas computacionais com grande quantidade

de elementos, a comunicação entre os processos envolvidos pode ser drasticamente reduzida

se a malha for decomposta de tal forma que o número de interfaces (elementos nas interfaces)

seja minimizado. Neste trabalho, o particionamento foi realizado por arestas, ou seja, as ares-

tas estão apenas em um processo, e as interfaces são dadas pelos nós. Neste trabalho foi de-

senvolvido um programa para gerar a estrutura de dados necessária para a utilização do

PARMETIS, que necessita da matriz de adjacências dos elementos que serão particionados

[Woll, 2004], veja apêndice A. Este programa paralelo lê uma malha serial, pois ainda não

temos disponível um gerador de malhas paralelo, sendo que cada processador lê apenas uma

parte da malha. Inicialmente o particionamento é feito simplesmente dividindo o número de

arestas pelo número de processadores, e cada processador lê uma parte desta divisão. Note

que desta forma é realizado um balanço de carga, mas não há controle sobre as interfaces en-

tre os processos envolvidos. A partir deste particionamento inicial, é feito uso extensivo do

MPI para gerar a matriz de adjacências em paralelo, pois várias trocas de informações são ne-

cessárias, e quando esta tarefa é realizada, o PARMETIS [Karypis et al., 2003] é chamado pa-

ra realizar o re-particionamento, agora com balanço de carga e minimização das interfaces.

A Fig.(4) mostra exemplos de um particionamento inicial (caracterizado por cores diferen-

tas para cada processo), a partir da divisão do número de elementos pelo número de processa-

dores alocados, evidenciando a grande quantidade de elementos com interfaces com outros

processadores. Este tipo de particionamento, apesar de ser considerado um particionamento

válido, pois possibilita a obtenção da solução via computação paralela, é extremamente inefi-

ciente, visto que o número de interfaces é elevado. E desta forma, como pelo menos a cada

76

passo de tempo informações devem ser trocadas entre processadores que compartilham os

mesmo nós, este custo computacional torna proibitivo este tipo de particionamento.

Figura 4. Exemplo do particionamento inicial e da comunicação entre processadores que

compartilham nós da malha.

A Fig.(5) mostra particionamentos obtidos via aplicação do PARMETIS, que a partir de

uma matriz de adjacências de um grafo de qualquer orientação, realiza o reparticionamento de

forma a balancear a carga para cada processador envolvido e ao mesmo tempo minimizando

as interfaces entre processadores, que está descrito em mais detalhes no Apêndice A. A matriz

de adjacências tem cada linha associada a um nó da malha, desta forma as informações sobre

o nó 1 estão na linha 1 da matriz e assim sucessivamente, sendo que os valores 0 e 1 represen-

tam a conectividade do nó da linha i com o nó da coluna j, ou seja, se o nó i tem conectividade

com a coluna j, então a posição ija da matriz tem entrada 1, caso contrário tem entrada 0, e

por definição se i=j, então a conectividade é 0. Quando tem-se um grafo não orientado, a ma-

triz de adjacências é simétrica. Obviamente neste caso o custo computacional com a comuni-

cação é muito inferior aos exemplos da Fig.(4).

Particionamento inicial da malha.

Compartilhamento de elementos.

Nó I Nó I

Nó K

Nó K

Nó J Nó J

77

Figura 5. Exemplo de particionamento via PARMETIS.

Alguns exemplos de particionamento de malhas são mostrados nas Figs. (6-12), a seguir,

considerando particionamentos realizados simplesmente dividindo o número de elementos pe-

la quantidade de processadores seguindo a seqüência de dados, e particionamentos através do

PARMETIS [Karypis et al., 2003]. É importante ressaltar que foi desenvolvido um programa

que realiza o primeiro particionamento, e a partir deste, monta a matriz de adjacências para

ser utilizada no PARMETIS [Karypis et al., 2003], sendo que todas as etapas são realizadas

em paralelo, ou seja, em nenhum momento a malha computacional é carregada totalmente em

apenas um nó.

Figura 6. Malha sobre o perfil NACA0012 com particionamento aleatório.

Grafo.

0 1 1 1 0 0

1 0 0 1 1 1

1 0 0 1 1 0

1 1 1 0 1 0

0 1 1 1 0 1

0 1 0 0 1 0

2

Matriz de adjacências.

1

3 4

6

5

78

Figura 7. Malha sobre o perfil NACA0012 reparticionada via PARMETIS.

Figura 8. Detalhe da malha sobre o perfil NACA0012 reparticionada via PARMETIS.

79

Figura 9. Malha sobre a geometria da artéria carótida 2D com particionamento aleatório.

Figura 10. Malha sobre a geometria da artéria carótida 2D reparticionada via PARMETIS.

80

Figura 11. Paticionamento aleatório sobre malha 3D.

Figura 12. Repaticionamento de malha 3D via PARMETIS.

As figuras Fig.(11) e Fig.(12), mostram uma malha tridimensional com um particiona-

mento aleatório, e a mesma malha reparticionada utilizando o PARMETIS, respectivamente,

81

onde fica evidente que as interfaces entre os processos foi reduzida, gerando grande economia

no tempo computacional que é desperdiçado nas comunicações entre os processos.

4.2.2 Gerenciamento da Comunicação Paralela via MPI

Na utilização de computadores paralelos com memória distribuída, é necessário que haja

troca de informações entre os processadores durante o processamento dos dados envolvidos

na simulação. Neste trabalho, toda a comunicação é gerenciada pelo MPI (Message Passing

Interface) [Pacheco, P. S., 1995 e 1997], que consiste em um conjunto de pacotes de funções

úteis e necessárias durante o processamento paralelo de informações em programas desenvol-

vidos utilizando linguagens de programação, tais como Fortran, C e C++. Uma grande vanta-

gem da utilização do MPI é que vários pacotes de ferramentas paralelas utilizam o MPI como

padrão de comunicação, tais como o PETSc [Balay et al., 2005] e o PARMETIS [Karypis et

al., 2003]. Neste trabalho o MPI foi utilizado de forma extensiva no particionamento da malha

e na etapa de mapeamento dos dados no simulador.

A Fig.(13) mostra um exemplo de uma malha que foi particionada em quatro subdomí-

nios, e cada subdomínio foi atribuído a um processador. Mesmo que as interfaces tenham sido

minimizadas, ainda existe a necessidade de troca de informações entre processadores. E neste

1 2

3 4

Figura 13. Exemplo de particionamento de malha com 4 processadores.

82

caso, como mostra a Fig.(14), existem interfaces entre os processadores 1 e 2 (I12), 1 e 3

(I13), 2 e 4 (I24), 3 e 4 (I34) 3 e 1, 2, 3 e 4 (I1234). Desta forma, algum tratamento especial

deve ser dado aos nós que estão nestas interfaces. Neste trabalho, o nó é considerado local em

apenas um processador e remoto em todos os outros que o contém, e sendo assim, ao final de

cada iteração temporal, todos os processadores com nós remotos recuperam as informações

atualizadas destes nós junto ao processador proprietário daquele nó.

4.2.3 Utilização do Toolkit PETSc

Também foi utilizado o PETSc (Portable Extensible Tollkit for Scientific Computation)

[Balay et al., 2005 ], que consiste em um conjunto de ferramentas que auxiliam em diversas

etapas do processamento de informações. O PETSc é essencialmente útil para a obtenção de

soluções numéricas de equações diferenciais parciais quando utiliza-se computação de alto

desempenho, pois consiste de rotinas que auxiliam na definição de estruturas de dados e estru-

turas matriciais em ambientes que envolvem computação em grande escala. Uma outra vanta-

gem é que o PETSc utiliza o MPI para viabilizar suas estruturas paralelas. O PETSc possui

aplicações desenvolvidas em Fortran, C e C++, tais como rotinas de álgebra matricial, neces-

sárias na simulação de escoamentos fluidos.

Processor 1

Processor 3

Processor 4

Processor 2

I: 1-2-3-4 I:1-3 I: 2-4

I: 1-2

I: 3-4

Figura 14. Esquematização da comunicação entre processadores.

83

A Fig.(15) mostra um vetor típico definido em ambiente paralelo utilizando-se o PETSc.

Na definição de um vetor paralelo, o usuário determina em quais processadores este vetor será

alocado e quantas posições do vetor estarão definidas localmente em cada processador. Desta

forma, como mostra a Fig.(15), o vetor U, que possui n posições, estará contido nos processa-

dores P1, P2, P3 e P4, sendo que cada processador armazena uma parte deste vetor. Sendo as-

sim, quando um processador necessita de uma informação em uma posição que não está sob

seu domínio local, é necessário que haja troca de informações entre os processadores.

Figura 15. Exemplo de variável vetorial paralela.

Para evitar que o custo destas trocas de informações seja elevado, após o particionamento

de malha, é importante que seja realizada também uma renumeração dos nós em cada proce-

sasdor de forma a melhorar o desempenho na inserção e obtenção de dados das matrizes e ve-

tores paralelos, visto que é muito interessante que as partições possuam o máximo possível de

nós locais àquele processo, minimizando a necessidade de comunicação. Esta renumeração é

realizada através de um mapeamento dos nós em cada processador, a partir do nó zero, no

processador zero, até o nó (n-1) no último processador alocado. Isto fica evidenciado na

Fig.(16) abaixo.

U0 U1 U2 … Un-2 Un-1 Un

P1

P2

P3

P4

U = ⊂

84

Figura 16. Renumeração dos nós após o particionamento da malha.

A Fig.(16) mostra uma malha particionada para três processos, sendo que o processo 0

mantém locais os nós 1, 4, 7, 13, 14, 18 e 20. O processo 1 mantém locais os nós 5, 6, 9, 11,

12, 15 e 17. O processo 2 mantém locais os nós 2, 3, 8, 10, 16, 19 e 21. Obviamente, nos

elementos que estão nas interfaces das partições, há a necessidade de que todos os processos

destas interfaces mantenham os nós que estão compartilhando, sendo que em um processo ele

permanece como nó local e nos outros, como nó remoto. Na definição de vetor paralelo no

PETSc, a questão é que este vetor também é particionado e cada processo mantém uma parte

desse vetor, como mostra a parte central da Fig.(16). O problema acontece quando, por

exemplo no lado esquerdo da Fig.(16), para o processo 0, os nós 13, 14, 18 e 20, que

correspondem às posições 13, 14, 18 e 20 do vetor paralelo, estão locais para o processo e

remotos para o vetor, pois para armazenar, ou recuperar informações nestas posições é preciso

comunicação entre processos via MPI, ocasionando desperdício de tempo computacional. Isto

pode ser minimizado realizando uma renumeração dos nós em cada processo, seguindo uma

seqüência, a partir do processo 0. Por exemplo, veja o lado direito da Fig.(16), que apresenta a

mesma malha, com o mesmo particionamento, porém com uma numeração dos nós que

85

minimiza a troca de informações entre os processos, pois a maioria dos nós que estão locais

ao processo estão também em posições locais do vetor.

4.3 Programa NS_INCOMP

Neste trabalho foi desenvolvido um programa computacional para tratar escoamentos in-

compressíveis governados pelas equações de Navier-Stokes em domínios tridimensionais e u-

tilizando-se computadores paralelos com memória distribuída. A computação paralela se faz

necessária neste tipo de análise devido ao fato de que escoamentos fluidos em domínios tri-

dimensionais necessitam grande refinamento da malha discreta que representa o domínio con-

tínuo, caracterizando-se como um problema de grande porte, e desta forma proibitivo para

análises em computadores seqüenciais.

4.3.1 Aspectos Gerais

Este programa foi denominado NS-INCOMP3D e consiste em um conjunto estruturado de

rotinas escritas na linguagem Fortran 90 [Metcalf e Reid, 1996], aproveitando as vantagens

desta versão do Fortran, por exemplo, quanto ao uso de tipos estruturados na definição das va-

riáveis necessárias. A estruturação do programa permite a inclusão de forma facilitada de no-

vas variáveis dinamicamente alocáveis e até mesmo de novas rotinas, pois duas estruturas

denominadas “MeshDat” e “PhyDat” carregam todas as informações relevantes para a simu-

lação. A estrutura “MeshDat” contém as informações numéricas e a estrutura “PhyDat” con-

tém as informações físicas da análise. O programa foi desenvolvido para trabalhar com n

processadores, sendo n=1,2,3,...np.

A seguir apresentamos um fluxograma mostrando a associação do conjunto de rotinas de-

senvolvidas e a seguir, uma breve definição de cada rotina.

86

Figura 17. Fluxograma do programa NS_INCOMP3D.

87

A Fig.(17) mostra a seqüência lógica de funcionamento do programa NS_INCOMP3D,

onde estão apresentadas as diversas rotinas que foram desenvolvidas. A seguir descrevemos

as principais operações de cada rotina.

- Módulo malha: neste módulo são definidas todas as estruturas vetoriais e matriciais neces-

sárias ao longo do programa, assim como todas as variáveis relativas aos procedimentos nu-

méricos realizados, ou seja, todas as informações que dizem respeito a algum parâmetro

relativo à malha computacional, por exemplo, número de nós, número de arestas, vetor que

armazena as informações nodais de velocidade e pressão, vetor RHS (Right Hand Side), ma-

triz LHS (Left Hand Side) e outras estruturas.

- Módulo físico: neste módulo são definidas todas as variáveis que contém alguma informa-

ção sobre a física do escoamento analisado, como viscosidade, massa específica, tempo inici-

al, etc.

- Rotina openpre: nesta rotina são abertos os arquivos contendo as informações físicas e nu-

méricas do escoamento fluido analisado, e da malha computacional utilizada. Este procedi-

mento é realizado em paralelo, ou seja, cada processo abre os arquivos e armazena apenas as

informações que dizem respeito ao seu processo. Vale ressaltar que a malha já está particio-

nada e a estrutura de elementos não é necessária, pois o particionamento foi realizado por a-

restas. Os arquivos necessários são:

- Nome_do_problema.coe: arquivo de dados da malha (nós, elementos, arestas e coefi-

cientes das arestas)

- Nome_do_problema.inp: arquivo de dados físicos tais como, massa específica, vis-

cosidade dinâmica, velocidade de referência, condições de contorno e iniciais, tempo inicial,

tempo final, e parâmetros numéricos tais como tolerância para o Gauss-Seidel, tolerância para

o regime permanente, número máximo de passos de tempo e outros.

- Nome_do_problema.dat: arquivo que contém as informações sobre as arestas que são

locais a cada processador. Neste caso, a cada processador é atribuído um arquivo com apenas

as suas informações locais.

- Rotina mapping: esta rotina realiza um mapeamento das informações da malha computa-

cional de forma a minimizar a comunicação entre os processos envolvidos na simulação. Co-

mo os vetores e matrizes definidos em ambientes paralelos armazenam uma parte destas

88

estruturas em cada processo, é vantajoso que estas partes correspondam a exatamente os nós

deste processo. Sendo assim, os nós são renumerados em uma seqüência a partir de zero até n-

1, onde n é o número de nós da malha e do processo zero até k-1, onde k é o número de pro-

cessos envolvidos na simulação. Esta rotina está detalhada no Apêndice B.

- Rotina bcondit: esta rotina define as condições de contorno e iniciais do problema, e tam-

bém, a partir destas informações determina o número de graus de liberdade para as variáveis

de pressão e velocidade, necessárias para definir as dimensões do vetor RHS e da matriz LHS.

- Rotina restart: esta rotina é acionada quando trata-se de uma reinicialização de alguma aná-

lise, ou “restart” de algum problema, e desta forma, reinicializa todas as informações do esco-

amento a partir do último instante de tempo salvo em arquivos de reinicialização.

- Rotina precalculus: nesta rotina são calculadas algumas informações relevantes a cada pas-

so de tempo, tais como os parâmetros de estabilização para as equações de Conservação de

Quantidade de Movimento e Pressão, assim como também a velocidade convectiva que preci-

sa ser determinada por nó. Estas informações são necessárias a partir do primeiro passo de

tempo, por isso está nessa localização do programa, e novamente aparece após a rotina que

determina as variáveis do escoamento após cada iteração do processo tipo Gauss-Seidel.

- Rotina timestep: esta rotina obtém o passo de tempo global e os parâmetros de tempo locais

necessários para controlar a quantidade de estabilização necessária para as equações de pres-

são e quantidade de movimento. Estes parâmetros são importantes porque o passo de tempo

para fenômenos envolvendo escoamentos incompressíveis são geralmente muito pequenos,

diminuindo a capacidade de estabilização empregada.

- Rotina uvar: esta rotina resolve a equação da Quantidade de Movimento para a componente

da velocidade na direção x .

- Rotina vvar: esta rotina resolve a equação da Quantidade de Movimento para a componente

da velocidade na direção y .

- Rotina wvar: esta rotina resolve a equação da Quantidade de Movimento para a componen-

te da velocidade na direção z .

89

- Rotina pvar: esta rotina resolve a equação de Poisson para a pressão.

- Rotina psiuvar: esta rotina resolve a equação da projeção da convecção no espaço de Ele-

mentos Finitos na direção x .

- Rotina psivvar: esta rotina resolve a equação da projeção da convecção no espaço de Ele-

mentos Finitos na direção y .

- Rotina psiwvar: esta rotina resolve a equação da projeção da convecção no espaço de Ele-

mentos Finitos na direção z .

- Rotina csiuvar: esta rotina resolve a equação da projeção do gradiente de pressão no espaço

de Elementos Finitos na direção x .

- Rotina csivvar: esta rotina resolve a equação da projeção do gradiente de pressão no espaço

de Elementos Finitos na direção y .

- Rotina csiwvar: esta rotina resolve a equação da projeção do gradiente de pressão no espaço

de Elementos Finitos na direção z .

- Rotina converge: esta rotina determina as normas dos resíduos nas iterações Gauss-Seidel.

A convergência é determinada pelas variáveis de pressão e velocidade. Se o resíduo for menor

do que a tolerância definida, então avança para o próximo passo de tempo, senão, continua i-

terando.

- Rotina updating: esta rotina realiza a atualização das variáveis envolvidas, ou seja, veloci-

dade, pressão, projeção da convecção e projeção do gradiente de pressão. Note que esta atua-

lização somente é realizada após a convergência em cada passo de tempo.

- Rotina outputing: esta rotina cria os arquivos de saída de dados no formato VTK, para se-

rem visualizados no VisIt [https://wci.llnl.gov/codes/VisIt/] e/ou ParaView

[http://www.paraview.org/]. Ressaltamos aqui que como a estrutura de elementos não é man-

tida durante a simulação, foi desenvolvida uma ferramenta específica para que, a partir destes

90

arquivos (cada processo cria seu próprio arquivo) sejam criados os arquivos de entrada do Vi-

sIt e/ou ParaView que necessitam da estrutura de elementos para as interpolações necessárias

na visualização dos resultados.

As rotinas apresentadas acima constituem os principais blocos do programa desenvolvido.

Outras rotinas estão em desenvolvimento e serão incorporadas ao programa em momento o-

portuno. Dentre as rotinas em desenvolvimento, destacamos as rotinas de movimentação de

malha e cálculo das variáveis estruturais necessárias no tratamento de fenômenos de interação

fluido-estrutura. Vale ainda ressaltar que todas as rotinas são paralelas, ou seja, todos os pro-

cessos acessam as rotinas, porém cada processo realiza as operações necessárias sobre a sua

parte do domínio, obviamente havendo troca de informações sempre que necessário.

Ressaltamos ainda que os sistemas de equações originados para as variáveis de velocida-

de, de caráter não simétrico, e para as variáveis de projeção da convecção e do gradiente de

pressão no espaço de elementos finitos são resolvidos com a aplicação do método iterativo

GMRES (Generalized Minimum Residual) [Saad, 2003], que tem como objetivo minimizar a

norma do resíduo do sistema. Esta operação é realizada com a construção de uma base orto-

normal no subespaço de Krylov [Saad, 2003], porém como o algoritmo do GMRES padrão

pode ser impraticável computacionalmente, pois a cada iteração os produtos matriz por vetor

precisam ser armazenados, o que pode ocasionar problemas de falta de memória, quando é

necessário um grande número de iterações para resolver o sistema linear, utiliza-se a alterna-

tiva de limitar a dimensão do subespaço de Krylov. Essa alternativa recebe o nome de G-

MRES(m) e possui a desvantagem de não ter uma robustez como o GMRES tradicional, pois

no GMRES(m) a convergência não é garantida ou ainda, pode ser muito lenta. Diante disso,

neste trabalho utilizamos o subespaço de Krylov definido pelo PETSc, que define que a má-

xima dimensão do subespaço, não ultrapassa o parâmetro de “restart”, que tem valor 30 como

padrão, pois não possuímos experiência na limitação explícita do subespaço. Também em to-

dos os casos estudados, os parâmetros de convergência do GMRES foram utilizados como o

padrão do PETSc, que é baseada da norma 2l do resíduo, sendo que a convergência é baseada

em três quantidades: o decréscimo da norma do resíduo relativa a norma do RHS, 510rtol−= ,

a dimensão absoluta da norma do resíduo, 5010atol−= e o acréscimo relativo no resíduo,

510dtol = . Como pré-condicionamento destes sistemas, utilizamos um bloco Jacobi, sendo

um bloco por processo, e em cada bloco um pré-condicionamento do tipo ILU (Incomplete

Lower Upper Fatoration) [Saad, 2003]. Para o sistema de equações originado pela equação da

91

pressão, aproveitamos o fato de este sistema possuir uma matriz simétrica, positiva definida,

ou seja, possuindo todos os autovalores e pivôs positivos, sendo que neste caso uma boa esco-

lha é utilizar o Método dos Gradientes Conjugados [Saad, 2003] com o mesmo pré-

condicionador utilizado para a velocidade. O padrão de pré-condicionamento no PETSc [Ba-

lay et al., 2005] é pela esquerda, de forma a preservar o resíduo, sendo que para efetuá-lo pela

direita, é preciso especificar explicitamente. Neste trabalho, o pré-condicionamento padrão foi

utilizado. Em trabalhos futuros pretendemos aprofundar estudos nesta temática e investigar

novas alternativas e melhorias em termos de pré-condicionadores, e dos parâmetros de contro-

le mencionados, como por exemplo, a otimização de banda para efetuar o pré-processamento

ILU, via reordenação dos nós. Também ressaltamos que os sistemas gerados nas rotinas uvar,

vvar, wvar possuem características muito similares e poderiam ser resolvidos como um siste-

ma único, de forma monolítica, o mesmo acontecendo com as rotinas, csiu, csiv, csiw, e psiu,

psiv, psiw. Porém esta alternativa não foi adotada para mantermos um controle maior sobre

cada etapa do processamento, nesta fase de desenvolvimento do sistema computacional. No

futuro poderemos aproveitar deste fato para obter ganho de desempenho computacional.

92

Capítulo 5

Aplicações Numéricas

Nesta seção serão apresentados os resultados numéricos obtidos das simulações computa-

cionais realizadas ao longo deste trabalho. Primeiramente foi realizado um estudo consideran-

do-se problemas de difusão pura permanente e transiente e a seguir problemas de convecção-

difusão. Estes estudos preliminares serviram para testar a capacidade do programa desenvol-

vido em realizar simulações considerando múliplos processadores, verificando todas as roti-

nas implementadas, principalmente em relação aos mapeamentos e trocas de mensagens

necessárias ao longo da simulação, sendo que estes resultados são apresentados no Apêndice

D, e ressaltamos que foi considerado um fluido hipotético, onde as propriedades do escoa-

mento foram ajustadas para o número de Reynolds desejado. Posteriormente foram analisados

problemas de escoamentos fluidos incompressíveis, inicialmente analisando exemplos onde

existem regimes permanentes, desta forma buscou-se avaliar a capacidade da ferramenta de-

senvolvida em capturar o regime permanente, bem como os fenômenos físicos envolvidos e a

precisão dessa captura, inicialmente em 2D com implementação seqüencial [Antunes et al.,

2005, Antunes et al., 2006], e posteriormente em 3D com implementação paralela [Antunes,

et al., 2008]. Após vários testes envolvendo escoamentos que convergem para o regime per-

manente, analisamos o escoamento em torno de um cilindro, que configura um escoamento

que converge para um regime periódico para a faixa de número de Reynolds estudada, desta

93

forma buscando avaliar a ferramenta tanto na captura dos fenômenos físicos envolvidos quan-

to na precisão temporal, da formulação. Neste caso foi realizada uma análise em termos da

frequência de formação dos vórtices atrás do cilindro, e comparada com resultados encontra-

dos na literatura, demonstrando a capacidade da ferramenta em tratar estes fenômenos. Todos

os resultados numéricos foram obtidos utilizando-se os conceitos de computação paralela, em

malhas particionadas demonstrando a capacidade da ferramenta em tratar subdomínios de

forma adequada.

Ressaltamos que todas as malhas computacionais foram geradas no Gmsh

(www.geuz.org/gmsh/), que é um gerador de malhas do tipo GPL (General Public License),

desenvolvido por Cristophe Geuzaine da University of Liège e Jean-François Remacle da Ca-

tholic University of Louvain. Como o Gmsh é um gerador de malhas seqüencial, foi desen-

volvido um programa que a partir de uma malha seqüencial, realiza um particionamento

inicial e todos os mapeamentos necessários para montar uma matriz de adjacências das arestas

da malha e com o auxílio do PARMETIS [Karypis et al., 2003], otimiza-se este particiona-

mento em temos de balanço de carga e minimização das interfaces entre as partições. Este

particionador de malhas foi escrito em FORTRAN 90 [Metcalf e Reid, 1996], e em nenhum

momento mantém toda a malha em um único processador (ver Apêndices A e B). Ainda, to-

dos os resultados foram visualizados nos programas, VisIt [https://wci.llnl.gov/codes/VisIt/.]

que possui licença do tipo BSD (Berkeley Software Distribution), desenvolvido no LLNL

(Laurence Livermore National Laboratory) e Paraview [http://www.paraview.org/] que é um

visualizador com código fonte aberto, que foi desenvolvido em uma parceria da Kitware Inc.

e Los Alamos National Laboratory. Uma vantagem na utilização destes visualizadores é o fato

de ambos permitirem visualizações de dados paralelos e utilizarem a mesma formatação de

arquivo do tipo *.VTK.

5.1 Escoamento em Tubo

Neste exemplo nós estudamos o fluxo de um fluido viscoso e incompressível em um tubo

tridimenisional, sendo que a malha foi particionada para dois processos, como mostra a Fig.

(18), onde cada cor corresponde à parte do domínio que foi atribuída a cada processador. O

tubo possui 20m de comprimento e 2m de diâmetro, e a malha possui 67.706 elementos tetra-

édricos, 82.714 arestas e 12.164 nós.

94

Figura 18. Domínio de interesse e malha particionada.

As condições de contorno são:

- perfil de velocidade parabólico prescrito na entrada;

- pressão prescrita na saída: 0p = ;

- condição de não deslizamento nas paredes do canal:u = 0 .

Neste caso a simulação foi realizada com número de Reynolds 100, e maxU é a máxima ve-

locidade na entrada, ajustada para se obter o Reynolds desejado, e fazendo 1µ ρ= = , caracte-

rizando um fluido hipotéico.

As Figs.(19) e (20), a seguir, mostram o mapa de cores da pressão no regime permanente,

e o gráfico com solução numérica ao longo da direção axial mostrando a linearidade esperada,

já que o gradiente de pressão é constante nessa direção. Este resultado representa uma boa a-

proximação da solução exata [Antunes et al., 2005].

Figura 19. Mapa de cores do campo de pressão.

95

Figura 20. Corte longitudinalmente ao tubo e gráfico da pressão ao longo desta direção.

A Fig. (21) mostra o campo de velocidade, detalhando o perfil parabólico na saída do ca-

nal, e comparando a componente horizontal da velocidade ao longo do eixo y, onde fica evi-

dente o caráter parabólico do campo e a boa aproximação obtida.

Figura 21. Campo escalar de velocidade absoluta (a), (b) detalhe do perfil na saída.

Neste exemplo foi calculada a vazão mássica na entrada e na saída do canal, onde foi con-

firmada a conservação de massa, a menos de erros da ordem de 410− .

5.2 Backward Facing Step

Neste exemplo foi considerado um escoamento viscoso e incompressível em um canal

com um “backward facing step” [Armaly et al., 1983, Soto et al., 2004, Williams e Baker,

1997], como mostra a Fig. (22). A malha computacional utilizada neste exemplo possui

96

77.163 elementos tetraédricos, 96.506 arestas e 14.617 nós, e o número de partições utilizadas

foi 2, sendo que cada cor na Fig. (22) corresponde ao subdomínio que foi atribuído a cada

processador. A malha adotada é não uniforme conventrando elementos na região próxima ao

degrau.

Consideramos 1m a dimensão da entrada do canal, 0.94m a altura do degrau, e 31m o

comprimento total do canal.

Figura 22. Domínio de interesse.

As condições de contorno para este caso são:

- entrada: perfil de velocidade parabólico prescrito;

- saída: pressão prescrita 0=p e as “pseudo-tractions” prescritas como zero;

- condição de não deslizamento nas paredes do canal: u = 0 ;

- nas fronteiras frente e trás o tensor viscoso nas direções x e z foi prescrito como zero,

e a componente da velocidade na direção y foi prescrita como zero.

As simulações foram realizadas com números de Reynolds 100 e 1000, e maxU é a máxima

velocidade na entrada, ajustada para se obter o Reynolds desejado. Neste caso o número de

Reynolds é dado por [Soto et al., 2004, Williams e Baker, 1997]:

υ

hU 2)3/2(Re max= (5.1)

onde h é a dimensão característica do escoamento e υ é a viscosidade cinemática.

97

5.2.1 Número de Reynolds 100

A Fig. (23) mostra o campo de velocidade, detalhando o vórtice formado devido à presen-

ça do degrau.

(a) (b)

Figura 23. Campo de velocidade (a), e detalhe do campo na região de recirculação (b).

A Fig. (24) mostra o campo de pressão, onde podemos verificar na região de recirculação

uma queda de pressão, como esperado.

Figura 24. Mapa de cores do campo de pressão.

Neste exemplo comparamos com outros autores o comprimento do vórtice, sendo que nes-

te trabalho o vórtice teve comprimento de aproximadamente 2,6m, como mostra a Fig. (25).

Este valor está próximo dos valores encontrados na literatura (Armaly et al., 1983, Boivin et

al., 2000, Williams and Baker, 1997), que ficam em torno de 2,8m. Levando o fato de não ter

98

sido realizado um estudo de convergência de malha, os resultados obtidos são considerados

bons.

Figura 25. Detalhe dos campos de pressão e velocidade mostrando o comprimento do vórtice.

A Fig. (26) mostra a componente da velocidade na direção do escoamento, onde podemos

constatar a região do domínio onde esta componente da velocidade torna-se negativa. Esta é

uma maneira de auxiliar na estimativa do comprimento do vórtice.

Figura 26. Mapa de cores da componente da velocidade na direção do escoamento.

5.2.2 Número de Reynolds 1000

A Fig. (27) mostra o campo de velocidade, detalhando o segundo vórtice formado na parte

superior do domínio devido à presença do degrau e o alto número de Reynolds.

99

Figura 27. Campo de velocidade com regiões de recirculação.

A Fig. (28) mostra o campo de pressão, onde podemos verificar nas regiões de recircula-

ção uma queda de pressão, como esperado.

Figura 28. Mapa de cores do campo de pressão.

Neste exemplo tivemos como objetivo observar a captura dos fenômenos físicos esperados

no problema e que diz que a partir de Reynolds 1000, há a formação de um segundo vórtice

na parte superior do domínio.

100

5.3 Escoamento ao Redor de um Cilindro

Este exemplo envolve o fluxo de um fluido viscoso ao redor de um cilindro, que repre-

senta um exemplo muito utilizado [Costa, 2004, De Sampaio et al., 1993, De Sampaio, 1991,

Gresho et al., 1984, Li et al., 1991, Soto et al., 2001] para testar a capacidade da formulação

em capturar os vórtices periódicos que ocorrem atrás do cilindro em determinada faixa de

número de Reynolds. Neste caso o número de Reynolds é baseado no diâmetro do cilindro e

no perfil de velocidade prescrito na entrada. O domínio é composto por um retângulo de di-

mensões [ ] [ ]8,160,0 × , com o cilindro centrado no ponto [ ]4,4 , possuindo diâmetro 2. Neste

exemplo dois casos de interesse foram analisados. Primeiramente considerando número de

Reynolds 40, para verificar a existência de dois vórtices simétricos e estáticos posicionados

atrás do cilindro, e o caso com número de Reynolds 100 para verificar a ocorrência dos famo-

sos vórtices de Von Karmam, que ocorrem devido à instabilidade dos vórtices. Esse caso é

largamente utilizado para testar a acurácia no tempo do esquema numérico utilizado. O perío-

do das oscilações para 100Re = fica em torno de aproximadamente 6.0, como verificado na

literatura [Gresho et al., 1984, Li et al., 1991, Soto et al., 2001]. O número de Reynolds é a-

justado pela viscosidade dinâmica µ , sendo 1ρ =

As condições de contorno para este caso são:

- entrada: perfil de velocidade constante prescrito: 1u = , 0v w= = ;

- saída: pressão prescrita 0p = , e o tensor das tensões viscosas também fixado como

zero;

- condição de não deslizamento no cilindro:u = 0 ;

- nas fronteiras frente e trás o tensor viscoso nas direções x e z foi prescrito como zero,

e a componente da velocidade na direção y foi prescrita como zero;

- nas fronteiras de cima e embaixo o tensor viscoso nas direções x e y foi prescrito co-

mo zero, e a componente da velocidade na direção z foi prescrita como zero.

A Fig. (29) mostra a geometria e a malha utilizadas neste exemplo, que consiste em

69.828 elementos tetraédricos, 91.356 arestas e 14.592 nós, e cada cor correpondendo ao sub-

domínio que foi atribuído a cada procesador.

101

Figura 29. Malha computacional sobre cilindro tridimensional.

5.3.1 Número de Reynolds 40

Neste caso, o objetivo é obter a formação dos vórtices simétricos e estáticos atrás do cilin-

dro [Lv et al., 2006], o que foi conseguido com êxito, como mostra a Fig. (30) abaixo.

Figura 30. Campo de velocidade mostrando as linhas de corrente.

Neste exemplo podemos comparar também o tamanho dos vórtices atrás do cilindro, que

na literatura [Lv et al., 2006] encontram-se com uma taxa de aspecto de 2.35 diâmetros, e nes-

te trabalho encontramos 2,1 diâmetros, como mostra a Fig. (31), o que representa boa concor-

dância e indica qualidade dos resultados numéricos obtidos.

102

Figura 31. Linhas de Corrente mostrando os vórtices simétricos e estáticos.

A Fig. (32) mostra o campo de pressão evidenciando a queda de pressão na região de for-

mação dos vórtices.

Figura 31. Campo de pressão.

5.3.2 Número de Reynolds 100

Neste caso, o objetivo é obter os vórtices de von Karmam, ou esteira de vórtices, como é

conhecida a formação periódica de vórtices atrás do cilindro. A precisão temporal aqui é tes-

tada comparando-se o período da formação dos vórtices através da análise da componente da

velocidade ortogonal ao escoamento livre, em um ponto localizado atrás do cilindro. Uma a-

nálise da freqüência das oscilações pode ser obtida via FFT (Fast Fourier Transform), aplica-

da ao histórico da componente ortogonal da velocidade naquele ponto.

103

A Fig. (33) abaixo mostra o campo de velocidade e as linhas de corrente na esteira de vór-

tices.

Figura 33. Campo de velocidade e linhas de corrente.

A Fig. (34) abaixo mostra o campo de pressão com as isosuperfícies de pressão no instan-

te coincidente com a Fig. (33).

Figura 34. Campo de pressão e isosuperfícies de pressão.

A Fig. (35) mostra a queda de pressão na região de formação do vórtice, o que fica claro

observando-se o campo de velocidade.

104

Figura 36. Campos de velocidade e pressão.

A Fig. (37-a) abaixo mostra o histórico dos valores assumidos pela componente da veloci-

dade ortogonal ao escoamento livre em um ponto na esteira de vórtices, onde visualmente po-

demos observar e estimar o período da formação dos vórtices. A Fig. (37-b) mostra, através

de um “zoom” de um trecho da Fig. (37-a), a frequência predominante na esteira, cujos valo-

res estão de acordo com os encontrados na literatura.

(a) (b)

Figura 37. Comportamento da componente da velocidade na direção ortogonal ao escoamento

para um ponto na esteira de vórtices.

Visualmente pode-se observar, pela Fig. (37-b), que o período da esteira de vórtices está

muito próximo de 6 segundos, que representa o valor aproximado encontrado na literatura

[Brooks e Hughes, 1982, Li et al., 1991, Soto et al., 2001]. Este valor pode ser encontrado fa-

zendo-se uma análise da freqüência da formação de vórtices via FFT (Fast Fourier Transform)

analisando a componente da velocidade na direção ortogonal à direção do escoamento livre.

105

Esta análise foi realizada utilizando ferramentas disponíveis no MATLAB. A Fig. (38) mostra

o resultado obtido, cuja freqüência é 0,1675 Hz, logo o período é

sf

T 9701,51675,0

11≈== (5.2)

onde f é a frequência dada em Hertz, e T é o período dado em segundos.

Figura 38. Freqüência da formação de vórtices obtida via FFT.

5.4 Escoamento em Canal com Cavidade

Neste exemplo foi testada a capacidade da formulação em tratar domínios com estruturas

de pequena espessura. O modelo examinado é um canal tridimensional com uma cavidade

semicircular logo após uma placa no interior do canal [Lv et al., 2006]. Esta placa correspon-

de a uma membrana, que foi tratada como contorno sólido, logo, com condições de não desli-

zamento. O canal tem comprimento de 10L e largura de 1L, sendo L = 20mm. O raio da

cavidade semicircular é de 0,5L. A membrana tem comprimento de 0,5L e espessura de

0,5L/100 e faz um ângulo de 42,5º com a direção horizontal. As Figs. (39-41) abaixo mostram

a malha utilizada e neste caso evidenciamos a necessidade de um particionamento otimizado,

visando minimizar a comunicação entre os processadores. Neste caso foram utilizados dois

processadores, e a malha computacional possui 77.781 elementos tetraédricos, 96.225 arestas

e 14.372 nós, sendo que cada cor corresponde ao subdomínio que foi atribuído a cada proce-

sador.

106

Figura 39. Malha do canal com membrana e cavidade, com um particionamento realizado di-

vidindo-se o número de arestas pelo número de processadores.

Figura 40. Malha do canal com membrana e cavidade com particionamento otimizado mini-

mizando as interfaces entre os processadores e realizando balanço de carga.

107

Figura 41. Detalhe da malha do canal com membrana e cavidade com particionamento otimi-

zado mostrando a região da cavidade com grande concentração de nós.

Nas figuras anteriores fica evidente a redução na comunicação entre os processadores, vis-

to que as interfaces estão bastante reduzidas. Note no detalhe da malha reparticionada que

houve o balanceamento de carga e fez com que um processador tivesse seus nós concentrados

exatamente na região da cavidade, onde a malha é refinada, reduzindo significativamente as

interfaces entre os processadores.

As condições de contorno para este caso são:

- entrada: perfil de velocidade parabólico prescrito;

- saída: pressão prescrita 0p = e o tensor viscoso é prescrito como zero;

- condição de não deslizamento nas paredes do canal incluindo a membrana e a cavida-

de:u = 0 ;

- nas fronteiras frente e trás o tensor viscoso nas direções x e z foi prescrito como zero,

e a componente da velocidade na direção y foi prescrita como zero.

As simulações foram realizadas com números de Reynolds 200, e maxU é a máxima velo-

cidade na entrada, ajustada para se obter o Reynolds desejado. Neste caso o número de Rey-

nolds é dado por:

maxReU h

υ= (5.3)

Novamente ajustando o número de Reynolds pela viscosidade dinâmica µ .

108

Nas Figs. (42-45), mostramos o campo de pressão obtido no regime permanente, que de

forma qualitativa mostram a queda de pressão na região após a membrana, como esperado e

também a região de alta pressão onde o fluido encontra a membrana. Esta região de baixa

pressão sugere a formação de vórtices. As figuras a seguir mostram os resultados obtidos, on-

de estamos interessados na captura dos fenômenos físicos, como a formação de vórtices.

Figura 42. Campo de pressão no domínio tridimensional.

Figura 43. Detalhe do campo de pressão mostrando a queda de pressão após a membrana.

109

Figura 44. Detalhe do campo de pressão com isosuperficies de pressão.

Figura 45. Detalhe do campo de pressão com a malha utilizada e as isosuperfíceis de pressão.

A Fig. (46) a seguir mostra a formação de um grande vórtice após a membrana, que ocorre

devido à queda de pressão nesta região. As linhas de corrente mostradas na Fig. (47) compro-

vam a formação deste vórtice.

110

Figura 46. Detalhe do campo vetorial de velocidades.

Figura 47. Detalhe do campo vetorial de velocidades com as linhas de corrente.

5.5 Escoamento em Canal com Expansão Brusca

Neste exemplo foi testada a capacidade da formulação em capturar a física do problema,

que consiste na formação de vórtices simétricos logo após a expansão, e também a precisão

numérica na captura do número de Reynolds a partir do qual ocorre uma região de instabili-

dade que gera uma assimetria na formação dos vórtices. Esses vórtices permanecem simétri-

cos até um determinado número de Reynolds, e a partir desse número de Reynolds,

denominado crítico, a simetria é quebrada. A determinação desse número de Reynolds crítico

também é o objetivo desta análise. O modelo examinado é um canal tridimensional com uma

expansão brusca [Elias et al., 2005], como mostra a Fig. (48), que representa a geometria do

111

problema com as dimensões utilizadas. As Figs. (49) e (50) mostram a malha utilizada, consi-

derando um particionamento inicial e um particionamento otimizado, respectivamente, e cada

cor correspondendo ao subdomínio que foi atribuído a cada processador. Neste caso também

foram utilizados dois processadores, e a malha computacional possui 77.781 elementos tetra-

édricos, 96.225 arestas e 14.372 nós.

Figura 48. Geometria e dimensões do domínio de interesse.

Figura 49. Malha computacional com particionamento inicial.

112

Figura 50. Malha computacional com particionamento otimizado.

Para a determinação do número de Reynolds crítico para gerar instabilidade na simetria

dos vórtices, foi analisada a componente vertical do vetor velocidade em um ponto na região

de formação dos vórtices, que localiza-se conforme mostra a Fig. (51). A quebra da simetria é

determinada quando a componente da velocidade na direção ortogonal ao escoamento apre-

senta valores diferentes de zero.

Figura 51. Geometria e localização do ponto de análise.

As condições de contorno para este caso são:

- entrada: perfil de velocidade constante prescrito;

113

- saída: pressão prescrita 0p = e as “pseudo-tractions” prescritas como zero;

- condição de não deslizamento nas paredes do canal: u = 0 ;

- nas fronteiras frente e trás o tensor viscoso nas direções x e z foi prescrito como zero,

e a componente da velocidade na direção y foi prescrita como zero.

As Figs. (52-58) mostram os campos de velocidade e pressão considerando números de

Reynolds 40, 50, 80 e 130, respectivamente, onde podemos notar a simetria dos vórtices

quando Reynolds é 40 e 50, e a assimetria quando Reynolds é 80 e 130. E o número de Rey-

nolds é baseado na velocidade média prescrita na entrada, na dimensão do degrau da expan-

são, na massa específica do fluido e na viscosidade do fluido. Neste caso, o número de

Reynolds é controlado pela viscosidade, sendo que todas as outras variáveis são prescritas

como unitárias.

Figura 52. Campo de velocidade mostrando a simetria dos vórtices, Re = 40.

Note na Fig. (52) que os vórtices são simétricos, e portanto, não há regiões de instabilida-

de no escoamento.

114

Figura 53. Campo de velocidade mostrando a simetria dos vórtices, Re = 50.

Na Fig. (53), os vórtices continuam simétricos com número de Reynolds 50, porém apre-

sentando-se mais alongados do que para Reynolds 40. Este alongamento continuará a medida

que o número de Reynolds for aumentando.

As Fig. (54-56) mostram o escoamento pra Reynolds 80, onde já é possível observar cla-

ramente a assimetria dos vórtices.

Figura 54. Campo de velocidade mostrando a assimetria dos vórtices, Re = 80.

115

Figura 55. Detalhe do comprimento dos vórtices assimétricos, Re = 80.

Figura 56. Detalhe da assimetria dos vórtices e do domínio tridimensional, Re = 80.

As Figs. (57) e (58) mostram a assimetria dos vórtices principais, assim como a queda de

pressão na região de formação desses vórtices.

116

Figura 57. Campos de velocidade e pressão, Re = 130.

Figura 58. Campo de velocidade e pressão, Re = 130, mostrando a grande assimetria dos vór-

tices.

Podemos observar pelas Figs. (52-58) que com o aumento do número de Reynolds, ocorre

o alongamento dos vórtices simétricos formados após a expansão, e que este alongamento ge-

ra uma região de instabilidade na extremidade final dos vórtices, fazendo com que um dos

vórtices sobreponha sua energia sobre o outro causando a assimetria observada no regime

permanente, para altos números de Reynolds.

A Fig. (59) mostra o comportamento da componente vertical da velocidade para diferentes

números de Reynolds, onde podemos observar que o número de Reynolds crítico fica entre 56

e 58, sendo que estudos posteriores precisam ser efetuados para determinar com mais exatidão

o número de Reynolds crítico, porém podemos constatar que este resultado apresenta boa

concordância com os resultados obtidos por Elias et al., 2005, que obteve número de Rey-

117

nolds crítico 56, e com Battaglia et al., 1997 e Hawa e Rusak, 2000, que obtiveram número de

Reynolds crítico de 53,8, considerando-se fluidos Newtonianos e o número de Reynolds base-

ado na velocidade média na entrada. A quebra da simetria está claramente caracterizada na

Fig. (59). Desta forma, os resultados obtidos verificam além da captura dos fenômenos físicos

envolvidos, a boa precisão dos resultados numéricos obtidos quando comparados com os re-

sultados encontrados na literatura [Battaglia et al., 1997, Elias et al., 2005, Hawa e Rusak,

2000].

Figura 59. Componente vertical da velocidade no ponto P (especificado na Fig. (51)), para di-

ferentes números de Reynolds.

5.6 Estudos de Desempenho do Programa Paralelo

Em termos de desempenho de um programa paralelo, vários testes podem ser efetuados

para verificar a eficência do programa em termos de redução de custo computacional. Dentre

estes testes, destaca-se a análise que consiste, em linhas gerais, em determinar o tempo de

CPU à medida que aumentamos o número de processadores envolvidos na simulação, quando

comparado com àquele obtido com um único processador. Esta análise mostra como o uso de

técnicas de programação utilizando conceitos de computação paralela com memória distribuí-

da e/ou compartilhada pode aumentar o desempenho de aplicações científicas que demandam

grande esforço computacional.

118

O “speed-up” é determinado aplicando-se uma fórmula que relaciona os tempos de pro-

cessamento para um problema de referência com o mesmo problema aumentando o número

de processadores [Da Silva, 2008]. Em geral, o tempo de referência está relacionado com o

tempo de processamento considerando-se apenas um processador, ou seja, o problema se-

quencial. Neste trabalho não foi analisado o caso com apenas um processador, devido ao fato

de que não foi possível pré-processar a malha computacional em apenas uma CPU, por falta

de memória disponível. E também é sabido que a comparação com o resolvedor serial deve

ser realizada considerando-se o melhor resolvedor serial para o problema em análise, o que

também não dispomos. Sendo assim, nosso tempo de processamento de referência poderia ser

adotado como o tempo obtido com dois processadores, para algumas análises, e também usa-

remos uma extrapolação linear para determinar uma aproximação do custo computacional re-

lacionado à utilização de um único processador. A expressão matemática utilizada para

determinar o “speed-up” está mostrada na Eq. (5.4).

( )( , )

( , )refT n

S n pT n p

= (5.4)

onde n representa a carga de trabalho, p é o número de processadores envolvidos na simula-

ção, ( , )S n p é o “speed-up” determinado, ( )seqT n é o tempo de processamento utilizado do

algoritmo sequencial, utilizado como referência, e ( , )T n p é o tempo de processamento obti-

do com p processadores.

Para a realização dos estudos de desempenho utilizamos o problema da cavidade cúbica,

extensivamente utilizado para a verificação de programas de simulação de escoamentos flui-

dos [Ghia et al., 1982]. Neste problema os contornos da cavidade são paredes impermeáveis

fixas, e no contorno superior tem-se velocidade tangencial prescrita associada à movimenta-

ção da tampa, ou cavidade. A geometria e a malha obtida são encontradas na Fig. (59), onde

todas as dimensões são unitárias, ou seja, 1L = .

119

Figura 60. Domínio de interesse.

As condições de contorno para este caso são:

- tampa: perfil de velocidade unitário prescrito na direção x.

- esquerda e direita: tensor viscoso prescrito como zero nas direções x e z, e componen-

tes da velocidade na direção y prescrita como zero.

- frente, trás e fundo: condição de não deslizamento.

Neste exemplo foram obtidos resultados considerando Re = 1, e os resultados foram

comparados com os resultados encontrados na literatura (Ghia et al., 1982). A Fig. (61) mos-

tra o campo de velocidade com as linhas de corrente, evidenciando o vórtice formado, e as

Fig. (62) mostra a componente horizontal da velocidade, através das cores, e a Fig. (63) mos-

tra a componente vertical da velocidade, através das cores.

Figura 61. Campo de velocidade e linhas de corrente, Re = 1 .

120

Figura 62. Campo de velocidades e componente da velocidade na direção x.

Figura 63. Campo de velocidades e componente da velocidade na direção y.

Este exemplo foi resolvido considerando-se uma malha computacional com 1.149.802

nós, 7.124.082 elementos tetraédricos e 8.398.525 arestas, e esta malha foi particionada para

2, 4, 8, 12, 16, 20 e 32 processadores. As figuras a seguir mostram o domínio particionado em

8, 16 e 32 subdomínios.

121

Figura 64. Domínio particionado em 8 subdomínios.

Figura 65. Detalhe interno do domínio particionado em 8 subdomínios.

Figura 66. Domínio particionado em 16 subdomínios.

122

Figura 67. Detalhe interno do domínio particionado em 16 subdomínios.

Figura 68. Domínio particionado em 32 subdomínios.

Figura 69. Detalhe interno do domínio particionado em 32 subdomínios.

123

Figura 70. Campos de pressão e velocidade obtidos com 32 subdomínios, Re = 1.

Figura 71. Campo s de velocidade e pressão com 32 subdomínios evidenciando o vórtice cen-

tral, Re =1.

124

Figura 72. Campo de velocidade mostrando o vórtice central e alguns dos 32 subdomínios.

Os gráficos das Figs. (73) e (74) mostram os resultados obtidos em termos de desempe-

nho computacional, onde podemos verificar que com 8 processadores, o ganho de desempe-

nho não se comporta mais de forma linear, considerando as duas etapas da simulação. Isto

ocorre devido ao fato que a quantidade de nós por processador diminui com a quantidade de

processadores e a comunicação entre processadores aumenta, gerando custo computacional

nas trocas de mensagens. Ressaltamos que os tempos de processamento foram obridos consi-

derando-se duas etapas da simulação. A primeira etapa consiste no tempo necessário para a

aquisição dos dados geométricos e físicos, via arquivo de leitura, da alocação das variáveis

necessárias e dos mapeamentos realizados para minimizar a necessidade de comunicação pa-

ralela. Esta etapa consiste nas rotinas OpenPre, Mapping, Bcondit, e primeira passagem pela

rotina PreCalculus, como pode ser verificado na Fig. (17). A segunda etapa da simulação con-

siste no conjunto de rotinas envolvidas nas iterações que determinam o avanço no tempo, que

necessita intenso processamento de dados e grande quantidade de comunicação paralela. A

análise foi separada em etapas devido ao fato de que a primeira etapa é realizada apenas uma

única vez em cada simulação, e a segunda etapa é dependente do número de iterações neces-

sárias para a completa análise do problema estudado. A relevância da primeira etapa diminuiu

com a quantidade de iterações necessárias, sendo que o tempo necessário para o processamen-

to da primeira etapa ficou, neste exemplo, em torno de 5% do tempo necessário para uma ite-

ração completa no tempo. Este fato está de acordo com o esperado, visto que a segunda etapa

125

consiste na resolução de dez sistemas lineares, um para cada variável do problema, como po-

de ser observado na Fig. (17), que mostra o fluxograma geral do programa.

Figura 73. Redução do tempo de processamento com o aumento do número de processadores,

considerando a etapa 1.

126

Figura 74. Redução do tempo de processamento com o aumento do número de processadores,

considerando a etapa 2.

Podemos observar pelas Figs. (73) e (74) que com o aumento do número de processado-

res, diminui o tempo de processamento das duas etapas da simulação, sendo que esta redução

apresenta-se não linear, indicando que o custo das comunicações entre os processadores torna-

se cada vez maior.

Nas Figs. (75) e (76) mostramos a redução percentual do tempo de processamento em re-

lação ao tempo computado com dois processadores, quando comparamos os tempos da pri-

meira etapa, e em relação ao tempo extrapolado para um processador, quando comparamos os

tempos para a segunda etapa. Podemos observar que o incremento da redução diminui com a

quantidade de processadores, como esperado, visto que o custo da comunicação aumenta com

o número de processadores. Este custo de comunicação pode aumentar excessivamente, de tal

forma que o incremento da redução torne-se negativo.

127

4

59

8

85

12

89

16

90

20

91

32

93

0

20

40

60

80

100

Número deProcessadores

Redução Percentual

Figura 75. Redução percentual do tempo de processamento da primeira etapa em relação à

quantidade de processadores envolvidos na simulação.

2

50

4

73

8

88

12

91

16

93

20

94

32

96

0

20

40

60

80

100

Número deProcessadores

Redução Percentual

Figura 76. Redução percentual do tempo de processamento da segunda etapa em relação à

quantidade de processadores envolvidos na simulação.

Como neste trabalho não foi obtido o tempo de processamento considerando apenas um

processador, assumimos que a variação para dois e quatro processadores comporte-se de for-

ma linear e realizamos uma extrapolação linear para estimar o tempo de processamento par

um processador. Desta forma é possível obter o gráfico do “speed-up” obtido e compará-lo

com o “speed-up” ideal. Estes resultados estão na Fig. (77) para a segunda etapa de processa-

mento, que consiste na resolução dos sistemas lineares para todas as variáveis de análise. É

importante observar que o pré-condicionamento de sistemas lineares quando efetuados em sis-

temas distribuídos é realizado por subdomínio, o que implica no fato de que o pré-

condicionamento é função do número de subdomínios aos quais o domínio foi dividido, visto

que neste caso, cada subdomínio corresponte a um bloco ILU. E desta forma temos pré-

condicionamentos diferentes em cada simulação efetuada, fato que merece investigação no fu-

turo.

128

Figura 77. Curvas de “speed-up” para a segunda etapa de processamento.

Um parâmetro importante na avaliação de algoritmos paralelos é denominado escalabili-

dade, que é definido como a capacidade de manter constante a eficiência do algoritmo parale-

lo à medida que aumentamos o número de processadores envolvidos na análise e também a

carga de trabalho de cada processador. Escalabilidade, associado ao “solver”, será alvo de es-

tudos posteriores utilizando uma estratégia de aumentar o refinamento da malha computacio-

nal em cada subdomínio, de forma a manter a carga máxima de trabalho em cada processador,

visto que ainda não dispomos de um gerador de malhas paralelo.

Muitos estudos ainda devem ser efetuados na tentativa de obter ganho de desempenho

computacional, dentre estes estudos citamos o fato de que em computadores com múltipos

núcleos, a comunicação entre os núcleos de um mesmo nó é muito mais rápida do que a co-

municação entre núcleos que encontram-se em nós diferentes. Em vista disso, quando um

domínio é particionado em vários subdomínios, é interessante que subdomínios que apresen-

tem interesecção não nula sejam designados à núcleos que estejam em um mesmo nó para ti-

rar proveito desta comunicação mais rápida, e que sobdomínios que apresentem intersecção

129

nula, isto é, que não necessitam trocar informações entre si em nenhuma etapa do processa-

mento, sejam alocados preferencialmente em diferentes nós.

Ressaltamos ainda que todos os casos foram simulados no Núcleo de Atendimento em

Computação de Alto Desempenho (NACAD) da COPPE/UFRJ (http://www.nacad.ufrj.br/),

que é um laboratório especializado na aplicação de computação de alto desempenho a pro-

blemas de engenharia e ciências em geral, que está sob a coordenação dos professores Álvaro

Coutinho e Amit Bhaya, e que possui, entre outros, um “cluster” denominado Uranus do tipo

SGI Altix ICE 8200, com a seguinte configuração:

- 64 CPUs Quad Core Intel Xeon: 256 Cores;

- Memória: 512 Gbytes RAM (distribuída);

- Armazenamento em disco: SGI InfiniteStorage NAS (32 TBytes);

- Rede: Inifiniband e Gigabit;

- Sistema operacional: Suse Linux Enterprise Server + SGI ProPack;

- Compiladores: Intel e GNU (Fortran-90 e C/C++) com suporte OpenMP e MPI .

130

Capítulo 6

Conclusões e Sugestões para Trabalhos

Futuros

6.1 Conclusões

Uma formulação matemático-numérica, monolítica, implícita, e com propriedades de es-

tabilidade para a pressão e para a convecção, discretizada através da aplicação do Método dos

Elementos Finitos e o Método do Passo Fracional, baseado em uma estrutura de dados por a-

restas, e utilizando computadores com memória distribuída foi apresentada para resolver pro-

blemas governados pelas equações de Navier-Stokes escritas na sua forma incompressível,

nas variáveis primitivas e sobre domínios tridimensionais. Resultados numéricos foram obti-

dos, e comparados com resultados obtidos por outros autores, indicando de forma promissora

a eficiência da formulação para simular escoamentos laminares com diferentes números de

Reynolds. O caráter implícito da formulação é atrativo para análise de problemas transientes

onde o passo de tempo crítico do problema (limite de estabilidade da formulação explícita) é

ordens de magnitude menor do que o passo de tempo necessário para a obtenção de resultados

acurados, ou seja, o passo de tempo físico.

131

Os casos estudados neste trabalho compreendem diversas situações que podem ocorrer na

Dinâmica dos Fluidos Computacional, quando os fenômenos caracterizam escoamentos in-

compressíveis. Estudamos situações onde existe um regime permanente, obtendo dados que

permitem avaliar a ferramenta desenvolvida em termos de captura dos fenômenos físicos en-

volvidos e da própria solução em regime permanente. Estudamos também casos onde não há o

regime permanente, analisando a precisão temporal em situações onde, por exemplo, existe

um regime periódico, como a formação da esteira de vórtices no escoamento sobre um cilin-

dro circular. Todas as etapas deste trabalho estão descritas nos capítulos anteriores com deta-

lhes.

Este é um vasto campo de pesquisa e certamente são muitos os avanços a serem realizados

de maneira a tornar nossa ferramenta computacional mais versátil e eficaz. Contudo, ao final

deste trabalho obteve-se uma ferramenta capaz de tratar escoamentos governados pelas equa-

ções de Navier-Stokes na sua forma incompressível sobre domínios tridimensionais de manei-

ra eficiente e com baixo custo, viabilizando as análises de diferentes situações encontradas na

engenharia.

Ressaltamos ainda que este trabalho foi desenvolvido a partir da linha zero nos desenvol-

vimentos teórico e computacional, sendo que diversas etapas foram concluídas superando

grandes dificuldades. Dentre estas etapas, destacamos que a formulação matemática foi pri-

meiramente desenvolvida em 2D e implementada serialmente em MATLAB, onde constata-

mos a eficiência da formulação quanto à discretização espacial por Elementos Finitos e

Método do Passo Fracional e temporal utilizando “Euler Backward”. O desenvolvimento da

formulação numérica, através de uma estrutura de dados baseada nas arestas foi obtido, com o

desenvolvimento de um pré-condicionador para determinar os coeficientes das arestas. A par-

tir deste ponto foi desenvolvida a formulação matemática para domínios 3D, com a determi-

nação analítica dos coeficientes das arestas e desenvolvimento do pré-condicionador para

domínios 3D, que foi implementado pelo professor Ramiro Brito, e posteirormente por Rogé-

rio Soares. Neste ponto iniciaram-se os estudos de computação paralela, onde foi iniciado o

desenvolvimento de ferramentas para este tipo de arquitetura de computadores. Primeiramente

foi desenvolvido um programa para particionar um domínio 3D considerando-se as arestas

dos elementos tetraédricos da malha. Este programa utilizou extensivamente as capacidades

do MPI e do PETSc para montar uma matriz de adjacências das arestas, sem que, em nenhum

momento a malha computacional seja totalmente carregada por um único processador. E após

obter uma malha particionada iniciamos o desenvolvimento do programa de análise de esco-

amentos fluidos, que foi desenvolvido de forma que nenhuma etapa da simulação seja realiza-

132

da de maneira seqüencial. Novamente foram utilizadas em grande quantidade as ferramentas

do MPI e PETSc. As grandes dificuldades encontradas ao longo da realização deste trabalho

estão centradas na implementação paralela, pois como autor deste trabalho não possuía expe-

riência na área, e as questões relativas à computação paralela não são nada triviais. A cada di-

ficuldade encontrada havia uma extensa pesquisa para o desenvolvimento da solução. Sendo

assim, este trabalho constitui-se, antes de qualquer coisa, em um intenso aprendizado e en-

grandecimento científico e acadêmico do autor.

6.2 Trabalhos Futuros

O sistema computacional desenvolvido neste trabalho, como mencionado anteirormente,

nos possibilita a expandir nossa capacidade de análise, se considerarmos a inclusão de novas

flexibilidades, inicialmente com estudos mais aprofundados na investigação de resolvedores

para sistemas lineares não simétricos, de forma a determinar de maneira mais precisa os pa-

râmetros de convergência de métodos como o GMRES, e também na investigação de desem-

penho computacional, analisando curvas de “speed-up” para diferentes problemas de maneira

a se obter as melhores condições de processamento para cada caso. De forma mais detalhada,

citamos a seguir algumas linhas que serão desenvolvidas como seqüência deste trabalho:

Interação Fluido-Estrutura

Neste caso, a formulação matemática já foi desenvolvida para tratar fenômenos onde um

fluido interage com uma estrutura. A implementação computacional usando uma descrição do

tipo ALE (Arbitrary Lagrangean-Eulerian) já foi realizada em domínios 2D pelo autor deste

trabalho [Antunes et al., 2005], incluindo ferramentas de movimentação e atualização de ma-

lhas. Devido à forma estruturada como foi desenvolvido o algoritmo principal deste trabalho,

acreditamos que esta nova etapa será rapidamente concluída.

Este tema é de extremo interesse social e econômico, pois a simulação de escoamentos

incompressíveis sobre domínios tridimensionais, que caracteriza uma importante classe de a-

plicações em muitos campos da engenharia, tendo como exemplos, o escoamento sobre estru-

turas de grande altura, como prédios e torres, escoamentos internos em dutos flexíveis, como

“risers” utilizados na indústria do petróleo, escoamentos ao redor de estruturas aero e hidrodi-

133

nâmicas, permitindo a determinação de parâmetros que influenciam na integridade destas es-

truturas com antecedência. Nesta mesma linha de raciocínio, estruturas de grande altura estão

submetidas à ação de ventos que interagem com a estrutura gerando oscilações prejudiciais à

integridade da estrutura e desconfortável para os seres humanos, uma vez que o escoamento

fluido influencia o comportamento estrutural e vice-versa. Neste sentido, a simulação dos fe-

nômenos fluido e estrutural deve ser realizada de forma acoplada para modelarmos de forma

mais realista esta classe de problemas. Ainda neste sentido, é possível determinar as melhores

condições de funcionamento das estruturas conhecendo-se os efeitos dos escoamentos fluidos,

em termos de carregamentos, sobre estas estruturas. Neste caso necessita-se o desenvolvimen-

to das formulações matemáticas dos fenômenos fluido e estrutural, assim como o correto tra-

tamento da interface entre os campos e o tratamento da malha dinâmica. Em termos do

tratamento da interface entre os campos, ressaltamos que o desenvolvimento de ferramentas

específicas para este fim, possibilita o acoplamento não apenas de diferentes campos físicos,

mas também de diferentes programas computacionais que resolvam estes campos, por exem-

plo, é possível acoplar o nosso simulador de escoamentos fluidos com o ANSYS como simu-

lador do campo estrutural. Desta forma é extremamente importante esta etapa do

desenvolvimento futura. Também vale ressaltar que estes desenvolvimentos devem ser com-

patíveis para a simulação em computadores paralelos com memória distribuída. Desta forma

uma aplicação futura das ferramentas de alto desempenho desenvolvidas nesta tese será o es-

tudo de fenômenos fluido-estruturais acoplados.

Modelos de Turbulência

Incorporar ao sistema computacional desenvolvido a capacidade de analisar fenômenos

com elevados números de Reynolds também é extremamente importante, visto que muitos fe-

nômenos de interesse desenvolvem-se em regimes turbulentos. Em vista disso, alguns mode-

los de turbulência apresentam características muito similares às características da formulação

matemática desenvolvida neste trabalho, visto que os modelos de estabilização das equações

neste trabalho apresentam estreita relação com os modelos de turbulência de sub-escala [De

Sampaio et al., 2001, De Sampaio et al., 2007, Lesieur et al., 2005, Mathieu e Scott, 2000].

Os modelos de sub-escala ou Sub-Grid Scale (SGS) baseiam-se no fato de que apenas as me-

nores escalas devem ser modeladas, e que as grandes escalas são resolvidas pela malha com-

putacional. Neste sentido, a simulação numérica de grandes escalas (LES) é uma alternativa

134

que possui grande atrativo do ponto de vista computacional, pois aproveita-se do fato que os

maiores turbilhões dominam a física de qualquer escoamento turbulento, sendo os responsá-

veis pelo transporte da maior parte da quantidade de movimento e energia, especialmente nas

regiões de recirculação e nos escoamentos cisalhantes livres. Este objetivo é alcançado com a

aplicação de um filtro nas equações governantes do escoamento, que divide o campo de esco-

amento em grandes escalas, que representa a parte filtrada das equações, e pequenas escalas,

que correspondem às escalas submalha. Diversos modelos de escalas submalha podem ser tes-

tados, dentre eles citamos: Modelo de Smagorinsk, Combinação Linear e Modelo Dinâmico,

como os mais populares. Este tema também representa uma possível continuidade deste traba-

lho.

Fluidos Não-Newtonianos

Incluir modelos de fluidos não-newtonianos é interessante para podermos analisar fenô-

menos importantes, como por exemplo o escoamento sanguíneo, visto que atualmente uma

área em extrema expansão é a bioengenharia, e muitos avanços têm sido feitos nesta área.

Neste caso determinar expressões para a viscosidade de fluidos não-newtonianos é de interes-

se deste autor.

Adaptação de Malha

Mesmo que esta ferramenta seja funcional para ambientes paralelos possibilitando assim o

uso de malhas refinadas, é sempre interessante a capacidade de adaptar malhas e gerar com is-

so garantias de acurácia e ganho computacional na análise de problemas de grande porte.

Também é atrativa esta capacidade na captura dos fenômenos de pequena escala envolvidos,

visto que estes muitas vezes necessitam diferentes graus de refinamento em diferentes regiões

do domínio. Atualmente o grupo de pesquisa PADMEC (Processamento de Alto Desempenho

na Mecânica Computacional), ao qual o autor faz parte, possui integrantes desenvolvendo es-

tudos nesta área, e já com várias capacidades desenvolvidas, e que poderia vir a ser incorpo-

rada no nosso programa.

135

Otimização Estrutural com Interação Fluido-Estrutura

Também é colocado como projeto futuro a inclusão da capacidade de realização de otimi-

zação de forma, visto que a otimalidade em projetos é hoje uma crescente área de estudo, e

neste caminho, realizar esta otimização em fenômenos onde ocorrem iterações fluido-

estruturais é de grande importância em várias áreas da engenharia, como por exemplo, o pro-

jeto ótimo de estruturas que estão submetidas a carregamentos provenientes de escoamentos

fluidos, tais como na indústria aeronáutica e na indústria do petróleo gerando oscilações peri-

gosas para a sua integridade estrutural.

136

Referências Bibliográficas

Almeida, R. C., Galeão, A. C., 1996, An Adaptive Petrov-Galerkin Formulation for Compres-

sible Euler and Navier-Stokes Equations, Comput. Methods Appl. Mech. Engrg., vol 129,

p. 157-176.

Anderson, J. D., 1995, Computational Fluid-Dynamics, New York, McGraw-Hill, Inc., 547p.

Antunes, A. R. E., Lyra, P. R. M., Willmersdorf, R. B., Soares, R., 2008, Solução das

Equações de Navier-Stokes Incompressíveis 3D Via Método dos Elementos Finitos e

Fractional Step Utilizando Computadores Paralelos com Memória Distribuída, XXIX

Iberian Latin American Congress on Computational methods in Engineering, Maceió-AL,

Brasil.

Antunes, A. R. E., Lyra, P. R. M., Soares, R., 2008, Utilizando Computadores Paralelos com

Memória Distribuída e o Método dos Elementos Finitos em Domínios 3D, Encontro

Regional de matemática Aplicada, Natal-RN, Brasil.

Antunes, A. R. E, Lyra, P. R. M., Willmersdorf, R. B., Soares, R., 2008, Conceitos e

Aplicações de Computação de Alto Desempenho na Dinâmica dos FLuidos

Computacional, Seminário de Computação Paralela, Fortaleza-CE, Brasil.

Antunes, A. R. E., Lyra, P. R. M., Mendonça, A. P., 2006, Simulação Numérica do

Escoamento Sanguíneo na Bifurcação da Artéria Carótida Humana, Simpósio Mineiro de

Mecânica Computacional, Araxá-MG, Brasil.

Antunes, A. R. E., Lyra, P. R. M., Farah, D., 2005, O Método dos Elementos Finitos e o

Fractional Step com Estrutura de Dados por Arestas Aplicado às Equações de navier-

Stokes Incompressíveis, XXVI Iberian Latin American Congress on Computational

Methods in Engineering, Guarapari-ES, Brasil.

Antunes, A. R. E., Lyra, P. R. M., Willmersdorf, R. B., 2005, A Methodology and

Computational System for the Simulation of Fluid-Structure Interaction Problem, J. of the

Braz. Soc. Mech. Sci. & Eng., vol XXVII, No.3, p. 255-265.

Armaly, B. F., Durst, F., Pereira, J. C. F., Schönubg, B., 1983, Experimental and Theoretical

Investigation in Backward-Facing Step Flow, Journal Fluid Mech., vol 127, p. 473-496.

137

Babuska, I., 1070/1971, Error-bounds for Finite Element Method, Numer. Math., vol 16, p.

322-333.

Balay, S., Buschelman, K., Eijkhout, V., Gorpp, W., Kaushik, D., Knepley, M., McInnes, L.

C., Smith, B., Zhang, H., 2005, PETSc Users Manual, Argonne National Laboratory, 179p.

Barth, T. J., 1992, Aspects of Unstructured Grids and Finite-Volume Solvers for the Euler and

Navier-Stokes Equations, AGARD Report 787, p. 6.1-6.61.

Bathe, K., Zhang, H., 2002, A Flow-Condition-Based Interpolation Finite Element Procedure

for Incompressible Fluid Flow, Computers and Structures, vol 80, p. 1267-1277.

Battaglia, F., Tavener, S. J., Kulkarni, A. K., Merkle, C. L., 1997, Bifurcation of Low

Reynolds Number Flows in Symmetric Channels, AIAA J, p. 35:99.

Boivin, S., Cayré, F., Hérard, J., 2000, A Finite Volume Method to Solve Navier-Stokes

Equations for Incompressible Flows on Unstructured Meshes, Int. J. Therm. Sci, vol 39, p.

806-825.

Brazza, M., Chassaing, P. & Minh, H. H., 1986, Numerical Study and Physical Analysis of

the Pressure and Velocity Fields in the Near Wake of a Circular Cylinder, J. Fluid Mech,

vol 165, pp. 79-139.

Brezzi, F., 1974, On the Existence . Uniqueness and Approximation of Saddle Point Problems

Arising from Lagrangian Multipliers, RAIRO, Anal, Num., vol 8 (R2), p. 129-151.

Brooks, A. N. & Hughes, T. J. R., 1982, Streamline Upwind/Petrov-Galerkin Formulations

for Convection Dominated Flows with Particular Emphasis on the Incompressible Navier-

Stokes Equations, Computer Methods in Applied Mechanics and Engineering, vol 32, pp.

199-259.

Catabriga, L., 2000, Soluções Implícitas das Equações de Euler Empregando Estrutura de Da-

dos por Arestas, DSc Thesis, Universidade Federal do Rio de Janeiro, Departamento de

Engenharia Civil, Rio de Janeiro-RJ.

Catabriga, L., Coutinho, A. L. G. A., 2002, Implicit SUPG Solution of Euler Equations Using

Edge-Based Data Structures, Comput. Methods Appl. Mech Engrg, vol 191, p. 3477-3490.

138

Cebral, J., Yim, P., Löhner, R., Soto, O., Marcos, H., Choyke, P., 2001, New Methods for

Computational Fluid Dynamics of Carotid Artery from Magnetic Resonance Angiography,

in: Proceedings of SPIE Medical Imaging, vol 4321, paper No. 22, San Diego, California,

February.

Chang, W., Giraldo, F., Perot, B., 2002, Analysis of an Exact Fractional Step Method, J.

Computational Physics, vol 180, p. 183-199.

Chorin, A. J., 1967, A Numerical Method for Solving Incompressible Viscous Problem, J.

Comput. Phys., vol 2.

Chorin, A. J., 1968, Numerical Solution of the Navier-Stokes Equations, Mathematics of

Computation, vol 22, p. 745-762.

Codina, R., Blasco, J., 1997, A Finite Element Formulation for the Stokes Problem Allowing

Equal Velocity-Pressure Interpolation, Comput. Methods Appl. Mech. Engrg., vol 143, p.

373-391.

Codina, R., 2000, Stabilization of Incompressible and Convection Trough Orthogonal Sub-

scales in Finite Element Methods, Compt. Meth. Applied mechanics and Engineering, v.

190, p. 1579-1599.

Codina, R., 2001, Pressure Stability in Fractional Step Finite Element Methods for

Incompressible Flows, J. Computational Physics, vol 170, p. 112-140.

Codina, R., 2002, Stability Finite Element Approximation of Transient Flows Using

Orthogonal Subscales, Comp. Meth. Applied mechanics and Engineering, vol 191, p.

4295-4321.

Codina, R., Soto, O., 2003, Approximation of the Incompressible Navier-Stokes Equations

Using Orthogonal Subscale Stabilization and Pressure Segregation on Anisotropic Finite

Element Meshes, Comp. Methods Appl. Mech. Engrg, vol 193, p. 1403-1419.

Codina, R., Owen, H. C., Nthiarasu, P., Liu, C. B., 2006, Numerical Comparison of CBS and

SGS as Stabilization Techniques for the Incompressible Navier-Stokes Equations,

International Journal of Numerical Methods in Engineering, vol 66, p. 1672-1689.

139

Costa, G. K., 2004, Simulação Tridimensional de Escoamentos Compressíveis e Incompressí-

veis Através do Método dos Elementos Finitos, DSc Thesis, Universidade Federal de Per-

nambuco, Departamento de Engenharia Nuclear, Recife-PE.

Currie, I. G., 1974, Fundamental Mechanics of Fluids, New York, McGraw-Hill, Inc., 435p.

Da Silva, R. S., 2008, Simulação de Escoamento Bifásico Óleo-Água em Reservatórios de Pe-

tróleo usando Computadores Paralelos de Memória Distribuída, DSc Thesis, Universidade

Federal de Pernambuco, Departamento de Engenharia Civil, Recife-PE.

De Sampaio, P.A.B., 1991, A Petrov-Galerkin Formulation for the Incompressible Navier-

Stokes Equations using Equal Ordering for Velocity and Pressure, Int. J. Numer. Meth.

Engrg., vol. 31, pp. 1135-1149.

De Sampaio, P. A. B., Lyra, P. R. M., Morgan, K. & Weatherill, N. P., 1993, Petrov-Galerkin

Solutions of Incompressible Navier-Stokes Equations in Primitive Variables with

Adaptative Remeshing, Computer Methods in Applied Mechanics and Engineering, vol

106, pp. 143-178.

De Sampaio, P. A. B., Coutinho, A. L. G. A., Hallak, P. H., Pfeil, M. S., 2001, Numerically

Implicit Sub-grid Modeling of Turbulent Flows trough a Stabilized Finite Element Method,

Proceedings of the 6º US National Congress on Computational Mechanics, Dearborn-

Michigan, USA.

De Sampaio, P. A. B., Junior, M., A. G., Lapa, C. M. F., 2007, A CFD Approach to the

Atmospheric Dispersion of Radionuclides in the Vicinity of NPPs, Nuclear Engineering

and Design, Article In Press.

De Sampaio, P. A. B. & Coutinho, A. L. G. A., 1999, Simulation of Free and Forced

Convection Incompressible Flows an Adaptative Parallel/Vector Finite Element

Procedure, Int. J. Num. Meth. Engr., vol 29, pp. 289-309.

Donea, J., Giuliani, S., Laval, H., Quartapelle, L, 1982, Finite Element Solution of the Navier-

Stokes Equations by a Fractional Step Method, Comput. Methods Appl. Mech. Engrg, vol

30, p. 53-73.

Donea, J., Huerta, A., 2003, Finite Element Methods for Flow Problems, John Wiley & Sons,

1ª edition.

140

Elias, R. N., Martins, M. A. D., Coutinho, A. L. G. A., 2005, Parallel Edge-Based Solution of

Viscoplastic Flows with the SUPG/PSPG Formulation, Comp. Mech., Springer-Verlag,

Accepted in 17 September 2005.

Franca, L. P., Hughes, T. R. J., 1988, Two Classes of Finite Element Methods, Comput.

Methods Appl. Mech. Engrg., vol 69, p. 89-129.

Gmsh, online documentation, http://www.geuz.org/gmsh/.

Guermond, J. L., Minev, P., Shen, J., 2006, An Overview of Projection Methods for

Incompressible Flows, Comput. Methods Appl. Mech. Engrg, vol 195, p. 6011-6045.

Gresho, P. M., Chan, S. T., Lee, R. L. & Upson, C. D., 1984, A Modified Finite Element

Method for Solving the Time-Dependent, Incompressible Navier-Stokes Equations Part 2:

Application”, Int. J. Numer Methods Fluids, vol 4, pp. 619-640,.

Gresho, P. M., Sani, R. L., 1998, Incompressible Flow and the Finite Element Method, John

Wiley & Sons, vol 1.

Gresho, P. M., Sani, R. L., 1998, Incompressible Flow and the Finite Element Method, John

Wiley & Sons, vol 2.

Hawa, T., Rusak, Z., 2000, Viscous Flow in a Slightly asymmetric Channel with a Sudden

Expansion, Phys. Fluids, 12:2257.

Henriksen, M. O., Holmen, J., 2002, Algebraic Splitting for Incompressible Navier-Stokes

Equations, J. Computational Physics, vol 175, p. 438-453.

Hughes, T. J. R., Franca, L. P. & Ballestra, M., 1986, A New Finite Element Formulation for

Computational Fluid Dynamics: V. Circunventing the Babuska-Brezzi Condition: A

Stable Petrov-Galerkin Formulation for Stokes Problem Accommodating Equal-Ordering

Interpolation, Computer Methods in Applied Mechanics and Engineering, vol 59, p. 85-

99.

Hughes, T. J. R., 1987, The Finite Element Method, 1ª Edition, United States of America,

Prentice-Hall, Inc., 803p.

141

Hughes, T. R. J., Franca, L. P., Hulbert, G. M., 1989, A new Finite Element Formulation for

Computational Fluid mechanics: VIII. The Galerkin/Least-Squares metoh for Advective

Diffusive Equations, Comput. Methods Appl. Mech. Engrg., vol 73, p. 173-189.

Kalro, V., Tezduyar, T., 1998, 3D Computation of Unsteady Flow Past a Sphere with a

Parallel Finite Element Method, Comput. Methods Appl. Mech. Engrg., vol 151, p. 267-

276.

Karypis, G., Schloegel, K., Kumar, V., 2003, PARMETIS (Parallel Graph Partitioning and

Saprse Matrix Orderinging, Version 3.1, University of Minnesota, 29p.

Ladyzhenskaya, O. A., 1969, The Mathematical Theory of Viscous Incompressible Flow, 2nd

Ed., Gordon and Breach, New York, 1969.

Lesieur, M., Métais, O., Comte, P., 2005, Large-Eddy Simulations of Turbulence, Cambridge

University Press.

Li, J., Chambarel, M., Donneaud, M. & Martin, R., 1991, Numerical Study of Laminar Flow

Past One and two Circular Cylinders, Comput. Fluids, vol 19, pp. 155-170.

Löhner, R., 2001, Applied CFD Techniques - An Introduction Based on Finite Element

Methods, john Wiley & Sons Ltd..

Löhner, R., Galle, M., 2002, Minimization of Indirect Addressing for Edge-Based Field

Solvers, AIAA-02-0967.

Löhner, R., 2004, Multistage Explicit Advective Prediction for Projection-Type

Incompressible Flow Solvers, vol 195, p. 143-152.

Löhner, R., Yang, C., Cebral, J., Camelli, F., Soto, O., Waltz, J., 2006, Improving the Speed

and Accuracy of Projection-Type Incompressible Flow Solvers, vol 195, p. 3087-3109.

Lv, X, Zhao, Y., Huang, X. Y., Xia, G. H., Wang, Z. J., 2006, An Efficient

Parallel/Unstructured-Multigrid Preconditioned Implicit Method for Simulating 3D

Unsteady Compressible Flows Moving Objects, Journal of Computational Physics, vol

215, p 661-690.

Luo, H., Sharov, D., Baum, J., D., Löhner, R., 2003, Parallel Unstructured Grid GMRES+LU-

SGS method for turbulent flows, AIAA.

142

Mathieu, J., Scott, J., 2000, An Introduction to Turbulent Flow, Cambridge University Press.

Metcalf, M., Reid, J., 1996, FORTRAN 90/95 Explained, Oxford Science Publications, 345p.

Pacheco, P. S., 1995, MPI User’s Guide in FORTRAN .

Pacheco, P. S., 1997, Parallel Programming with MPI, Morgan Kaufmann Publishers, 419p.

Paraview, online documentation, http://www.paraview.org/.

Peraire, J., Vahdati, M., Morgan, K., 1993, Finite Element Multigrid Solution of Euler Flows

Past Installed Aero-Engines, J. Computational Mechanics, vol 11, p. 433-451

Perot, J. B., 1993, An Analysis of the Fractional Step Method, J. Computational Physics, vol

108, p. 51-58.

PETSc online documentation, Web page. http:/www-unix.mcs.anl.gov/petsc/petsc-

as/documentation.

Ramamurti, R., Löhner, R., 1996, A Parallel Implicit Incompressible Flow Solver Using

Unstructured Meshes; Computers and Fluids, vol 5, p. 119-132.

Saad, Yousef, 2003, Iterative Methods for Sparse Linear Systems, 2nd edition, SIAM.

Soto, O., Löhner, R., Cebral, J., 2001, An Implicit Monolithic Accurate Finite Element

Scheme for Incompressible Flow Problems. In: 15th AIAA Computational Fluid Dynamics

Conference, Anaheim, EUA, june.

Soto, O., Löhner. R., Cebral, J., Camelli, F. 2004, A Stabilized Edge-Based Implicit

Incompressible Flow Formulation, Comput. Methods Appl. Mech. Engrg., vol 193, p.

2139-2154.

Soulamani, A., Saad, Y., Rebaine, A., 2001, An Edge Based Stabilized Finite Element

Method for Solving Compressible Flows: Formulation and Parallel Implementation,

Comput. Methods, Appl. Mech. Engrg, vol 190, p. 6735-6761.

Temam, R., 1969, Sur l’approximation de la Solutin des Équaciones de Navier-Stokes par la

Méthode dês pás Fractionaires (I), Arch. Ration. Mech. Anal., vol 32.

143

Tezduyar, T., Aliabadi, S., Behr, M., Johnson, A., Kalro, V. & Litke, M., 1997, High

Performance Computing Techniques for Flow Simulations in: Parallel Solution Methods

in Computational Mechanics, Papadrakakis, M. (Ed.), John Wiley and Sons, London, UK.

VisIt, online documentation, https://wci.llnl.gov/codes/VisIt/.

White, F. C., 1991, Viscous Fluid Flow, 2ª Edition, Singapore, McGraw-Hill, Inc., 614p.

Williams, P. T., Baker, A. J., 1997, Numerical Simulations of Laminar Flow Over a 3D

Backward-Facing Step, Int. Journal Num. Meth. In Fluids, vol 24, p. 1159-1183.

Woll, C., 2004, A Modern Approach for the Computational Simulation of Flow in Porous

Media with Parallel Computing, Master Thesis, Technische Universität Braunschweig,

Institut für Thermodynamik, Braunschweig, Germain.

Zienkiewicz, O. C., Morgan, K., 1983, Finite Elements and Approximation, John Wiley &

Sons, Inc, 1ª edition.

Zienkiewicz, O. C., Taylor, R. L., 1988, The Finite Element Method: Basic Formulation and

Linear Problems, vol 1, MacGraw-Hill, 4ª edition.

Zienkiewicz, O. C. & Wu, J. 1991, Incompressible Without Tears-How to Avoid Restrictions

of Mixed Formulation, Int. J. Num. Meth. Eng., vol 32, pp 1189-1203.

Zienkiewicz, O. C., Mithiarasu, P., Codina, R., Ortiz, P., 1999, The Chacarteristic-Based-Split

Procedure: an Efficient and Accurate Algorithm for Fluid Problems, Int. Journal Num. Me-

th. In Fluids, vol 31, p. 359-392.

144

Apêndice A: Matriz de Adjacências

Neste apêndice, discutimos a criação e montagem da matriz de adjacências necessária para

se obter um particionamento otimizado de um grafo não orientado. Esta montagem é baseada

no trabalho de Woll, 2004.

O particionamento de grafos através do PARMETIS é baseado na matriz de adjacências

deste grafo, que é armazenada em uma estrutura de dados do tipo CSR (Compressed Sparse

Row). Com esta estrutura e a função MatCreateMPIAdj do PETSc torna-se relativamente fácil

particionar um grafo, pois a interface do PETSc para a chamada do PARMETIS é bastante

amigável. A dificuldade então está na montagem da matriz de adajcências, que neste trabalho

consiste em uma matriz de adjacências de arestas, ou seja, é necessário conhecer todas as a-

restas que têm conectividade não nula com cada aresta da malha. Sendo assim, cada aresta

precisa conhecer suas arestas vizinhas. Esta tarefa seria fácil se realizada de forma serial, po-

rém como estamos partindo sempre de um particionamento inicial, ou seja, cada processador

possui informações apenas sobre um subconjunto das arestas totais, são necessárias várias tro-

cas de mensagens entre os processadores para se obter tal estrutura de dados.

A criação da estrutura de dados do tipo aresta ao redor de aresta neste trabalho, seguiu os

seguintes passos:

Passo 1: Criação de uma lista de nós conectados às arestas locais – neste caso um “loop” é

realizado em todas as arestas locais a cada processador e são obtidos todos os nós que estão

em cada processador;

Passo 2: Determinação de nós com duplicidade de processadores – cada processador envia

sua lista de nós obtida no passo 1, de forma a determinar quais nós possuem cópias em mais

de um processador e em quais, simplesmente comparando com a sua lista local. Esta troca de

mensagem é realizada utilizando um esquema de comunicação do tipo “Round-Robin”, com

funções MPI_Send e MPI_Receive, ou com a função MPI_Allgatherv do MPI. Destacamos

145

que ao final desta tese, quando houve a necessidade de troca de informações de tamanho

grande entre processadores, o esquema Round-Robin não funcionou, sendo necessário uma al-

teração, onde substituímos este esquema pelas funções do MPI, MPI_Allgather e

MPI_Allgatherv.

Passo 3: Criação de uma lista local de arestas ao redor de nós – um “loop” é realizado nas

arestas e é criado um vetor para cada nó, de forma a armazenar as arestas que tem conectivi-

dade com cada nó. Esta estrutura é necessária, pois os nós que estiverem em interfaces entre

processadores estarão compartilhando as mesmas arestas em processadores diferentes.

Passo 4: Criação da lista completa de arestas ao redor de nós – cada processador envia a

lista criada no passo 3 para os outros processadores e inclui na sua lista local as arestas que

têm conectividade com os seus nós locais. Neste caso cada processador envia e recebe apenas

para os processadores com quem compartilha nós, sendo assim, utilizamos um tipo de comu-

nicação sem bloqueio, MPI_ISEND e MPI_IRECEIVE.

Passo 5: Criação da lista de arestas ao redor de arestas – é criado um vetor para armazenar

as conectividades de cada aresta e um vetor apontador destas conectividades. Neste caso é

preciso tomar cuidado com os índices das arestas, pois devem coincidir com as linhas da ma-

triz de adjacências que será criada. Em caso de programação em FORTRAN, deve começar

com índice zero, pois o PETSC inicia seus vetores e matrizes na linha e coluna zero.

Após a criação da matriz de adjacências, e através da chamada de algumas funções do

PETSc, o particionamento é obtido. Porém, o PETSc não move dados entre processadores, ele

apenas retorna um mapeamento, indicando que aresta deve ser alocada em cada processador.

Desta forma, comunicação adicional deve haver para realocar as arestas nos novos processa-

dores. Neste trabalho, a lista local obtida em cada processador é enviada a todos os processa-

dores e, estes, obtém quais arestas devem permanecer locais e quais serão remotas. Note que

como estamos realizando um particionamento por arestas, as interfaces são dadas pelos nós,

pois não temos uma aresta em mais de um processador. Isto significa que o “loop” nas arestas

necessário para montar o LHS e o RHS será dado sempre em estruturas locais em cada pro-

cessador. Isto é interessante pelo fato de que os coeficientes das arestas não precisam de có-

pias remotas. Antes de afetuar o particionamento propriamente dito, as seguintes etapas

devem ser efetuadas:

146

Passo 1: Mapeamento entre o particionamento antigo e o novo – cada processador envia

sua lista de arestas para os outros processadores, via esquema Round-Robin, e compara para

determinar quais arestas permanecerão locais a cada processador.

Passo 2: Mapeamento entre os nós locais e remotos – cada processador cria uma lista de

nós locais e envia para cada processador, via esquema Round-Robin, de forma a comparar

com sua lista local e descobrir com quais processadores haverá interface.

Neste ponto cada processador contém as informações sobre suas arestas e nós locais e so-

bre seus nós que possuem cópias remotas em outros processadores.

A seguir, mostramos as funções necessárias para realizar o particionamento:

! Create adjacency matrix

call MatCreateMPIAdj(PETSC_COMM_WORLD,rows,cols,ia_new,ja_new,

PETSC_NULL_INTEGER,Adj,ierr)

!

call MatAssemblyBegin(Adj,MAT_FINAL_ASSEMBLY,ierr)

!

call MatAssemblyEnd(Adj,MAT_FINAL_ASSEMBLY,ierr)

!

! View adjacency matrix

call MatView(Adj,PETSC_VIEWER_STDOUT_WORLD,ierr)

!

! Once the connectivity matrix has been created the following code will generate the

! renumbering required for the new partition

!

! Create partitions

call MatPartitioningCreate(PETSC_COMM_WORLD,part,ierr)

!

! Set partitions with adjacency matrix

call MatPartitioningSetAdjacency(part,Adj,ierr)

!

! Options for partitioning (PARMETIS)

call MatPartitioningSetFromOptions(part,ierr)

147

!

! Create mapping, partitioning is ready

call MatPartitioningApply(part,is_new,ierr)

!

! View partitioning

call MatPartitioningView(part,PETSC_VIEWER_STDOUT_WORLD,ierr)

!

! View mapping is_new -> mapping between partitions

call ISView(is_new,PETSC_VIEWER_STDOUT_WORLD,ierr)

!

! New global numbering of edges

call ISPartitioningToNumbering(is_new,isg_new,ierr)

!

! View mapping isg_new after renumbering

call ISView(isg_new,PETSC_VIEWER_STDOUT_WORLD,ierr)

!

! Extract indices from IS is_new

! Get the number of indices in the is_new

call ISGetLocalSize(is_new,n_index,ierr)

allocate(is_array(n_index)) ! allocate vector to save indices

do ip = 1, n_index ! initialize vector

is_array(ip) = mrank + ip

end do

!

! Get indices in the is_new

call ISGetIndices(is_new,is_array,i_is,ierr)

allocate(part_index(n_index))

do ip = 1,n_index

part_index(ip) = is_array(i_is + ip)

end do

!

! Restore index in is_new

call ISRestoreIndices(is_new,is_array,i_is,ierr)

!

148

A seguir apresentamos o esquema de comunicação do tipo Round –Robin

Figura 59. Esquema Round-Robin.

Na Fig. (59) apresentamos o esquema de comunicação do tipo Round-Robin, que consiste

em cada processador passar adiante a informação que possui, na ordem especificada. Neste

caso, por exemplo, o processador 0 envia sempre para o processador 1 e recebe sempre do

processador 3, porém as informações que recebe e envia são modificadas a cada passo de co-

municação. A comunicação ocorre n-1 vezes, onde n é o número de processadores envolvidos

na comunicação.

Segue abaixo, o algoritmo do esquema Round-Robin utilizado em um processo de troca

de informações de índices de nós.

! left process neighbors

if (mrank.EQ.0) then

left = nproc - 1

else

left = mrank - 1

end if

!

! right process neighbors

if (mrank.EQ. (nproc - 1)) then

right = 0

else

right = mrank + 1

149

end if

!

do iproc = 1,nproc-1

!

! WE NEED KNOWN WHAT REMOTE PROCESS IS TALKING WITH ME

! remote node

if (imarc1.EQ.0) then

recv_from = right + (iproc - 1)

else

recv_from = recv_from + 1

end if

!

if (recv_from.GE.nproc) then

recv_from = 0

imarc1 = 1

end if

!

! WE NEED KNOWN WHAT REMOTE NODE I AM TALKING TO

! remote node

if (imarc2.EQ.0) then

send_to = left - (iproc - 1)

else

send_to = send_to - 1

end if

!

if (send_to.LT.0) then

send_to = nproc

imarc2 = 1

end if

!

position = 1 ! position in buffer

do ip = 1,nploc

i1 = send_recv_nod(ip)

call MPI_PACK(i1,1,MPI_INTEGER,buffer_nod_out,len(buffer_nod_out),

150

position,MPI_COMM_WORLD,ierr)

end do

!

tag_node_list = 1 ! message identifier

!

! How many nodes I must send to

! nploc = size(send_recv_nod)

!

! Deallocate vector "send_recv_nod" (old entries)

deallocate(send_recv_nod)

!

! Send this size of node list to left process neighbors

call MPI_SEND(nploc,1,MPI_INTEGER,left,tag_node_list,MPI_COMM_WORLD,

ierr)

!

! Receive the size of node list from right process neighbors

call MPI_RECV(nploc,1,MPI_INTEGER,right,tag_node_list,

MPI_COMM_WORLD,status,ierr)

!

! Send node list to left process neighbors

call MPI_SEND(buffer_nod_out,abs(position),MPI_INTEGER,left,

tag_node_list,MPI_COMM_WORLD,ierr)

!

! Cleaning buffer

pos1 = 1

i3 = 0

do in = 1, position

call MPI_UNPACK(buffer_nod_out,abs(len(buffer_nod_out)),pos1,i3,1,

MPI_INTEGER,MPI_COMM_WORLD,ierr)

end do

!

! Test the memory size in buffer_nod

if (nploc.GT.len(buffer_nod_in)) then

write(*,*)'no sufficient memory in buffer - PROGRAM ABORTED'

151

stop

else

! Receive list from right process neighbors

call MPI_RECV(buffer_nod_in,abs(len(buffer_nod_in)),MPI_INTEGER,

right,tag_node_list,MPI_COMM_WORLD,status,ierr)

end if

!

! Message received from process "recv_from"

!

! Allocate vector "send_recv_nod" (new entries)

allocate(send_recv_nod(nploc))

!

! Receiving remote nodes and assembly in "send_recv_nod" vector

position = 1 ! position in buffer

do ip = 1,nploc

call MPI_UNPACK(buffer_nod_in,len(buffer_nod_in),position,i1,1,

MPI_INTEGER,MPI_COMM_WORLD,ierr)

send_recv_nod(ip) = i1

end do

!

! Cleaning buffer

MPI_UNPACK(buffer_nod_in,len(buffer_nod_in),pos1,i1,1,

MPI_INTEGER,MPI_COMM_WORLD,ierr)

!

! Looking for matching nodes (local and remote)

cproc = 0

do ip = 1,nploc

i1 = send_recv_nod(ip) ! remote node

marc = 0

indi1 = 1

do while ((indi1.LE.size(nodlc)).AND.(marc.EQ.0))

i2 = nodlc(indi1) ! local node

indi1 = indi1 + 1

if (i1.EQ.i2) then ! if remote node equal local node

152

icont = icont + 1

id_cnod_aux(icont) = i2

marc = 1

cproc = 1

end if

end do

end do

!

if (cproc.EQ.1) then

ecount = ecount + 1

id_proc_aux(ecount) = recv_from

id_apnt(ecount+1) = icont

end if

!

end do

!

153

Apêndice B: Mapeamento dos Nós

Neste apêndice, discutimos a rotina de mapeamento necessária para minimizar o endere-

çamento indireto na manipulação de vetores e matrizes paralelos. Nesta rotina, o mapeamento

é realizado a partir do processador zero até o processador n-1, numerando os nós em ordem

crescente a partir do índice zero. Desta forma se, por exemplo, a malha contiver 15 nós e for

particionada em três subdomínios, o processador 0 conterá os nós 0, 1, 2, 3 e 4, o processador

1 conterá os nós 5, 6, 7, 8 e 9, e o processador 2 conterá os nós, 10, 11, 12, 13 e 14. Note que

mesmo em FORTRAN, os nós estão numerados a partir do nó zero para facilitar o armazena-

mento nas estruturas vetoriais.

Como cada processador conhece seus nós locais e seus nós que pertencem a mais de um

processador, cada processador conta quantos nós são locais e os renumera, sendo que os nós

que tem cópias remotas são renumerados apenas uma única vez, sempre pelo processador de

menor ranking que possuir este nó. A seguir cada processador envia uma lista de seus nós que

tem cópias remotas, via esquema de comunicação Round-Robin, e ocorre uma atualização da

sua numeração em cada processador. Ao mesmo tempo é criado um mapeamento utilizando a

função AOCreatebasic e AOAppplicationToPetsc, do PETSc para gerenciar o mapeamento

criado. Desta forma é possível identificar cada nó pela sua numeração inicial e após o mape-

amento. A seguir apresentamos a parte da rotina que realiza o mapeamento:

! Test - Creating AO Mapping (local and global numbering mapping)

allocate(locvec(0:locnod-1))

allocate(glbvec(0:locnod-1))

icont = 0

do ip = 1, meshDat%numNodesLocal

if (meshDat%node%nodesIdentRenb(ip,1).NE.-1) then

! start in ident 0 -> PETSC specification

locvec(icont) = meshDat%node%nodesIdentRenb(ip,1)

154

! start in ident 0 -> PETSC specification

glbvec(icont) = meshDat%node%nodesIdent(ip,1) - 1

icont = icont + 1

end if

end do

write(10,*)'number of local nodes is',locnod

do ip = 1,locnod

write(10,'(2i8)')locvec(ip-1),glbvec(ip-1)

end do

meshDat%numNodesRenb = locnod

!

! Create mapping (needed for scalability of matrix and vector)

call AOCreateBasic(PETSC_COMM_WORLD,locnod,glbvec,locvec,ao_mapp,ierr)

call AOView(ao_mapp,PETSC_VIEWER_STDOUT_WORLD,ierr)

!

deallocate(locvec,glbvec)

allocate(locvec(0:1))

locnod = 2

! Test mapping

write(10,*)'PETSC Mapping EDGE'

do ip = 1, meshDat%numEdgesLocal

i1 = meshDat%edge%edgesNodes(ip,1) ! first edge node

i2 = meshDat%edge%edgesNodes(ip,2) ! second edge node

locvec(0) = i1 - 1

locvec(1) = i2 - 1

call AOApplicationToPetsc(ao_mapp,locnod,locvec,ierr)

write(10,*)' global ident ------ local ident'

write(10,*) i1-1, i2-1,' ',locvec(0), locvec(1)

meshDat%edge%edgesNodesRenb(ip,1) = locvec(0)

meshDat%edge%edgesNodesRenb(ip,2) = locvec(1)

end do

!

deallocate(locvec)

allocate(locvec(0:0))

155

locnod = 1

write(10,*)'PETSC Mapping NODE'

do ip = 1, meshDat%numNodesLocal

i1 = meshDat%node%nodesIdentRenb(ip,1) ! flag node

if (i1.EQ.-1) then

locvec(0) = meshDat%node%nodesIdent(ip,1) - 1

call AOApplicationToPetsc(ao_mapp,locnod,locvec,ierr)

meshDat%node%nodesIdentRenb(ip,1) = locvec(0)

write(10,*)' global ident ------ local ident'

write(10,*) meshDat%node%nodesIdent(ip,1),' ',locvec(0)

end if

end do

!

! Mapping edges nodes in their positions in nodes vector

do is = 1,meshDat%numEdgesLocal

i1 = meshDat%edge%edgesNodesRenb(is,1) ! edge node i1 (local numbering)

i2 = meshDat%edge%edgesNodesRenb(is,2) ! edge node i2 (local numbering)

! Find node identification in vector

j1 = 0

idno = -1

do while (i1.NE.idno)

j1 = j1 + 1

idno = meshDat%node%nodesIdentRenb(j1,1) ! node identification

end do

! Find node identification in vector

j2 = 0

idno = -1

do while (i2.NE.idno)

j2 = j2 + 1

idno = meshDat%node%nodesIdentRenb(j2,1) ! node identification

end do

!

meshDat%edge%edgesNodesRenb(is,3) = j1 ! position in nodes vector

meshDat%edge%edgesNodesRenb(is,4) = j2 ! position in nodes vector

156

end do

!

157

Apêndice C: Comparação entre Estru-

turas de Dados por Elementos, Aresta

Pura e Aresta CSR (Compressed Sparse

Row)

Neste trabalho, testamos também e eficiência computacional, em termos de tempo de pro-

cessamento, para diferentes tipos de estrutura de dados. Esta etapa foi importante para defi-

nirmos a forma de armazenamento dos dados necessários seria utilizada ao longo de todo o

trabalho.

Classicamente o Método dos Elementos Finitos necessita de uma estrutura de elementos,

baseada na topologia da malha computacional, e sendo assim, em cada passo iterativo da si-

mulação, um “loop” nos elementos é realizado de forma a se obter a contribuição de cada e-

lemento para os nós que o compõem. Este tipo de estrutura de dados, baseada nos elementos,

necessita que seja obtida uma matriz para cada elemento, da ordem do número de vértices do

elemento. Esta matriz pode ser obtida em cada passo, gerando custo computacional, ou ser ob-

tida no início da simulação e armazenada para cada elemento, gerando custo de armazena-

mento e custo de tempo computacional devido ao acesso aos dados de cada matriz de

elemento. Alternativamente a este tipo de estrutura de dados, foram propostas diferentes for-

mas de armazenar as informações necessárias à solução de equações diferenciais parciais via

métodos numéricos [Löhner, 2001, Löhner e Galle, 2002], entre elas, vamos analisar uma es-

trutura de dados baseada nas arestas, denominada pura, onde o “loop” é realizado diretamente

nas arestas, e uma estrutura baseada nas arestas indiretamente, armazenando os nós e as suas

respectivas conectividades nodais, denominada CSR (Compressed Sparse Row).

Para realizar esta comparação entre estruturas de dados, utilizamos um programa de simu-

lação de escoamentos do tipo convecção-difusão transiente em domínios bidimensionais, que

158

foi desenvolvido no MATLAB. Os custos computacionais foram comparados apenas em ter-

mos de tempo de simulação, medindo o tempo necessário para a realização da simulação

completa. A Fig. (C.1) mostra os tempos obtidos para as diferentes estruturas de dados utili-

zadas.

1 1,01

1,8

0

0,2

0,4

0,6

0,8

1

1,2

1,4

1,6

1,8

Tempo

Adimension

al

Aresta Pura Aresta CSR Elemento

Tipo de Estrutura de Dados

Figura C.1. Comparação entre diferentes tipos de estruturas de dados.

A Fig. (C.1) mostra os resultados obtidos para solução numérica via Método dos Elemen-

tos Finitos do tipo Petrov-Galerkin de uma equação diferencial parcial do tipo convecção-

difusão transiente, com número de Peclet de 100, com termos de estabilização e sobre um

domínio bidimensional retangular discretizado com uma malha computacional com 94 ele-

mentos triangulares, 74 pontos nodais e 167 arestas, como mostra a Fig. (C.2). Os resultados

obtidos mostram a eficiência computacional em termos de tempo de simulação da estrutura de

arestas pura em relação às estruturas de elementos e arestas CSR. O mesmo exemplo foi ana-

lisado em malhas com diferentes graus de refinamento, apresentando resultados semelhantes.

Estes resultados estão de acordo com os resultados obtidos por Löhner, 2001 e Löhner e Gal-

le, 2002, porém como estamos utilizando um domínio bidimensional estes efeitos são minimi-

zados em termos de endereçamento indireto, quando comparados com domínios

tridimensionais discretizados com tetraedros, como mostra Elias et al., 2005. Sendo assim,

nesse trabalho optamos por utilizar uma estrutura de dados baseada puramente nas arestas. Es-

ta estrutura apresentou vantagens em relação à paralelização, pois desta forma, as arestas es-

tão em apenas um subdomínio, e consequentemente, os coeficientes relativos às arestas não

precisam manter cópias remotas em outros processadores, o que não aconteceria se utilizás-

159

semos uma estrutura de arestas CSR, que resulta em “loop” nos nós, pois nós remotos são ne-

cessários, visto que as interfaces são definidas pelos nós.

Figura C.2. Malha utilizada na comparação entre estruturas de dados.

As Figs. (C.3,C.4,C.5) mostram os resultados numéricos obtidos para um problema de

dispersão de poluentes com estabilização do termo de transporte convectivo, com controle da

difusão numérica aplicada, onde uma concentração inicial é prescrita na entrada e fluxo nulo é

prescrito nas fronteiras superior e inferior da Fig. (C.2).

Figura C.3. Dispersão de poluentes, com Peclet 100.

Figura C.4. Dispersão de poluentes, com Peclet 100 mostrando corte ao longo do eixo longi-

tudinal.

160

Figura C.5. Gráfico da concentração ao longo do corte mostrado na Fig. (C.4).

Também realizamos testes com operador de captura de gradientes elevados, comparando

os resultados obtidos com as aproximações de Galerkin, com estabilização caracterizando um

método “upwind”, e com operador definindo a quantidade de estabilização necessária em cada

região do domínio, como mostra a Fig. (C6).

Figura C6. Comparação entre os métodos de Galerkin, Upwind e estabilização com captura de

gradientes acentuados.

Observando a Fig. (C6) podemos verificar que, como esperado, o método de Galerkin a-

presenta oscilações espúrias na presença de fenômenos com convecção dominante, necessi-

tando de termos de estabilização. Porém a inclusão de termos de estabilização sem controle da

161

ordem da aproximação leva a um excesso de difusão, resultando em uma aproximação do tipo

“upwind’, de primeira ordem em todo o domínio. Sendo assim, necessita-se ainda de algum

parâmetro que controle a ordem da aproximação, mantendo segunda ordem em regiões suaves

e baixando para primeira ordem nas regiões de gradientes elevados, evitando assim oscilações

no campo analisado e melhorando a captura dos fenômenos físicos estudados.

162

Apêndice D Resultados Preliminares

Obtidos para Difusão Pura e Convecção-

Difusão Transiente

Neste apêndice mostramos os resultados obtidos para problemas de difusão pura e con-

vecção difusão, ambos transientes, e que serviram para testar os mapeamentos e as trocas de

mensagens necessárias ao longo do processo de simulação. Para tal, analisamos exemplos

simples, mas que necessitam todos os passos do processo de simulação computacional de fe-

nômenos da Dinâmica dos Fluidos Computacional utilizando conceitos de computação parale-

la. Inicialmente estudamos um caso simples de um problema de difusão pura, com termo

transiente, caracterizado por um tubo com condições prescritas em ambos os lados e fluxo nu-

lo nas paredes do tubo. Os resultados podem ser conferidos nas Fig. (D1-D4), com o campo

escalar e a solução ao longo da linha central, longitudinalmente ao tubo. Estes resultados a-

presentam-se de acordo com o esperado, mostrando uma solução linear no regime estacioná-

rio.

Figura D1. Difusão Pura Transiente, passo 1.

163

Figura D2. Difusão Pura Transiente, passo 10.

Figura D3. Difusão Pura Transiente, passo 20.

Figura D4. Difusão Pura Transiente, passo 30.

164

Também foi analisado um problema regido por uma equação de convecção e difusão e em

regime transiente, considerando um número de Peclet (número adimensional definido como a

razão entre os efeitos convectivos e difusivos no escoamento) igual a 2. O número de Peclet

foi considerado baixo para evitar oscilações espúrias no campo escalar, visto que o objetivo

neste momento não é analisar os termos de estabilização espacial, e sim testar a capacidade do

programa desenvolvido em lidar com múltiplos domínios. Os resultados apresentam-se de a-

cordo com o esperado, ou seja, sem oscilações, devido ao baixo número de Peclet e com difu-

são ao longo do domínio. As Figs. (D5-D10) mostram a solução no regime transiente, com o

pulso precorrendo o domínio, devido ao termo de convecção, e difundindo-se devido ao termo

de segunda ordem, e a solução ao longo da linha central na direção longitudinal ao tubo.

Figura D5. Convecção-Difusão Transiente, passo 1.

Figura D6. Convecção-Difusão Transiente, passo 10.

165

Figura D7. Convecção-Difusão Transiente, passo 20.

Figura D8. Convecção-Difusão Transiente, passo 30.

Figura D9. Convecção-Difusão Transiente, passo 40.

166

Figura D10. Convecção-Difusão Transiente, passo 57.