Download PDF
ads:
ESTRUTURAS DE DADOS POR ARESTAS PARA A SIMULAÇÃO PARALELA
DE ESCOAMENTOS INCOMPRESSÍVEIS PELO MÉTODO ESTABILIZADO DE
ELEMENTOS FINITOS.
Renato Nascimento Elias
TESE SUBMETIDA AO CORPO DOCENTE DA COORDENAÇÃO DOS
PROGRAMAS DE PÓS-GRADUAÇÃO DE ENGENHARIA DA UNIVERSIDADE
FEDERAL DO RIO DE JANEIRO COMO PARTE DOS REQUISITOS
NECESSÁRIOS PARA A OBTENÇÃO DO GRAU DE DOUTOR EM CIÊNCIAS
EM ENGENHARIA CIVIL.
Aprovada por:
Prof. Alvaro Luiz Gayoso de Azeredo Coutinho, D.Sc.
Prof. Luiz Landau, D.Sc.
Prof. Nelson Francisco Favilla Ebecken, D.Sc.
Prof. Abimael Fernando Dourado Loula, D.Sc.
Dr. Paulo Augusto Berquo de Sampaio, Ph.D.
RIO DE JANEIRO, RJ – BRASIL
DEZEMBRO DE 2007
ads:
Livros Grátis
http://www.livrosgratis.com.br
Milhares de livros grátis para download.
ii
ELIAS, RENATO NASCIMENTO
Estruturas de Dados por Arestas para a
Simulação Paralela de Escoamentos
Incompressíveis pelo Método Estabilizado de
Elementos Finitos. [Rio de Janeiro] 2007.
XIII, 213, p. 29,7 cm (COPPE/UFRJ,
D.Sc., Engenharia Civil, 2007).
Tese – Universidade Federal do Rio de
Janeiro, COPPE.
1. Elementos Finitos
2. Estruturas de dados por arestas
3. Escoamento incompressível
4. Computação paralela
I. COPPE/UFRJ II. Título (série)
ads:
iii
À minha querida esposa Elizandra e ao
pedacinho de nós dois que está por nascer.
iv
Agradecimentos
Gostaria de oferecer meus sinceros votos de agradecimento a todos àqueles que direta
ou indireta contribuíram para a realização deste trabalho, em especial agradeço a:
− Ao Prof. orientador, e grande amigo, Alvaro Luiz Gayoso de Azeredo Coutinho
que investiu de forma incondicional na realização deste trabalho, por sua
amizade, compreensão, incentivo, dedicação e apoio, por todas suas qualidades
como pesquisador, professor e exemplo de profissional da área acadêmica;
− Ao fiel companheiro e irmão Marcos André Duarte Martins, pelas eternas
discussões sobre os mais diversos assuntos, pela EdgePack, pelo ouvido amigo e
conselheiro nos momentos de dificuldade. Grande parte desse trabalho foi
incentivado e direcionado por este exemplar profissional e grande amigo;
− A Elizandra Elias, minha eterna companheira e defensora, por sua compreensão
nos momentos de ausência, por seu amor, carinho e companheirismo não só nos
momentos prósperos como nos difíceis. Este trabalho é o resultado da nossa
busca por uma vida melhor, agora, ao lado do fruto de nossa união que está por
nascer...
− Aos profissionais e amigos do Núcleo de Atendimento em Computação de Alto
Desempenho (NACAD): Mara (quem nos mantém atualizados dos nossos
compromissos e a salvo de todos os males da burocracia), Albino e Myrian
(gurus de administração, redes e Unix e que aturaram minhas eternas
reclamações sobre computadores, sistemas operacionais, filas, privilégios de
acesso, etc...), Orlando (o cara mais procurado do Programa de Eng. Civil,
quiçá, da COPPE inteira e que tem sempre uma solução na sua bolsa preta
mágica para manter qualquer computador funcionando), D. Helenice (a senhora
mais jovem que conheço e que nos trata como filhos sem nunca esquecer do
cafezinho e do chá na hora certa);
− Aos amigos do “Laboratório-irmão” de Métodos Computacionais em
Engenharia (LAMCE): Landau, José Alves, Mônica e Sérgio;
v
− À Denis de Araujo Filgueiras de Souza pelas frutíferas discussões acerca de
nossos trabalhos, com destaque por sua ajuda no desenvolvimento da solução de
transporte presente neste trabalho;
− À Agência Nacional do Petróleo (ANP/PRH-02) pelo apoio financeiro
indispensável para a realização deste trabalho;
− À Deus pela oportunidade de poder contar com tão bons amigos na realização
deste trabalho.
vi
Resumo da Tese apresentada à COPPE/UFRJ como parte dos requisitos necessários
para a obtenção do grau de Doutor em Ciências (D.Sc.)
ESTRUTURAS DE DADOS POR ARESTAS PARA A SIMULAÇÃO PARALELA
DE ESCOAMENTOS INCOMPRESSÍVEIS PELO MÉTODO ESTABILIZADO DE
ELEMENTOS FINITOS.
Renato Nascimento Elias
Dezembro/2007
Orientador: Alvaro Luiz Gayoso de Azeredo Coutinho
Programa: Engenharia Civil
A solução de problemas envolvendo o acoplamento do escoamento incompressível
com o transporte advectivo-difusivo é de grande interesse por permitir que uma ampla
gama de fenômenos de interesse sejam estudados. Com os recentes avanços nas técnicas
de estabilização o método dos elementos finitos tem ganhado popularidade também na
solução de problemas de dinâmica de fluidos computacional. Além disto, a constante
busca por soluções cada vez mais realistas e ricas em detalhes têm requerido o
desenvolvimento de procedimentos de solução acoplados, tridimensionais e transientes
com natural aumento dos custos computacionais.
Este trabalho apresenta o desenvolvimento do software EdgeCFD, software este
baseado na formulação estabilizada SUPG/PSPG de elementos finitos para a solução de
problemas de escoamento incompressível acoplada, ou não, à solução do transporte
advectivo-difusivo de uma grandeza escalar pelo método SUPG/
YZ
. Diversas
técnicas, algoritmos e métodos foram explorados na busca pela solução eficiente
cobrindo as mais diversas arquiteturas de processamento de alto desempenho.
vii
Abstract of Thesis presented to COPPE/UFRJ as a partial fulfillment of the
requirements for the degree of Doctor of Science (D.Sc.)
EDGE-BASED DATA STRUCTURES FOR THE PARALLEL STABILIZED FINITE
ELEMENT SIMULATION OF INCOMPRESSIBLE FLOWS.
Renato Nascimento Elias
December/2007
Advisor: Alvaro Luiz Gayoso de Azeredo Coutinho
Department: Civil Engineering
The coupled solution for the incompressible fluid flow and advective-diffusive
transport problems is important since it allows that a wide range of practical interest
phenomena are treated. With the recent advances of stabilization techniques the finite
element method has became popular also for the solution of computational fluid
dynamics problems. Moreover, the constant search for more realistic and high detailed
solutions has required the development of coupled, three dimensional and transient
solution procedures which naturally demands for more computational resources.
This work presents the development of the EdgeCFD software which is based on
the SUPG/PSPG finite element formulation for the incompressible fluid flow solution,
coupled, or not, with the advective-diffusive scalar transport using the SUPG/
YZ
method. Several techniques, algorithms and methods were explored for the efficient
solution covering the most different high performance computing architectures.
viii
Índice Analítico
Capítulo 1....................................................................................................................1
Introdução................................................................................................................1
1.1. Motivação......................................................................................................1
1.2. Objetivos.......................................................................................................5
1.3. Estrutura do Texto.......................................................................................6
Capítulo 2....................................................................................................................8
Equações Governantes............................................................................................8
2.1. Escoamento incompressível..........................................................................8
2.2. Transporte advectivo-difusivo....................................................................10
2.3. Fluidos não-Newtonianos...........................................................................11
2.4. Escoamentos gravitacionais.......................................................................14
2.5. Escoamentos turbulentos...........................................................................16
2.5.1. Modelagem de turbulência...................................................................18
2.5.2. Método da simulação dos vórtices de larga escala...............................20
Capítulo 3..................................................................................................................22
Formulação de Elementos Finitos.........................................................................22
3.1. Equações de Navier-Stokes........................................................................22
3.2. Transporte advectivo-difusivo....................................................................25
3.3. Integração no tempo...................................................................................26
3.3.1. Algoritmo preditor multicorretor.........................................................27
3.4. Solução não-linear......................................................................................29
3.4.1. Método de Newton inexato...................................................................30
3.5. Solução linear.............................................................................................35
3.5.1. Método dos resíduos mínimos generalizados (GMRES)......................35
3.6. Pré-condicionamento..................................................................................40
3.6.1. Pré-condicionamento diagonal.............................................................41
3.6.2. Pré-condicionamento bloco diagonal nodal..........................................42
3.7. Controle adaptativo de passo de tempo......................................................43
3.8. Acoplamento escoamento incompressível – transporte.............................47
Capítulo 4..................................................................................................................49
Escoamentos de superfície-livre............................................................................49
4.1. Tratamento da superfície-livre...................................................................49
ix
4.1.1. Método do volume-de-fluido.................................................................52
4.1.2. Método das curvas de nível..................................................................53
4.2. Correção de massa......................................................................................56
4.3. Integração de volume para as fases...........................................................57
4.4. Cálculo de funções-distância......................................................................57
4.5. Desativação Dinâmica em Paralelo............................................................59
Capítulo 5..................................................................................................................62
Estruturas de dados por arestas...........................................................................62
5.1. Grafos, malhas e matrizes esparsas...........................................................63
5.2. Incidências de elemento e aresta................................................................65
5.3. Construção das matrizes de aresta............................................................67
5.3.1. Elementos tetraédricos lineares...........................................................69
5.4. Produto matriz-vetor por aresta.................................................................70
5.5. Custo computacional...................................................................................74
5.6. Técnicas de agrupamento e ordenação de dados........................................76
Capítulo 6..................................................................................................................81
Técnicas computacionais.......................................................................................81
6.1. Computação paralela..................................................................................81
6.1.1. Troca de mensagens.............................................................................83
6.1.2. Multitarefa...........................................................................................89
6.1.3. Paralelismo Misto ou Híbrido..............................................................92
6.2. Vetorização e pipelining .............................................................................94
6.3. Pré-processamento......................................................................................97
6.4. Pós-processamento......................................................................................98
Capítulo 7................................................................................................................103
Conclusões...........................................................................................................103
Trabalhos futuros................................................................................................106
Referências..............................................................................................................107
Apêndices................................................................................................................117
Apêndice A...............................................................................................................118
Escoamento de fluidos visco-plásticos.................................................................118
Apêndice B...............................................................................................................136
Tratamento de dados por arestas........................................................................136
Apêndice C...............................................................................................................167
Cálculo de funções-distância para a simulação utilizando level-sets.................167
x
Apêndice D..............................................................................................................184
Escoamento de superfície-livre............................................................................184
xi
Índice de Figuras
Figura 2.1 - Representação da reologia dos fluidos considerados............................13
Figura 2.2 – Modelo de Bingham “bi-viscoso”..........................................................14
Figura 2.3 – Transporte de corante em escoamento de transição entre os regimes
laminar e turbulento (a) Ilustração de vórtices de diferentes escalas (b)
Flutuação de velocidade (Figura obtida em http://dspace.mit.edu).................18
Figura 3.1 - Pré-condicionador Bloco-Diagonal Nodal.............................................43
Figura 3.2 – Seleção do passo de tempo como um problema de controle retro-
alimentado........................................................................................................44
Figura 4.1 – Formas de tratamento da superfície-livre (a) acompanhamento da
interface (b) captura da interface.....................................................................50
Figura 4.2 – Função de marcação para o método do volume-de-fluidos (VOF).......52
Figura 4.3 – Difusão da função degrau do método VOF..........................................53
Figura 4.4 – Função de marcação para o método de curvas de nível (LS)...............54
Figura 4.5 – Desativação dinâmica em paralelo utilizando dois processadores......61
Figura 5.1 – Grafo de uma matriz 5 × 5 e sua matriz correspondente....................64
Figura 5.2 – Relação entre malha, grafos e matrizes...............................................65
Figura 5.3 – Construção das incidências de elemento e aresta................................66
Figura 5.4 – Contribuições dos elementos para a matriz da aresta a
1
....................68
Figura 5.5 – Desmembramento da matriz de rigidez do elemento tetraédrico em
suas matrizes de aresta correspondentes........................................................69
Figura 5.6 – Produto matriz-vetor aresta-por-aresta...............................................72
Figura 5.7 – Tetraedro linear - Custos computacionais para as estruturas aresta-
por-aresta (EDS), elemento-por-elemento (EBE) e compressed sparse row
(CSR).................................................................................................................75
Figura 5.8 – Exemplo de ordenações........................................................................78
Figura 5.9 – Técnicas de reordenação realizadas pela EdgePack............................80
Figura 6.1 – Paralelismo por troca de mensagens aplicado a malha computacional
..........................................................................................................................84
Figura 6.2 – Conceito de volume de comunicação e balanço de carga aplicados a
uma malha não-estruturada............................................................................85
Figura 6.3 – Malha original......................................................................................87
Figura 6.4 – Malha particionada mantendo a numeração original.........................88
xii
Figura 6.5 – Malha particionada e renumerada......................................................88
Figura 6.6 – Modelo de paralelismo por multitarefas ou memória compartilhada.90
Figura 6.7 – Laço com dependência de memória......................................................91
Figura 6.8 – Laço sem dependência de memória......................................................92
Figura 6.9 – Tratamento de dados para paralelismo...............................................93
Figura 6.10 – Produto matriz-vetor híbrido.............................................................94
Figura 6.11 – Operação de adição de dois vetores (a) escalar (b) vetorial...............95
Figura 6.12 – Instruções sendo efetuadas sem o uso de pipeline.............................96
Figura 6.13 – Instruções sendo efetuadas com o uso de pipeline.............................96
Figura 6.14 – Visualização remota em paralelo.....................................................102
xiii
Índice de Algorítmos
Algoritmo 3.1 – Integração no tempo das equações de Navier-Stokes
incompressíveis.................................................................................................27
Algoritmo 3.2 – Método de Newton inexato..............................................................34
Algoritmo 3.3 – Método dos resíduos mínimos generalizados (GMRES).................38
Algoritmo 3.4 – Escolha adaptativa de passo de tempo por controlador PID..........46
Algoritmo 3.5 – Laço principal do software desenvolvido mostrando o acoplamento
entre os procedimentos de solução para o escoamento e o transporte............48
Algoritmo 5.1 – Construção das matrizes de aresta.................................................70
Algoritmo 5.2 – Produto matriz-vetor aresta por aresta..........................................71
Algoritmo 5.3 – Produto matriz-vetor aresta por aresta no esquema Gizjen..........73
Algoritmo 6.1 – Pré-processamento..........................................................................97
1
Capítulo 1
Introdução
1.1. Motivação
A simulação computacional de escoamentos incompressíveis tem se tornado de suma
importância no auxílio à compreensão de diversos fenômenos de interesse econômico,
industrial, acadêmico e do nosso dia-a-dia. De forma bastante sucinta, os escoamentos
podem ser classificados como incompressíveis quando a variação da massa específica
do fluido é desprezível [20]. Contudo deve-se observar que, na prática, nenhum fluido é
de fato incompressível, mas para fins de engenharia tal suposição permite que uma
ampla gama de fenômenos possa ser estudada dentro da precisão desejada [5].
As equações que descrevem o movimento dos fluidos, batizadas em homenagem
aos trabalhos do engenheiro e físico francês Claude Louis Marie Henri NAVIER (1785-
1836) [56] e do matemático e físico irlandês George Gabriel STOKES (1819-1903)
[57], estabelecem que as mudanças na quantidade de movimento de um volume de
fluido infinitesimal estão relacionadas à soma da dissipação de forças viscosas (tensões
cisalhantes), da variação de pressão (tensões normais), da força gravitacional, dentre
outras forças que atuam no interior do mesmo. O que, de fato, constitui-se como uma
aplicação da segunda lei de Newton.
A abrangência e versatilidade desse conjunto de equações é tão vasto que elas são
utilizadas para descrever diversos fenômenos físicos que vão desde modelagem
2
climática, correntes oceânicas, escoamento no interior de tubos, escoamento ao redor de
aerofólios e asas até o movimento de estrelas nas galáxias. Uma conseqüência direta de
tamanha versatilidade se reflete em inúmeras aplicações práticas, tais como: previsões
meteorológicas; projeto de carros, aviões, navios; estudo do escoamento de sangue em
veias e artérias; projeto de usinas de energia; análise dos efeitos da poluição dentre
outras inúmeras aplicações que poderiam ser extensamente listadas.
As equações de Navier-Stokes também são de grande interesse do ponto de vista
matemático. Talvez por sua versatilidade e aplicabilidade exista o consenso entre os
matemáticos de que uma solução analítica geral e em três dimensões simplesmente não
exista, ou, se existir, esta será finita e suave, ou seja, isenta de singularidades e
descontinuidades. Este desafio é conhecido como prova da existência e suavidade das
equações de Navier-Stokes e constitui um dos sete problemas mais importantes da
matemática atual, sendo oferecido um prêmio de um milhão de dólares, pelo Clay
Mathematics Institute
1
, para quem apresentar uma solução ou contra-prova.
Uma conseqüência da possível inexistência de uma solução analítica para as
equações de Navier-Stokes, e da necessidade de soluções para os problemas e
aplicações deste conjunto de equações, foi o impulso para o desenvolvimento de
métodos numéricos capazes de prover soluções aproximadas. Tais métodos numéricos
freqüentemente lançam mão do uso intensivo de computadores para viabilizar que
soluções cada vez mais precisas sejam obtidas. A este ramo da mecânica de fluidos deu-
se a designação de dinâmica de fluidos computacional, do inglês, Computational Fluid
Dynamics ou simplesmente CFD [82]. A dinâmica de fluidos computacional atualmente
se constitui uma área bastante prolífica e que torna possível o estudo de fenômenos que
antes eram impossíveis ou difíceis de serem analisados. Conseqüentemente, o
desenvolvimento deste segmento tem permitido a substituição de estudos experimentais,
muitas vezes dispendiosos ou impossíveis de serem realizados, por simulações
numéricas. A dinâmica de fluidos computacional é uma área intrinsecamente
1
O Clay Mathematics Institute é uma fundação privada, sem fins lucrativos e localizada em Cambridge,
Massachusetts, EUA. Este instituto, fundado em 1998, destina-se à divulgação e disseminação de
conhecimentos da matemática oferecendo prêmios e bolsas para descobertas inovadoras no ramo da
matemática .
3
interdisciplinar que requer o conhecimento técnico de engenheiros, as observações e
modelos das ciências naturais, os métodos numéricos de matemáticos e técnicas e
algoritmos eficientes provenientes da ciência da computação [27].
O método dos elementos finitos constitui uma ferramenta numérica poderosa que
possui como principais atributos a sua habilidade intrínseca em lidar com geometrias
complexas, o tratamento consistente de condições de contorno e a possibilidade de ser
programado de uma forma flexível e genérica [10]. Segundo FELLIPA [17], sua
aplicação em mecânica dos fluidos data dos anos 50 na indústria aerospacial. Durante
muito tempo este método ficou restrito ao tratamento de problemas de mecânica dos
sólidos, onde o tratamento do contínuo só se faz necessário em casos bastante
específicos [89].
A aproximação de elementos finitos padrão baseia-se na formulação de Galerkin do
método dos resíduos ponderados. Esta formulação se mostra bastante eficaz para a
solução de problemas elípticos e parabólicos, entretanto, para soluções hiperbólicas o
método de Galerkin apresenta soluções poluídas por oscilações e que, para serem
tratadas, necessitam de refinamento severo da malha nas regiões onde estas oscilações
se apresentam ou do tempo, o que inviabiliza o uso do método [10]. Tal dificuldade
motivou o aparecimento das técnicas de estabilização que eliminam tais requisitos. O
método de Petrov-Galerkin propõe o controle destas oscilações através da modificação
dos espaços das funções de peso. Em 1982, BROOKS e HUGHES [6] propõem o uso
de funções de peso descontínuas com a formulação Streamline Upwind Petrov-Galerkin
(SUPG) abrindo um novo horizonte na utilização do método dos elementos finitos em
mecânica dos fluidos. HUGHES et al em [33] propõem uma formulação que permite o
uso de funções de interpolação de mesma ordem para problemas de escoamentos de
Stokes. Em [9] de SAMPAIO apresenta uma formulação Petrov-Galerkin generalizada
para escoamentos incompressíveis, utilizando também funções de interpolação de
mesma ordem, porém nesta o parâmetro de estabilização está relacionado à escolha do
passo de tempo local. Mais tarde, em [80], TEZDUYAR estende a formulação proposta
por HUGHES et al [33] para o tratamento de escoamentos incompressíveis no método
conhecido por Streamline Upwind Petrov-Galerkin/Pressure Stabilizing Petrov-
Galerkin (SUPG/PSPG).
4
Um fator comum a todos os métodos que se dedicam a solução de escoamentos tem
sido a necessidade do acoplamento e tratamento de problemas envolvendo múltiplas
físicas, além disso, a solução está geralmente relacionada à qualidade e grau de
refinamento da discretização. Conseqüentemente, o desenvolvimento de algoritmos e
métodos que utilizem eficientemente os recursos de processamento de alto desempenho
se torna imperativo para viabilizar soluções com níveis de detalhamento cada vez mais
realistas.
A estrutura de dados por arestas tem provado, em diversos trabalhos, sua habilidade
em acelerar processos iterativos de solução. VENKATAKRISHNAN E MAVRIPLIS
[86] apresentaram os primeiros testes de desempenho utilizando estruturas de dados por
arestas para a solução implícita das equações de Navier-Stokes compressíveis utilizando
o método dos volumes finitos. CATABRIGA [7], mais tarde, empregou com sucesso
esta estrutura de dados na solução das equações de Euler em 2D utilizando triângulos
lineares e formulação SUPG com operador de captura de descontinuidade. MARTINS
[49] aplicou a estrutura de dados por arestas para elementos tetraédricos e hexaédricos
na solução de problemas de plasticidade. Em [41] KRAFT et al empregam a estrutura
de dados por arestas na solução das equações de Navier-Stokes incompressíveis. Neste
trabalho, os autores adotam uma formulação de elementos finitos segregada que utiliza
passo de tempo local para simetrização do sistema de equações resultante e para
utilização de funções de interpolação de mesma ordem. Recentemente RIBEIRO e
COUTINHO [62] apresentaram diversas comparações com estruturas de dados
convencionais e mostraram que a estrutura de dados por arestas possui desempenho
comparável à estrutura Compressed Sparse Row (CSR) para a maioria dos problemas e
chega a ser superior para problemas envolvendo quatro ou mais graus de liberdade por
nó.
Neste trabalho desenvolveu-se um software destinado a solução de problemas de
grande porte envolvendo escoamentos incompressíveis acoplados, ou não, ao transporte
de uma grandeza escalar qualquer. Este software, registrado sobre o nome de
5
EdgeCFD
2
, emprega o método dos elementos finitos estabilizado como método
numérico de discretização. Na busca por um software eficiente e robusto, capaz de
resolver problemas de interesse prático, foram associados algoritmos e métodos tal
como: o método de Newton-inexato [14], [15], controle adaptativo de passo de tempo
por controlador PID [83], [84], desativação dinâmica para problemas de transporte [44],
métodos de correção de massa, dentre outros. Além disso, técnicas de processamento de
alto desempenho, tais como: paralelismo híbrido, vetorização e pipelining
implementados em linguagem Fortran90 atribuíram flexibilidade no uso do software nas
mais diversas plataformas.
1.2. Objetivos
Este trabalho se desenvolveu com foco multidisciplinar e com a preocupação contínua
de sintetizar, na forma de um software, diversas técnicas que representam o estado da
arte em simulação computacional de escoamentos incompressíveis acoplados, ou não,
ao transporte de grandezas escalares (temperatura, concentração, marcador, etc...)
utilizando métodos de elementos finitos estabilizados. Algumas dessas técnicas são
completamente novas enquanto outras ainda não haviam sido reunidas e aplicadas ao
tratamento de escoamentos incompressíveis e/ou transporte advectivo-difusivo. Atenção
especial foi dada aos aspectos de implementação e otimização do software, resultante
desta pesquisa, para que este pudesse ser efetivamente utilizado na análise de casos reais
que requeiram tanto recurso computacional quanto for necessário para se atingir a
qualidade de resultado esperado. Para este fim, técnicas computacionais como:
paralelismo híbrido, reordenação de dados, pós-processamento e visualização de alto
desempenho, dentre outras, foram reunidas, avaliadas e consideradas. As escolhas
seguiram as tendências atuais de outros grupos de pesquisa além de considerar as
experiências pregressas de pesquisadores que desenvolveram algumas das técnicas que
foram adotadas. Aspectos de qualidade e engenharia de software, tais como:
documentação, controle de versões e adoção da linguagem Fortran como padrão de
2
Edge é uma referência à estrutura de dados utilizada pelo software e CFD por se tratar de uma
implementação tão genérica quanto o termo CFD representa e que sugere que no futuro haverá extensões
do EdgeCFD para o tratamento de outros tipos de escoamentos além do incompressível.
6
desenvolvimento também foram considerados para que objetivos como: portabilidade,
flexibilidade, facilidade de uso e extensões futuras fossem alcançados.
1.3. Estrutura do Texto
A estrutura deste texto se apresenta dividida em duas partes. Na primeira parte os
assuntos tratados neste trabalho são sucintamente discutidos. A segunda parte,
apresentada sobre a forma de Apêndices, contém os quatro artigos mais relevantes e
publicados em periódicos indexados representando marcos de aceitação das hipóteses
assumidas ao longo desta pesquisa. Estes artigos estão anexados em seus formatos
originais em inglês. É importante observar que por se tratarem de trabalhos científicos
completos, estes artigos reúnem as referências mais relevantes assim como os
problemas empregados na validação das implementações que estão sendo propostas e,
desta forma, problemas de validação e testes foram omitidos da primeira parte e devem
ser consultados a partir dos correspondentes Apêndices. Em contrapartida, os tópicos
que não são tratados nestes artigos se encontram detalhados ao longo da primeira parte.
Na primeira parte a estrutura do texto consiste em: no Capítulo 2 são apresentadas
as equações governantes para o a solução de problemas incompressíveis e transporte de
uma grandeza escalar juntamente com as considerações requeridas para o tratamento de
escoamentos de fluidos não-Newtonianos, escoamentos de correntes gravitacionais
regidos pela hipótese de Boussinesq e as premissas utilizadas no tratamento de
escoamentos turbulento pelo método da simulação dos vórtices de larga escala (LES).
No Capítulo 3 são apresentadas as formulações estabilizadas SUPG/PSPG de elementos
finitos para o tratamento do escoamentos incompressíveis utilizando funções de
interpolação de mesma ordem e a formulação SUPG provida do operador
YZ
de
captura de descontinuidades. Neste mesmo capítulo são apresentados: o método de
integração no tempo pelo algoritmo preditor-multicorretor, o método de Newton-inexato
para a solução não-linear das equações de Navier-Stokes, o método dos resíduos
mínimos generalizados (GMRES), utilizado para a solução do problema linear e a
técnica de controle de adaptação de passo de tempo com o uso de controlador do tipo
proporcional-integral-derivativo (PID). Ao final deste capítulo é apresentado o esquema
de acoplamento entre as soluções de escoamento e transporte. No Capítulo 4 a solução
7
de problemas de superfície-livre, pelos métodos do volume-de-fluido e curvas de nível,
é apresentado dentro do contexto de suas implementações neste trabalho. Neste mesmo
capítulo são descritos: os algoritmos de integração de volume e correção de quantidade
de matéria, adotados na melhoria da qualidade da solução, e o algoritmo de desativação
dinâmica, utilizado para acelerar a solução de transporte da função de marcação. O
Capítulo 5 é devotado à descrição da estrutura de dados por arestas a partir de seus
conceitos básicos, construção da matriz de rigidez do problema e implementação do
produto matriz-vetor. Ao final deste capítulo são apresentados os custos computacionais
e algumas das técnicas de ordenação e agrupamento de dados empregados para fins de
otimização. No Capítulo 6 são apresentadas as principais técnicas computacionais
empregadas na implementação do software resultado deste trabalho. Tais técnicas
incluem: paralelismo por troca de mensagens, paralelismo por multitarefas, vetorização
e pipelining. Neste capítulo também são descritos o pré-processador utilizado na
preparação dos dados e as técnicas de pós-processamento e visualização de grandes
volumes de dados. Finalmente, no Capítulo 7 são feitas as conclusões e sugestões para
trabalhos futuros.
8
Capítulo 2
Equações Governantes
Nesta seção são apresentadas as equações de Navier-Stokes desenvolvidas para o
tratamento de escoamentos incompressíveis e seu caso particular aplicado ao transporte
advectivo-difusivo de uma grandeza escalar qualquer. Será tratado também o
acoplamento entre o escoamento e o transporte de um escalar aplicado à solução de
algumas classes de problemas, tais como: transporte e dispersão de poluentes e energia e
tratamento de superfície-livre.
2.1. Escoamento incompressível
Considera-se um fluido viscoso, incompressível, ocupando o domínio
ndim

¡
, de
contorno
, onde
ndim
é o número de dimensões no espaço. A velocidade e a pressão,
u e p, são governadas, em regime transiente, pelos seguintes balanços de quantidade de
movimento e massa respectivamente
0,
em
f
t
t





u
uuf
(2.1)
0,
em
f
t



u0∇
(2.2)
9
onde
é a massa específica do fluido, assumida constante para escoamento
incompressível,
t
é o tempo,
f
o vetor de forças externas ou de forças de corpo e
é
o tensor de tensões
,
pp

uIT
(2.3)
onde
p
representa a pressão hidrostática,
I
é o tensor identidade e
T
é o tensor de
tensão desviadora ou de tensões cisalhantes.
No caso de um fluido Newtoniano, o tensor de tensão desviadora relaciona-se à taxa
de deformação do fluido através de uma constante de proporcionalidade que representa
a difusão da quantidade de movimento experimentada pelo fluido ao escoar [59], assim,
a definição da tensão desviadora é expressa por
2
Tu
ε
(2.4)
onde
é a constante de proporcionalidade conhecida como viscosidade cinemática e
ε
é o tensor taxa de deformações ou
1
2
T
uuuε
(2.5)
As equações (2.1) e (2.2) são conhecidas como equações de Navier-Stokes para
escoamento incompressível. Estas estão sujeitas às seguintes condições de contorno de
Dirichlet e Neumann
,
,
emonde1,...,
emonde1,...,
gd
hd
dndim
dndim


ug
nh
(2.6)
onde
g
e
h
são subconjuntos do contorno Γ e a condição inicial,
0
,0
uxu
(2.7)
10
De um ponto de vista numérico, o sistema de equações formado pelas equações
(2.1) e (2.2) apresenta algumas dificuldades. Primeiramente a equação de quantidade de
movimento é não-linear devido à presença do termo de convecção. Além disso, esta
equação possui seu comportamento dominado pelo termo de difusão, em escoamentos
altamente viscosos ou em baixa velocidade, ou pelo termo de convecção, em
escoamentos pouco viscosos ou em alta velocidade.
Neste trabalho será aplicada a formulação estabilizada de elementos finitos,
proposta inicialmente por TEZDUYAR [80] para as variáveis primitivas velocidade e
pressão.
2.2. Transporte advectivo-difusivo
O transporte advectivo-difusivo pode ser entendido como um caso particular das
equações de Navier-Stokes onde uma grandeza tal como, temperatura, fração mássica,
traçador ou outra função escalar qualquer é deslocada por um campo de velocidade
conhecido ou computado [26]. O princípio de conservação da grandeza de interesse
pode ser declarado através da seguinte equação diferencial parcial,
0
t


uK∇∇
(2.8)
onde
é a grandeza escalar a ser transportada,
K
é o tensor difusividade e
u
é o
campo de velocidades fornecido ou calculado a partir das equações de Navier-Stokes
discutida na seção anterior. O termo
u
∇
representa a parcela advectiva
3
responsável
pelo transporte em larga-escala da grandeza
enquanto o termo

K
∇
representa a parcela difusiva
4
e que se refere ao espalhamento sofrido pela grandeza. No
3
Em alguns casos o termo convecção é utilizado no lugar de advecção, entretanto, cabe ressaltar que o
termo convecção refere-se ao transporte de uma forma mais abrangente e seria mais propriamente
empregado para se referir à soma dos efeitos advectivos e difusivos [34].
4
O termo difusão provém do latim de diffundere e que significa espalhar, derramar [34].
11
contexto do transporte da concentração de uma espécie química, o termo difusão refere-
se ao movimento espontâneo das moléculas de regiões de mais alta concentração para
regiões menos concentradas [34].
2.3. Fluidos não-Newtonianos
Os fluidos que não obedecem a equação (2.4) são denominados fluidos não-
Newtonianos. A principal característica dos fluidos não-Newtonianos, reside na
dependência da viscosidade com outros parâmetros do escoamento, como a taxa de
cisalhamento ou até mesmo da história de deformação do fluido como se verifica em
alguns fluidos denominados tixotrópicos e reopéticos [20]. Os fluidos não-Newtonianos
adotados neste trabalho foram tratados através dos modelos de Lei de Potência (Power
Law) e Bingham, por representarem bem uma ampla gama de materiais comuns na
indústria, tal como: tintas, lamas de perfuração de poços, pastas e resinas poliméricas
[20].
Nas definições a seguir a equação (2.4) foi modificada de forma a refletir o caráter
não-Newtoniano do fluido tornando sua definição mais ampla, desta forma:
2

Tu
ε
&
(2.9)
onde
&
é a taxa de cisalhamento do fluido e

&
passa a ser chamada de viscosidade
aparente de fluidos não-Newtonianos.
A partir da equação (2.9) as leis constitutivas aqui consideradas passam a ser
definidas por:
− Fluido Newtoniano:
const


&
(2.10)
− Fluido de Lei de Potência (Power Law):
12
1
00
1
000
n
n
kse
kse



&&&
&
&&&
(2.11)
onde
k
denota o índice de consistência,
0
é uma viscosidade nominal,
n
o expoente
da Lei de Potência e
0
&
é a taxa de cisalhamento limite.
− Fluido de Bingham clássico:
0
Y


σ
&
&
(2.12)
na expressão acima
Y
σ
é a tensão limite para fluidos de Bingham e
0
é a viscosidade
plástica.
A Figura 2.1 mostra a variação da tensão cisalhante experimentada por alguns
fluidos em função da taxa de cisalhamento. Nota-se a relação linear existente nos
fluidos Newtonianos, onde a inclinação da reta representa a viscosidade. Os fluidos que
seguem a Lei de Potência apresentam dois comportamentos distintos. Nos casos em que
a viscosidade diminui com o aumento da taxa de deformação
1
n
o fluido é
chamado de pseudo-plástico (torna-se delgado com o aumento das tensões tangenciais).
A maioria dos fluidos não-Newtonianos enquadra-se neste grupo e os exemplos incluem
soluções poliméricas e suspensões coloidais. Os fluidos nos quais a viscosidade
aumenta com o aumento da taxa de deformação
1
n
são chamados de dilatantes
(tornam-se mais espessos com o aumento das tensões tangenciais) [20]. As suspensões
de amido e areia são exemplos de fluidos dilatantes.
Um “fluido” que se comporta como um sólido até que uma tensão limite seja
excedida, e subseqüentemente apresenta uma relação linear entre a tensão e a taxa de
cisalhamento é denominado de fluido de Bingham.
13
Lei de Potência (pseudo plástico) n < 1
Lei de Potência (dilatante) n > 1
Newtoniano
Taxa de Cisalhamento ( γ )
Tensão Cisalhante ( T )
Bingham
Clássico
.
Figura 2.1 - Representação da reologia dos fluidos considerados
Cuidados especiais devem ser tomados quando se deseja simular o escoamento de
fluidos de Bingham conforme definido pela equação (2.12). Dentro das regiões do
domínio onde a tensão desviadora apresenta-se abaixo da tensão limite
Y
σ
, o material
comporta-se como um corpo rígido e a viscosidade tende ao infinito, i.e.


&
.
Para contornar tal dificuldade, o fluido de Bingham será tratado por um modelo “Bi-
viscoso” conforme sugerido em [23] e [1], desta forma:
0
0
0
YY
r
Y
r
r
se
se






σσ
σ
&
&
&
&
(2.13)
onde
r
é uma viscosidade Newtoniana normalmente escolhida para ser algumas
ordens de magnitude maior que a viscosidade plástica, com o objetivo de reproduzir as
características de um fluido de Bingham. A Figura 2.2 abaixo apresenta o modelo de
Bingham “bi-viscoso” definido pela equação (2.13).
14
Taxa de Cisalhamento ( γ )
Tensão Cisalhante (
T
)
Bingham
"bi-viscoso"
0
µ
µ
r
µ
r
− µ
0
σ
Y
Y
σ
.
Figura 2.2 – Modelo de Bingham “bi-viscoso”.
2.4. Escoamentos gravitacionais
Na dinâmica de fluidos a aproximação de Boussinesq é freqüentemente empregada na
modelagem de escoamentos gravitacionais ou sob o efeito de flutuação. A hipótese de
Boussinesq assume que as diferenças de massa específica são pequenas o suficiente para
serem desprezadas a não ser naqueles termos onde há o efeito da aceleração da
gravidade. Conseqüentemente, a aproximação de Boussinesq considera que os efeitos
inerciais são pequenos quando comparados com as diferenças de peso específico
experimentada pelo fluido.
Escoamentos deste tipo são bastante comuns. Na natureza, correntes marinhas,
deslocamento de frentes atmosféricas, estratificação de sedimentos, formação da brisa
do mar, deslocamento de nuvens piroclásticas e avalances são alguns exemplos. Na
indústria os escoamentos gravitacionais podem se manifestar na forma da dispersão de
nuvens de gases densos, descarga de gases de chaminés, aproveitamento de ventilação
natural, projetos de sistemas de refrigeração/aquecimentos, dentre outros inúmeros
15
exemplos. A aproximação de Boussinesq permite que vários desses fenômenos sejam
calculados com uma boa precisão simplificando a formulação física e matemática [27].
A principal simplificação surge quando a diferença entre as massa específicas é inferior
a 5%. Nestes casos, a aproximação de Boussinesq considera para a parte inercial
somente a massa específica de referência, associando o efeito da diferença de massa
específica ao termo gravitacional. A introdução da aproximação de Boussinesq à
Equação (2.1) se dá a partir da modificação do termo associado às forças de corpo, ou
seja,

fg
%
(2.14)
onde
g
é o vetor aceleração da gravidade,
é a massa específica de referência e
%
é a
massa específica ligeiramente superior ou inferior à de referência.
Pequenas diferenças de massa específica estão freqüentemente associadas a
diferenças de temperatura sofridas pelo fluido. Isto provoca o surgimento de forças de
flutuação que podem ser modeladas a partir da hipótese de Boussinesq. Neste caso, a
Equação (2.14) torna-se
()
T
fg
(2.15)
onde
()
T
é a relação de dependência da massa específica do fluido com a temperatura
T
, definida por
()1(,)
TTtT




x
(2.16)
onde
x
é o vetor posição,
T
é uma temperatura de referência e
é o parâmetro
conhecido por coeficiente de expansão característico do material e definido por
1
T

(2.17)
16
Note que existe um acoplamento entre o campo de velocidades gerado pela diferença
nas massas específicas e obtido a partir da solução das Equações (2.1) e (2.2) e a
distribuição de temperatura calculada a partir da Equação (2.8). Conseqüentemente, o
único fator responsável pelo escoamento é o gradiente de temperatura. Este tipo de
escoamento é conhecido como convecção-natural [27].
2.5. Escoamentos turbulentos
A maior parte dos escoamentos presentes na natureza, assim como aqueles de interesse
na engenharia são, de fato, turbulentos. Daí vem a importância em simulá-los
numericamente. Embora escoamentos turbulentos possam ser facilmente observados no
dia-a-dia, não existe um consenso acerca da definição de turbulência [65]. Entretanto,
diversas características são aceitas para descrever um escoamento turbulento, tais como
[65], [18]:
− Instabilidade: Escoamentos turbulentos são transientes e podem se manifestar
de maneira periódica ou caótica dependendo do número de Reynolds do
escoamento;
− Vorticidade: Os escoamentos turbulentos ocorrem com a formação de
redemoinhos (vórtices) de diversos tamanhos. Os maiores podem ser da ordem
de magnitude do domínio do problema enquanto os menores podem ser tão
pequenos quanto poucos milímetros. Os vórtices maiores carregam a maior parte
da energia do escoamento. Esta energia é constantemente transferida para os
vórtices de menores escalas. Este processo de transferência de energia é
conhecido por cascata de energia. Evidências comprovam que energia também é
transferida no sentido oposto; dos vórtices menores para os maiores; porém, em
escala bem menos reduzida. Este fenômeno é conhecido como cascata inversa
(do inglês, inverse cascade) ou simplesmente backscattering [27];
− Tridimensionalidade: Escoamentos turbulentos são tridimensionais e
conseqüentemente com vórtices interagindo nas três direções;
17
− Diferentes escalas: As flutuações das grandezas presentes em escoamentos
turbulentos ocorrem em uma ampla gama de escalas de comprimento e tempo.
Desta característica surge a necessidade de introduzir a noção de escala para se
descrever quantitativamente o fenômeno de turbulência [65].
− Difusividade: A taxa de difusão, ou seja, a velocidade com a qual, inércia,
massa e calor são espalhados é substancialmente maior em escoamentos
turbulentos. A razão para isto é a interação existente entre os vórtices de
diferentes tamanhos que contribuem no espalhamento dessas grandezas [27].
− Não-linearidade: Pequenas perturbações em um escoamento turbulento podem
ser amplificadas ao ponto de se tornarem grandes perturbações as quais podem
permanecer estáveis por um longo período de tempo [27].
A Figura 2.3a ilustra o transporte de um corante em uma corrente de transição entre
os regimes laminar e turbulento. Na região laminar o fluido segue paralelamente à linha
de corrente conforme ilustrado pela linha reta do corante. Nesta região as forças
viscosas são preponderantes em relação às forças de inércia [27]. Na região turbulenta
ocorre a formação de vórtices de diferentes escalas (tamanhos). A partir deste ponto o
caminho percorrido pelo corante passa a ser influenciado pelas linhas de corrente; no
sentido preponderante do escoamento; juntamente com seus vórtices. Os vórtices
maiores geralmente possuem mais energia e são responsáveis pela parcela mais efetiva
do transporte enquanto os vórtices menores, e menos energéticos, misturam e espalham
o filamento de corante [18]. As interações entre os vórtices de diferentes escalas geram
flutuações no campo de velocidades conforme sugerido na Figura 2.3b para o ponto A
da Figura 2.3a, entretanto, segundo SAGAUT [65], sabe-se que os escoamentos
turbulentos apresentam estruturas coerentes e com certo grau de organização que se
manifesta na formação de estruturas e aglomerados de vórtices com aproximadamente o
mesmo tamanho e forma e que se repetem com certa periodicidade.
18
(a)
(b)
Figura 2.3 – Transporte de corante em escoamento de transição entre os regimes laminar
e turbulento (a) Ilustração de vórtices de diferentes escalas (b) Flutuação de velocidade
(Figura obtida em http://dspace.mit.edu)
2.5.1. Modelagem de turbulência
As equações de Navier-Stokes transientes e em três dimensões descrevem
matematicamente os escoamentos de qualquer natureza, sejam eles turbulentos ou não
[27]. Conseqüentemente, poder-se-ia modelar escoamentos turbulentos através da
solução numérica das equações de Navier-Stokes em uma malha computacional tão fina
quanto fosse o menor vórtice formado. De fato, esta possibilidade existe e é chamada de
simulação numérica direta, do inglês, Direct Numerical Simulation ou simplesmente
DNS. Entretanto, segundo KOLMOGOROV [40], o tamanho dos menores vórtices em
um escoamento turbulento depende da viscosidade cinemática
do fluido e é
proporcional a
14
3
()
e , onde
e
é a taxa de dissipação da energia cinética. Esta relação
é conhecida como comprimento de escala de Kolmogorov. Uma conseqüência desta
relação é que a resolução espacial, determinada pelo comprimento da malha
computacional
h
, necessária para capturar o menor vórtice do escoamento deve ser
menor ou igual a este valor, ou seja,
14
3
he
(2.18)
como
3
'
euL
, onde
'
u
é a raiz quadrada da velocidade média do escoamento e
L
é
o comprimento característico do escoamento [65].
19
Substituindo-se a expressão do número de Reynolds na Equação (2.18), obtém-se uma
relação para a resolução espacial mínima para se modelar a turbulência pelo método
DNS
94
3
Re
N
(2.19)
onde
N
é o número de pontos da malha.
Supondo uma aplicação típica de engenharia envolvendo o escoamento de ar, o qual
possui viscosidade da ordem de 10
-6
, seria necessária uma resolução espacial da ordem
de 10
13
. Por conseguinte, embora o poder de processamento dos computadores tenha
avançado consideravelmente nas últimas décadas, a simulação de escoamentos
turbulentos através de DNS ainda está restrita a escoamentos com Reynolds moderados
utilizando geometrias simples.
De modo a contornar os custos computacionais proibitivos de simulações DNS em
aplicações reais, utilizando geometrias complexas e números de Reynolds elevados, a
modelagem da turbulência pode ser feita tal que somente algumas escalas do
escoamento sejam capturadas. Não obstante, as grandezas relacionadas às macro escalas
do escoamento são, em geral, mais importantes do ponto de vista de engenharia [27].
Conseqüentemente, os métodos de modelagem de turbulência possuem o objetivo
comum de reduzir o custo computacional através da solução dos comportamentos
médios do escoamento, geralmente associados às escalas maiores, representando de
alguma forma, os efeitos das escalas menores [65]. O principio desses métodos passa,
primeiramente, pela separação de escalas presentes no escoamento e identificação das
escalas de interesse.
Conforme ilustrado na Figura 2.3b, as grandezas relacionadas às equações de
Navier-Stokes podem ser separadas em um comportamento médio e uma componente
flutuante, geralmente associada às escalas menores do escoamento [27], ou seja,
'

uuu
(2.20)
20
onde
u
é a velocidade do escoamento em um dado instante,
u
é a velocidade média
resolvida e
'
u
é a velocidade flutuante, geralmente, o foco dos modelos de turbulência.
A seguir será descrito o método de modelagem de turbulência adotado neste trabalho.
2.5.2. Método da simulação dos vórtices de larga escala
Neste trabalho adotou-se a modelagem de turbulência pelo método da simulação
dos vórtices de larga escala, do inglês, Large Eddy Simulation, ou simplesmente LES.
Em LES, o procedimento de modelagem consiste em resolver as equações de Navier-
Stokes para as escalas maiores e aproximar os efeitos das escalas menores a partir de
um modelo matemático. A solução das escalas maiores é obtida a partir da aplicação de
um filtro passa-baixa (low-pass), caracterizado, no contexto de métodos discretos, pela
menor escala de comprimento representada pela malha computacional, ou seja, o
comprimento do elemento ou célula [48], [65]. Para filtros homogêneos a aplicação da
filtragem possui propriedade comutativa com as derivadas das equações de Navier-
Stokes, assim, as Equações (2.1) e (2.2), podem ser reescritas na forma filtrada por:
t
p


u
uuTT
I∇
(2.21)

u0
∇
(2.22)
onde a barra superior refere-se à filtragem das grandezas. Desta forma,
T
representa as
tensões cisalhantes filtradas enquanto
T
é a tensão cisalhante residual, também
conhecida por tensão de sub-malha e definida, a partir da notação de Newton (onde
índices repetidos significam somatórias), por
ijijij
Tuuuu

(2.23)
A tensão residual é originada pelo processo de filtragem e reflete o efeito das escalas
excluídas (escalas não-resolvidas) sobre as escalas capturadas e resolvidas pela malha
computacional [48]. O traço do tensor de tensão residual é incorporado à pressão
filtrada enquanto a parcela desviadora ou cisalhante é obtida, em LES, a partir de um
21
modelo matemático. Neste trabalho adotou-se o modelo de Smagorinsky [70] definido a
partir da seguinte relação para a tensão cisalhante residual
2
1
2
3
d
s
ChTTTI
εε
(2.24)
onde
h
é o comprimento do filtro associado a malha computacional, geralmente
calculado, em elementos finitos, como sendo a raiz cúbica do volume do elemento,
s
C
é a constante de Smagorinsky que recebe valores que variam entre 0,1 e 0,2 e
ε
é o
tensor taxa de deformação filtrado
1
2
T
uuuε
(2.25)
e
ε
é a sua norma. Note que por similaridade o termo
2
s
Ch
ε
na Equação (2.24)
assemelha-se à viscosidade da Equação (2.4), desta forma, este termo é geralmente
chamado, em inglês, de eddy viscosity, que em uma tradução livre seria a viscosidade
existente entre os vórtices do escoamento turbulento.
22
Capítulo 3
Formulação de Elementos Finitos
Neste capítulo são apresentadas as formulações estabilizadas de elementos finitos para o
escoamentos incompressíveis e o transporte advectivo-difusivo de uma grandeza escalar
qualquer. A semi-discretização de elementos finitos para as equações de Navier-Stokes
incompressível tridimensional adota a formulação SUPG/PSPG conforme descrita por
TEZDUYAR [80]. Esta formulação, além de evitar oscilações em problemas fortemente
convectivos, devido à parcela SUPG, permite que funções de interpolação de mesma
ordem sejam utilizadas para a velocidade e pressão. Para o problema de transporte
utiliza-se uma associação da estabilização SUPG ao termo de captura de
descontinuidades
YZ
. Ambos os problemas semi-discretos são integrados no tempo
por diferenças finitas com o uso do algoritmo preditor-multicorretor. Para o problema
não-linear associado ao escoamento o método de Newton-inexato é utilizado como laço
de multicorreção. Os sistemas de equações lineares resultantes são tratados pelo método
dos resíduos mínimos generalizados (GMRES). Ao final deste capítulo são apresentadas
as técnicas de pré-condicionamento, controle adaptativo de passo de tempo e
acoplamento entre o escoamento e o transporte.
3.1. Equações de Navier-Stokes
Considere a discretização do domínio
pela sua subdivisão em elementos
e
,
1,2,...,
enel
, onde
nel
é o número de elementos utilizados nesta subdivisão, tal que
23
ee
ij

I e
e

U
. Associados a esta discretização, considere a definição dos
seguintes espaços para as funções de interpolação de velocidade e pressão
hhh1hhh
{|[],em}
ndim
g
SHg
u
uuu&
(3.1)
hhh1hh
{|[],0em}
ndim
g
VH

u
www&
(3.2)
hhhh1h
{|}
pp
SVqqH

(3.3)
onde
ndim
é o número de dimensões no espaço e
1h
()
H
é um espaço de funções de
dimensão finita de derivada quadrado integrável sobre
. Este espaço é formado pelo
uso, no domínio do elemento, de polinômios de primeira ordem. A formulação
estabilizada proposta por TEZDUYAR [80] para as equações (2.1) e (2.2) pode ser
escrita como: Encontrar
hh
S
u
u
e
hh
p
pS
tal que,
hh
V

u
w
e
hh
p
qV
 ,
1
1
:,
,
0
1
el
e
el
e
h
h
hhhhhhhhh
n
h
hhhhhhhe
SUPGPSPG
e
n
hhe
LSIC
e
dpddqd
t
qpd
t

















u
wuufwuwhu
u
uwuuuf
wu
∇εσ∇
∇∇∇∇σ
∇∇
(3.4)
Na equação acima as primeiras quatro integrais correspondem aos termos da formulação
padrão de Galerkin, enquanto as integrais remanescentes são termos de estabilização
acrescentados à formulação. Note que os termos de estabilização são somatórias
realizadas no nível de elemento. A primeira somatória diz respeito à estabilização
Streamline-Upwind/Petrov-Galerkin (SUPG) utilizada para evitar oscilações espúrias
em problemas fortemente convectivos e introduzida por BROOKS e HUGHES em [6].
A segunda somatória é o termo de estabilização Pressure Stabilizing/Petrov Galerkin
(PSPG), proposto inicialmente por HUGHES et al [33] para escoamentos de Stokes e
posteriormente generalizado para escoamentos incompressíveis por TEZDUYAR [80].
Esta estabilização permite que elementos de mesma ordem sejam utilizados satisfazendo
24
a deficiência de posto existente no acoplamento entre os balanços de inércia e massa das
equações de Navier-Stokes. Finalmente, a terceira somatória se deve a estabilização
LSIC (Least Squares Incompressibility Constraint) para escoamentos com números de
Reynolds elevados. Os termos de estabilização SUPG e PSPG, utilizados neste trabalho,
foram obtidos seguindo a referência [79], conseqüentemente:
12
2
2
2
2
4
9
h
SUPGPSPG
h
h













u
(3.5)
onde
h
u
é o vetor de velocidade local e
é a viscosidade cinemática. Na equação (3.5)
o termo de estabilização LSIC é calculado, segundo BEHR et al [4], por:
2
h
LSIC
h
u
(3.6)
A discretização espacial da equação (3.4) dá origem ao seguinte sistema semi-discreto
de equações não-lineares a ser resolvido em cada passo de tempo.
()()()()()


MMaNuNuKKLuGGpff
(3.7)
()
T


GuMaNuKuGpee
(3.8)
onde
u
é o vetor de incógnitas nodais de
h
u
,
a
é a derivada temporal de
u
e
p
é o
vetor de incógnitas nodais de
h
p
. As matrizes
M
,
K
e
G
correspondem aos termos de
derivada temporal (matriz de massa), tensões cisalhantes (viscosa) e tensões normais
(pressão) respectivamente, os vetores
()
Nu
são as derivadas dos termos convectivos,
L
é a matriz de estabilização LSIC e os termos
f
e
e
são as contribuições para os
vetores de força. Os índices subscritos
,
e
dizem respeito aos termos de
estabilização SUPG, PSPG e LSIC respectivamente.
25
3.2. Transporte advectivo-difusivo
Para o problema de transporte advectivo-difusivo expresso pela equação (2.8), considere
os seguintes espaços de aproximação de elementos finitos para as funções de teste e
peso respectivamente
hhh1hhh
{|,em}
g
SHg

(3.9)
hhh1hh
{|,0em}
g
VwwHw

(3.10)
a solução estabilizada de elementos finitos consiste em: encontrar
hh
S
, tal que,
hh
wV
,
1
1
h
el
e
el
e
h
hhhhhhh
n
h
hhhhh
SUPG
e
n
hhh
e
wdwdwhd
t
wd
t
wdwfd













uK
uuK
∇∇∇
∇∇∇∇
∇∇
(3.11)
onde as primeiras duas integrais representam a formulação de Galerkin da Equação
(3.11), a primeira somatória no nível de elemento corresponde à estabilização SUPG,
adicionada a formulação pelos mesmos fins discutidos para as equações de Navier-
Stokes. A segunda somatória corresponde ao termo de captura de descontinuidades
introduzido à formulação na forma de difusão não-linear e é utilizado para prevenir a
ocorrência de oscilações da solução em regiões de gradientes elevados [71]. Neste
trabalho, o parâmetro
se baseia no resíduo da equação (3.11). Este termo foi
recentemente introduzido na referência [3] sob o nome de estabilização
YZ
e sua
definição é dada por
26
1
2
2
3
11
1
(
())
2
h
heh
e
i
i
j
h
R
x






(3.12)
onde
()
eh
i
R
é o resíduo da Equação (3.11):
()·
h
ehh
h
R
t

 
h
u K∇
(3.13)
É importante observar que se
1
e o escalar de referência
1
o termo de captura
de descontinuidades
YZ
se torna o Consistent Approximated Upwind ou CAU
proposto por GALEÃO e DO CARMO em [22]. O parâmetro de estabilização
SUPG
é
calculado por uma relação similar a Equação (3.5)
A discretização espacial da equação (3.11) dá origem ao seguinte sistema semi-
discreto de equações não-lineares,
()()()



MMNNKKKff
&
(3.14)
onde
é a grandeza escalar a ser transportada,
&
é sua derivada temporal, as matrizes
M
,
N
e
K
correspondem aos termos de derivada temporal (matriz de massa),
advecção e difusão respectivamente.
()
K
diz respeito à linearização do termo de
captura de descontinuidades e o subscrito
refere-se aos termos de estabilização
SUPG.
3.3. Integração no tempo
A discretização no espaço das equações descritas nas seções anteriores resulta em um
sistema de equações diferenciais ordinárias (EDO) de primeira ordem em relação ao
tempo, sendo este procedimento conhecido como semi-discretização [10]. A
discretização das equações diferenciais parciais (EDP) completa-se com a aplicação de
27
um método de integração no tempo para as EDO’s resultantes do método dos elementos
finitos [10]. Neste trabalho, a dimensão do tempo foi discretizada por diferenças finitas
centradas originando no algoritmo preditor multicorretor conforme descrito a seguir.
3.3.1. Algoritmo preditor multicorretor
Para construção de um procedimento de integração no tempo das Equações (3.7), (3.8) e
(3.14) considere que
max
n
é o número máximo de passos de tempo,
max
i
é o número
máximo de iterações de multicorreção e
é o parâmetro que controla a aproximação da
derivada temporal
a
. Seguindo o procedimento descrito por FRANCA e FREY em [21]
obtém-se o seguinte algoritmo de integração no tempo para as equações de Navier-
Stokes:
Algoritmo 3.1 – Integração no tempo das equações de Navier-Stokes incompressíveis
Passo 1: Inicie
0
0
u
,
0
0
p
e
0
0
a
Passo 2: Atribua
0
n
e
0
i
Passo 3: Fase de predição:
()
()()
1
()
1
()
()
1
(1)
0
i
ii
nnn
i
n
i
i
nn
t

uua
a
pp
(3.15)
Passo 4: Montar as matrizes
*
M
,
*
G
,
*
()
T
G
,
L
e
G
()()
*
11
*
*
()(),
,
()()
ii
nn
TT
t
t









N
N
MMMuuKL
uu
GGG
GMNuKG
(3.16)
28
Passo 5: Montar os resíduos
r
e
s
()()()()()()
111111
()()
11
()()()()()()
111111
()()
()
iiiiii
nnnnnn
ii
nn
iiiiii
T
nnnnnn












rffMMaNuNu
KLuGG
seeGuMaNuKuG
p
p
(3.17)
Passo 6: Resolver o sistema incremental:
()()
**
11
*()()
11
()
ii
nn
Tii
nn














MGar
GG
sp
(3.18)
Passo 7: Fase de correção
(1)()()
111
(1)()()
111
(1)()()
111
iii
nnn
iii
nnn
iii
nnn
t






uua
aaa
ppp
(3.19)
Passo 8: Se
max
ii
, então
1
ii

e retorna ao Passo 4, senão,
0
i
e:
max
max
max
()
()
11
()
()
11
()
()
11
i
i
nn
i
i
nn
i
i
nn



uu
aa
pp
(3.20)
Passo 9: Se
max
nn
, então
1
nn

e retorna ao Passo 3, senão PARE.
OBSERVAÇÕES:
− Para elementos lineares as matrizes viscosas
K
e
K
podem ser desprezadas
pela existência de derivadas de segunda ordem;
− A linearização dos termos de convecção nas Equações (3.7) e (3.8) torna os
passos correspondentes à multi-correção no algoritmo acima (passos 4 a 7)
29
similares ao método de Newton. Neste trabalho, o método de Newton inexato
(detalhado na próxima seção) foi utilizado em substituição a estes passos;
− Neste trabalho foi adotada uma estratégia de controle de adaptação de passo de
tempo por controlador Proporcional-Integral-Derivativo (PID) para escolha do
valor de
t
, também descrita adiante;
− O sistema de equações lineares correspondente ao Passo 6 é não-simétrico.
Neste trabalho as variáveis
a
e
foram resolvidas de forma acoplada com
o método dos resíduos mínimos generalizados (GMRES);
Para a integração no tempo das equações de transporte o algoritmo é semelhante ao
descrito no Algoritmo 3.1, bastando substituir
u
por
,
a
por
&
, omitir os termos
referentes à pressão e substituir as Equações 3.16-3.18, por:
()
*
1
(),
i
n
t





K
MMMNNKK
(3.21)
()()()()()()
111111
()()
iiiiii
nnnnnn






rffMMNNKKK
&
(3.22)
()()
*
11
ii
nn


Mar
(3.23)
3.4. Solução não-linear
Na solução do sistema de equações descrito na Equação (3.18) diversos métodos de
solução não-lineares podem ser utilizados, dentre os quais, podem-se citar: o método de
Picard, também conhecido como método de ponto fixo ou iterações sucessivas, o
método de Broyden e a família dos métodos de Newton como estratégias
freqüentemente utilizadas [19]. No presente trabalho, optou-se pela utilização do
método de Newton inexato com base nos resultados obtidos no estudo conduzido por
ELIAS et al [14], [15].
30
3.4.1. Método de Newton inexato
Para descrição do método de Newton-inexato, primeiramente, considere o sistema de
equações não-lineares descrito pela Equação (3.18) e reescrito de maneira compacta
NdF
(3.24)
A cada iteração, dada uma solução aproximada
()
i
d
calcula-se a correção
()
i
d
Δ
, tal que
()()ii

NddF
Δ
(3.25)
onde
()
i
d
Δ
pode ser obtido a partir da seguinte expansão em série de Taylor em torno de
()
i
d
()
()()()()
i
iiii

d
N
NddNdd
d
ΔΔ
(3.26)
substituindo-se (3.26) em (3.25), e fazendo
()()
ii
rFNd o resíduo da iteração
não-linear, obtém-se:
()
()()
i
ii
d
N
dr
d
Δ
(3.27)
onde a atualização da solução é feita por:
(1)()()
iii

ddd
Δ
(3.28)
No procedimento descrito acima, a derivada do lado esquerdo da equação (3.27) é
comumente chamada de Jacobiano de
N
e corresponde aos termos

Nu
da Equação
(3.16). Para simplificar a notação, esta matriz será chamada de
J
, ou seja,
31
()
()
()
i
i
d
N
Jd
d
(3.29)
O método descrito pelas equações (3.27) e (3.28) corresponde ao método de
Newton-clássico e converge quadraticamente desde que sejam observados os seguintes
aspectos:
− A solução inicial
(0)
d
deve estar próxima o suficiente da solução do sistema,
caso contrário, o método diverge, a não ser que sejam tomadas providências no
sentido de contornar tal problema;
− A inversa de
J
existe em todas as iterações necessárias até a convergência ser
atingida;
− O Jacobiano
J
seja obtido de forma exata;
− A solução do sistema de equações lineares (Equação (3.27)) seja obtida por um
método analítico ou direto dentro da precisão da máquina;
Caso os termos correspondentes às derivadas que formam a matriz de rigidez
tangente não sejam avaliados e os valores provenientes de iterações anteriores sejam
utilizados na linearização, o método de Newton recai no método de iterações sucessivas
(IS) também conhecido como método de Picard [19].
Os diversos métodos tipo Newton são baseados na expansão em série de Taylor da
equação (3.26) próxima da solução. Sendo assim, o cálculo das soluções lineares no
início das iterações não-lineares não necessita de grande precisão [60]. Esta necessidade
aumenta na medida em que as iterações não-lineares vão progredindo em direção à
solução final no equilíbrio. Portanto, a idéia do método de Newton Inexato, também
conhecido por Newton Truncado, é:
Minimizar o esforço empregado para a solução dos sistemas de equações lineares
quando a direção de busca não-linear estiver longe da solução de equilíbrio.
32
A análise matemática do efeito da aplicação desta técnica ao comportamento da
convergência da solução calculada pelo método não-linear pode ser encontrada em [12],
[13] e [38].
A taxa de convergência assintótica e o esforço computacional no método de Newton
Inexato são controlados por um parâmetro calculado,
i
, o qual especifica o critério de
terminação truncado requerido para a solução linear dentro da iteração não-linear
i
[38]. Um método de determinação do critério de terminação truncada é sugerido por
KELLEY [38] onde a medida do grau de aproximação da iteração não-linear
i
em
relação à solução é dada por
2
()
()
2
(1)
i
i
A
i

r
r
(3.30)
onde
(0,1]
e
g
é a norma Euclideana. Para a implementação do algoritmo, se
()
i
A
é uniformemente limitada longe de 1, então, fazendo
()
()
i
i
A

para
0
i
produz
uma convergência q-quadrática [38]. Dessa forma, a tolerância utilizada na solução do
sistema linear é eficientemente explorada. No intuito de se determinar o valor da
tolerância na primeira iteração, bem como nas demais, mantendo o limite da escolha
sempre longe de 1, tem-se
max
()
()
max
,0,
min,,0.
i
i
B
A
i
i

(3.31)
em (3.31), o parâmetro
max
é um limite superior da seqüência
()
i
,
max
1,2,...,
ii
.
Os valores sugeridos por EISENSTAT [12] são
0,99
e
max
0,9999
.
Pode ocorrer, entretanto, que
()
i
B
seja pequeno demais para uma ou mais iterações
enquanto a iteração não-linear ainda esteja longe da solução. Para evitar este tipo de
33
comportamento, KELLEY [38] sugere que se
(1)
i
for suficientemente grande, o valor
seguinte, η
i
não pode diminuir muito mais que um fator de
(1)
i
, ou seja,
max
()()(1)
2
max
()(1)(1)
22
max
,0,
min,,0,()0,1,
min,max,(),0,()0,1.
iii
CA
iii
A
i
i
i





(3.32)
sendo o valor
0,1
arbitrário.
Há casos em que a iteração não-linear reduz
()
i
r
a valores muito pequenos, além
do nível de precisão desejado, produzindo, portanto, elevado custo no cálculo da
solução linear, sendo esta mais precisa que o necessário. Este esforço pode ser reduzido
ao necessário, se a norma do resíduo da iteração não-linear corrente for comparada com
a norma do resíduo determinante do final das iterações não-lineares
(0)
res
 r
(3.33)
e limitando
()
i
por uma constante múltipla de
()
i
r
. KELLEY [38] sugere
()
()
max
min,max,0.5
i
i
Ci
 r
(3.34)
O cálculo de
conforme a equação (3.30) resulta em um algoritmo adaptativo para
a solução do problema não-linear. Ou seja, quando a iteração não-linear está longe da
solução,
()
i
r
é grande e, portanto, pouco esforço é necessário para se obter uma
solução que satisfaça à tolerância corrente. Conforme as iterações não-lineares vão se
aproximando da solução, o valor de
()
i
r
vai diminuindo e, portanto, o esforço para o
cálculo da solução é maior, pois a precisão necessária aumenta.
34
O Algoritmo 3.2 apresenta os passos principais correspondentes ao método de
Newton inexato. Para o algoritmo completo, consulte ELIAS [14].
Algoritmo 3.2 – Método de Newton inexato
Passo 1: Encontrar um
[0,1)
i
a partir das Equações 3.30-32
Passo 2: Efetuar controle de oversolving conforme Equações (3.33) e (3.34)
Passo 3: Encontrar a solução para
d
de modo que a seguinte condição seja satisfeita
()()(1)()()
()
iiiii

rJddr
(3.35)
Passo 4: Atualização da solução:
(1)()()
iii

ddd
(3.36)
Passo 5: Se
(1)(0)
toler
i
rr e
max
ii
retorne ao Passo 1, senão, PARE.
OBSERVAÇÕES:
− Neste trabalho, o método dos resíduos mínimos generalizados (GMRES) foi
utilizado para satisfazer a condição requerida no Passo 3 (Equação (3.35));
− O critério de convergência baseado no resíduo não-linear relativo,
(1)(0)
i
rr
, conforme descrito no Passo 5, pode não ser suficiente para
representar bem o estado de convergência do sistema. Segundo [19], o critério
de convergência também deve levar em consideração o incremento de passo
relativo, definido por:
()()
ii
dd, pois ambos tendem a zero quando o
problema não-linear converge.
35
3.5. Solução linear
Os métodos de solução de sistema de equações nos quais a matriz de coeficientes é
fatorada são chamados de métodos diretos [64]. Tais métodos se tornam
computacionalmente caros para matrizes esparsas, especialmente, em problemas onde a
matriz de coeficientes é não simétrica [64].
Os métodos iterativos são uma alternativa viável aos métodos diretos. Tais métodos
geram uma seqüência de aproximações da solução e, essencialmente, utilizam a matriz
de rigidez somente em multiplicações com vetores [64]. Neste trabalho o método dos
resíduos mínimos generalizados (GMRES) foi utilizado como solucionador para o
tratamento dos sistemas de equações lineares dos problemas de escoamento e transporte.
3.5.1. Método dos resíduos mínimos generalizados (GMRES)
Considere os sistemas de equações (3.18) (3.35) representados de forma compacta por:
Axb
(3.37)
onde
1
,...,
T
neq
xxx é o vetor de variáveis nodais sendo
neq
o número de
incógnitas;
A
uma matriz esparsa não-simétrica
neqneq
e
b
é o vetor de termos
independentes, também de ordem
neq
.
A matriz
A
e o vetor
b
são construídos a partir da montagem das contribuições dos
elementos da malha computacional, isto é,
1
nel
e
e
AA
A
(3.38)
1
nel
e
e
bb
A
(3.39)
36
O método iterativo GMRES, do inglês, Generalized Minimum RESidual, [64], [68],
é um método iterativo particularmente atraente para a solução de sistemas de equações
não-simétricos. Esse método tem por objetivo básico minimizar a norma residual do
sistema (3.37).
Considerando
0
x
como solução inicial, uma solução aproximada irá apresentar a
forma
0
xz
onde
z
é um vetor do espaço de Krylov definido como
21
0000
,,,...,
k
span
rrrr
ΑΑΑ= , sendo
00

rbAx
e
k
a dimensão de
. O
método GMRES determina
z
tal que a norma do resíduo seja mínima, isto é,
0
xz
é
a solução de (3.37) se
0

bAxz
é mínima, onde
é uma norma, geralmente,
euclideana.
Uma base ortonormal de
pode ser obtida pelo processo modificado de Gram-
Schmidt [68]:
0
1
0
r
u
r
Para
1,...,
ik
faça:
1
ii
uAu
%
Para
1,...,
ji
faça:
1,1
,
ijij

uu
%%
111,
iiijj


uuu
%%
Fim do
j
1
1
1
i
i
u
u
u
%
%
Fim do
i
Seja
12
,,...,
kk



Uuuu
. Pode-se então, mostrar que
1
kkk
AUUH
onde
k
H
é a matriz superior de Hessenberg de ordem
1
kk

,
37
2,13,1,11,1
23,2,21,2
3
,1
1,
1
0
000
kk
kk
k
kk
kkk
k















u
u
H
u
u
L
%
L
%
OMM
MMOM
%
MMO
%
L
(3.40)
Seja
1
k
jj
j
yu
z e
0
,0,...,0
T
er com
1
k
termos. Observe que
01
k
rUe
, portanto, tem-se:
00
1
0
1
k
jj
j
k
kk
k
xzyu




bArA
rAUy
UeHy
eHy
(3.41)
Assim o problema de minimização em questão pode ser reescrito por:
0
minmin
k
zk
xz

y
bAeHy
¡
κ
(3.42)
A matriz
k
H
tem uma estrutura quase triangular superior, necessitando-se somente
eliminar a diagonal contendo
231
,,...,
i
uuu
%%%
. A transformação de
k
H
em uma
matriz triangular superior pode ser feita pelo algoritmo QR, usando em cada passo o
processo de rotação de Givens para obter um novo problema de minimização. Maiores
detalhes podem ser vistos em [68]. Assim, obtém-se:
min
k
k
y
eHy
¡
(3.43)
onde
k
H
é uma matriz triangular superior. O sistema pode então, ser resolvido
facilmente por retro-substituição e um novo resíduo é calculado com a aproximação
encontrada. Se a convergência não for atingida dentro de um ciclo (formado pelo
38
número de vetores na base de Krylov), a solução mais recente será usada para iniciar a
próxima iteração. Resumindo, uma iteração do Método GMRES consiste de:
− Formação de uma base ortogonal para o espaço de Krylov pelo processo
modificado de Gram-Schmidt;
− Triangularização da Matriz de Hessenberg pelo algoritmo QR;
− Retro-substituição para obter a solução aproximada.
Com base nestas considerações descreve-se um algoritmo padrão do método
GMRES. Para tal, definem-se os seguintes parâmetros:
ε
tol
→ tolerância desejada;
k → número de vetores de Krylov;
l
max
→ número máximo de iterações.
Algoritmo 3.3 – Método dos resíduos mínimos generalizados (GMRES)
Dados A, b, k, ε
tol
e l
max
(Inicializações)
tol
b

0
x
(ciclo GMRES)
Para
max
1,...,
ll
faça
11
1


ubLAUx
1
e
u
1
1
1
u
u
u
(iterações GMRES)
Para
1,...,
ik
faça
11
1
ii

uLAUu
39
(Ortogonalização modificada de Gram-Schmidt)
Para
1,...,
ji
faça
1,1
,
ijij

uu
111,
iiijj


uuu
Fim do
j
(Fim da ortogonalização modificada de Gram-Schmidt)
1,11,1
,...,,
T
i
iiii


hu
1
1
1
i
i
i
u
u
u
(Algoritmo QR)
Para
1,...,1
ji

faça


11
ii
jjjj
ii
jj
jj
hcsh
sc
hh















Fim do
j


22
1
ii
ii
rhh

i
i
i
h
c
r
1
i
i
i
h
s
r
i
i
hr
1
0
i
i
h
1
iii
ese

iii
ece
(Fim do algoritmo QR)
Teste de convergência: Se
1
itol
e
, fim do
i
(Fim da iteração GMRES)
Resolver
y
:
40


11
111
11
1
11
11
0
00
ii
ii
ii
ii
i
ii
i
hhh
ye
ye
hh
ye
h



























L
OMM
MM
MO
L
Cálculo da Solução:
1
i
jj
j
yu

xx
Teste de Convergência: Se
1
i
e
, fim do
l
(Fim do ciclo GMRES)
x
é a solução aproximada
3.6. Pré-condicionamento
Segundo SAAD [64], a falta de robustez dos métodos iterativos é uma fraqueza em
relação aos métodos diretos. Tal característica muitas vezes impede que os métodos
iterativos sejam amplamente utilizados em aplicações industriais, mesmo sendo eles
apropriados para a solução de grandes sistemas de equações. Uma forma de melhorar a
robustez e eficiência dos métodos iterativos é conhecida por pré-condicionamento. O
pré-condicionamento consiste na transformação do sistema de equações lineares
original em outro de mesma solução, porém, mais bem condicionado e
conseqüentemente de solução mais fácil [64]. Por exemplo, considere a matriz
A
positiva definida e seu pré-condicionador
M
. O pré-condicionador
M
é uma matriz
que aproxima
A
de alguma forma. Assumindo-se que
M
também é positiva definida, o
único requisito para
M
é que o custo computacional da solução de
Mxb
seja baixo.
Isto se deve ao fato de que esta matriz de pré-condicionamento será utilizada no auxílio
à solução do sistema de equações original em cada iteração não-linear para cada passo
de tempo discutido nas seções anteriores. Assim, o sistema pré-condicionado poderia
ser escrito da seguinte forma:
11

MAxMb
(3.44)
41
ou ainda
11
,


AMxbxMu
(3.45)
As formas como o pré-condicionamento é comumente aplicado aos sistemas de
equações são relacionadas abaixo:
11

LAxLb
(3.46)
111

LAUUxLb
(3.47)
11

AUUxLb
(3.48)
onde (3.46) é conhecido como pré-condicionador à esquerda, (3.47) é conhecido como
pré-condicionador à esquerda e a direita e finalmente (3.48) é conhecido como pré-
condicionador à direita.
L
e
U
podem representar simples decomposições ou
complicadas transformações de
A
. Segundo SAAD em [64], encontrar um bom pré-
condicionador para resolver um sistema linear é freqüentemente visto como uma
combinação de arte e ciência. Resultados teóricos são raros e algumas técnicas oferecem
resultados surpreendentes para algumas aplicações, e para outras nem tanto. A
eficiência da operação matriz-vetor, mais uma vez é a maior preocupação, pois nela
baseia-se a transformação do sistema. A seguir são apresentadas duas técnicas de pré-
condicionamento, particularmente convenientes para uso com estruturas de dados por
arestas e, utilizadas neste trabalho como alternativa para acelerar a convergência do
método GMRES citado na seção anterior.
3.6.1. Pré-condicionamento diagonal
Uma estratégia bastante simples de pré-condicionamento consiste no escalonamento das
linhas da matriz
A
por um pré-condicionador à esquerda, ou das suas colunas por um
pré-condicionador à direita [52]. Desta forma, dada uma matriz
11
,...,
nn
diagddD
42
como sendo um pré-condicionador à esquerda ou à direita, existem várias possibilidades
para a obtenção de
D
, por exemplo:
1
ii
d
a
(3.49)
1
ii
ii
d
a
(3.50)
1
1
ii
n
ij
j
d
a
(3.51)
12
2
1
1
ii
n
ij
j
d
a


(3.52)
1,...,
1
max
ii
ij
jn
d
a
(3.53)
O pré-condicionamento diagonal apresenta algumas vantagens óbvias sobre outras
técnicas, tal como, fácil implementação e baixo custo de armazenamento da matriz de
pré-condicionamento. Em contrapartida, não é esperada grande aceleração do
solucionador linear, por ser a matriz
D
uma aproximação bastante “ingênua” da inversa
de
A
. Note que as formas (3.49), (3.51), (3.52) e (3.53) são usadas à esquerda ou à
direita conforme as equações (3.46) ou (3.48), enquanto que a forma (3.50) é usada à
esquerda e a direita de acordo com a equação (3.47).
3.6.2. Pré-condicionamento bloco diagonal nodal
Seja uma matriz
B
da forma:
block
BA
(3.54)
43
onde o operador
block
extrai de
A
as matrizes bloco-diagonais nodais de ordem
mm
, onde
m
é o número de graus de liberdade por nó (
4
m
) para as equações
de Navier-Stokes incompressíveis em três dimensões, conforme mostrado na Figura 3.1.
Figura 3.1 - Pré-condicionador Bloco-Diagonal Nodal
Para os sistemas em consideração,
B
deve ser positiva-definida e de fácil inversão
de forma que o sistema (3.37) transformado fique:
11

BAxBb
(3.55)
O bloco-diagonal nodal representa uma técnica de pré-condicionamento de fácil
implementação, com custo de armazenamento e computacional maior que no pré-
condicionador diagonal, porém com uma maior eficiência na aceleração da solução
iterativa de sistemas de equações.
3.7. Controle adaptativo de passo de tempo
A solução do sistema apresentado pela Equação (3.18) é computacionalmente cara e
envolve a solução de um sistema não simétrico a cada iteração de multi-correção dentro
de cada passo de tempo. Segundo LOHNER [47], embora os métodos implícitos de
integração no tempo, como o adotado no presente trabalho, sejam considerados
incondicionalmente estáveis em relação ao passo de tempo
t
, o uso de valores
arbitrariamente grandes é raro em aplicações práticas. Problemas de interesse da
engenharia geralmente apresentam um alto grau de não-linearidade e as técnicas
utilizadas na linearização normalmente interferem no critério de estabilidade dos
métodos implícitos [47]. Portanto, uma alternativa que permita que valores de passo de
B =
44
tempo grandes sejam utilizados sem comprometer a precisão e estabilidade do
procedimento de solução torna-se particularmente atraente.
O controle adaptativo de passo de tempo é uma técnica desenvolvida com o intuito
de minimizar o custo computacional na obtenção de uma solução aproximada
preservando-se a qualidade da solução. Tal técnica baseia-se na escolha do passo de
tempo mantendo-se o controle a partir de uma métrica de erro [83], [84]. VALLI [84]
analisou o uso de controladores tipo PID (Proporcional-Integral-Derivativo) aplicado a
problemas envolvendo o acoplamento entre escoamento incompressível e transporte.
VALLI [84] focou a solução de problemas bidimensionais utilizando formulação
penalizada e realizou um estudo detalhado das respostas do controlador frente aos
diversos parâmetros empregados pelo mesmo. No presente estudo, o trabalho
desenvolvido em [83] e [84] foi estendido para o tratamento da mesma classe de
problemas, entretanto, utilizando a formulação estabilizada SUPG/PSPG para o
escoamento incompressível e SUPG/
YZ
para o transporte em problemas
tridimensionais de larga escala e em paralelo.
Para utilização do controlador PID, considere o problema de controle do passo de
tempo como sendo um problema de controle automático retro-alimentado conforme
sugerido por GUSTAFSSON [28] e ilustrado na Figura 3.2.
Figura 3.2 – Seleção do passo de tempo como um problema de controle retro-
alimentado
No contexto utilizado neste trabalho, o processo é representado pelo método de
integração no tempo.
Controlador Processo
tol
∆t
erro
retroalimentação
45
Em [83] VALLI et al propõem duas estratégias de controle. A primeira estratégia
baseia-se na taxa de mudança das grandezas de interesse (velocidade, pressão e escalar)
e a segunda estratégia monitora a taxa de mudança da energia cinética. Neste trabalho,
seguindo as conclusões apresentadas por VALLI et al [83], foi adotado o controle pela
taxa de mudança das grandezas de interesse, definido por
PIP
kkk
DDD
k(2k)k
12
111
prev
nnn
tt
eee













(3.56)
onde
max(,,)
np
eeee
u
(3.57)
*
*
1
,
nn
n
e
ee
tol

u
uu
u
uu
u
(3.58)
*
*
1
,
p
nn
pp
p
n
e
pp
ee
tol
p
 (3.59)
*
*
1
,
nn
n
e
ee
tol


 (3.60)
e
t
é o novo passo de tempo obtido,
prev
t
é o passo de tempo de iterações
anteriores,
P
k
,
I
k
e
D
k
são os ganhos do controlador correspondente às ações
proporcional, integral e derivativa respectivamente e
tol
u
,
p
tol
e
tol
são tolerâncias
fornecidas ao algoritmo para a taxa normalizada de mudança (Equações 3.58-60) da
velocidade, pressão e escalar respectivamente. Alternativamente a condição de
estabilidade de Courant-Friedrichs-Levy (CFL), pode ser utilizada como parâmetro de
controle ao invés do passo de tempo propriamente dito. Neste caso, considerando-se que
a condição de estabilidade determina que
1.0
CFL
[89], o passo de tempo pode ser
obtido a partir das seguintes relações:
46
e
e
t
h
CFLt
h

u
u
(3.61)
onde o
t
considerado será o menor valor obtido para todos os elementos da malha.
Com base nas definições para o controlador PID, descritas pelas Equações 3.56-60,
o algoritmo para a estratégia de controle de adaptação de passo de tempo, descrito em
[83], se resume ao Algoritmo 3.4.
Algoritmo 3.4 – Escolha adaptativa de passo de tempo por controlador PID.
Dados de entrada
n
u
,
1
n
u
,
n
p
,
1
n
p
,
n
,
1
n
,
t
,
t
,
n
,
nrej
,
,
P
k
,
I
k
,
D
k
,
tol
u
,
p
tol
,
tol
,
min
t
e
max
t
Passo 1: Se
1
n
, então inicie:
2
1.0
n
e
,
1
1.0
n
e
,
min
prev
tt

,
0
nrej
Passo 2: Calcule
n
e
a partir das Equações 3.57-60
Passo 3: Obtenha um novo passo de tempo ou rejeite uma solução ruim:
(rejeita passo de tempo)
SE
min
((1.0)e())
n
ett

ENTÃO:
1
1
1
1,
,
,
,
1
1
nn
nn
nn
n
nrejnrej
pp
ttt
nn
fatore




uu
Se
(0,8)
fator
então
0,8
fator
47
min
2
max(,)
prevprev
tfatortt
ttt


(calcula novo passo de tempo)
SENÃO:
Calcule o novo passo de tempo pela Equação (3.56)
min
max
21
1
max(,)
min(,)
prev
nn
nn
ttt
ttt
tt
ee
ee




FIM SE
3.8. Acoplamento escoamento incompressível – transporte
A solução de alguns problemas requer o acoplamento entre os procedimentos de solução
para o escoamento e o transporte enquanto outros necessitam de somente um destes
procedimentos. Neste trabalho os sistemas de equações referentes ao escoamento e ao
transporte são tratados em estágios separados. O Algoritmo 3.5 ilustra a estrutura
principal do software desenvolvido. Note que os procedimentos de solução para o
escoamento e o transporte são módulos independentes e podem ser utilizados
separadamente. O acoplamento entre os problemas ocorre naturalmente a partir das
variáveis em comum (velocidade, pressão e escalar). Um conjunto de chaves de controle
(flags) é utilizado para controlar o fluxo do programa, assim como, alocação e reuso de
memória. Desta forma, o usuário pode, por exemplo, fornecer um campo de velocidades
obtido analiticamente (ou calculado por outro software), e realizar uma simulação do
transporte de uma grandeza qualquer, considerando que não haja maiores interferências
pela falta do acoplamento.
48
Algoritmo 3.5 – Laço principal do software desenvolvido mostrando o acoplamento
entre os procedimentos de solução para o escoamento e o transporte.
FAZER de
0
t
até
max
t
SE (Resolve Escoamento) ENTÃO:
FAZER até convergir ou
max
ii
:
Resolve escoamento incompressível utilizando:
§ Método não-linear: Newton-Inexato
§ Método linear: GMRES
§ Pré-condicionador: bloco diagonal nodal
§ Estrutura de dados: aresta-por-aresta
FIM FAZER.
FIM SE.
SE (Resolve Transporte) ENTÃO:
Fazer até convergir ou
max
ii
Resolve transporte utilizando:
§ Método não-linear: Newton simples (multi-correções)
§ Método linear: GMRES
§ Pré-condicionador: diagonal
§ Estrutura de dados: aresta-por-aresta
FIM FAZER.
FIM SE.
SE (Utilizar controlador PID) ENTÃO:
Determina novo passo de tempo
FIM SE.
FIM FAZER.
4
9
Capítulo 4
Escoamentos de superfície-livre
Neste capítulo é apresentado o método de simulação de escoamentos de superfície-livre
desenvolvido no presente trabalho. O método consiste no transporte Euleriano de uma
função de marcação utilizada na definição das regiões preenchidas pelos diferentes
fluidos, geralmente, água e ar. Técnicas de correção de massa e redução do custo
computacional também foram desenvolvidas para melhorar a qualidade da solução sem
afetar a eficiência computacional.
4.1. Tratamento da superfície-livre
Os escoamentos nos quais a região de interesse é a interface entre dois fluidos são
considerados escoamentos de superfície-livre [39]. Tais tipos de escoamento são
importantes em diversos problemas de hidrodinâmica. O movimento de fluidos no
interior de tanques (fenômeno conhecido em inglês como sloshing), ondas se quebrando
contra navios, decks, áreas costeiras, plataformas, inundação de conveses, danos
ocasionados por tsunamis são exemplos de problemas que se enquadram neste tipo de
escoamento [39]. O maior desafio na simulação computacional de escoamentos de
superfície-livre é a determinação da posição da interface que separa os fluidos. Existe
uma ampla gama de métodos destinados à solução de problemas de superfície-livre.
Tais métodos freqüentemente são classificados de acordo com a forma em que o
problema é tratado. Segundo TEZDUYAR [77], métodos no qual a malha
computacional acompanha o movimento da superfície-livre sendo, portanto,
50
Lagrangeanos, são classificados como métodos de acompanhamento de interface, do
inglês, interface tracking methods. Alternativamente, os métodos que fazem uso de uma
função auxiliar para marcar a localização dos fluidos, sendo, portanto, métodos
Eulerianos, são conhecidos por métodos de captura de interface, do inglês, interface
capturing methods [77]. A Figura 4.1 apresenta esquematicamente estas duas formas de
tratamento da superfície-livre através de métodos discretos.
(a)
(b)
Figura 4.1 – Formas de tratamento da superfície-livre (a) acompanhamento da interface
(b) captura da interface.
Observa-se que no método de acompanhamento de interface (Figura 4.1a) a malha deve
se acomodar ao posicionamento da interface e um procedimento de atualização e/ou
reconstrução da malha deve ser adicionado ao procedimento de solução enquanto que
no método de captura de interface (Figura 4.1b) a interface é implicitamente definida
por uma função de marcação com o custo adicional de se resolver uma equação de
transporte para a função de marcação. A Tabela 4.1 sintetiza as principais virtudes e
dificuldades associadas aos métodos de tratamento de superfície-livre citados
anteriormente.
51
Tabela 4.1 – Dificuldades e qualidades de acordo com a forma
de tratamento da superfície-livre
Acompanhamento
(Lagrangeano) - Figura 4.1a
Captura
(Euleriano) - Figura 4.1b
Virtudes § Precisão
§ Utiliza malhas grossas
§ Malha fixa (não necessita
atualização da malha);
§ Flexível (fácil implementação);
§ Permite grandes deformações,
quebras da interface, etc...
Dificuldades § Restrito a pequenas
deformações;
§ Necessidade de mover a
malha;
§ Surgimento de elementos
distorcidos (necessidade de
recriar a malha)
§ Necessidade de resolver uma
equação de transporte;
§ Necessidade de malhas bastante
refinadas (pelo menos na região
onde a interface se localiza)
Dada a facilidade de implementação, bem como a flexibilidade de uso na solução de
problemas envolvendo ondas se quebrando e fragmentando, neste trabalho adotou-se
exclusivamente o método de captura de interface para o tratamento da superfície-livre
conforme detalhado no Apêndice D. É importante observar que, devido as
características de implementação do método de captura de interfaces, qualquer software
que contemple a solução das equações de escoamento acopladas a uma equação de
transporte possui, intrinsecamente, os requisitos necessários para resolver problemas de
superfície-livre.
Os métodos de captura de interface, conforme citados anteriormente, somente
necessitam de uma função auxiliar para se determinar a posição dos fluidos envolvidos.
Assumindo-se que não há dissolução entre os fluidos, a parcela de difusão da Equação
(2.8) pode ser desprezada e a equação de transporte passa a ter caráter hiperbólico ou
puramente advectivo, ou seja, a Equação (2.8) se reduz simplesmente a:
0
t

u ∇
(4.1)
onde
é a função de marcação e
u
é o campo de velocidades proveniente da solução
da equação (2.1). Normalmente, os métodos Eulerianos de superfície-livre se
52
distinguem de acordo com a função de marcação escolhida. Nas próximas seções as
funções de marcação adotadas neste trabalho são apresentadas.
4.1.1. Método do volume-de-fluido
O método conhecido por Volume-de-Fluido foi inicialmente proposto por HIRT e
NICHOLS [30] para a simulação de escoamentos de superfície-livre utilizando malhas
estruturadas. Neste método a função de marcação determina a fração mássica ou
volumétrica dos fluidos envolvidos, ou seja, a função de marcação irá determinar se
uma dada região do domínio estará 100%
1,0
preenchida de um fluido, ou não
(
0
) e desta forma a interface estará determinada nos pontos do domínio de fração
igual a 50%
0,5
. A Figura 4.2 apresenta, à direita, a função de marcação para o
método VOF com as regiões dos fluidos correspondentes, à esquerda.
Figura 4.2 – Função de marcação para o método do volume-de-fluidos (VOF).
No método VOF, as propriedades do fluido necessárias para a solução do escoamento
são interpoladas a partir das seguintes expressões:
10
1


(4.2)
10
1


(4.3)
53
Note que a função de marcação possui caráter descontínuo e tende a apresentar
picos (overshoots) e/ou depressões (undershoots) para soluções que utilizam somente a
parcela de estabilização SUPG para a solução da Equação (4.1). Uma prática que
remedia tal comportamento é a utilização do termo de captura de descontinuidades
descrito na seção 3.2. Entretanto, mesmo com a introdução deste termo, a solução pode
sofrer de difusão numérica conforme ilustrado na Figura 4.3. Na prática, isto se reflete
em perda de massa e descaracterização da função de marcação. Neste trabalho foram
estudadas diversas combinações dos termos de estabilização para a equação de
transporte associados a algoritmos de correção de massa, além de transformações da
função de marcação e truncamento dos valores da mesma entre os valores 0 e 1. Os
resultados desta análise comparativa são apresentados no artigo do Apêndice D.
Figura 4.3 – Difusão da função degrau do método VOF
4.1.2. Método das curvas de nível
O método das curvas de nível, do inglês, Level Sets ou simplesmente LS, possui uma
vasta gama de aplicações em computação gráfica e foi apresentado, no contexto da
simulação de escoamentos de superfície-livre, devido aos trabalhos de SETHIAN [67] e
OSHER [58]. Neste método, a função de marcação determina a menor distância de um
dado ponto da malha até a interface. Assim sendo, a interface entre os fluidos estará
implicitamente determinada pelos valores da função de marcação iguais a 0 (
0
),
ou seja, a curva de nível 0. Conseqüentemente a função de marcação utilizada é do tipo
função distância. Geralmente esta função é dotada de sinal para distinguir as regiões do
domínio que pertencem a cada fluido, assim, os valores de
0
podem ser
54
relacionados ao Fluido A, enquanto
0
indica a região preenchida pelo Fluido B. A
Figura 4.4 apresenta a função de marcação para o método de curvas de nível com sua
respectiva correspondência às regiões do domínio à esquerda. Note a sobreposição entre
as funções de marcação do método VOF em relação ao método de curvas de nível.
Figura 4.4 – Função de marcação para o método de curvas de nível (LS)
Uma característica interessante deste tipo de função de marcação é a existência de um
gradiente constante de propagação da função distância. Isto se reflete na velocidade de
propagação das curvas de nível a partir do nível 0 e é freqüentemente utilizada em
aplicações que vão além da simulação de superfície-livre, como por exemplo:
determinação da menor distância em ambientes congestionados, detecção de colisão,
reconstrução de imagens, planejamento de rotas de fuga, etc... [58], [67]. Na prática,
esta característica pode ser prontamente utilizada, em escoamentos de superfície-livre,
na modelagem de tensão superficial do fluido a partir da curvatura da função distância,
definida por
()k



(4.4)
assim, a força de corpo devido à tensão superficial do fluido resulta em
()()
kH
W

f
(4.5)
55
onde
H
é a função de Heaviside
00
120
10
se
Hse
se


(4.6)
e
W
é o número de Weber definido por
2
0
t
LU
W
s
(4.7)
onde
0
é a massa específica do fluido,
L
e
U
são o comprimento e a velocidade
característicos respectivamente e
t
s
é a tensão superficial.
Para fins de simulação de superfície-livre basta que a função distância esteja bem
constituída somente em uma região próxima da interface, normalmente 2 ou 3
comprimentos de em cada direção [58], [67].
No método das curvas de nível, a interpolação das propriedades do fluido para a
solução do escoamento se dá a partir das equações
10
1
HH


(4.8)
10
1
HH


(4.9)
Segundo SETHIAN [67], a função de Heaviside representada pela Equação (4.6)
apresenta resultados ruins por assumir que a interface possui espessura nula e
conseqüentemente a transição entre os fluidos se dá de forma abrupta. Uma alternativa à
Equação (4.6) seria atribuir uma espessura para a interface onde a transição de um
fluido para o outro fosse mais suave, assim:
56
0se
11
1se
2
1se
Hsen












(4.10)
onde
é a espessura da interface, geralmente assumida com valor igual a
h
, onde
h
é o comprimento da malha, e
1
permite que transições abruptas através da
interface sejam adequadamente suavizadas [67].
A principal dificuldade associada ao método de curvas de nível consiste na
construção de uma função distância para geometrias arbitrárias em malhas não
estruturadas. Além disso, o caráter da função distância é perdido conforme esta vai
sendo transportada o que exige que a mesma seja periodicamente reconstruída. É
importante observar que os métodos de curvas de nível eram tradicionalmente aplicados
em associação com métodos de malhas cartesianas e em computação gráfica (ver
referências no artigo do Apêndice C). Uma conseqüência disto é que somente nos
últimos anos começaram a surgir trabalhos na área de métodos baseados em malhas
não-estruturadas fazendo uso deste método e, mesmo assim, freqüentemente estes
trabalhos apresentam soluções somente para geometrias simples, ou seja, o problema de
se obter funções distância de forma eficiente, em malhas não-estruturadas e para
geometrias complexas ainda é um problema de pesquisa em aberto. Neste trabalho foi
proposto um método inteiramente pioneiro para a realização de tal cálculo. O método
proposto está sucintamente descrito na seção 4.4 e os resultados contendo a discussão
detalhada são apresentados no artigo do Apêndice C.
4.2. Correção de massa
A maior dificuldade para um bom método numérico de tratamento de superfície-livre
reside na habilidade de se preservar a massa das espécies. Nos métodos tipo VOF, a
perda de massa está associada a fatores como: difusão numérica da função degrau
utilizada como marcador, campo de velocidade numérico não atender ao critério de
divergência nula, ondulações na advecção da função de marcação. Nos métodos de
curva de nível, a perda de massa geralmente ocorre por perda da propriedade da função
57
distância a qual deve ser reconstruída periodicamente. No presente estudo, o algoritmo
de preservação de massa adotado por LOHNER em [46] foi empregado para garantir a
preservação global de massa. Neste algoritmo, as perdas ou ganhos materiais são
obtidos a partir da comparação entre a quantidade inicialmente existente, adicionada ou
subtraída da quantidade de matéria injetada ou retirada a cada passo de tempo. A
quantidade de matéria residual é então distribuída pelos elementos pertencentes à
interface entre os dois fluidos a partir da seguinte relação
n
v

u
(4.11)
A quantidade calculada a partir da Equação (4.11) garante que a massa seja
preservada principalmente nas regiões onde a interface se move com maior velocidade
enquanto mantém as regiões de baixa mobilidade da solução praticamente intactas.
Detalhes acerca da implementação em paralelo deste algoritmo bem como testes de
eficiência na preservação de matéria são discutidos no artigo do Apêndice D (seção
3.1).
4.3. Integração de volume para as fases
O sucesso do algoritmo de correção de massa descrito na seção anterior depende
diretamente da precisão com que se determina a massa/volume das fases envolvidas.
Para este tipo de cálculo, não basta integrar o volume dos elementos baseados na média
no baricentro das frações mássicas nodais. Neste estudo, foi criado um método baseado
em iso-contornos para a determinação precisa do volume das fases. Este método está
detalhadamente descrito no artigo do Apêndice D (seção 3.1).
4.4. Cálculo de funções-distância
O método de captura de interface, utilizando funções marcadoras do tipo curvas de nível
(level sets), necessitam da geração e reconstrução periódica de campos de distância
onde a menor distância entre um ponto do domínio é obtida em relação a uma superfície
58
(em 3D) ou linha (em 2D) representando a interface entre dois fluidos. O algoritmo
mais simples para construção de tais campos seria um algoritmo de força bruta, onde
para cada nó fora da interface calcula-se a distância em relação a todos os nós da
superfície e armazena-se o menor valor. Entretanto, este algoritmo para problemas
tridimensionais em malhas não-estruturadas possui um custo computacional inviável.
Existem poucos métodos para cálculo de funções distância disponíveis na literatura e os
que existem são de implementação complexa requerendo a utilização de estruturas de
dados em árvores binárias e listas. No presente trabalho foi dada prioridade para o
desenvolvimento de um algoritmo inteiramente novo destinado ao cálculo de funções
distância em malhas não-estruturadas utilizando o próprio método dos elementos finitos
como base.
A geração de funções distância pode ser resumida na solução da equação de
Eikonal:
1
TF

(4.12)
onde
F
é a velocidade de propagação da frente de cálculo e
T
é o tempo de viagem da
frente. A solução da Equação (4.12) pode ser interpretada como um balão sendo inflado
onde dado um ponto
p
qualquer no espaço, o tempo de viagem
()
T
p
é o tempo que a
superfície do balão leva para alcançar este ponto enquanto inflado. Note que se
1
F
,
a frente (superfície do balão) se move com velocidade unitária e o tempo de chegada ao
ponto
p
equivale a distância percorrida pela frente desde o instante
0
t
e a Equação
(4.12) torna-se
1
d

(4.13)
O método proposto surgiu da observação de que, para um elemento finito triangular,
dados dois nós conhecidos, o terceiro pode ser calculado resolvendo-se a Equação
(4.13) no nível do elemento. A equação resultante é uma equação do segundo grau e o
esforço computacional desta solução consiste em resolver-se uma equação do segundo
grau para cada elemento que possuir um nó de incógnita. Além do cálculo também foi
proposto um algoritmo para percorrer a malha de modo que os elementos fossem sendo
59
adicionados e eliminados do cálculo conforme a disponibilidade de informações.
Detalhes acerca deste estudo, para malhas compostas por tetraedros, estão disponíveis
no artigo do Apêndice C.
4.5. Desativação Dinâmica em Paralelo
O algoritmo de desativação dinâmica é outro dispositivo importante na redução do custo
computacional da solução de problemas puramente advectivos, como o transporte da
função de marcação para escoamentos de superfície-livre. O algoritmo de desativação
dinâmica desenvolvido no presente trabalho é uma extensão da desativação dinâmica
apresentada por LOHNER e CAMELI [44] para o transporte de contaminantes, todavia,
aplicada ao tratamento em paralelo de problemas de superfície-livre. Em [73] SOUZA
et al utilizam uma técnica semelhante aplicada à problemas de transporte. Esta técnica
baseia-se na que foi inicialmente proposta por TEZDUYAR e LIOU [78] e consiste na
seleção de partes do domínio que podem ser resolvidas utilizando integração temporal
implícita e explícita. Conseqüentemente, a solução de um sistema de equações só se faz
necessária para os elementos considerados como implícitos enquanto a partição
explícita pode ser resolvida utilizando soluções de passos anteriores [73]. Esta técnica
foi chamada por seus autores de partição implícita/explícita.
Na solução de problemas de superfície-livre a função de marcação deve ser
transportada por uma equação hiperbólica (Equação (4.1)). Com o algoritmo de
desativação dinâmica os graus de liberdade ativos para o sistema de equações estarão
restritos somente àquelas entidades de malha que circundam a região da interface. Em
[44], LOHNER e CAMELI aplicaram a seguinte equação para selecionar as entidades
ativas no domínio:
e
m


(4.14)
onde
e
é a norma Euclideana do gradiente da solução por elemento,
m
é o
gradiente médio da solução calculada para todo o domínio e
0,1


é um parâmetro
60
que controla a sensibilidade na seleção dos elementos. Neste trabalho o critério de
seleção utilizando a Equação (4.14) foi mantido somente para problemas de transporte
advectivo-difusivo, entretanto, para problemas de superfície-livre este critério foi
substituído por um método de seleção explícita onde os elementos pertencentes à
interface entre os dois fluidos são automaticamente selecionados como ativos. Além da
seleção efetuada na primeira passada do algoritmo um número adicional de camadas é
selecionado para garantir que a solução seja mantida dentro da partição ativa a cada
passo de tempo. Posteriormente, a seleção das entidades ativas (nós, equações, arestas e
elementos) é utilizada como restrição para os laços mais computacionalmente intensos
que são, normalmente, a construção da matriz de rigidez, vetor de resíduos e produto
matriz-vetor. Assim, enquanto o laço original descrito abaixo,
do ie=1,ntent
! ... cálculos intensos utilizando a entidade ie...
enddo
seria executado ntent vezes, onde ntent é o número de entidade para todo o domínio,
com a desativação dinâmica o laço é substituído por:
do ie=1,nae
ia = EntidadeAtiva(ie)
! ... cálculos intensos utilizando a entidade ia
enddo
e passa a ser executado nae vezes, sendo nae o número de entidade ativas (geralmente
muito menor que ntent).
A Figura 4.5 ilustra o funcionamento do algoritmo de desativação dinâmica em uma
simulação de escoamento de superfície-livre paralela utilizando dois processadores.
61
Figura 4.5 – Desativação dinâmica em paralelo utilizando dois processadores
Note que a desativação dinâmica operando com paralelismo distribuído pode gerar
desbalanceamento de carga entre os processos. Entretanto, cabe salientar que os ganhos
obtidos na desativação das partes irrelevantes do problema poderiam ser degradados
caso houvesse a tentativa de recuperação do balanço de esforço computacional, devido
às operações envolvidas.
É importante observar também que, no contexto de superfície-livre, esta mesma
técnica poderia ser estendida para restringir os cálculos referentes à solução das
equações de Navier-Stokes somente na região de interesse, normalmente a fase líquida.
Entretanto, para que isto seja possível, algumas técnicas de projeção da solução da fase
líquida para a fase gasosa, como a proposta por LOHNER et al em [46], são requeridas
para manter a consistência nos cálculos.
Os testes de desempenho deste algoritmo assim como outros detalhes a respeito da
implementação desta técnica estão disponíveis no artigo do Apêndice D, seção 3.2.
62
Capítulo 5
Estruturas de dados por arestas
Neste capítulo é apresentada a estrutura de dados utilizada no armazenamento dos
coeficientes da matriz oriunda da aplicação do método dos elementos finitos. Como
visto nos capítulos anteriores, a solução de elementos finitos para as equações de
Navier-Stokes incompressíveis e de advecção-difusão, via integração no tempo
implícita, resulta na solução de um problema não-linear para cada passo de tempo. O
método de Newton inexato é utilizado como técnica de aceleração da solução não-
linear, entretanto, embora esta técnica reduza o esforço computacional através do
controle da tolerância linear, deve-se considerar que um número maior de iterações não-
lineares é requerido para se atingir o mesmo critério de parada que em um método de
Newton exato atingiria com poucas iterações de custo computacional bastante elevado
[15], [14]. Para a série de problemas lineares gerados é utilizado o método iterativo
GMRES empregando pré-condicionadores bloco diagonal e diagonal para o escoamento
e o transporte respectivamente. Conforme visto anteriormente, o método GMRES
necessita realizar diversas operações de multiplicação matriz-vetor durante o
procedimento de solução. Geralmente, grande parte do esforço computacional
consumido por uma simulação se concentra na realização desta operação.
Conseqüentemente, a redução dos custos computacionais associados à realização de
produtos matriz-vetor torna-se algo desejável e imprescindível na busca por
implementações de procedimentos eficientes. A redução de custos computacionais
referentes ao consumo de memória, número de vezes em que as informações necessitam
ser buscadas e/ou gravadas na memória (endereçamentos indiretos) e quantidade de
operações de ponto flutuante (do inglês, floating point operations ou simplesmente flop)
63
estão diretamente relacionadas com a forma com que a matriz é armazenada e,
conseqüentemente, utilizada durante os acessos. Nas próximas seções a estrutura de
dados aresta-por-aresta é discutida como uma alternativa eficiente na redução destes
custos computacionais.
5.1. Grafos, malhas e matrizes esparsas
Uma particularidade dos métodos discretos é que eles originam matrizes com alto grau
de esparsidade, onde esparsidade está definida, vagamente, em relação ao número
reduzido de elementos não-nulos presentes na matriz [64]. A forma mais ingênua e
didática de armazenamento de uma matriz seria na sua forma bidimensional, onde cada
coeficiente não-nulo é armazenado na posição correspondente aos índices da linha e
coluna deste arranjo. Entretanto, uma grande quantidade de memória seria inutilizada
para se armazenar poucos coeficientes não-nulos, o que não se justifica para matrizes
esparsas. Assim sendo, a utilização de estruturas de dados que tirem proveito da
esparsidade das matrizes torna-se importante na construção de métodos eficientes.
A teoria de grafos é freqüentemente utilizada no estudo e definição de estruturas de
dados. Um grafo é definido com base em dois conjuntos [64], um conjunto de vértices
12
,,...,
n
Vvvv
(5.1)
e um conjunto de arestas
E
que define os pares
(,)
ij
vv
, onde
i
v
,
j
v
são elementos de
V
, isto é,
EVV

. Conseqüentemente, o grafo
(,)
GVE
é uma forma de
representar relações binárias entre objetos do conjunto
V
e pode ser utilizado na
representação de matrizes esparsas onde os vértices de
V
representam as incógnitas da
matriz
A
. Desta forma, haverá uma aresta unindo o nó
i
a incógnita
j
para cada
coeficiente de
0
ij
a
conforme ilustrado na Figura 5.1, onde os símbolos ×
representam termos não-nulos. A teoria de grafos ainda classifica as formas como os
vértices do conjunto
V
se relacionam. Desta forma, para matrizes com perfil simétrico
de coeficientes não-nulos os grafos são classificados como grafos não-direcionados (ou
64
seja, não importa a direção em que a aresta ligando os vértices
i
v
e
j
v
foi construída).
Nestes casos para todo
0
ij
a
haverá um termo não-nulo, igual ou não, na posição
ji
a
conforme mostra a Figura 5.1.
Figura 5.1 – Grafo de uma matriz 5 × 5 e sua matriz correspondente.
Uma conseqüência dos conceitos apresentados até o momento é que existe uma relação
direta entre malhas computacionais, matrizes e grafos [11], [37]. De fato, as malhas
computacionais podem ser representadas na forma de grafos. Para isso, os conceitos de
grafo nodal e grafo dual devem ser apresentados primeiramente. Em cálculos realizados
percorrendo-se os nós da malha, haverá um vértice relacionado a cada nó da malha e
uma aresta do grafo correspondente a cada aresta da malha e este grafo é classificado
como grafo nodal. Nos casos em que os cálculos são realizados percorrendo-se os
elementos da malha, haverá um vértice correspondente a cada elemento e uma aresta do
grafo correspondente aos elementos vizinhos (arestas em 2D e faces em 3D) [11], [37].
Para topologias simples, como as de triângulos e tetraedros, o grafo nodal corresponde à
própria malha formada por estes elementos
5
[64]. A Figura 5.2 ilustra as relações entre
os grafos nodal (b) e dual (c), assim como, a matriz (d) de coeficientes não nulos
correspondente ao grafo nodal de uma malha de triângulos (a). Na Figura 5.2 a
numeração dos elementos foi substituída por letras para facilitar a compreensão.
5
No contexto de elementos finitos, a correspondência entre malha e grafo nodal está relacionada somente
aos elementos triangulares e tetraédricos que utilizam funções de interpolação lineares.
1
2
4
3
5
12345
1
2
3
4
5
















65
(a) Malha
(b) Grafo nodal
(c) Grafo dual
123456
1
2
3
4
5
6



















(d) Matriz de coeficientes não-nulos
Figura 5.2 – Relação entre malha, grafos e matrizes.
A estrutura de dados aresta-por-aresta tira proveito das relações ilustradas na Figura 5.2.
Na próxima seção é discutida a utilização desta estrutura de dados, utilizada no escopo
deste trabalho, para reduzir os custos computacionais relacionados às operações de
multiplicação matriz-vetor mencionadas anteriormente.
5.2. Incidências de elemento e aresta
As informações relacionais apresentadas até o momento são implementadas na forma de
arranjos de incidência ou conectividade (do inglês, connectivity). Tais arranjos podem
se apresentar de diversas formas de acordo com a conveniência de seu uso. Para a
utilização da estrutura de dados por arestas, neste trabalho foram adotados os seguintes
arranjos de incidência:
IEN(nnoel,nel): conectividade dos nós que formam o elemento;
IEDGE(2,nare) conectividade dos nós que formam a aresta;
66
IEEI(nae,nel): conectividade das arestas que formam o elemento;
onde nnoel é o número de nós por elemento, nel é o número de elementos, nare é o
número de arestas e nae é o número de arestas por elemento. Para exemplificar a
construção destes arranjos, considere a Figura 5.3.
Figura 5.3 – Construção das incidências de elemento e aresta.
os arranjos IEN, IEDGE e IEEI para a malha da Figura 5.3 correspondem a
IEN:
ie
1
2
3
ien(1,ie)
:
1
4
5
ien(2,ie)
:
4
2
2
ien(3,ie)
:
3
3
4
IEDGE:
ia
1
2
3
4
5
6
7
iedge(1,ia)
:
1
3
4
2
3
4
4
iedge(2,ia)
:
4
1
3
5
2
5
2
IEEI:
ie
1 2 3
ieei(1,ie)
1 6 7
ieei(2,ie)
3 -4 -5
ieei(3,ie)
2 -7 -3
Note que no arranjo IEEI foi convencionado atribuir valor negativo para as arestas que
não possuem incidência nodal correspondente àquela apresentada pelo elemento. Na
prática os arranjos IEN e IEEI são utilizados na construção das matrizes de aresta
67
enquanto o arranjo IEDGE é empregado na construção do arranjo de localização das
equações (arranjo LM) a serem utilizadas nas operações de multiplicação matriz-vetor.
Estas operações serão apresentadas nas próximas seções.
É importante salientar que as informações referentes aos arranjos IEDGE e IEEI
não são normalmente disponibilizadas pelos softwares de geração de malha.
Conseqüentemente, estas informações devem ser construídas para que a estrutura de
dados possa ser utilizada. Uma forma eficiente de obter estes arranjos envolve a
aplicação de tabelas hash e listas encadeadas. Neste trabalho estas informações, assim
como outros tratamentos específicos à estrutura de dados por aresta, foram obtidos
através do uso da biblioteca EdgePack desenvolvida por MARTINS [49], [50].
5.3. Construção das matrizes de aresta
Tomando-se como exemplo o elemento “A” da Figura 5.2a, nota-se que este elemento é
formado pelas arestas 1-3, 3-6, 6-1 do grafo nodal (Figura 5.2b). Desta forma, as
matrizes de aresta podem ser construídas a partir do desmembramento dos coeficientes
da matriz de rigidez do elemento
6
correspondente às arestas que o compõem [49].
Conseqüentemente, a matriz de aresta para um grau de liberdade por nó, será, para
qualquer elemento (já que qualquer aresta invariavelmente é criada a partir da união de
dois vértices do grafo), dada por
ee
ijij
e
ij
ee
ijij
T





kk
kk
(5.2)
onde
e
ij
k
é o coeficiente correspondente à aresta
ij
do elemento
e
. Note que para mais
de um grau de liberdade,
e
ij
k
corresponde a blocos de matrizes de dimensão
nglngl
,
onde
ngl
é o número de graus de liberdade por nó.
6
O termo matriz de rigidez de elemento segue a convenção adotada no método dos elementos finitos e
sua definição formal está disponível em HUGHES [32] pag. 40.
68
Na Figura 5.2a, a aresta 3-6 é compartilhada pelos elementos “A” e “B”,
diferentemente dos elementos, as matrizes de arestas são, na sua maioria, formadas pela
contribuição de mais de um elemento (uma exceção a esta regra se dá para as arestas do
contorno do domínio como, por exemplo, a aresta 1-3 da Figura 5.2a), assim, torna-se
conveniente formular as matrizes das arestas com base nos seus elementos concorrentes,
ou seja
11
11
nn
ee
ijij
a
ee
ij
nn
ee
ijij
ee
T













kk
kk
(5.3)
onde
n
é o número total de elementos que compartilham a aresta
ij
e
e
ij
k
é a submatriz
do elemento
e
correspondente aos graus de liberdade dos nós
i
e
j
da aresta
a
. A
Figura 5.4 ilustra a contribuição de 6 elementos para os coeficientes de uma única aresta
(aresta a
1
destacada em vermelho).
Figura 5.4 – Contribuições dos elementos para a matriz da aresta a
1
69
5.3.1. Elementos tetraédricos lineares
Para elementos tetraédricos lineares a formulação estabilizada do presente trabalho dá
origem a matrizes não simétricas, conseqüentemente, o desmembramento em suas
matrizes de aresta é conforme ilustrado na Figura 5.5 adiante.
11
22
3
3
13
31
23
32
24
42
12
2
4
43
3
44
4
1
1
41
e









k
k
k
k
k
k
k
k
k
k
k
k
K
kk
k
k
Matriz de rigidez do elemento
onde
ij
nglnnoelnglnnoel




k
11
22
12
21
e
A





k
T
k
k
k
Matriz da aresta
A
22
33
23
32
e
B





k
T
k
k
k
Matriz da aresta
B
11
33
13
31
e
C





k
T
k
k
k
Matriz da aresta
C
14
41
11
44
e
D





k
T
k
k
k
Matriz da aresta
D
24
42
22
44
e
E





k
T
k
k
k
Matriz da aresta
E
34
43
33
44
e
F





k
T
k
k
k
Matriz da aresta
F
Figura 5.5 – Desmembramento da matriz de rigidez do elemento tetraédrico em suas
matrizes de aresta correspondentes.
Na Figura 5.5 observe que a construção da matriz de rigidez do elemento, assim como
das matrizes de aresta correspondentes, segue as incidências nodais para o elemento e
aresta respectivamente. A montagem das contribuições de todos os elementos que
compartilham de uma mesma aresta se dá conforme apresentado na Equação (5.3). Note
que as arestas possuem sentido convencionado para suas incidências (indicado na
Figura 5.5 por setas), desta forma, a matriz da aresta 1-3, por exemplo, corresponde à
transposta da matriz da aresta 3-1. Esta informação é incorporada ao arranjo IEEI
através do sinal da aresta conforme apresentado na seção 5.2. Além disso, observe que
70
os termos da diagonal da matriz de rigidez do elemento representam a parcela de
contribuição deste elemento para o nó correspondente. Tal característica pode ser
explorada na construção do pré-condicionador bloco diagonal nodal (ou diagonal
dependendo-se do sistema a ser tratado) conforme apresentado na seção 3.6. Na prática,
estes termos são armazenados separadamente das matrizes de aresta para fins de
economia de memória conforme será mostrado mais adiante. O Algoritmo 5.1 apresenta
o laço básico de construção das matrizes de aresta.
Algoritmo 5.1 – Construção das matrizes de aresta
Para cada elemento fazer:
Passo 1: Recuperar incidência nodal do elemento (arranjo ien)
Passo 2: Calcular coeficientes da matriz de rigidez de elemento
Passo 3: Recuperar incidência de aresta do elemento (arranjo ieei)
Passo 4: acumular contribuições para cada aresta
Fim fazer
Note que no Algoritmo 5.1 os passos 1 e 2 são comuns a qualquer programa de
elementos finitos [32].
5.4. Produto matriz-vetor por aresta
As operações de multiplicação matriz-vetor são freqüentemente requisitadas durante o
processo de solução linear iterativa conforme visto na seção 3.5.1. Uma das principais
influências no uso da estrutura de dados por arestas se dá na aceleração desta operação
[49]. Nesta estrutura de dados, o produto matriz-vetor é realizado localmente no nível
das arestas, seguindo um procedimento similar ao que se observa para a estrutura de
dados elemento-por-elemento
7
, e os resultados locais são, então, combinados ao global.
Conseqüentemente, nestas duas estruturas de dados a matriz de rigidez global nunca é
explicitamente formada e o efeito global só é obtido após a realização do produto
7
Na estrutura de dados elemento-por-elemento as matrizes de rigidez de elemento são armazenadas
individualmente. Maiores detalhes acerca desta estrutura de dados podem ser obtidos em HUGHES [32].
71
matriz-vetor. O produto matriz-vetor aresta por aresta é obtido a partir da seguinte
expressão
1
nare
aa
a

zApAp
(5.4)
onde nare é o número de arestas da malha,
z
é o vetor global que armazena a solução,
a
p
são as componentes de
p
restritas aos graus de liberdade da aresta e
a
A
é a matriz
da aresta
a
armazenada em um arranjo
(2,)
nglnglnare

A
para as equações de
Navier-Stokes e
(2,)
nare
A
para as equações de transporte. Nestes arranjos não estão
sendo considerados os termos bloco-diagonal e diagonal, pois os mesmos são
armazenados separadamente em arranjos
(,)
nglnnos
B
e
()
nnos
B
para o escoamento e
transporte respectivamente. Note que basta alocar espaço de armazenamento para os
arranjos correspondentes ao problema de escoamento e reutilizar a área para o problema
de transporte. O Algoritmo 5.2 apresenta o produto matriz-vetor aresta por aresta para
um grau de liberdade por nó
Algoritmo 5.2 – Produto matriz-vetor aresta por aresta
Para cada aresta
a
fazer
Passo 1: Localize as componentes de
a
p
no vetor global
p
Passo 2: Realize o produto matriz-vetor local:
aaa
zAp
Passo 3: Espalhe o resultado local
a
z
no vetor global
z
Fim do laço
No Algoritmo 5.2 a localização dos coeficientes de
a
p
, invocada no Passo 1, é feita
com o uso do arranjo de localização LM(2×ngl, nare). Este arranjo, presente na maioria
das implementações de elementos finitos [32], mapeia os graus de liberdade da aresta no
sistema de equações global e possui relação direta com a incidência da aresta (arranjo
IEDGE). Por convenção, se LM(i, j) for igual a zero, o grau de liberdade é eliminado da
solução do sistema. A Figura 5.6 ilustra o produto matriz vetor aresta por aresta.
72
Figura 5.6 – Produto matriz-vetor aresta-por-aresta.
Conforme citado anteriormente, neste trabalho o armazenamento dos termos da
diagonal é feito de forma separada dos termos fora da diagonal. Estes termos são
armazenados em um arranjo global que combina a influência de todos os elementos (e
conseqüentemente das arestas) incidentes sobre um mesmo nó. O motivo da utilização
de tal esquema reside no aumento da eficiência do produto matriz vetor que, agora,
passa a ser realizado em mais estágios do que aquele apresentado no Algoritmo 5.2. Tal
técnica foi proposta inicialmente por van GIJZEN [85] para os termos da diagonal
principal e que aqui foi estendido para o bloco diagonal, a exemplo do que foi feito por
CATABRIGA [7] na solução de problemas de escoamentos compressíveis. Assim, o
produto matriz-vetor expresso pela Equação (5.4) pode ser reescrito como:
1
nare
aa
a

zBpAp
%
(5.5)
onde
B
é o arranjo bloco diagonal nodal global para problemas de escoamento e
diagonal para problemas de transporte e
a
A
%
é a matriz de aresta desprezando-se o termo
bloco diagonal (ou diagonal) local, ou seja
73
11121112
21222221
00
00
aaa
aaaa
aaaa






AAB
%
(5.6)
Note que a primeira parcela da Equação (5.5) consiste em uma multiplicação global em
nível nodal enquanto a segunda parcela é similar ao produto matriz-vetor convencional,
em nível de aresta (Equação (5.4)) operando-se somente nos termos fora da diagonal.
Note que a segunda parcela, geralmente, requer um esforço computacional maior, pois,
em malhas não estruturadas, o número de arestas que compartilham de um mesmo nó
pode variar resultando um número maior de arestas do que o de nós. Assim sendo, os
efeitos desta modificação normalmente se traduzem em:
− Menor consumo de memória no armazenamento dos coeficientes da matriz de
rigidez, pois os termos da diagonal que antes seriam armazenados por aresta
agora são armazenados por nó;
− Redução do número de operações de ponto flutuante no produto matriz vetor por
arestas, pois a cada iteração desta operação deixa-se de operar nos termos da
diagonal, operando-os mais tarde em um laço nodal;
− Redução do número de endereçamentos indiretos por motivo similar ao
mencionado anteriormente.
O Algoritmo 5.3 apresenta as modificações necessárias para a implementação do
produto matriz-vetor no esquema Gizjen.
Algoritmo 5.3 – Produto matriz-vetor aresta por aresta no esquema Gizjen.
Para cada aresta
a
fazer
Passo 1: Localize as componentes de
a
p
no vetor global
p
Passo 2: Realize o produto da matriz-vetor local:
aaa
zAp
Passo 3: Espalhe o resultado local
a
z
no vetor global
z
Fim do laço
Para cada nó
n
fazer
74
Passo 4: Multiplique
n
B
por
n
p
Passo 5: Adicione o resultado a
n
z
, ou seja:
nnnn

zzBp
Fim do laço
5.5. Custo computacional
O principal ganho na utilização da estrutura de dados por aresta se reflete na redução
dos custos computacionais referentes ao uso de memória, número de operações de ponto
flutuante e quantidade de endereçamentos indiretos. Ganhos adicionais podem ser
obtidos com a aplicação de técnicas de ordenação como as que serão vistas na próxima
seção.
Em [62] RIBEIRO e COUTINHO fazem uma comparação detalhada dos custos
computacionais relacionados a três estruturas de dados freqüentemente utilizadas em
métodos discretos, incluindo-se a estrutura aresta-por-aresta, conforme adotada neste
trabalho. A Figura 5.7 apresenta a compilação dos custos computacionais de
armazenamento (Figura 5.7a), número de operações de ponto flutuante (Figura 5.7b) e
quantidade de endereçamento indireto (Figura 5.7c) descritas na referência [62]. Os
valores apresentados na Figura 5.7 consideram o caso regular onde os coeficientes da
diagonal são armazenados junto dos termos fora da diagonal.
75
10
100
1000
10000
1 2 3 4 5 6
ndof
bytes
EBE
EDS
CSR
(a) consumo de memória
10
100
1000
10000
1 2 3 4 5 6
ndof
flops
EBE
EDS
CSR
(b) operações de ponto flutuante
10
100
1000
1 2 3 4 5 6
ndof
i/a
EBE
EDS
CSR
(c) endereçamentos indiretos
Figura 5.7 – Tetraedro linear - Custos computacionais para as estruturas aresta-por-
aresta (EDS), elemento-por-elemento (EBE) e compressed sparse row (CSR).
Nota-se, pela Figura 5.7, que a estrutura de dados aresta-por-aresta (EDS), em geral, se
mostra superior à estrutura elemento-por-elemento (EBE). Para problemas com mais de
3 graus de liberdade por nó a estrutura EDS é a mais eficiente de todas as outras listadas
na comparação ao utilizar um menor espaço de armazenamento e necessitando de uma
quantidade menor de endereçamentos indiretos do que a estrutura Compressed Sparse
Row (CSR). Cabe ressaltar que para solução de escoamentos incompressíveis, conforme
adotada neste trabalho, considera-se o acoplamento entre velocidade e pressão e
conseqüentemente resolve-se 4 graus de liberdade simultaneamente. Para manter uma
única estrutura de dados no código, resolve-se o transporte usando-se tambem a
estrutura de dados aresta-por-aresta, apesar de esta ser menos eficiente do que a CSR
para um grau de liberdade por nó.
76
5.6. Técnicas de agrupamento e ordenação de dados
Os ganhos obtidos com a adoção da estrutura de dados aresta-por-aresta podem ser
estendidos ainda mais com a utilização de técnicas de agrupamento e ordenação de
dados que tirem proveito da arquitetura onde o software irá realizar os cálculos. Em
particular, o esquema de ordenação ótimo deve explorar a hierarquia de memória da
máquina de forma eficiente e, conseqüentemente, levar em consideração o fluxo de
acesso aos dados (Apêndice B). Diversas são as possibilidades na aplicação desta
técnica e, em geral, não existe um esquema ótimo e único de dispor os dados que
contemple, de maneira eficiente, todas as arquiteturas de computadores existentes
(Apêndice B).
Alguns esquemas são de simples compreensão e válidos como estratégia de
otimização mínima. Considere, por exemplo, as ordenações apresentadas na Figura 5.8.
Nesta figura, os esquemas de ordenação em relação às arestas e nós da malha são
apresentados a esquerda enquanto os padrões de acesso à memória são mostrados a
direita para a matriz global e vetores equivalentes. Tais padrões de acesso correspondem
aos laços de aresta que são freqüentemente executados durante a realização de produtos
matriz-vetor. Note que o acesso aos termos diagonais foi omitido tornando a ilustração
similar ao que ocorre na primeira parte do produto matriz-vetor para o esquema Gizjen,
descrito no Algoritmo 5.3. Para melhor compreensão deste exemplo considere os
seguintes aspectos:
− Os dados envolvidos nos cálculos estão armazenados em áreas contíguas de
memória, caso típico para armazenamento de dados em arrays
8
;
− O processador armazena os dados em memórias velozes contidas no interior do
núcleo de processamento, conhecidas como registradores e memória cache [11],
[74], este armazenamento temporário é realizado com intuito de acelerar o
acesso aos dados;
8
O uso de listas encadeadas, estruturas compostas e listas de objetos representam exemplos de esquemas
de armazenamento de dados em áreas não contíguas de memória e que devem ser evitados para
processamento intensivo [74].
77
− O processador se antecipa na leitura dos dados esperando que os mesmos sejam
utilizados em ciclos seguintes, mecanismo este conhecido por pré-carga
(prefetch em inglês [11], [74]).
Neste exemplo considere também o cenário em que o processador possui a capacidade
de realizar a pré-carga dos coeficientes de 3 nós e que a área consumida pela pré-carga
correspondente ao tamanho da memória cache do processador. A pré-carga e a área da
memória cache estão representadas na Figura 5.8 por retângulos tracejados em
vermelho. Note que a velocidade de execução das iterações estará relacionada à taxa de
acerto na pré-carga, representada pela área do retângulo, juntamente com a reutilização
da memória cache, representada pelo deslocamento do mesmo. Conseqüentemente,
todas as vezes que ocorre um erro na pré-carga (acessos fora da área do retângulo) o
processador deverá retornar à memória principal para buscar a informação necessária,
neste caso ocorre um erro de cache, do inglês, cache miss [11], [74]. A ocorrência de
erros de cache consome vários ciclos de processamento e degrada de forma substancial
o desempenho [11]. Note também, que para malhas fixas, um erro de cache tende a se
propagar durante toda a simulação já que o padrão de acesso à memória está relacionado
ao grafo da malha e conseqüentemente a ordem em que os nós e arestas se relacionam.
Com base nestas premissas, pode-se concluir, do exemplo apresentado na Figura 5.8,
que:
− Figura 5.8a: A malha original apresenta ordenação ruim por apresentar 3 erros
de cache (iterações 1, 5 e 9) e baixa taxa de reutilização da cache;
− Figura 5.8b: A ordenação de aresta permite a reutilização de dados da cache para
algumas iterações, entretanto, apresenta 3 erros de cache (iterações 2, 3 e 8);
− Figura 5.8c: A ordenação nodal não apresenta nenhum erro de cache, entretanto
a taxa de reutilização é baixa;
− Figura 5.8d: A associação das técnicas de ordenação das arestas e dos nós
alcança o nível de acesso mais eficiente com perfeito aproveitamento das pré-
cargas e com boa taxa de reutilização da cache.
78
(a) Malha original
(b) Ordenação de aresta
(c) Ordenação nodal
(d) Ordenação de aresta e nodal
Figura 5.8 – Exemplo de ordenações
O estudo apresentado por COUTINHO et al no Apêndice B propõe um método
heurístico de determinação do arranjo ideal de acordo com a arquitetura em questão. Na
metodologia proposta por estes autores, as operações de custo computacional mais
12345
25
79
241
,,,,
783
94
123
86
5136
456
1
2
3
4
5
6


















iterações
6789
,,,,
















644444444444444474444444444444448
12345
25
214
5163
,,,,
4689
38
1234
7
97
56
1
2
3
4
5
6


















iterações
6789
,,,,
















644444444444444474444444444444448
12345
12
134
2356
,,,,
457
123456
1
2
8
65 79
89
3
4
6

















iterações
6789
,,,,
















644444444444444474444444444444448
12345
12
98
143
,,,,
976
84
123456
1
2
3
4
5
6
75
2365

















iterações
6789
,,,,
















644444444444444474444444444444448
79
elevado, em uma simulação de elementos finitos, são ensaiadas um número determinado
de vezes para diversas combinações de esquemas de ordenação e o melhor esquema é,
então, determinado. Normalmente as operações mais intensivas são os laços de
construção da matriz de rigidez e vetor de resíduos juntamente com o produto matriz-
vetor. No presente trabalho, as técnicas descritas nas referências [8] (Apêndice B) e [50]
são aplicadas aos dados da simulação em uma fase de pré-processamento através do uso
do software EdgePack [50]. As técnicas de ordenação e agrupamento desenvolvidas
nestas referências são apresentadas na Figura 5.9. Nesta figura, os esquemas de
ordenação são apresentados na forma de gráficos onde as linhas verdes e azuis
representam a numeração dos nós
i
e
j
de uma aresta respectivamente. Na Figura 5.9a
é apresentada a numeração original (similar a Figura 5.8a), sem qualquer tipo de
reordenação e, geralmente, retornada pelo software de geração de malhas e as demais
Figuras 5.9 (b)-(e) demonstram o perfil da numeração dos nós de acordo com as
seguintes características:
(b) Minimização de endereçamentos indiretos. Estes esquema é recomendado para
problemas que tratam muitos graus de liberdade simultaneamente como o que
foi tratado neste trabalho para o escoamento incompressível;
(c) Aumento da localidade de dados levando-se em consideração a dependência de
memória. Neste esquema, o tratamento da dependência de memória se faz
necessária para uso em máquinas vetoriais e/ou para uso de processamento
paralelo por multitarefas (OpenMP). Estes requisitos serão abordados no
próximo capítulo;
(d) Aumento da localidade de dados desprezando-se a dependência de dados. Este
caso é o mesmo ilustrado na Figura 5.8d;
(e) Aumento da localidade de dados otimizando-se o uso dos registradores da CPU;
80
0
20
40
60
80
100
0 100 200 300 400 500 600 700
x 1000
x 1000
Edges
Nodes
Node i
Node j
0
20
40
60
80
100
0 100 200 300 400 500 600 700
x 1000
x 1000
Edges
Nodes
Node i
Node j
(a)
0
20
40
60
80
100
0 100 200 300 400 500 600 700
x 1000
x 1000
Edges
Nodes
Node i
Node j
(b)
0
20
40
60
80
100
0 100 200 300 400 500 600 700
x 1000
x 1000
Edges
Nodes
Node i
Node j
0
20
40
60
80
100
0 100 200 300 400 500 600 700
x 1000
x 1000
Edges
Nodes
Node i
Node j
(c)
0
20
40
60
80
100
0 100 200 300 400 500 600 700
x 1000
x 1000
Edges
Nodes
Node i
Node j
(d)
0
20
40
60
80
100
0 100 200 300 400 500 600 700
x 1000
x 1000
Edges
Nodes
Node i
Node j
0
20
40
60
80
100
0 100 200 300 400 500 600 700
x 1000
x 1000
Edges
Nodes
Node i
Node j
(e)
Figura 5.9 – Técnicas de reordenação realizadas pela EdgePack.
81
Capítulo 6
Técnicas computacionais
A qualidade dos resultados obtidos com o método dos elementos finitos está
normalmente relacionada ao número de graus de liberdade utilizados na análise.
Conseqüentemente, o uso eficiente de processamento de alto desempenho (PAD), do
inglês, High Performance Computing ou HPC, torna-se fundamental para viabilizar
cálculos cada vez mais precisos e ricos em detalhes obtidos em um período de análise
cada vez menor. A seguir serão apresentadas as técnicas de processamento de alto
desempenho e as características do software desenvolvido neste trabalho. Tais técnicas
abordam paralelismo distribuído, compartilhado e híbrido, vetorização, pipelining e,
finalizando, com as técnicas de pós-processamento dos resultados das simulações.
6.1. Computação paralela
Um computador paralelo é, por definição, um computador contendo mais de um
processador [11]. No passado, cada processador de um sistema multi-processado era
disposto isoladamente, entretanto, recentemente, com a introdução dos processadores
multicores, vários processadores residem em uma mesma unidade lógica. Existem
diversos tipos de sistemas paralelos, estes normalmente são diferenciados de acordo
com a interconexão entre os processadores (também conhecidos por elementos de
processamento ou PEs) e memória [24]. Um modelo amplamente difundido, a
Taxonomia de Flynn, classifica os computadores paralelos (e seriais) de acordo com o
número de instruções executadas ao mesmo tempo, assim, sistemas classificados como
82
SIMD (Single Instruction/Multiple Data) aplicam uma mesma instrução a vários dados
simultaneamente enquanto os sistemas MIMD (Multiple Instruction/Multiple Data)
aplicam instruções distintas a dados distintos. Outra forma de classificar computadores
paralelos se baseia na arquitetura de memória. Nesta classificação, computadores de
memória compartilhada são aqueles em que diversos processadores acessam e utilizam
a mesma área de memória. Dentro desta mesma classificação, ainda podem ser
definidos os sistemas UMA (Uniform Memory Access) onde o tempo de acesso é o
mesmo para todos os processadores que compõem o sistema ou NUMA (Non-Uniform
Memory Access) onde o tempo de acesso dependerá de fatores tal como a interconexão
entre os processadores. Outra classe de computadores paralelos que tem ganhado
destaque nos últimos anos é composta pelos sistemas de memória distribuída onde cada
processador possui sua própria memória, sendo estes sistemas geralmente chamados de
clusters [11]. Sistemas paralelos com milhares de processadores são classificados como
sistema massivamente paralelos, desta forma, ainda existe a designação de
processadores paralelos de larga e pequena escala. Tal designação depende da
característica do processador, assim, máquinas paralelas formadas por diversas
máquinas seriais geralmente são classificadas como máquinas de pequena escala. As
máquinas de processamento paralelo ainda podem ser classificadas como simétricas ou
assimétricas de acordo com a existência de um único tipo de processador ou não [24].
A presente pesquisa foi desenvolvida visando contemplar todas as arquiteturas
paralelas citadas no parágrafo anterior. O objetivo foi obter um simulador de elementos
finitos portável, flexível e eficiente para uma ampla gama de arquiteturas paralelas e
seriais. Embora seja um objetivo difícil de ser plenamente alcançado, neste trabalho,
este problema foi abordado a partir da adoção dos modelos de programação paralela
normalmente utilizados por estes sistemas juntamente com o tratamento de dados de
acordo com a arquitetura em questão.
Dentre os diversos modelos de programação para arquiteturas paralelas, dois se
tornaram dominantes nas últimas décadas, são eles: o paralelismo por troca de
mensagens (message passing) e por multitarefa (multitheading). Estes dois modelos se
diferenciam na forma como os segmentos de execução do programa compartilham e
sincronizam dados entre si [11]. Ambos possuem suas vantagens e desvantagens que
dependerão da arquitetura em questão. Embora o modelo de multitarefa seja
83
particularmente vantajoso para arquiteturas de memória compartilhada, os
computadores paralelos atuais têm sido construídos favorecendo o uso de memória
distribuída devido ao menor custo de construção e manutenção.
Nas próximas seções estes modelos de programação paralela são apresentados, no
contexto de métodos discretos, e baseados em dois padrões amplamente utilizados na
implementação destes modelos, o MPI (Message Passing Interface) e o OpenMP (Open
Multi-Processing).
6.1.1. Troca de mensagens
O modelo de programação paralela por troca de mensagens tem sido o mais utilizado
pelos sistemas de computação paralela de grande porte (o multithread é mais comum
em pequenos sistemas devido ao custo de construção de máquinas paralelas de grande
porte utilizando memória compartilhada). Neste modelo de programação os cálculos são
realizados por vários processos que se comunicam através da chamada de rotinas de
envio e recepção de mensagens. A comunicação entre os processos se realiza de forma
cooperativa onde um dado enviado só será recebido se a chamada da rotina de recepção
correspondente também tiver sido efetuada por outro processo [11]. Este modelo possui
duas grandes virtudes, sendo a primeira a portabilidade. Praticamente qualquer
computador paralelo suporta eficientemente o paralelismo por troca de mensagens,
independentemente do seu hardware. O mesmo não ocorre no modelo multitarefas.
Entretanto, é importante observar que tal portabilidade só pôde ser alcançada após o
desenvolvimento do padrão, Message Passing Interface, ou simplesmente MPI [54] que
sintetizou em uma biblioteca o conjunto de funcionalidades necessárias para
implementação deste modelo de paralelismo. A segunda virtude reside no controle
explícito do armazenamento, localização e utilização dos dados por parte do
programador. Tal virtude algumas vezes é vista como desvantagem por implicar em
aumento na dificuldade de programação.
No contexto de métodos discretos, o paralelismo por troca de mensagens pode ser
implementado a partir da divisão do problema original em diversos problemas menores
que serão tratados, cada um, por um processador ou computador distinto. Esta
84
metodologia, além de intuitiva, é uma forma amplamente empregada em códigos de
elementos finitos. TEZDUYAR et al apresentam em [76] diversas referências que
reportam casos de sucesso na utilização desta metodologia com estruturas de dados
elemento-por-elemento (EBE) e Compressed Sparse Row (CSR). Em [35], [36] JOHAN
et al utilizam técnica similar em um sistema massivamente paralelo. Segundo RIBEIRO
e FERREIRA [63] as estruturas de dados que armazenam somente informações locais,
tal como o esquema elemento-por-elemento e aresta-por-aresta, são particularmente
favoráveis para implementação neste tipo de paralelismo pois não exigem operações de
combinação para montagem da matriz de rigidez, o que geralmente requer a troca de um
grande volume de informações. A Figura 6.1 ilustra como geralmente é feita a divisão
de um problema de elementos finitos em outros de menor porte.
(a) malha original
(b) malha particionada
Figura 6.1 – Paralelismo por troca de mensagens aplicado a malha computacional
Na Figura 6.1 o domínio da malha original (a) é particionado em dois sub-domínios
menores (b) os quais serão tratados de forma quase independente em cada processo.
Nota-se que a única parcela de informação em comum entre estes dois sub-domínios
está representada na Figura 6.1b por linhas tracejadas. A esta parcela de informação
comum costuma-se denominar interface de comunicação e é neste trecho que os
processos devem trocar informações entre si através do envio e recebimento de
mensagens.
Um aspecto importante a ser considerado, e que deriva da forma como a interface
de comunicação se constitui, consiste no volume de comunicação. Quanto menor for o
85
volume de informações a ser enviado/recebido por troca de mensagem, menor será o
tempo consumido pelos processos nesta tarefa e conseqüentemente menos dependentes
eles serão para efetuar seus próprios cálculos. Na prática isso se reflete no desempenho
sendo, portanto, um aspecto muito importante a ser considerado. Outro ponto
fundamental em paralelismo por troca de mensagens diz respeito à quantidade de
cálculos ou esforço computacional a ser efetuado em cada processo. Para sistemas
simétricos, aqueles formados por somente um tipo de processador, é desejável que o
esforço computacional seja igualmente distribuído entre os processos. A distribuição do
esforço computacional costuma-se designar por balanço de carga. A Figura 6.2 ilustra
os conceitos de volume de comunicação e balanço de carga aplicados a uma malha não-
estruturada.
(a) Carga balanceada e grande volume de comunicação
(b) Carga desbalanceada e baixo volume de comunicação
(c) Carga balanceada e baixo volume de comunicação
Figura 6.2 – Conceito de volume de comunicação e balanço de carga aplicados a uma
malha não-estruturada
Na Figura 6.2, associando-se a carga computacional ao número de elementos
triangulares por partição e a quantidade de mensagens trocadas entre os processos
(volume de comunicação) ao número de nós de interface de comunicação pode-se
constatar que:
− Embora a partição representada na Figura 6.2a esteja balanceada, a mesma
apresenta um grande volume de comunicação, ou seja, irá demandar uma grande
quantidade de troca de mensagens entre os dois processos;
86
− A partição representada na Figura 6.2b possui um baixo volume de
comunicação, entretanto, está desbalanceada por ter mais elementos em uma
partição do que na outra. É importante salientar que embora tal característica
seja prejudicial em sistemas simétricos, a mesma pode ser desejável para
sistemas assimétricos, pois neste caso as partições podem ser dimensionadas de
acordo com a capacidade de processamento de cada processador.
− A Figura 6.2c apresenta o caso mais comum e desejável para sistemas simétricos
onde o volume de comunicação é baixo e as partições balanceadas.
Com base no que foi apresentado surge a seguinte questão: “− Como obter um bom
particionamento de dados de acordo com as características de hardware do sistema a
ser utilizado?”
Este assunto foge ao escopo deste trabalho, tendo em vista que diversas pesquisas
anteriores apresentam soluções eficientes para este problema. Para uma revisão
detalhada acerca deste assunto sugere-se o Capítulo 18 da referência [11]. Neste
trabalho, foi utilizada a biblioteca Metis [37] para realizar o particionamento de dados.
Este procedimento ocorre em uma fase de pré-processamento conforme descrito na
seção 6.3
9
. A Metis possui como principais características:
− Produz partições balanceadas ou não de acordo com a escolha do usuário;
− Bastante eficiente. Malhas com milhões de elementos são particionadas em
poucos segundos mesmo em computadores de pequeno porte;
− Desenvolvida na linguagem C ANSI com interface para utilização com outras
linguagens de programação, incluindo Fortran;
− Portátil;
− Fácil utilização;
9
Cabe salientar que para procedimentos que utilizem malhas adaptativas o (re)particionamento de dados
deve ser integrado ao procedimento de solução para o balanceamento de carga ser recuperado sempre que
necessário.
87
− Particiona grafos ou estruturas de dados equivalentes tais como: matrizes e
malhas homogêneas (contendo somente um tipo de elemento) ou não;
− Possui algoritmos de reordenação e redução de largura de banda de matrizes
esparsas.
Na implementação do paralelismo por troca de mensagens no método dos elementos
finitos é comum a utilização de estruturas que facilitem o mapeamento dos sub-
domínios nos processos. No presente trabalho, após o particionamento dos dados da
malha computacional, estes são renumerados individualmente para cada partição. Desta
forma, cada sub-malha é mapeada e tratada pelo mesmo código computacional como se
fossem problemas independentes. As Figuras 6.3 a 6.5 ilustram o processo de
particionamento e uma possibilidade de renumeração. Nestas Figuras a numeração dos
elementos foi substituída por letras para facilitar a compreensão. Os nós coloridos em
azul representam nós internos da partição enquanto os nós vermelhos são nós da
interface de comunicação.
Figura 6.3 – Malha original
88
Figura 6.4 – Malha particionada mantendo a numeração original
Figura 6.5 – Malha particionada e renumerada
A Figura 6.3 apresenta uma malha contendo 16 nós e 20 elementos. Note que após a
criação de duas partições a numeração global torna-se descontínua devido à distribuição
das entidades (Figura 6.4). Para restabelecer a continuidade da numeração deve-se criar
um esquema de numeração que respeite a quantidade local de entidades em cada
partição. A Figura 6.5 apresenta uma possível renumeração onde os nós da interface de
comunicação (nós vermelhos) são inseridos ao final da lista de nós internos (nós azuis).
Note que, na prática, os nós de interface devem ser sempre utilizados em cálculos de
grandezas globais (para a malha original), assim sendo, estes nós, e conseqüentemente
as equações associadas a eles, devem ser de fácil localização. Embora a renumeração
sugerida pela Figura 6.5 apresente um esquema de fácil implementação, o mesmo não
preserva os agrupamentos e ordenações que foram abordados na seção 5.6. Uma
possível solução, e que foi adotada neste trabalho, consiste na criação de índices de
89
localização para os nós de interface de modo que os valores associados a estes nós
sejam facilmente recuperados e montados de uma forma consistente em um buffer de
comunicação.
6.1.2. Multitarefa
Conforme visto na seção anterior, o modelo de programação paralela por troca de
mensagens assume que diferentes componentes de um programa possuem, cada um,
espaços de memória independentes e que dados são explicitamente movidos, copiados e
sincronizados através desses espaços de memória com o uso de operações de envio e
recebimento de mensagens. O modelo de programação por multitarefas, também
conhecido por modelo de memória-compartilhada, é exatamente o oposto no sentido de
que um único programa se subdivide em diversas unidades (linhas de execução) que
concorrem na utilização do mesmo espaço de memória [11]. Conseqüentemente, neste
modelo de programação não existe a necessidade de operações especiais de cópia de
informações, já que diferentes linhas de execução lêem e gravam na memória através do
uso de atribuições de variáveis em comum. A Figura 6.6 ilustra a forma como o
paralelismo por multitarefas é implementado em um programa. Nesta figura, o esforço
computacional está representado a partir da espessura da linha contínua de cor preta.
Nota-se que a linha de execução do programa pode se dividir em diversas regiões
paralelas onde o esforço computacional se distribui igualmente, ou não, pelas múltiplas
linhas de execução (daí o nome multithread).
Como a comunicação nos programas de memória compartilhada é implícita, em
geral, a responsabilidade de atualizar as informações armazenadas na memória é
inteiramente do hardware [11]. Uma conseqüência deste fator se reflete na dificuldade e
nos custos de implementação de máquinas de memória compartilhada, pois, os
processadores realizam cópias da memória convencional em memórias de velocidades
mais altas (cache e registradores). Algumas vezes, essas memórias de alta velocidade
são dispostas em diversos níveis e é de responsabilidade do hardware, juntamente com o
compilador, gerar e manter tais cópias coerentes [11].
90
Figura 6.6 – Modelo de paralelismo por multitarefas ou memória compartilhada
A implementação do paralelismo de memória-compartilhada em computação de
alto desempenho faz uso de um conjunto de API’s
10
(Application Program Interface)
construídas sobre o padrão OpenMP (Open Multi-Processing). O OpenMP consiste de
diversas diretivas de compilação
11
, rotinas e variáveis de ambiente que controlam o
comportamento de um programa em tempo de execução. Pela forma de sua
implementação, a partir da inserção de diretivas de compilação no código-fonte, o
mesmo programa pode ser compilado e executado de forma serial sem maiores
dificuldades.
Este modelo de programação apresenta como principal virtude a facilidade de
implementação, principalmente em programas já existentes, entretanto, é de
responsabilidade do programador tratar as possíveis dependências de memória e
situações de concorrência entre as linhas de execução (race condition). A Figura 6.7
ilustra uma situação bastante comum ao se utilizar o modelo de programação de
memória compartilhada em métodos discretos, tal como o método dos elementos finitos.
Nesta Figura, o laço à esquerda foi paralelizado com o uso da diretiva OpenMP
“C$OMP PARALLEL DO”. Supondo que este laço seja efetuado em duas ou mais
linhas de execução, cada uma delas irá se encarregar de efetuar os cálculos referentes a
um dos cinco elementos ilustrados no lado direito da Figura 6.7. Como a área de
memória representada pelo vetor “x” é compartilhada pelas linhas de execução, logo,
haverá uma situação de concorrência pelo uso e atualização das informações referentes
10
Conjunto de rotinas que são disponibilizadas para uso por um programa.
11
Construção sintática que controla a forma como o código-fonte deve ser tratado pelo compilador.
Tarefa mestre
(serial)
Tarefa mestre
(serial)
Regiões paralelas
91
ao nó três (destacado em vermelho). Este exemplo caracteriza uma dependência de
memória implicitamente imposta pela malha computacional.
Figura 6.7 – Laço com dependência de memória
Uma possível solução para o problema descrito na Figura 6.7 consiste em agrupar os
elementos em blocos disjuntos de forma que dentro de um mesmo bloco não haja
entidades compartilhadas. Esta operação recebe a designação de blocagem ou coloração
de malha, do inglês, mesh coloring e o algoritmo envolvido neste procedimento está
detalhado em SHAKIB et al [68]. A Figura 6.8 mostra um exemplo de laço que pode ser
seguramente efetuado em memória compartilhada após a aplicação da coloração de
malha. Neste exemplo, várias linhas de execução podem atuar simultaneamente dentro
de uma mesma cor, verde por exemplo, sem que haja concorrência pela atualização de
áreas compartilhadas.
Outra virtude do modelo de paralelismo por multitarefas parte da premissa de que
processadores podem se comunicar mais rapidamente através de um canal comum,
viabilizado pela memória, do que em qualquer outro meio de comunicação, como por
exemplo, através de redes conforme normalmente é utilizados em grandes sistemas
distribuidos [54].
C$OMP PARALLEL DO
do i=1,nel
! recupere os nós do elemento
x(no) = x(no) + a
enddo
C$OMP END PARALLEL DO
92
Figura 6.8 – Laço sem dependência de memória
6.1.3. Paralelismo Misto ou Híbrido
Devido à recente proliferação dos sistemas de memória distribuída compartilhada (do
inglês, Distributed Shared-Memory ou simplesmente DSM)
12
tem havido uma constante
busca pelo melhor aproveitamento dos dois mundos. Embora o modelo baseado em
troca de mensagens disponibilize uma maneira eficiente de comunicação entre vários
núcleos de máquinas, este não tira proveito da velocidade do sistema de memória
compartilhada, quando este é mais eficiente do que a comunicação convencional através
de redes. Atualmente existem iniciativas de software que facilitam a utilização e
conversão de programas já existentes, paralelos ou não, para tirarem proveito de
plataformas DSM (ver, por exemplo, [53]), entretanto, cabe observar que softwares que
foram projetados para usar ambos os modos de programação paralela já são,
intrinsecamente, preparados para uso em sistemas DSM.
Em métodos discretos, o paralelismo híbrido pode ser implementado tratando-se os
dados para os dois modelos de paralelismo conforme citado nas seções anteriores.
Conseqüentemente, os dados devem ser particionados e cada partição colorida em
blocos de entidades disjuntas conforme ilustrado na Figura 6.9.
12
A interligação de máquinas multi-processadas através de uma rede seria um exemplo de memória
distribuída compartilhada
ielm = 0
do icor = 1, ncores
nvec = ielblk(icor)
C$OMP PARALLEL DO
C$OMP& FIRSTPRIVATE(NVEC)
do i = ielm+1, ielm+nvec
! recupere os nós do elemento
x(no) = x(no) + a
enddo
C$OMP END PARALLEL DO
ielm = ielm+nvec
enddo
93
(a) Troca de mensagens
(b) Multitarefa
(c) Híbrido
Figura 6.9 – Tratamento de dados para paralelismo.
Na Figura 6.9a a malha original foi particionada para ser utilizada em 4 processos
paralelizados por troca de mensagens (MPI) enquanto que na Figura 6.9b os elementos
da malha foram agrupados em blocos de elementos disjuntos afim de remover as
dependências de memória que impedem o paralelismo em sistemas de memória
compartilhada. Finalmente, na Figura 6.9c os dados são preparados para utilização em
plataformas híbridas que combinem os dois modelos de paralelismo. A Figura 6.10 a
seguir ilustra como o produto matriz-vetor (elemento-por-elemento ou aresta-por-aresta)
pode ser efetuado com base nos dois modelos de paralelismo citados anteriormente em
uma plataforma DSM.
94
Figura 6.10 – Produto matriz-vetor híbrido
6.2. Vetorização e pipelining
A vetorização é um artifício utilizado em computadores para converter instruções, que
antes eram efetuadas termo a termo, para que possam ser aplicadas simultaneamente a
um agrupamento de dados dispostos na forma de um vetor (série de valores adjacentes
na memória). Considere, por exemplo, a operação de adição de dois vetores “a” e “b”
conforme ilustrado na Figura 6.11. Nesta Figura, o resultado da operação de adição será
armazenado em um terceiro vetor “c”. A linha tracejada representa o dado sobre o qual
a instrução está sendo aplicada. Se para operar cada termo dos vetores “a” e “b” um
processador escalar consome 1 ciclo de processamento, serão necessários 3 ciclos para
completar a operação para todos os termos do vetor (Figura 6.11a). Em um processador
vetorial a mesma operação poderia ser efetuada simultaneamente para vários termos dos
vetores a cada ciclo de processamento. Supondo-se um processador vetorial que consiga
aplicar a instrução de soma à 3 posições dos vetores “a” e “b” simultaneamente, seria
necessário somente 1 ciclo para realizar a mesma operação que em um processador
escalar consome-se três.
...Multitarefa...
iside = 0
DO iblk = 1, nedblk
nvec = ia_edblk(iblk)
!dir$ ivdep
!$OMP PARALLEL DO
DO ka = iside+1, iside+nvec, 1
...calculos do MATVEC...
ENDDO
!$OMP END PARALLEL DO
ENDDO
...Troca de mensagens...
#ifdef MPICODE
call MPI_AllReduce
#endif
95
a(1)a(2)a(3)
+
b(1)b(2)b(3)
=
c(1)c(2)c(3)
a(1)a(2)a(3)
+
b(1)b(2)b(3)
=
c(1)c(2)c(3)
a(1)a(2)a(3)
+
b(1)b(2)b(3)
=
c(1)c(2)c(3)
ciclo 1ciclo 2ciclo 3
(a)
a(1)a(2)a(3)
+ + +
b(1)b(2)b(3)
= = =
c(1)c(2)c(3)
ciclo 1
(b)
Figura 6.11 – Operação de adição de dois vetores (a) escalar (b) vetorial
O pipelining é outro artifício utilizado nas arquiteturas de computadores para aumentar
o numero de instruções executadas por tempo (métrica que em inglês é conhecida por
throughput [74]). No pipelining as operações são entrelaçadas de forma a aproveitar os
ciclos do processador mais eficientemente. Para ilustrar como o pipelining funciona,
considere que para realizar uma dada instrução sejam necessários os seguintes estágios:
− F: Fetch; recuperação de um dado ou instrução;
− D: Decode: decodificação da instrução;
− E: Execute: execução da instrução;
− WB: Write back: escrita do resultado da instrução na memória
Sem utilizar pipeline a execução de 4 instruções pelo processador seria similar ao
ilustrado na Figura 6.12, onde para ser iniciada uma nova instrução deve-se aguardar o
final de outra.
96
F
D
E
WB
F
D
E
WB
F
D
E
WB
F
D
E
WB
TEMPO
Figura 6.12 – Instruções sendo efetuadas sem o uso de pipeline.
Com o uso do pipeline as 4 instruções podem ser executadas de forma entrelaçada, ou
seja, uma instrução pode ser iniciada sem que seja necessário aguardar o término da
última. A Figura 6.13 ilustra este mecanismo.
F
D
E
WB
F D E WB
F
D
E
WB
F
D
E
WB
INSTRUÇÃO
TEMPO
Figura 6.13 – Instruções sendo efetuadas com o uso de pipeline.
Embora a vetorização e o pipelining sejam artifícios que aumentam a eficiência no
uso dos processadores, as duas técnicas sofrem do mesmo problema de dependência de
memória exemplificado para o paralelismo de memória compartilhada (seção 6.1.2).
Note, por exemplo, que se na Figura 6.11 o resultado da soma dos vetores “a” e “b”
fosse armazenado em um desses dois vetores, o cálculo incorreria em erro caso a
operação fosse efetuada da forma vetorial. O mesmo se aplica ao caso do pipeline, onde
o resultado de uma instrução é utilizado em instruções que já tenham sido iniciadas e
que ainda não foram completadas. Da mesma forma como foi descrito para o problema
de paralelismo de memória compartilhada, a dependência de memória para utilização
dessas duas técnicas pode ser removida com a aplicação da técnica de coloração de
malha [68].
97
6.3. Pré-processamento
Conforme normalmente ocorre em softwares de simulação, neste trabalho foi
desenvolvida uma etapa de pré-processamento onde são realizadas as operações de
preparação dos dados para o procedimento de solução. Tais operações são apresentadas
no Algoritmo 6.1.
Algoritmo 6.1 – Pré-processamento.
1. Construção do modelo computacional representando a geometria e o domínio
do problema;
2. Criação da malha computacional de tetraedros lineares;
3. Definição e aplicação de condições de contorno;
4. Particionamento do domínio (caso seja requerida uma simulação em paralelo
utilizando o MPI);
5. Para cada domínio (ou partição paralela) fazer:
5.1 Extrair incidência de arestas;
5.2 Extrair incidência dos elementos pelas arestas;
5.3 Agrupar e aplicar reordenações requeridas pelo usuário ou
automaticamente definidas pela EdgePack;
5.4 Gerar arquivos para o procedimento de solução.
OBSERVAÇÕES:
− Paras os passos 1 a 3 são utilizados softwares de terceiros. Em especial, para
facilitar a definição de novos problemas foi desenvolvido um procedimento
genérico de criação de modelos utilizando o software comercial Ansys [1]. Este
procedimento consiste na utilização deste software, na realização criação de
modelos, geração de malhas e aplicação de condições de contorno e iniciais, de
forma tão natural quanto se estivesse criando um problema para o próprio
Ansys. Depois de concluída a etapa 3, os dados são exportados através do uso de
uma macro especificamente desenvolvida para este fim (macro
Ansys2EdgeCFD). Os arquivos gerados por esta macro são, então, enviados ao
98
software de pré-processamento, também parte integrante do sistema
desenvolvido neste trabalho, o qual realiza as etapas subseqüentes (etapas 4 e 5);
− No passo 4, a criação das partições paralelas é opcional e requerida somente para
as simulações utilizando paralelismo por troca de mensagens. O número de
partições a ser gerado é dado de entrada para o usuário e o particionamento é
gerado com o uso da biblioteca Metis [37] conforme mencionado na seção 6.1.1;
− Os passos 5.1 a 5.3 são realizados pelo software EdgePack. Uma biblioteca
desenvolvida por MARTINS et al [50] para a construção e tratamento de
estruturas de dados por arestas.
6.4. Pós-processamento
A interpretação dos dados gerados durante ou após uma simulação é um aspecto
fundamental da dinâmica de fluidos computacional [27]. Tal aspecto é grandemente
auxiliado com o uso de técnicas de visualização. O conceito de visualização refere-se à
transformação de dados e/ou informações em figuras e gráficos, ou seja, a visualização
explora o principal sentido do ser humano, a visão, e, conseqüentemente, expande a sua
capacidade de raciocinar sobre o que está sendo observado [66]. Com o aumento na
capacidade de processamento dos computadores, o volume de dados gerados por uma
simulação do escoamento de fluidos tridimensional e transiente pode facilmente
ultrapassar a capacidade de memória de uma única máquina de visualização. Nestes
casos, o uso de técnicas de pós-processamento paralelo se torna imperativo [66]. Além
disso, os recentes avanços nos algoritmos de visualização permitem que o mesmo dado
seja explorado de diversas formas diferentes facilitando a sua compreensão.
Outro ponto que merece destaque se refere à forma como as informações são
processadas. Freqüentemente o computador que realiza os cálculos não é o mesmo
utilizado na geração das imagens para visualização. Geralmente, estes computadores
estão montados em grandes centros de computação e o acesso é feito através de redes
enquanto as máquinas destinadas ao pós-processamento são normalmente em menor
número e dotadas de placas gráficas com grande poder de processamento. Esta
99
separação geográfica acaba trazendo à tona o problema do tráfego de grandes volumes
de informação através de redes com capacidades de transmissão limitadas.
No presente trabalho não houve preocupação com o estudo e desenvolvimento de
novas ferramentas de visualização, tendo em vista que existem bibliotecas e softwares
gratuitos e código-fonte aberto que estão alinhados na solução dos problemas citados
nos parágrafos anteriores. Dentre estas bibliotecas a mais proeminente e difundida
atualmente é o Visualization Toolkit ou simplesmente VTK [66]. Formalmente, o VTK
é um sistema orientado a objetos de código-fonte aberto e destinado à computação
gráfica, visualização e processamento de imagens [81]. Na prática o VTK constitui-se
de um amplo acervo contendo centenas de classes, desenvolvidas na linguagem C
++
, que
auxiliam na construção de aplicações gráficas e de processamento de imagens. Dentre
as diversas aplicações gráficas baseadas no VTK os pós-processadores ParaView
(http://www.paraview.org) e VisIt (http://www.llnl.gov/visit/) são exemplos de
softwares gratuitos com grande aplicabilidade na visualização de simulações de grande
porte, foco deste trabalho. Dentre os recursos de particular interesse, disponibilizados
por estes softwares, destacam-se os seguintes:
− Pós-processamento em paralelo: Através deste recurso, cada processador ou
máquina paralela trata da sua porção de informações e constrói uma parcela da
imagem final o que permite a visualização de simulações de grande porte;
− Pós-processamento remoto: Estes softwares podem ser utilizados para pós-
processamento remoto onde o sistema que realiza a simulação está
geograficamente distante do sistema de visualização. Nestes casos, um recurso
que constrói a imagem, sem a necessidade do tráfego de grandes volumes de
dados, pode ser usado para acelerar o processo de visualização. Este mesmo
recurso, chamado em inglês de offscreen rendering, permite que o processo de
geração de imagens seja efetuado por software ao invés de placas gráficas. Isto
permite que mesmo máquinas com baixo poder de processamento gráfico
possam ser utilizadas no pós-processamento de grandes volumes de
informações;
100
− Portabilidade: Ambos os softwares são distribuídos na forma de binários ou
código-fonte e são testados para compilação nos sistemas operacionais mais
difundidos em processamento de alto desempenho;
− Segurança: O ParaView utiliza um sistema reverso de comunicação entre as
máquinas servidoras de dados/imagens e as máquinas clientes. Desta forma, a
conexão pode ultrapassar, de dentro para fora, os sistemas de proteção (firewall)
normalmente utilizados em centros de pesquisa que mantém grandes máquinas
paralelas. O VisIt utiliza conexão encriptada entre máquinas cliente e servidoras
com o uso do protocolo SSH (Secure Shell);
− Extensibilidade: Ambos os softwares são de código-fonte aberto e gratuitos e
possuem grandes comunidades de desenvolvedores e colaboradores. Isto permite
que futuras modificações e/ou adaptações do software original sejam realizadas
com pouco investimento de homens/hora. Não obstante, por serem softwares
baseados no VTK, eles carregam todos os recursos de visualização
disponibilizados pela biblioteca, ou seja, mesmo que o software não possua
algum recurso de imediato, o mesmo pode ser adicionado caso já exista dentro
do VTK ou desenvolvido caso se trate de algo de interesse comum;
− Estado-da-arte em visualização: O VTK é dotado de diversos algoritmos que
representam o estado-da-arte em visualização, desta forma, algoritmos como
renderização volumétrica, estereoscopia, cavernas, renderização em escala
ampliada (paredes de visualização), dentre outros, são intrinsecamente
contemplados nestes softwares.
Dentre os diversos formatos de arquivos suportados pelo VTK, e conseqüentemente
pelos softwares ParaView e VisIt, optou-se pela utilização do formato Ensight Gold por
possuir as seguintes características:
− Simplicidade de implementação;
− Flexibilidade. O formato Ensight Gold é implicitamente suportado por qualquer
software baseado no VTK, além do próprio software Ensight, expandindo as
possibilidades de visualização;
101
− Suporte a arquivos paralelos através do formato Server-of-Server (SOS);
− Formato de arquivo incremental, ou seja, os resultados para cada passo de tempo
são gravados em arquivos separados, o que permite que resultados parciais
sejam visualizados antes do término da simulação;
− O formato Ensight Gold permite que sejam escritos arquivos em formato ASCII,
favoráveis a depuração de problemas pequenos, ou formato binário, que são
apropriados para armazenar grandes volumes de dados.
A Figura 6.14 mostra um caso típico de visualização remota em paralelo onde uma
simulação é executada por 2 processadores de uma máquina paralela remota e os
resultados são visualizados em uma máquina de pequeno porte através da internet.
Neste caso, o trabalho de construção da imagem (renderização) pode ser, ou não,
realizado por processadores de cálculo da máquina paralela (offscreen rendering) e o
resultado final é enviado para a máquina cliente através da internet. Note que este
procedimento evita que grandes volumes de informação trafeguem por redes de
capacidade de comunicação limitada. Outras configurações de visualização, diferentes
da mostrada na Figura 6.14, podem ser montadas de acordo com os recursos disponíveis
ou a necessidade. Note que, embora o domínio seja particionado a imagem final é de um
modelo único. Isto ocorre, pois o software de visualização é provido de mecanismos que
unificam e costuram o modelo paralelo em uma única representação. Para uma revisão
de métodos avançados de visualização as referências [29], [81] são sugeridas.
102
Figura 6.14 – Visualização remota em paralelo
103
Capítulo 7
Conclusões
Diversas técnicas foram reunidas tendo como base a formulação estabilizada
SUPG/PSPG de elementos finitos para a solução de problemas de grande porte
envolvendo escoamentos incompressíveis Newtonianos e não-Newtonianos. Em
consonância com a simulação do escoamento, a solução do problema de transporte
advectivo-difusivo foi desenvolvida, permitindo que uma ampla gama de fenômenos
fosse estudada e abrindo-se fronteira para o tratamento de novos problemas de interesse.
Tendo como foco contínuo a solução de problemas de grande porte, especial
atenção foi devotada à eficiência computacional dos métodos e algoritmos utilizados. A
estrutura de dados por arestas, associada a técnicas de reordenação, foi aplicada para
aumentar a eficiência no processo de solução, resultando no uso eficiente da memória
cache (de acordo com as características do processador), menor consumo de memória
(quando comparado a técnicas tradicionais tais como EBE e CSR) e redução da
quantidade de endereçamentos indiretos. Associado à estrutura de dados por arestas, foi
utilizado o método de Newton inexato no procedimento de solução não-linear o qual
permite que os valores de tolerância linear para o método GMRES sejam adaptados de
acordo com a proximidade da solução desejada assim como pelo histórico dos resíduos.
Para a integração no tempo o algoritmo preditor-multicorretor foi utilizado juntamente
com o algoritmo de controle de adaptação de passo de tempo através de controlador do
tipo PID, o qual seleciona o passo de tempo de acordo com a tolerância ao erro ajustada
para o controlador.
104
Especificamente para o problema de transporte, o algoritmo de desativação
dinâmica foi desenvolvido e aplicado para concentrar o esforço computacional somente
nas regiões onde a solução sofre alterações relevantes. Embora esta técnica tenha sido
desenvolvida para problemas de superfície-livre a mesma pode ser estendida para
problemas de transporte de massa e/ou energia ou até mesmo para reduzir o esforço
computacional na solução das equações de Navier-Stokes. Ainda para o problema de
transporte, foi desenvolvido um algoritmo de correção de massa para problemas de
superfície-livre que praticamente anula as perdas sofridas ao se transportar a função de
marcação.
Com base no acoplamento entre as soluções do escoamento e transporte foram
tratados problemas de superfície-livre e escoamentos gravitacionais dentro da validade
do modelo de Boussinesq. Para problemas de superfície-livre, os métodos de volume-
de-fluido (VOF) e o método de curvas de nível (LS) foram adotados para o
acompanhamento Euleriano da evolução da interface. A principal distinção entre estes
métodos reside na função de marcação utilizada para determinar a localização da
interface. Para o método VOF, a associação da formulação SUPG com o termo de
captura de descontinuidade
YZ
e o algoritmo de correção de massa permitiu que
problemas complexos de superfície-livre, envolvendo o impacto e fragmentação da
interface, fossem simulados com sucesso. Um dos principais desafios na utilização do
método de curvas de nível para simulação de problemas de superfície-livre consiste na
criação e manutenção de uma função de marcação do tipo função distância em malhas
não estruturadas e geometrias complexas. Neste sentido, foi desenvolvido um algoritmo
especificamente destinado a este fim. Este algoritmo, além de fácil implementação é
eficiente e produz resultados precisos conforme apresentados na publicação do
Apêndice C.
Todos os algoritmos e métodos desenvolvidos para o processo de solução são
paralelos e preparados para tirar proveito dos principais artifícios empregados em
processamento de alto desempenho, dos quais citam-se: paralelismo de memória
distribuída, paralelismo de memória compartilhada, vetorização e pipelining. Estes
artifícios podem ser utilizados de forma isolada ou combinada de acordo com a
arquitetura disponível. Tal flexibilidade é alcançada a partir do tratamento adequado dos
dados que vão desde o particionamento balanceado da malha, provido pela biblioteca
105
Metis [37], remoção da dependência de dados, realizada pelo algoritmo de coloração de
malha [68], até a ordenação dos dados tirando-se proveito das propriedades do
hardware.
De forma a viabilizar a solução de problemas de grande porte, métodos e softwares
de visualização apropriados foram explorados permitindo que simulações paralelas
fossem acompanhadas mesmo à distância através da internet.
Este trabalho resultou no desenvolvimento do software registrado sobre o nome de
EdgeCFD, software este desenvolvido em linguagem de programação Fortran90 e
formado pelas bibliotecas Metis [37] (partitionamento de dados), EdgePack [50]
(extração, ordenação e agrupamento de dados) e dos programas EdgeCFDPre
(preparação de arquivos para o solver a partir da utilização das bibliotecas Metis e
EdgePack) e EdgeCFDSolver (procedimento de solução onde a maior parte desta
pesquisa se concentrou).
106
Trabalhos futuros
Durante a realização desta pesquisa surgiram as seguintes sugestões como extensão para
trabalhos futuros:
− Modelagem de turbulência livre de parâmetros e modelos adicionais a partir do
método variacional multiscala segundo as referências [2], [25], [31] e [65];
− Avaliação de outros métodos para o transporte de choques e descontinuidades
que preservem as características da função de marcação empregada em
problemas de superfície-livre. Nesta avaliação sugerem-se os métodos do
Transporte de Fluxo Corrigido, do inglês, Flux Corrected Transport ou
simplesmente FCT [42], [45], [88] e o método de Interpolação-Cúbica com
Coordenadas de Volume/Área, do inglês, Cubic-Interpolated with Volume/Area
ou simplesmente CIVA [75];
− Desenvolvimento de metodologia para o tratamento de movimento de malha
pelo método Lagrangeano e Euleriano Arbitrário (ALE) [51];
− Estudo e implementação de métodos de refinamento adaptativo de malha;
− Estudo e desenvolvimento de métodos de contorno imersos no qual a geometria
do problema é implicitamente inserida em um domínio de geometria simples
(cubo ou paralelepípedo, por exemplo) [43];
− Aprimorar a implementação paralela existente de forma a contemplar
arquiteturas massivamente paralelas;
− Desenvolvimento de pré-condicionadores paralelos utilizando estrutura de dados
por aresta;
− Desenvolvimento de interface gráfica, manuais de usuário e cursos de
treinamento.
107
Referências
[1] ANSYS INC., Ansys Theory Reference, cp. 7 – Fluid Flow, URL:
http://www.ansys.com.
[2] AKKERMAN, I., BAZILEVS, Y., CALO, V. M., COTTRELL, J. A., HUGHES,
T. J. R., HULSHOFF, S., “The Role of Continuity in Residual-Based
Variational Multiscale Modeling of Turbulence”, Computational Mechanics,
2007, DOI 10.1007/s00466-007-0193-7.
[3] BAZILEVS, Y., CALO, V. M., TEZDUYAR, T. E. e HUGHES, T. J. R.,
YZ
“Discontinuity Capturing for Advection-Dominated Processes with Application
to Arterial Drug Delivery”, International Journal for Numerical Methods in
Fluids, 54:593-608, 2007.
[4] BEHR, M. A., FRANCA, L. P., e TEZDUYAR, T. E., "Stabilized Finite Element
Methods for the Velocity-Pressure-Stress Formulation of Incompressible
Flows", Computer Methods in Applied Mechanics and Engineering, 104:31-48,
1993.
[5] BIRD, R. B., STEWART, W. E. AND LIGHTFOOT, E. N., Transport
Phenomena, John Wiley & Sons, 1960.
[6] BROOKS, A. N. and HUGHES, T. J. R., “Streamline Upwind/Petrov-Galerkin
Formulations for Convection Dominated Flows with Particular Emphasis on
the Incompressible Navier-Stokes Equation”, Comput. Methods Appl. Mech.
Engrg. 32, pp. 199-259, 1982.
[7] CATABRIGA, L., Soluções Implícitas das Equações de Euler Empregando
Estrutura de Dados por Arestas, Tese de Doutorado - COPPE/UFRJ, Maio
2000.
108
[8] COUTINHO, A. L. G. A., MARTINS, M. A., SYDENSTRICKER, R. M., ELIAS,
R. N., “Performance Comparison of Data-Reordering for Sparse Matrix-Vector
Multiplication in Edge-Based Unstructured Grid Computations”, Int. J. for
Num. Meth. Engrg., 66:431-460, 2006.
[9] de SAMPAIO, P. A. B., “A Petrov-Galerkin Formulation for the Incompressible
Navier-Stokes Equations using Equal Order Interpolation for Velocity and
Pressure”, Int. J. for Num. Meth. Engrg., 31:1135-1149, 1991.
[10] DONEA, J. e HUERTA, A., Finite Element Methods for Flow Problems, John
Wiley & Sons, ISBN 0-471-49666-9.
[11] DONGARRA, J., FOSTER, I., FOX, G., GROPP, W., KENNEDY, K.,
TORCZON, L. e WHITE, A., Sourcebook of Parallel Computing, Morgan
Kaufmann Publishers, 2003, São Francisco, EUA, ISBN: 1-55860-871-0
[12] EISENSTAT, S. C. AND WALKER, H. F., “Choosing the Forcing Terms in
Inexact Newton Method”, SIAM. J. Sci Comput. 17(1):16–32, 1996.
[13] EISENSTAT, S. C. AND WALKER, H. F., “Globally Convergent Inexact New-
ton Methods”, SIAM J. Optim., 4: 393-422, 1994.
[14] ELIAS R. N. Métodos Tipo-Newton Inexatos para a Solução de Problemas não-
Lineares Resultantes da Formulação SUPG/PSPG das Equações de Navier-
Stokes Incompressíveis em Regime Permanente, Dissertação de MSc,
Departamento de Engenharia Civil, COPPE/UFRJ, 2003.
[15] ELIAS, R. N., COUTINHO, A. L. G. A. e MARTINS, M. A. D., “Inexact Newton-
Type Methods for the Solution of Steady Incompressible Viscoplastic Flows
with the SUPG/PSPG Finite Element Formulation”, Comput. Methods Appl.
Mech. Engrg., 195:3145-3167, 2006.
[16] ELIAS, R. N., COUTINHO, A. L. G. A. e MARTINS, M. A. D., “Parallel Edge-
Based Inexact-Newton Solution of Steady Incompressible 3D Navier-Stokes
Equations”, Euro-par, 2005, Lisboa, Portugal, publicado no Lecture Notes in
Computer Science 3648.
109
[17] FELIPPA, C. A., “A Historical Outline of Matrix Structural Analysis: A Play in
Three Acts”, Comput. Struct, 79(14):673-705, 2001.
[18] FERZIGER, J. H. e PERIC, M., Computational Methods for Fluid Dynamics,
Springer 3
rd.
edition, 2001.
[19] FLUENT INC., Theory Manual FIDAP 8, December 1998.
[20] FOX, R. W. AND McDONALD, A. T., Introdução à Mecânica dos Fluidos, LTC,
5a. ed., 1998.
[21] FRANCA, L. P. e FREY, S. F., Stabilized Finite Element Methods: II. The
Incompressible Navier-Stokes Equations, Comput. Methods Appl. Mech.
Engrg, 99:209-233, 1992.
[22] GALEÃO, A. C., DO CARMO, E. G. D. , “A Consistent Approximate Upwind
Petrov-Galerkin Method for Convection-Dominated Problems”. Comput.
Methods Appl. Mech. Engrg, 68(1):83-95, 1988.
[23] GARTLING, D.K., "Finite Element Methods for Non-Newtonian Flows", report
SAND92-0886, CFD Dept., Sandia National Laboratories, Albuquerque, NM
1992.
[24] GRAMA, A., GUPTA, A., KARYPIS, G., and KUMAR, V., Introduction to
Parallel Computing, 2003 ISBN 0-201-64865-2.
[25] GRAVEMEIER, V., LENZ, S., WALL, W. A., “Towards A Taxonomy For
Multiscale Methods in Computational Mechanics: Building Blocks of Existing
Methods”, Comput Mech, 2007, DOI 10.1007/s00466-007-0185-7
[26] GRESHO, P. M. e SANI, R. L., Incompressible Flow and the Finite Element
Method: Advection-Diffusion and Isothermal Laminar Flow, John Wiley &
Sons Ltd, 1999, ISBN 0-471-96789-0.
[27] GRIEBEL, M., DORNSEIFER, T. AND NEUNHOEFFER T., Numerical
Simulation in Fluid Dynamics – A Practical Introduction, SIAM, Philadelphia,
1998.
110
[28] GUSTAFSSON, K., “Control Theoretic Techniques for Stepsize Selection in
Implicit Runge-Kutta Methods”, ACM Transactions on Mathematical
Software, 20:496-517, 1994.
[29] HANSEN, C. e JOHNSON, C. R., The Visualization Handbook, Elsevier, 2004.
ISBN: 978-0-12-387582-2
[30] HIRT, C. W. e NICHOLS, B. D., “Volume-of-Fluid (VOF) Method for the
Dynamics of Free Boundaries”, Journal of Computational Physics, 39:201-
225, 1981.
[31] HUGHES, T. J. R., “Multiscale phenomena: Green’s functions, the Dirichlet-to-
Neumann formulation, subgrid scale models, bubbles and the origins of
stabilized methods”. Comput Methods Appl Mech Eng, 127:387–401, 1995;
[32] HUGHES, T. J. R., The Finite Element Method - Linear Static and Dynamic Finite
Element Analysis, Prentice-Hall, Inc., 1987.
[33] HUGHES, T. J. R., FRANCA, L. P. e BALESTRA, M., “A New Finite-Element
Formulation for Computational Fluid Dynamics: V. Circumventing the
Babuska-Brezzi Condition: A Stable Petrov-Galerkin Formulation of the
Stokes Problem Accomodating Equal-Order Interpolations”, Comput. Methods
Appl. Mech. Engrg, 59:85-99, 1986.
[34] INCROPERA, F., P.; WITT, D. P. de, Fundamentals of Heat and Mass Transfer,
3rd Ed., John Wiley & Sons, 1990. ISBN 0-471-51729-1
[35] JOHAN, Z., HUGHES, T. J. R., MATHUR, K., K. e JOHNSSON, S. L., “A Data
Parallel Finite Element Method for Computational Fluid Dynamics on the
Connection Machine System”, Comput. Meth. Appl. Mech. Engrg, 99:113-134,
1992.
[36] JOHAN, Z., MATHUR, K. K., JOHNSSON, S. L. e HUGHES, T. J. R, “An
Efficient Communications Strategy for Finite Element Methods on the
Connection Machine CM-5 System”, Comput. Meth. Appl. Mech. Engrg,
(113)363-387, 1994.
111
[37] KARYPIS, G., KUMAR, V., Metis 4.0: Unstructured Graph Partitioning and
Sparse Matrix Ordering System. Technical report. Department of Computer
Science, University of Minnesota, Mineapolis, EUA, 1998.
http://glaros.dtc.umn.edu/gkhome/metis/metis/overview.
[38] KELLEY, C. T., Iterative Methods for Linear and Nonlinear Equations, Frontiers
in applied mathematics, SIAM, Philadelphia, 1995.
[39] KLEEFSMAN, K. M. T., Water Impact Loading on Offshore Structure: A
Numerical Study, PhD thesis, University of Gröningen, Gröningen, 2005.
[40] KOLMOGOROV, A., “The Equations of Turbulent Motion in an Incompressible
Fluid”, Izv. Sci. USSR Phys., 6:56-58, 1942.
[41] KRAFT, R. A., COUTINHO, A. L. G. A. e de SAMPAIO, P. A. B., “Edge-Based
Data Structures for the Incompressible Navier-Stokes Equations with Heat
Transfer”, Int. J. Num. Meth. Fluids, 53:1473-1494, 2007.
[42] KUZMIN, D., MÖLLER, M. e TUREK, S., “Multidimensional FEM-FCT
Schemes for Arbitrary Time-Stepping”, Int. J. Numer. Meth. Fluids, 42:265–
295, 2003.
[43] LOHNER, R., BAUM, J. D., MESTREAU, E., SHAROV, D., CHARMAN, C. e
PELESSONE, D., “Adaptive Embedded Unstructures Grid Methods”, Int. J.
Numer. Meth. Engrg, 60:641-660, 2004.
[44] LOHNER, R., CAMELLI F. “Dynamic Deactivation for Advection-Dominated
Contaminant Transport”, Comm. Num. Meth. Engrg, 20:639-646, 2004.
[45] LOHNER, R., MORGAN, K, PERAIRE, J. e VAHDATI, M., “Finite Element
Flux-Corrected Transport (FEM-FCT) for the Euler and Navier-Stokes
Equations”, Int. J. Num. Meth. Fluids, 7:1093-1109, 1987.
[46] LOHNER, R, YANG, C. e OÑATE, E. “On the Simulation of Flows with Violent
Free-Surface Motion”, Comput. Methods Appl. Mech. Engrg, 195:5597-5620,
2006.
112
[47] LOHNER, R., Applied CFD Techniques: An Introduction Based on Finite Element
Methods, John Wiley & Sons, LTD., 2001, ISBN: 0-471-49843-2.
[48] MARTÍNEZ, A. E. T. e KENNETH, J., “On the Interaction Between Dynamic
Model Dissipation and Numerical Dissipation due to Streamline
Upwind/Petrov-Galerkin Stabilization”, Comput. Methods Appl. Mech. Engrg,
194:1225-1248, 2005.
[49] MARTINS, M. A. D., Estrutura de Dados por Arestas para a Solução de Pro-
blemas de Plasticidade Computacional, Tese de Doutorado, Programa de
Engenharia Civil, COPPE/UFRJ, Rio de Janeiro, 2001.
[50] MARTINS, M., ELIAS R., COUTINHO, A., “EdgePack: A Parallel Vertex and
Node Reordering Package for Optimizing Edge-Based Computations in
Unstructured Grids”, 7th. High Performance Computing Science –
VECPAR’06, LNCS 4395, 2006, Rio de Janeiro, Brazil.
[51] MASUD, A, HUGHES, T.J.R., “A Space-Time Galerkin/Least-Squares Finite
Element Formulation of The Navier Stokes Equations for Moving Domain
Problems”, Comput. Methods Appl. Mech. Engrg., 146:91-126, 1997
[52] MEISTER, A. AND VÖMEL, C., “Efficient Preconditioning of Linear Systems
Arising from the Discretization of Hyperbolic Conservation Laws”, Com-
putacional Fluid Dynamics and Data Analysis, 10 – July 1999
[53] MENDES, R., WHATELY, L., BENTES, C. e AMORIM, C. L., “Implementation
and Preliminary Evaluation of a New Software DSM System for Clusters”, in
Proceedings of the VI Workshop on High Performance Computer Systems, Rio
de Janeiro, 2005.
[54] Message Passing Interface Forum. MPI: A message-passing interface standard,
International Journal of Supercomputer Applications and High Performance
Computing, 8(3/4):165-414, 1994. (Special issue on MPI, also available at
ftp://www.netlib.org/mpi/mpi-report.ps)
113
[55] NAGRATH, S., Adaptive Stabilized Finite Element Analysis of Multi-Phase Flows
Using Level Set Approach, PhD thesis, Rensselaer Polytechnic Institute, New
York, USA, 2004.
[56] O'CONNOR, J. J.; EDMUND F. R. “Claude-Louis Navier”. MacTutor History of
Mathematics archive, URL: http://www-history.mcs.st-
andrews.ac.uk/Biographies/Navier.html
[57] O'CONNOR, J. J; EDMUND F. R. “George Gabriel Stokes”. MacTutor History of
Mathematics archive, URL: http://www-history.mcs.st-
andrews.ac.uk/Biographies/Stokes.html
[58] OSHER, S, FEDKIW, R, Level Set Methods and Dynamic Implicit Surfaces,
Applied mathematical sciences, 153, Springer Verlag, 1st edition, 2002.
[59] OWENS, R. G. e PHILLIPS, T. N., Computational Rheology, Imperial College
Press, 2002, ISBN: 1-86094-186-9.
[60] PAPADRAKAKIS, M., Solving Large-Scale Problems in Mechanics: The De-
velopment and Application of Computational Solution Procedures, John Wiley
& Sons, England, 1993.
[61] POLYANIN, A. D., KUTEPOV, A. M., VYAZMIN, A. V., and KAZENIN, D. A.,
Hydrodynamics, Mass and Heat Transfer in Chemical Engineering, Taylor &
Francis, London, 2002. ISBN 0-415-27237-8
[62] RIBEIRO, F. L. B. e COUTINHO, A. L. G. A., “Comparison between Element,
Edge and Compressed Storage Schemes for Iterative Solutions in Finite
Element Analysis”, Int. J. Numer. Methods in Engrg. 36(4):659-685, 2005.
[63] RIBEIRO, F. L. B e FERREIRA, I. A., “Parallel Implementation of the Finite
Element Method Using Compressed Data Structures”, Comput. Mech., 41:31-
48, 2007.
[64] SAAD, Y., Iterative Methods for Sparse Linear Systems, Boston, PWS Publish-ing
Company, 1996
114
[65] SAGAUT, P., DECK, S. e TERRACOL, M., Multiscale and Multiresolution
Approaches in Turbulence, Imperial College Press, 2006, ISBN 1-86094-650-
X
[66] SCHROEDER, W., MARTIN, K. e LORENSEN, B., The Visualization Toolkit: An
Object-Oriented Approach to 3D Graphics, 3
rd
edition, Kitware, Inc.
publishers, ISBN 1-930934-12-2, 2004.
[67] SETHIAN J. A., Level Set Methods and Fast Marching Methods – Evolving
Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision,
and Materials Science, Cambridge Monographs on Applied and Computational
Mathematics, 2nd edition, 1999.
[68] SHAKIB, F., HUGHES, T. J. R. AND JOHAN, Z., “A Multi-element Group
Preconditioned GMRES Algorithm for Nonsymmetric Systems Arising in
Finite Element Analysis”, Comput. Methods Appl. Mech. Engrg., 75:415-456,
1989.
[69] SIMPSON, J., Gravity Currents, 2nd. Ed., Cambridge University Press, 1997.
[70] SMAGORINSKY, J., “General Circulation Experiments with Primitive Equations,
Part I: The Basic Experiment”, Monthly Wea. Rev. 91:99-152, 1963.
[71] SOUZA, D. A. F., Algoritmo Adaptativo Implícito/Explícito por Arestas para a
Solução de Problemas de Transporte Tridimensionais, dissertação de MSc,
Programa de Engenharia Civil, COPPE/UFRJ, 2002.
[72] SOUZA, D. A. F.; COUTINHO, A. L. G. A.; TEZDUYAR, T. E.. “New-Shock-
Capturing Operators for Edge-Based SUPG Finite Element Computations of
3D Euler Equations”. In: ASME International Mechanical Engineering
Congress and Exhibition, 2006, Chicago. Proceedings of ASME International
Mechanical Engineering Congress and Exhibition, 2006
[73] SOUZA, D. A. F., MARTINS, M. A. D., COUTINHO, A. L. G. A., “Edge-based
Adaptive Implicit/Explicit Finite Element Procedures for Three Dimensional
Transport Problems”, Comm. Num. Meth. Engrg, 2005, 21:545-552.
115
[74] STEFAN, G. E ADOLFY, H., Performance Optimization of Numerical Intensive
Codes, SIAM, 2001
[75] TANAKA, N., “Development of a Highly Accurate Interpolation Method for
Mesh-Free Flow Simulations. I. Integration of Gridless, Particle and CIP
Methods”, Int. J. Num. Meth. Fluids, 30:957-976, 1999
[76] TEZDUYAR, T. E., http://www.mems.rice.edu/TAFSM/PUB_PRE/tezduyar.html
[77] TEZDUYAR, T. E., “Interface-Tracking And Interface-Capturing Techniques For
Finite Element Computation of Moving Boundaries And Interfaces”, Comput.
Methods Appl. Mech. Engrg, 195:2983-3000, 2006.
[78] TEZDUYAR, T. E., LIOU, J., “Adaptive Implicit-Explicit Finite Element
Algorithms for Fluid Mechanics Problems”, Comput. Methods Appl. Mech.
Engrg, 78:165-179, 1990.
[79] TEZDUYAR, T. E., MITTAL, S., RAY, S. E., e SHIH, R., “Incompressible Flow
Computations with Stabilized Bilinear and Linear Equal-Order-Interpolation
Velocity-Pressure Elements”. Comput. Methods Appl. Mech. Engrg.; 95:221-
242, 1992.
[80] TEZDUYAR, T. E., “Stabilized Finite Element Formulations for Incompressible
Flow Computations”, Advances in Applied Mechanics, 28:1-44, 1992.
[81] The VTK user’s guide: Install, Use and Extend the Visualization Toolkit, Kitware,
Inc. publishers, ISBN 1-930934-13-0.
[82] WESSELING, P. Principles of Computational Fluid Dynamics, Springer, 2001.
[83] VALLI, A. M. P., CAREY, G. F. e COUTINHO, A. L. G. A., “Control Strategies
for Timestep Selection in Finite Element Simulation of Incompressible Flows
and Coupled Reaction-Convection-Diffusion Processes”, Int. J. Num. Meth.
Fluids, 47:201-231, 2005.
[84] VALLI, A. M. P., Estratégias de Controle para a Seleção de Passo de Tempo para
Análise de Escoamentos Incompressíveis Acoplados com Transporte de Calor
116
e Massa, Tese de doutorado, programa de engenharia civil, COPPE/UFRJ,
2001.
[85] VAN GIJZEN, M. B., “Large Scale Finite Element Computations with GMRES-
like Methods on a CRAY Y-MP, Appl. Numer. Math., 19:51-62, 1995.
[86] VENKATAKRISHNAN, V. MAVRIPLIS, D. J., “Implicit Solvers on
Unstructured Meshes”, J. Comp. Physics, 105:83-91, 1993.
[87] WESSELING, P., Principles of Computational Fluid Dynamics, Springer, 2001.
[88] ZALESAK, S. T.. “Fully multidimensional flux-corrected transport algorithms for
fluids”. Journal of Computational Physics; 31:335-362, 1979
[89] ZIENKIEWICZ, O. C., e TAYLOR, R. L., The Finite Element Method - Volume 3:
Fluid Dynamics, Butterworth & Heinemann, 5a. edição, 2000, ISBN 0-7506-
5050-8.
117
Apêndices
118
Apêndice A
Escoamento de fluidos visco-plásticos
Este trabalho foi originalmente apresentado no FEF05 - Thirteenth Conference on
Finite Elements for Flow Problems, realizado no período de 04 a 06 de abril de 2005 em
Swansea, País de Gales, Reino Unido e posteriormente convidado para publicação na
Computational Mechanics.
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
150
Apêndice B
Tratamento de dados por arestas
Trabalho submetido ao International Journal for Numerical Methods in Engineering
após apresentação no VECPAR’06 – Seventh International Meeting on High
Performance Computing for Computational Science, realizado no período de 10 a 13 de
julho de 2006 no Rio de Janeiro, Brasil com o título: “EdgePack: A Parallel Vertex and
Node Reordering Package for Optimizing Edge-based Computations in Unstructured
Grids”.
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
Apêndice C
Cálculo de funções-distância para a
simulação utilizando level-sets
Trabalho destinado ao uso com o método das curvas de nível e desenvolvido para
solucionar o cálculo de funções distâncias de geometrias complexas em malhas não
estruturadas. O método e algoritmo foram submetidos para aprovação do International
Journal for Numerical Methods in Engineering em 2006 e aprovado para publicação em
2007.
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
Apêndice D
Escoamento de superfície-livre
Trabalho originalmente apresentado no 7th World Congress on Computational
Mechanics realizado em Los Angeles, EUA no período de 16 a 22 de julho de 2006 e
posteriormente convidado para publicação no International Journal for Numerical
Methods in Fluids.
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
Livros Grátis
( http://www.livrosgratis.com.br )
Milhares de Livros para Download:
Baixar livros de Administração
Baixar livros de Agronomia
Baixar livros de Arquitetura
Baixar livros de Artes
Baixar livros de Astronomia
Baixar livros de Biologia Geral
Baixar livros de Ciência da Computação
Baixar livros de Ciência da Informação
Baixar livros de Ciência Política
Baixar livros de Ciências da Saúde
Baixar livros de Comunicação
Baixar livros do Conselho Nacional de Educação - CNE
Baixar livros de Defesa civil
Baixar livros de Direito
Baixar livros de Direitos humanos
Baixar livros de Economia
Baixar livros de Economia Doméstica
Baixar livros de Educação
Baixar livros de Educação - Trânsito
Baixar livros de Educação Física
Baixar livros de Engenharia Aeroespacial
Baixar livros de Farmácia
Baixar livros de Filosofia
Baixar livros de Física
Baixar livros de Geociências
Baixar livros de Geografia
Baixar livros de História
Baixar livros de Línguas
Baixar livros de Literatura
Baixar livros de Literatura de Cordel
Baixar livros de Literatura Infantil
Baixar livros de Matemática
Baixar livros de Medicina
Baixar livros de Medicina Veterinária
Baixar livros de Meio Ambiente
Baixar livros de Meteorologia
Baixar Monografias e TCC
Baixar livros Multidisciplinar
Baixar livros de Música
Baixar livros de Psicologia
Baixar livros de Química
Baixar livros de Saúde Coletiva
Baixar livros de Serviço Social
Baixar livros de Sociologia
Baixar livros de Teologia
Baixar livros de Trabalho
Baixar livros de Turismo