ads:
Fabio Issao Nakamura
Animação Interativa de Fluido baseada em Partículas pelo
Método SPH
Dissertação de Mestrado
Dissertação apresentada como requisito parcial para
obtenção do título de Mestre pelo Programa de Pós-
Graduação em Informática da PUC-Rio.
Orientador: Prof. Waldemar Celes Filho
Rio de Janeiro, março de 2007
PUC-Rio - Certificação Digital Nº 0511019/CA
ads:
Livros Grátis
http://www.livrosgratis.com.br
Milhares de livros grátis para download.
Fabio Issao Nakamura
Animação Interativa de Fluido baseada em Partículas pelo
Método SPH
Dissertação apresentada como requisito parcial para
obtenção do título de Mestre pelo Programa de Pós-
Graduação em Informática da PUC-Rio. Aprovada pela
Comissão Examinadora abaixo assinada.
Prof. Waldemar Celes Filho
Orientador
Departamento de Informática – PUC-Rio
Prof. Marcelo de Andrade Dreux
Departamento de Engenharia Mecânica – PUC-Rio
Prof. Bruno Feijó
Departamento de Informática – PUC-Rio
Prof. José Eugênio Leal
Coordenador(a) Setorial do Centro Técnico Científico - PUC-Rio
Rio de Janeiro, 30 de março de 2007
PUC-Rio - Certificação Digital Nº 0511019/CA
ads:
Todos os direitos reservados. É proibida a reprodução total
ou parcial do trabalho sem autorização da universidade, do
autor e do orientador.
Fabio Issao Nakamura
Graduou-se em Engenharia da Computação pela PUC-Rio
em 2004. Durante o mestrado foi bolsista CAPES e
pesquisador no laboratório de Computação Gráfica,
Tecgraf, em projetos financiados principalmente pela
Petrobrás
Ficha Catalográfica
Nakamura, Fabio Issao
Animação interativa de fluido baseada em partículas
pelo método SPH / Fabio Issao Nakamura ; orientador:
Waldemar Celes Filho. – 2007.
64 f. : il. ; 30 cm
Dissertação (Mestrado em Informática)–Pontifícia
Universidade Católica do Rio de Janeiro, Rio de Janeiro,
2007.
Inclui bibliografia
1. Informática – Teses. 2. Animação. 3. Simulação de
fluidos. 4. SPH. I. Celes Filho, Waldemar. II. Pontifícia
Universidade Católica do Rio de Janeiro. Departamento de
Informática. III. Título.
CDD: 004
PUC-Rio - Certificação Digital Nº 0511019/CA
Agradecimentos
Ao meu orientador prof. Waldemar Celes e ao prof. Marcelo de Andrade Dreux,
que sempre me apoiaram e me orientaram para a conclusão deste trabalho.
Agradeço à minha família em especial à minha mãe, que sempre acreditou no
meu potencial e sempre me deu forças para continuar.
Um obrigado especial à Patrícia, minha namorada, por estar sempre ao meu lado
nos momentos mais difíceis em que eu mais precisei.
À CAPES e à PUC-Rio, pelos auxílios concedidos.
PUC-Rio - Certificação Digital Nº 0511019/CA
Resumo
Nakamura, Fabio Issao; Celes, Waldemar. Animação Interativade Fluido
baseada em Partículas pelo Método SPH. Rio de Janeiro, 2007. 64p.
Dissertação de Mestrado - Departamento de Informática, Pontifícia
Universidade Católica do Rio de Janeiro.
Neste trabalho foi feito um estudo investigativo sobre animação de fluidos
utilizando sistemas de partículas. Baseado nas propostas apresentadas por Muller
et al., esta dissertação objetiva investigar e compreender o uso do método
Lagrangeano baseado em partículas, conhecido como Smoothed Particle
Hydrodynamics (SPH), para simulação de fluidos. A validação do método foi feita
através da implementação de uma biblioteca capaz de animar fluidos a taxas
interativas. Para testar a eficácia e eficiência do método, a biblioteca desenvolvida
permite a instanciação de diferentes configurações, incluindo o tratamento de
colisões do fluido com obstáculos, o tratamento da interação entre dois fluidos
distintos e o tratamento de forças externas exercidas pelo usuário via um
mecanismo de interação.
Palavras-chave
Animação; Simulação de Fluidos; SPH;
PUC-Rio - Certificação Digital Nº 0511019/CA
Abstract
Nakamura, Fabio Issao; Celes, Waldemar. Fluid Interactive Animation
Based on Particle System using SPH Method. Rio de Janeiro, 2007. 64p.
MSc. Dissertation - Departamento de Informática, Pontifícia Universidade
Católica do Rio de Janeiro.
This work investigates the use of particle-based system for fluid animation.
Based on proposals presented by Müller et al., the goal of this dissertation is to
investigate and fully understand the use of a Lagrangian method known as
Smoothed Particle Hydrodynamics (SPH) for fluid simulations. A library has been
implemented in order to validate the method for fluid animation at interactive rate.
To demonstrate the method effectiveness and efficiency, the resulting library
allows the instantiation of different configurations, including the treatment of
fluid-obstacle collisions, interaction between two distinct fluids, and fluid-user
interaction.
Keywords
Animation; Fluid Simulation; SPH;
PUC-Rio - Certificação Digital Nº 0511019/CA
Sumário
1 Introdução 13
2 Trabalhos Relacionados 15
3 Método SPH 19
3.1. Formalismo Básico 19
3.2. Aproximação por Partículas 21
3.3. Propriedades da função de suavização 23
4 SPH em Animação de Fluido 26
4.1. Modelagem de Fluidos com Partículas 26
4.1.1. Cálculo da Massa específica 28
4.1.2. Força de Pressão 28
4.1.3. Força de Viscosidade 29
4.1.4. Tensão Superficial 30
4.1.5. Forças Externas 31
4.1.6. Interação entre Fluidos 31
4.2. Funções de Suavização 33
5 Implementação 36
5.1. Estruturação da biblioteca 36
5.2. Evolução do Sistema 38
5.3. Discussão 40
6 Resultados 43
6.1. Propriedades de Fluidos 43
6.1.1. Volume 43
6.1.2. Viscosidade 44
6.1.3. Tensão Superficial 45
PUC-Rio - Certificação Digital Nº 0511019/CA
6.2. Copo de Água 47
6.3. Interação entre Fluidos 49
6.4. Quebra de Barragem 50
6.5. Descarga de Água 54
6.6. Custo Computacional 56
7 Conclusão 61
8 Bibliografia 63
PUC-Rio - Certificação Digital Nº 0511019/CA
Lista de figuras
Figura 1 Trabalho de Müller et al. [24]. O enchimento do copo com a água é
simulado a uma taxa de 5 quadros por segundo. Imagens extraídas de [24]. 17
Figura 2 Alguns exemplos obtidos no trabalho de Müller et al. [26], demonstrando
a interação entre fluidos. A esquerda fluidos com massas específicas
diferentes são misturados, ao meio a simulação de um “lava lamp” e à direita
a iteração entre 3 fluidos, água, ar e a chama. Imagens extraídas de [26]. 17
Figura 3 As partículas claras possuem massa específica de repouso X. As
partículas escuras possuem massa específica de repouso 2X. As setas
indicam o gradiente de massa específica e conseqüentemente o gradiente de
pressão resultante na interface dos fluidos. 32
Figura 4 Gráfico da função de suavização W
poly6
e suas derivadas usada por
Müller et al [24]. O eixo das ordenadas representa o valor da função de
suavização e o eixo das abscissas representa a distância entre os pontos. O
gráfico mais à esquerda representa a função de suavização. O gráfico do
meio representa o gradiente de W
poly6
e o mais à direita representa o
Laplaciano de W
poly6
. 33
Figura 5 Gráfico da função de suavização W
Spike
e suas derivadas usada por
Müller et al [24]. O eixo das ordenadas representa o valor da função de
suavização e o eixo das abscissas representa a distância entre os pontos.O
gráfico mais à esquerda representa a função de suavização.O gráfico do meio
representa o gradiente de W
Spike
e o mais à direita representa o Laplaciano de
W
Spike
. 34
Figura 6 Gráfico da função de suavização W
Viscosity
e suas derivadas usada por
Müller et al [24]. O eixo das ordenadas representa o valor da função de
suavização e o eixo das abscissas representa a distância entre os pontos.O
gráfico mais à esquerda representa a função de suavização.O gráfico do meio
representa o gradiente de W
Viscosity
e o mais à direita representa o Laplaciano
de W
Viscosity
. 35
Figura 7 Estruturação da biblioteca desenvolvida. 36
Figura 8 A célula cinza clara possui a partícula de interesse. A busca para
PUC-Rio - Certificação Digital Nº 0511019/CA
determinar quais partículas estarão interagindo com a partícula de interesse
será feita na própria célula e nas células em cinza escuro. 37
Figura 9 Uma caixa enchida com um fluido representado por 8000 partículas. A
linha reta (em branco) foi utilizada para a amostragem dos valores da norma
do vetor normal à superfície e de Cs. 40
Figura 10 O gráfico com os valores da norma do vetor normal à superfície
amostrados em diferentes pontos do fluido. O eixo das abscissas do gráfico
representa o eixo X. 41
Figura 11 O gráfico com os valores do campo de cor (Cs) amostrados em
diferentes pontos do fluido. O eixo das abscissas do gráfico representa o eixo
X. 41
Figura 12 Comparação entre o volume dos fluidos. O fluido da imagem à
esquerda possui volume de 50 litros e o fluido da imagem à direita possui
volume de 10 litros. 44
Figura 13 Configuração inicial do fluido utilizado no experimento. 44
Figura 14 Diferença de comportamento entre o fluido menos viscoso (imagem à
esquerda) e o fluido mais viscoso (imagem à direita). 45
Figura 15 A imagem à esquerda ilustra um fluido com o valor do coeficiente de
tensão superficial 0 e à direita o mesmo fluido com o valor do coeficiente de
tensão superficial 1. 46
Figura 16 A imagem de cima, ilustra um fluido com o valor do coeficiente de
tensão superficial 0 e a imagem de baixo o mesmo fluido com o valor do
coeficiente de tensão superficial 1. 46
Figura 17 Três instantes da simulação. As imagens da esquerda mostram os
resultados obtidos no trabalho desenvolvido. As imagens da direita foram
retiradas do trabalho de Müller et al.[24]. 48
Figura 18 Comparação entre os resultados obtidos e os resultados obtidos em [26].
As imagens da esquerda foram obtidas através da biblioteca desenvolvida, e
as imagens da direita foram retiradas de [26]. 49
Figura 19 Instante inicial da simulação. Um fluido está confinado por uma
barragem. A imagem mais acima foi retirada de [17]. A imagem do meio
mostra a visualização ortográfica da simulação e a imagem mais abaixo
mostra a vista perspectiva da simulação. 51
PUC-Rio - Certificação Digital Nº 0511019/CA
Figura 20 Instante após a abertura da barragem em que o fluido encontra-se com o
obstáculo (rampa). A imagem mais acima foi retirada de [17]. A imagem do
meio mostra a visualização ortográfica da simulação e a imagem mais abaixo
mostra a vista perspectiva da simulação. 52
Figura 21 O fluido interagindo com a rampa (obstáculo) formando uma onda. A
imagem mais acima foi retirada de [17]. A imagem mais acima foi retirada
de [17]. A imagem do meio mostra a visualização ortográfica da simulação e
a imagem mais abaixo mostra a vista perspectiva da simulação. 53
Figura 22 Instante que a comporta é aberta para vazão da água. A imagem da
esquerda foi retirada de [17], a imagem do meio e a imagem da direita são a
vista ortográfica e perspectiva, respectivamente, da simulação obtida neste
trabalho. 54
Figura 23 Neste instante devido à forte força de pressão as partículas de água são
expelidas pelo compartimento aberto formando uma cavidade. A imagem da
esquerda foi retirada de [17], a imagem do meio e a imagem da direita são a
vista ortográfica e perspectiva, respectivamente, da simulação obtida neste
trabalho. 55
Figura 24 Diferença entre os resultados obtidos em [17] e os resultados obtidos
neste trabalho. A imagem da esquerda foi retirada de [17], a imagem do meio
e a imagem da direita são a vista ortográfica e perspectiva, respectivamente,
da simulação obtida neste trabalho. 55
Figura 25 Gráfico dos dados apresentados na Tabela 2. 57
Figura 26 Gargalo da simulação com 1000 partículas. 58
Figura 27 Identificação do gargalo com 4000 partículas. 59
Figura 28 Identificação do gargalo com 8000 partículas. 59
PUC-Rio - Certificação Digital Nº 0511019/CA
Lista de tabelas
Tabela 1 Os principais atributos das partículas de um fluido presentes neste
trabalho. Os últimos 3 atributos são necessários para o caso de interação
entre fluidos (Seção 4.1.6). 38
Tabela 2 Valores obtidos na comparação entre o número de partículas e o tempo
médio de simulação por passo. 56
PUC-Rio - Certificação Digital Nº 0511019/CA
1
Introdução
Na natureza, podemos encontrar fluidos em todos os lugares, seja na forma
líquida, seja na forma gasosa. Muitos fenômenos físicos podem ser explicados e
estudados através da dinâmica de fluidos. Entre eles podemos citar: escoamento
de líquidos e gases, movimento das marés, estudo de fenômenos atmosféricos, etc.
Dinâmica de fluidos estuda os efeitos das forças aplicadas ao fluido.
Fluidos respeitam as leis de conservação de massa, de quantidade de movimento
(momento angular e momento linear) e de energia. As leis de conservação que
regem o comportamento dos fluidos são descritas por equações diferenciais.
Muitas vezes não é possível encontrar uma solução analítica para resolver tais
equações. Por isso, Dinâmica de Fluidos Computacional (CFD – Computational
Fluid Dynamics) visa encontrar, através de métodos numéricos uma aproximação
satisfatória capaz de descrever com certa precisão o comportamento de fluidos.
Mesmo com os atuais computadores, ainda é praticamente impossível se
obter uma simulação com alto grau de precisão em tempo real, devido à
complexidade das equações e dos sofisticados métodos numéricos necessários
para tratá-las. Por outro lado, se somente for necessário uma aproximação do real
comportamento de um fluido para efeito de animação, é possível conseguir
resultados satisfatórios em tempo real, sem a necessidade de máquinas especiais
ou modelos numéricos complexos. Para tais animações existe uma variedade de
aplicações, entre elas: simulações médicas, jogos, ambientes virtuais, animações
cinematográficas. Modelos físicos mais simplificados podem também ser usados
como protótipos para validar simulações mais realistas.
O objetivo deste trabalho é investigar a simulação de fluidos em aplicação
para animação. Baseado nos trabalhos de Müller et al. [24, 26], esta dissertação
objetiva investigar e compreender o uso da abordagem Lagrangeana baseada em
partículas, conhecida como Smoothed Particle Hydrodynamics – SPH, para
animação de fluido em tempo interativo.
PUC-Rio - Certificação Digital Nº 0511019/CA
Introdução 14
Para isso foi desenvolvida uma biblioteca para simulação de fluidos,
reproduzindo os trabalhos descritos em [24, 26]. Este desenvolvimento permitiu
avaliar e validar o uso do método SPH em aplicação de tempo real, para animação
de fluidos baseada em modelos físicos. As principais contribuições desta
dissertação são:
• Compreender, implementar e validar o uso do método SPH para
animação de fluidos como descrito em [24, 26].
• Experimentar e avaliar diferentes animações conseguidas em
diferentes configurações.
• Analisar os recursos computacionais necessários de uma
implementação do método SPH em CPU para aplicação de
animação em tempo real.
Esta dissertação está organizada da seguinte maneira. No capítulo 2 alguns
trabalhos relacionados serão discutidos. O capítulo 3 descreve os fundamentos
básicos do método SPH. O capítulo 4 descreve como o método SPH pode ser
aplicado para solucionar a simulação de fluidos utilizando partículas. O capítulo 5
apresenta uma breve descrição da estruturação da biblioteca desenvolvida. O
capítulo 6 descreve os testes realizados e os resultados obtidos. Finalmente, no
capítulo 7, a conclusão e sugestões para trabalhos futuros são apresentadas.
PUC-Rio - Certificação Digital Nº 0511019/CA
2
Trabalhos Relacionados
Podemos destacar dois principais modelos usados em simulação
computacional de fluidos, são eles: o modelo Euleriano baseado em grid e o
modelo Lagrangeano baseado em partículas. O modelo Euleriano baseado em grid
divide (discretizam) o domínio do problema em células que são fixas no espaço
durante toda a simulação. As grandezas (massa, massa específica, energia, etc)
são simuladas através do fluxo que passa pelas células. Já o modelo Lagrangeano,
baseado em partículas, não necessita da discretização do domínio do problema em
células. No modelo Lagrangeano baseado em partículas, as grandezas (massa,
massa específica, energia, etc) são representadas por partículas e tais partículas
movem-se junto com o fluido, ao contrário do modelo Euleriano baseado em grid,
que permanece estático durante toda a simulação. Ambos os modelos apresentam
vantagens e desvantagens.
Em relação aos modelos Eulerianos, podemos citar alguns trabalhos que
tiveram maior importância. Em 1990, Kass e Miller [14] propuseram um novo
método para animação de água e ondas. O método era baseado em uma versão
simplificada das equações de água rasa em 2D, junto com um campo de altura
para representar a superfície do fluido. Em 1996, Foster e Metaxas [10, 11] foram
os primeiros a proporem uma maneira capaz de resolver a Equação de Navier-
Stokes em três dimensões, utilizando um grid regular. O trabalho também
apresenta um algoritmo para extração da superfície livre do fluido. Tratamentos
de fronteira contra corpos sólidos e contra o ar, também faziam parte da
simulação. Apesar do trabalho de Foster e Metaxas ser capaz de animar fluidos
com alto grau de precisão, a simulação não consegue animar em tempo real,
mesmo com os atuais computadores, devido à necessidade de um passo de
simulação muito pequeno para garantir a estabilidade.
Em 1999, Stam [32] propôs um esquema semi Lagrangeano capaz de
contornar o problema do passo de simulação encontrado em Foster e Metaxas [10,
11]. O modelo apresentado por Stam [32] permitia simulações estáveis mesmo
PUC-Rio - Certificação Digital Nº 0511019/CA
Trabalhos Relacionados 16
com passos de simulação grandes ao custo de forças dissipativas artificiais. No
entanto, o modelo não trata fluidos de superfície livre como água, devido à
dificuldade de se extrair a superfície do fluido. Enright e Fedkiw [8] estenderam o
trabalho de Stam [32] para lidar com líquidos usando modelos híbridos com
partículas para auxiliar a extração da superfície do fluido. Carlson et al.[3]
desenvolveram um sistema Euleriano capaz de animar fluidos líquidos viscosos e
objetos capazes de simular o derretimento de objetos.
Com o aumento do poder computacional dos processadores gráficos
(GPUs), diversos trabalhos [13, 33] estão sendo propostos para usar GPUs na
animação de fluidos, obtendo excelentes resultados e tempo de animação
interativo.
Em 1977, Lucy [20] e Gingold e Monaghan [21] desenvolveram o método
conhecido como Smoothed Particle Hydrodynamics (SPH). A priori o método foi
desenvolvido para lidar com problemas astrofísicos. Reeves [31], em 1983,
introduziu na comunidade de computação gráfica sistemas de partículas para
modelar uma classe de objetos fuzzy. Sistema de partículas também é usado em
diversos tipos de efeitos, como explosões, chuva e neve. Miller e Pearce [21]
juntaram em seus trabalhos idéias de dinâmica moleculares com um sistema de
partículas interativo (onde uma partícula interage com outras partículas) para
modelagem de materiais deformáveis e materiais líquidos como lava, lama, óleo,
entre outros.
Desbrun e Gascuel [5], foram pioneiros ao aplicarem o método SPH na área
de computação gráfica, a princípio, para simular substâncias altamente
deformáveis.
Premože et al. [29], utilizando o método Lagrangeano Moving Particles
Semi Implicit (MPS) [16], uma variante do método SPH, obtiveram resultados
bastante realistas na simulação de fluidos incompressíveis.
Müller et al. [24], em 2003 apresentaram uma simulação e visualização de
fluidos, baseado em partículas, em tempo interativo. Müller et al. [24] incluíram
tensão superficial ao fluido, baseado nas idéias de Morris [23] (Figura 1).
PUC-Rio - Certificação Digital Nº 0511019/CA
Trabalhos Relacionados 17
Figura 1 Trabalho de Müller et al. [24]. O enchimento do copo com a água é simulado a
uma taxa de 5 quadros por segundo. Imagens extraídas de [24].
Em 2004 Müller et al. [25] desenvolveram um trabalho focado na
interação entre fluidos e corpos sólidos deformáveis. O trabalho desenvolvido
permite além da simulação do fluido em tempo interativo, a simulação das forças
de repulsão, adesão e atrito entre o fluido e os corpos sólidos.
Clavet et al. [4] propuseram um novo método para simulação de fluidos
visco elásticos baseado em partículas. No trabalho desenvolvido, Clavet et al. [4]
forçam tanto a incompressibilidade do fluido quanto a não formação de
aglomerações de partículas através do método chamado de double density
relaxation. O método apresentado é robusto, estável e capaz de simular fluidos em
diversos cenários aà taxas interativas.
Müller et al. [26] apresentaram uma proposta para simulação e interação
entre fluidos. O trabalho estende o trabalho de Müller et al. [24] acrescentando
forças de interface entre fluidos distintos. A proposta apresentada é capaz de
simular a interação de fluidos de diferentes características, permitindo por
exemplo a simulação de um objeto “lava lamp” (Figura 2).
Figura 2 Alguns exemplos obtidos no trabalho de Müller et al. [26], demonstrando a
interação entre fluidos. A esquerda fluidos com massas específicas diferentes são
misturados, ao meio a simulação de um “lava lamp” e à direita a iteração entre 3 fluidos,
água, ar e a chama. Imagens extraídas de [26].
PUC-Rio - Certificação Digital Nº 0511019/CA
Trabalhos Relacionados 18
Em 2007 Paiva et al.[1] apresentaram um trabalho propondo uma técnica
capaz de animar o derretimento de objetos utilizando simulação de fluidos não
newtonianos, baseado em partículas. A técnica consiste em modelar objetos
através da transição de fluidos não newtonianos altamente viscosos para líquidos
pouco viscosos. Para isto, Paiva et al.[1] utilizam uma variação do modelo
lagrangeano SPH.
A biblioteca desenvolvida neste trabalho, segue as propostas feitas por
Müller et al. [24, 26]. Os resultados obtidos no trabalho de Müller et al. [24, 26],
nos levou a crer que seria uma ótima forma de validar o uso do método SPH para
simulação e animação de fluidos em tempo interativo.
PUC-Rio - Certificação Digital Nº 0511019/CA
3
Método SPH
SPH (Smoothed Particle Hydrodynamics) é um método numérico
Lagrangeano baseado em partículas. O método foi desenvolvido em 1977 por
Gingold e Monaghan [12] e Lucy [20] para o tratamento de problemas
astrofísicos. Por possuir uma formulação simples e por ser baseado em partículas,
o método é eficiente o suficiente para permitir a simulação em tempo real de
certos fenômenos físicos (dependendo do grau de realismo) sendo geral o
suficiente para ser aplicado a vários tipos de problemas como por exemplo:
dinâmica de fluidos, deformações de materiais, ondas de choque, fluxos de lava,
etc.
Os fundamentos do método SPH estão na teoria da interpolação. O método
utiliza funções de suavização chamados de kernels. Para se determinar o valor de
uma grandeza qualquer num ponto
P
r
no espaço, o método, através das funções
de suavização, interpola valores amostrados dentro da vizinhança do ponto
P
r
.
Nas seções seguintes, iremos apresentar uma formulação matemática do método.
A formulação apresentada segue à exposta por Liu e Liu [17].
3.1.
Formalismo Básico
Qualquer função ƒ definida em um domínio de interesse (domínio de
suporte), representando uma propriedade física ou uma massa específica, pode ser
aproximada através dos valores de um conjunto de amostras aleatórias ponderadas
por uma função de suavização com raio finito.
Dada uma função ƒ definida sobre todo o domínio Ω podemos expressar f
como:
(
)
xdxxxfxf
r
r
r
r
r
′′
−
′
=
∫
Ω
)()(
δ
(3.1)
PUC-Rio - Certificação Digital Nº 0511019/CA
Método SPH 20
onde Ω é o volume da integral que contém x
r
, e )( xx
′
−
r
r
δ
é a função delta de
Dirac definida por :
xxse
xxse
xx
′
≠
′
=
→
=
′
−
rr
r
r
rr
0
1
)(
δ
(3.2)
A Equação 3.1 implica que uma função ƒ pode ser representada através de
uma integral. Esta integral é chamada de representação integral da função ƒ. Se ƒ
for definida e contínua em Ω, sua representação integral é garantidamente exata.
Se substituirmos a função delta de Dirac por uma função de suavização do
tipo
(
)
hxxW ,
′
−
r
r
, a aproximação de ƒ na forma integral será dada por:
∫
Ω
′′
−
′
= xdhxxWxfxf
r
r
r
r
r
),()()( (3.3)
Uma vez que a função de suavização W não é exatamente igual a função
delta de Dirac, a Equação 3.3 representa uma aproximação da Equação 3.1. Na
função de suavização, o termo h representa o comprimento (raio) de suavização
(smoothing length) que define a área (ou volume) de influência da função de
suavização.
Em diversas situações, há a necessidade de se aproximar a derivada da
função f ao invés da própria função f . Aplicando a representação integral à
aproximação da derivada da função f obtemos:
∫
Ω
′′
−
′
∇=∇ xdhxxWxfxf
r
r
r
r
r
r
r
),()]([)( (3.4)
Aplicando a integração por partes, temos:
∫∫
ΩΩ
′′
−∇
′
−
′′
−
′
∇=∇ xdhxxWxfxdhxxWxfxf
r
r
r
r
r
r
r
r
r
r
r
r
),()()],()([)( (3.5)
Usando o teorema da divergência (Teorema de Gauss), podemos transformar a
primeira integral do lado direito da Equação 3.5 em uma integral de superfície.
Com isso obtemos:
∫∫
Ω
′′
−∇
′
−
′
−
′
=∇ xdhxxWxfdSnhxxWxfxf
S
r
r
r
r
r
r
r
r
r
r
r
),()(),()()( (3.6)
PUC-Rio - Certificação Digital Nº 0511019/CA
Método SPH 21
onde
n
r
é o vetor normal à superfície. Se o domínio de suporte estiver contido
dentro do domínio do problema (domínio onde as equações parciais diferenciais
estão contidas), então a integral de superfície é zero e podemos reescrever a
Equação 3.6 da seguinte forma:
∫
Ω
′′
−∇
′
−=∇ xdhxxWxfxf
r
r
r
r
r
r
r
),()()( (3.7)
Ou seja, a aproximação da derivada da função f é a aproximação integral da
função f usando-se o gradiente da função de suavização
(
)
hxxW ,
′
−∇
r
r
r
, ao invés
de simplesmente a função de suavização
(
)
hxxW ,
′
−
r
r
. Pode-se mostrar também
que o mesmo processo é válido para a segunda derivada de f (Laplaciano de f):
∫
Ω
′′
−∇
′
=∇ xdhxxWxfxf
r
r
r
r
r
),()()(
22
(3.8)
3.2.
Aproximação por Partículas
No método SPH, o sistema é representado por um número finito de
partículas que possuem massa e ocupam um volume do espaço. A representação
integral da função ƒ (Equação 3.3) pode ser escrita de forma discreta, aplicando o
somatório sobre todas as partículas dentro do domínio de suporte. Este processo
chama-se aproximação por partículas.
Se o volume infinitesimal xd
′
r
da representação integral de f na posição da j-
ésima partícula for substituído pelo volume ocupado pela partícula,
j
V
∆
, teremos:
∑
∆−=
j
jj
VhxxWxfxf ),()()(
rrrr
(3.9)
Se substituirmos
j
V
∆
por
j
j
m
ρ
, onde
j
m e
j
ρ
são a massa e a massa específica da
partícula j, respectivamente, podemos reescrever a Equação 3.9 como:
∑
−=
j
jj
j
j
hxxWxf
m
xf ),()()(
rrrr
ρ
(3.10)
PUC-Rio - Certificação Digital Nº 0511019/CA
Método SPH 22
Ou seja, podemos aproximar o valor da função f na posição
x
r
através dos valores
amostrados na vizinhança de x
r
, ponderados por uma função de suavização.
Da mesma forma que podemos aproximar a derivada da função f, através
de sua representação integral, podemos aproximar a derivada de f através da
aproximação por partículas. Se aplicarmos o mesmo raciocínio que aplicamos
para determinar a Equação 3.7, obteremos a seguinte aproximação por partículas
para a derivada de f.
∑
−∇−=∇
j
jjj
j
j
hxxWxf
m
xf ),()()(
rr
r
rr
r
ρ
(3.11)
Note que o gradiente ),( hxxW
j
r
r
−
∇
é referente à partícula j. Como queremos o
valor da aproximação da derivada em relação à posição x
r
, a Equação 3.11 passa a
ser:
∑
−∇=∇
j
jj
j
j
hxxWxf
m
xf ),()()(
rr
r
rr
r
ρ
(3.12)
O mesmo desenvolvimento pode ser aplicado à aproximação do
Laplaciano da função f na posição
x
r
, resultando em:
∑
−∇=∇
j
jj
j
j
hxxWxf
m
xf ),()()(
22
rrrr
ρ
(3.13)
As Equações 3.12 e 3.13 fazem com que o método SPH seja simples e ao
mesmo tempo eficiente. Em inúmeras situações ocorre a necessidade de se
calcular tanto o gradiente quanto o Laplaciano da função f. Com o método SPH,
podemos obter os valores aproximados tanto do gradiente quanto do Laplaciano
de f sem ter que derivar a função f. Para isto basta derivar a função de suavização
que, na maioria dos casos, é muito mais simples que a função f.
Resumindo, podemos aproximar o valor de um campo escalar ou vetorial
A, na posição x
r
, através da seguinte equação:
∑
−=
j
jj
j
j
i
hxxWxA
m
xA ),()()(
rrrr
ρ
(3.14)
PUC-Rio - Certificação Digital Nº 0511019/CA
Método SPH 23
A derivada e o laplaciano do campo escalar ou vetorial será dado pelas
Equações 3.15 e 3.16, respectivamente.
∑
−∇=∇
j
jj
j
j
i
hxxWxA
m
xA ),()()(
rrr
r
r
r
ρ
(3.15)
∑
−∇=∇
j
jj
j
j
i
hxxWxA
m
xA ),()()(
22
rrrr
ρ
(3.16)
3.3.
Propriedades da função de suavização
A escolha das funções de suavização no método SPH é primordial para
estabilidade, velocidade e precisão do método. Diversas funções de suavização
são propostas na literatura [17, 24]. Uma função para ser usada como função de
suavização no método SPH deve satisfazer certas propriedades. Iremos destacar as
principais propriedades presentes nas funções de suavização encontradas na
literatura.
1. As funções de suavização devem ser normalizadas dentro do domínio
de suporte, ou seja:
∫
Ω
=
′′
− 1),( xdhxxW
r
r
(3.17)
A integral da função de suavização dentro do domínio de suporte deve
ser unitária para garantir que a interpolação seja correta.
2. As funções de suavização devem possuir um domínio de suporte
compacto, ou seja:
hxxxxW
κ
>
′
−
∀
=
′
−
||,0)(
r
r
r
r
(3.18)
O domínio de suporte define a área (diferente de zero) efetiva de
atuação da função de suavização onde k é uma constante.
3. Dependendo do tipo de problema a ser modelado, o valor da função de
suavização )( xxW
′
−
r
r
deve ser maior ou igual a zero para qualquer
ponto x
r
dentro do domínio de suporte. Esta propriedade não é
necessária para garantir convergência da função de suavização. No
PUC-Rio - Certificação Digital Nº 0511019/CA
Método SPH 24
entanto, se o problema for a modelagem de um fenômeno físico,
muitas vezes o valor de )( xxW
′
−
r
r
sendo menor que zero não tem
significado físico. Por exemplo, em problemas hidrodinâmicos, não há
sentido ter valores de massa específica negativos.
4. O valor de )( xxW
′
−
r
r
no ponto x
′
r
deve decrescer à medida que x
′
r
se
afasta de
x
r
. Ou seja, as partículas mais próximas de
x
r
devem ter
maior influência do que aquelas mais afastadas de x
r
.
5. As funções de suavização devem satisfazer à função delta de Dirac
quando o raio de influência h tender a zero, ou seja:
)(),(lim
0
xxhxxW
h
′
−
=
′
−
>−
r
r
r
r
δ
(3.19)
A Equação 3.3 garante que o valor aproximado da função no ponto x
r
se aproxima do valor real da função no ponto
x
r
. Ou
seja,
)()( xfxf
r
r
=
.
6. As funções de suavização dever ser funções pares. Esta propriedade
também é conhecida como propriedade de simetria. Partículas
afastadas da mesma distância de
x
r
devem contribuir da mesma forma,
mesmo estando em diferentes posições.
7. As funções de suavização devem ser suficientemente suaves. Para
garantir uma melhor aproximação as funções de suavização e suas
derivadas devem ser contínuas.
Se as condições acima forem respeitadas, então pode-se garantir funções de
suavização com erro de interpolação de segunda ordem. Se )(xf
r
for
diferenciável, podemos usar a expansão de Taylor sobre )(xf
′
r
em torno de x
r
,
obtendo:
∫
Ω
′′
−−
′
+−
′′
+= xdhxxWxxrxxxfxfxf
r
r
r
r
r
r
r
r
r
r
),()])(())(()([)(
2
(3.20)
onde o termo
r
representa o resíduo da série de Taylor. Podemos reescrever a
Equação 3.20 como:
PUC-Rio - Certificação Digital Nº 0511019/CA
Método SPH 25
∫
∫
Ω
Ω
+
′′
−−
′′
+
′′
−=
)(),()()(
),()()(
2
hrxdhxxWxxxf
xdhxxWxfxf
rrrrrr
r
r
r
r
r
(3.21)
Usando a Equação 3.17 e o fato da função de suavização ser uma função par (o
termo de primeira ordem desaparece), teremos:
)()()(
2
hrxfxf +=
r
r
(3.22)
Isto é, o erro de interpolação será de segunda ordem.
PUC-Rio - Certificação Digital Nº 0511019/CA
4
SPH em Animação de Fluido
O comportamento de um fluido é governado por 3 leis físicas de
conservação, são elas: conservação de massa, conservação de momento e
conservação de energia. Para simular o comportamento de um fluido devemos
resolver essas três equações. Porém, dificilmente essas equações apresentam
soluções analíticas. Por isso, usam-se métodos numéricos para obter uma
aproximação satisfatória. Neste capítulo, iremos mostrar, baseado nos trabalhos de
Müller et al. [24, 26], como utilizar o método SPH para simular, através de
partículas, o comportamento de fluidos incompressíveis e isotérmicos.
4.1.
Modelagem de Fluidos com Partículas
Para modelar o comportamento de um fluido incompressível e isotérmico,
devemos resolver o problema de conservação de massa e de conservação de
momento. Pelo fato do fluido ser isotérmico, a conservação de energia é garantida
automaticamente e não precisa ser modelada.
A incompressibilidade do fluido garante que a massa específica do fluido
não irá variar com o tempo (durante a simulação), já que a relação entre o volume
ocupado pelo fluido e sua massa será constante. A conservação de massa de um
fluido é dada pela chamada Equação da Continuidade (Equação 4.1):
0).( =∇+
∂
∂
v
t
r
r
ρ
ρ
(4.1)
onde
ρ
e v
r
são a massa específica e a velocidade do fluido, respectivamente.
Como a massa específica não varia com o tempo, podemos reescrever a Equação
4.1 da seguinte forma:
0).( =∇ v
r
r
ρ
(4.2)
Pelo fato do fluido ser modelado por partículas, não é necessário resolver a
Equação
4.2, já que a conservação de massa está automaticamente garantida. Isto
PUC-Rio - Certificação Digital Nº 0511019/CA
SPH em Animação de Fluido 27
ocorre porque o número de partículas usadas para representar o fluido é constante
e cada partícula possui massa também constante.
Por último a conservação de momento é dada pela Equação de Navier-
Stokes:
vgPvv
t
v
r
r
r
r
r
r
r
r
2
).( ∇++∇−=∇+
∂
∂
µρρ
(4.3)
A Equação 4.3 representa uma forma simplificada da Equação de Navier–Stokes
para fluidos newtonianos e incompressíveis. Novamente o uso de partículas
simplifica o problema. O lado esquerdo da Equação 4.3 pode ser substituído por
dt
d
ν
r
(derivada material), uma vez que as partículas movem-se junto com o fluido,
não sendo necessário o termo convectivo
vv
r
r
r
∇.
. Ou seja, a variação do campo de
velocidade é simplesmente a taxa de variação da velocidade das partículas no
tempo. Podemos reescrever a Equação 4.3 como sendo:
)(
1
2
vgP
dt
vd
r
r
r
r
r
∇++∇−=
µρ
ρ
(4.4)
O lado direito da Equação 4.3 representa três campos de força, são eles:
campo de pressão )( P∇−
r
, gravidade )( g
r
ρ
e força de viscosidade )(
2
v
r
∇
µ
, onde
µ
representa a viscosidade do fluido. A soma desses três campos de força
determina a força resultante aplicada ao fluido. Esta força resultante determina a
variação do momento do fluido. Com isso, a aceleração da i-ésima partícula que
representa o fluido é dado por:
i
ii
i
dt
d
ρ
fv
a
r
r
r
== (4.5)
onde
i
a
r
,
dt
d
i
ν
r
,
i
f
r
e
i
ρ
são a aceleração da i-ésima partícula, a variação da
velocidade em relação ao tempo da i-ésima partícula, a força resultante aplicada à
i-ésima partícula e a massa específica da i-ésima partícula, respectivamente.
Nas seções seguintes será mostrado como as forças atuantes (lado direito
da Equação 4.4) no fluido podem ser modeladas com o método SPH.
PUC-Rio - Certificação Digital Nº 0511019/CA
SPH em Animação de Fluido 28
4.1.1.
Cálculo da Massa específica
Para o uso do método SPH, é necessário o cálculo das massas específicas
amostradas nas posições das partículas. Aplicando o método SPH para o cálculo
da massa específica, obtemos o seguinte resultado:
),()( hxxWmx
ji
j
j
j
ji
rrr
−=
∑
ρ
ρ
ρ
(4.6)
que corresponde a:
),()( hxxWmx
ji
j
ji
r
r
r
−=
∑
ρ
(4.7)
Ou seja, a massa específica de cada partícula i é dada pelo somatório das massas
das partículas na vizinhança da partícula i, ponderadas pela função de suavização.
A Equação 4.7 também é chamada de summation density.
4.1.2.
Força de Pressão
O termo )( P∇−
r
presente na Equação 4.3 representa a variação do campo de
pressão de um fluido. Esta variação no campo de pressão gera um campo de força
de pressão. Aplicando o método SPH (Equação 3.15) ao termo )( P∇−
r
obtém-se:
∑
−∇−=∇−=
j
ji
j
j
ji
i
pres
hxxW
P
mxP ),()(
rr
r
r
r
r
ρ
f
(4.8)
A força gerada pela Equação 4.8 não é simétrica. Se interagirmos somente duas
partículas, i e j, a força de pressão atuante em i será determinada somente pela
pressão amostrada na posição da partícula j e vice versa. Uma vez que a pressão
dada nas posições das partículas i e j podem ser diferentes, teremos forças
assimétricas:
),( hxxW
P
ji
j
j
pres
i
rr
r
r
−∇−=
ρ
f
),( hxxW
P
ij
i
i
pres
j
rr
r
r
−∇−=
ρ
f
PUC-Rio - Certificação Digital Nº 0511019/CA
SPH em Animação de Fluido 29
A falta de simetria é um problema inerente ao método de SPH. Para
contornar tal problema e tornar esta força simétrica (garantindo a conservação de
momento e a 3ª lei de Newton de ação e reação), Müller et al. [24] propuseram a
seguinte solução:
∑
−∇
+
−=∇−=
j
ji
j
ij
ji
pres
i
hxxW
PP
mxP ),(
2
)(
rr
r
r
r
r
ρ
f (4.9)
A Equação 4.9 resolve o problema de simetria presente na Equação 4.8 uma vez
que é usada a média aritmética das pressões das partículas interagindo.
O cálculo do valor da pressão de cada partícula é necessário para se
determinar o gradiente de pressão e, conseqüentemente, o campo de força gerado
por ele. Usamos a mesma proposta encontrada em Müller et al. [24]. A pressão
das partículas do fluido pode ser determinada através da Equação de Estado dos
gases:
)(
ii
kP
ρ
=
(4.10)
onde k representa a constante de gás dependente da temperatura do fluido. Müller
et al. [24] utilizam uma versão modificada da Equação 4.10, proposta por
Desbrun [5]:
)(
0
ρ
ρ
−
=
ii
kP (4.11)
onde
0
ρ
representa a massa específica de repouso (rest density) do fluido e é dado
como parâmetro de simulação. Uma vez que a força de pressão depende do
gradiente de pressão, a subtração do termo
0
ρ
não afeta o campo de força de
pressão, já que não afeta o gradiente do campo de pressão.
4.1.3.
Força de Viscosidade
A força de viscosidade é dada pelo termo )(
2
v
r
∇
µ
da Equação 4.3.
Aplicando o método SPH ao termo )(
2
v
r
∇
µ
, obtém-se
∑
−∇=∇=
j
ji
j
j
ji
visc
i
hxxW
v
mxv ),()(
22
rr
r
rr
r
ρ
µµ
f (4.12)
PUC-Rio - Certificação Digital Nº 0511019/CA
SPH em Animação de Fluido 30
Novamente é gerada uma força não simétrica, já que as velocidades de cada
partícula são diferentes. Müller et al. [24] propuseram o uso da seguinte equação
para resolver o problema de simetria.
∑
−∇
−
=∇=
j
ji
j
ij
ji
visc
i
hxxW
vv
mxv ),()(
22
rr
r
r
rr
r
ρ
µµ
f (4.13)
O fato da Equação 4.13 depender somente da diferença das velocidades, e não do
valor absoluto das velocidades, naturalmente garante a conservação do momento e
o princípio de ação e reação.
4.1.4.
Tensão Superficial
Apesar de não aparecer na equação de Navier–Stokes, Müller et al. [24]
modelam as forças de tensão superficial atuantes sobre o fluido. A tensão
superficial surge devido a um desbalanceamento de forças. As moléculas do
fluido estão sujeitas a uma força de atração de suas moléculas vizinhas. Dentro do
fluido essas forças intermoleculares são iguais em todas as direções, balanceando
uma às outras. Porém, as forças atuando sobre as moléculas da superfície do
fluido não são balanceadas. A força resultante gera o que chamamos de tensão
superficial. Devido a tensão superficial, que é tangente à superfície do fluido,
surge uma força que é proporcional a curvatura do fluido, atuando na direção
inversa à normal à superfície, isto é, apontando para dentro do fluido.
Implementamos o modelo apresentado no trabalho de Müller et al.[24] que
por sua vez foi baseado nas idéias de Morris [23]. Müller et al.[24] descrevem
uma grandeza chamada de campo de cor. O campo de cor é um campo onde o
valor amostrado em cada partícula vale 1 e em qualquer outro lugar do espaço
vale 0. Aplicando o método SPH no campo de cor, obtemos o seguinte:
),(
1
)( hxxWmxCs
ji
j
j
ji
rrr
−=
∑
ρ
(4.14)
A normal da superfície do fluido é dada pelo gradiente do campo de cor:
),(
1
)( hxxWmxCsn
ji
j
j
ji
rr
r
r
r
r
−∇=∇=
∑
ρ
(4.15)
PUC-Rio - Certificação Digital Nº 0511019/CA
SPH em Animação de Fluido 31
A direção da normal está orientada para dentro do fluido, onde a curvatura da
superfície é dada pelo Laplaciano do campo de cor:
||
2
n
Cs
r
∇
−=
κ
(4.16)
o sinal de menos é necessário para garantir curvatura positiva e superfícies
convexas.
Müller et al.[24] definem então a força gerada pela tensão superficial
atuante sobre a superfície do fluido como sendo:
||
2
n
n
Csnf
surf
r
r
r
r
∇==
σκσ
(4.17)
onde
σ
é uma constante que controla o quanto de tensão superficial se deseja e é
dado como parâmetro da simulação.
4.1.5.
Forças Externas
Forças externas como gravidade, por exemplo, são aplicadas diretamente
às partículas, sem a necessidade de usar o método SPH. O problema de colisão
pode ser tratado como força externa, ou seja, a resposta à colisão é aplicada
diretamente às partículas.
4.1.6.
Interação entre Fluidos
Em [26] Müller et al. apresentam uma forma de estender a simulação de
fluidos apresentado em [24] para lidar com diferentes fluidos. No trabalho de
Müller et al. [26], os autores propõem um novo modelo de simulação capaz de
simular mudanças de fase de um fluido, interação entre ar e água, interação entre
dois fluidos e até mesmo a simulação de um “lava lamp”. Neste trabalho
abordamos apenas os aspectos necessários para a interação entre dois fluidos
(líquidos) seguindo as idéias apresentadas em [26].
Para a interação entre dois fluidos, os autores de [26] propõem a mudança
de alguns atributos do fluido considerados “globais”, isto é, iguais para todas as
partículas, para atributos que podem variar de partícula para partícula. A primeira
PUC-Rio - Certificação Digital Nº 0511019/CA
SPH em Animação de Fluido 32
modificação sugerida em [26], é permitir que as partículas possam ter valores de
coeficiente de viscosidade, constante de gás e massa específica de repouso
diferentes uma das outras. A segunda modificação é a mudança do cálculo da
força de viscosidade, para levar em conta os coeficientes de viscosidades por
partículas. Müller et al. [26] propõem a seguinte equação para força de
viscosidade:
∑
−∇
−
+
=∇
j
ji
j
ij
j
ji
i
hxxW
vv
mxv ),(
2
)(
22
rr
r
r
rr
ρ
µ
µ
µ
(4.18)
Usando a Equação 4.18 o cálculo da força de viscosidade passa a considerar a
média dos coeficientes de viscosidade das partículas interagindo.
Com essas duas modificações, é possível simular o efeito de empuxo entre
dois fluidos. O empuxo é simulado com o uso por partícula da massa específica de
repouso. A massa específica de repouso (termo
0
ρ
na Equação 4.11) é responsável
por atrair ou repelir partículas a fim de chegar a um estado de equilíbrio. Quando
dois fluidos com diferentes massas específicas são misturados, surge um gradiente
de massa específica entre as partículas e conseqüentemente um gradiente de
pressão na interface dos fluidos (Figura 3).
Figura 3 As partículas claras possuem massa específica de repouso X. As partículas
escuras possuem massa específica de repouso 2X. As setas indicam o gradiente de
massa específica e conseqüentemente o gradiente de pressão resultante na interface
dos fluidos.
O gradiente de pressão na interface dos fluidos será responsável por fazer o
fluido menos denso emergir sobre o fluido mais denso. Os autores de [26]
propõem uma nova força chamada de força de empuxo artificial. Esta força é
necessária para a correta simulação do empuxo. O método SPH, por sua natureza,
somente obtém resultados satisfatórios se o domínio for bem amostrado. Quando
Interface entre os fluidos
PUC-Rio - Certificação Digital Nº 0511019/CA
SPH em Animação de Fluido 33
dois fluidos com massas específicas de repouso diferentes estão interagindo, é
possível que somente algumas poucas partículas do fluido menos denso estejam
dentro do outro fluido (mais denso). Ao se calcular a massa específica dessas
partículas isoladas, o valor da massa específica será mais alto do que deveria ser,
pois no cálculo da massa específica (Equação 4.7), a contribuição será muito
maior das partículas do outro fluido. Conseqüentemente essas partículas isoladas
ou emergirão de forma lenta ou nem mesmo conseguirão emergir do fluido mais
denso. Para contornar este problema de amostragem, os autores propõem a força
artificial de empuxo dada por:
gbf
empuxo
r
)(
0
ρρ
−= (4.19)
onde b é uma constante de controle dada como parâmetro e g
r
é a
gravidade.
4.2.
Funções de Suavização
Neste trabalho utilizamos as mesmas funções de suavização propostas por
Müller et al. [24]. Todas as funções de suavização propostas por Müller et al. [24]
obedecem às condições levantadas no Capítulo 3, Seção 3.3.
A primeira função de suavização proposta por Müller et al. [24] é:
≤≤
−
=
′
−
contráriocaso
hr
rh
h
hxxW
poly
0
0
)(
64
315
),(
322
9
6
π
rr
(4.20)
Figura 4 Gráfico da função de suavização W
poly6
e suas derivadas usada por Müller et al
[24]. O eixo das ordenadas representa o valor da função de suavização e o eixo das
abscissas representa a distância entre os pontos. O gráfico mais à esquerda representa
a função de suavização. O gráfico do meio representa o gradiente de W
poly6
e o mais à
direita representa o Laplaciano de W
poly6
.
PUC-Rio - Certificação Digital Nº 0511019/CA
SPH em Animação de Fluido 34
A escolha de Müller et al. [24] por esta função de suavização se deve ao
fato de somente
2
r
(quadrado da distância) aparecer na equação, eliminando a
necessidade do cálculo de raiz quadrada presente no cálculo da distância entre os
pontos. Isto faz com que a função de suavização seja avaliada eficientemente.
Müller et al.[24] utilizam a função de suavização W
poly6
para o cálculo da massa
específica das partículas e da tensão superficial.
Para o cálculo da força de pressão )( P∇−
r
, se fosse utilizada a função
W
poly6
, iria ocorrer o efeito de aglomeração das partículas do fluido, especialmente
em situações de alta pressão. Esse efeito se daria pelo fato de que à medida que as
partículas se aproximassem, o gradiente de W
poly6
se aproximaria de zero (Figura
4) e, conseqüentemente, a força de repulsão também se aproximaria de zero.
Müller et al.[24] propuseram uma segunda função de suavização, para contornar o
problema de aglomeração das partículas. A segunda função de suavização
proposta por Müller et al.[24] foi proposta originalmente por Desbrun [5] e é
usada somente no cálculo da força de pressão.
≤≤
−
=
contráriocaso
hr
rh
h
W
spike
0
0
)(
15
3
6
π
(4.21)
Figura 5 Gráfico da função de suavização W
Spike
e suas derivadas usada por Müller et al
[24]. O eixo das ordenadas representa o valor da função de suavização e o eixo das
abscissas representa a distância entre os pontos.O gráfico mais à esquerda representa a
função de suavização.O gráfico do meio representa o gradiente de W
Spike
e o mais à
direita representa o Laplaciano de W
Spike
.
Ao se usar a Equação 4.21, a força de repulsão entre as partículas do fluido
aumenta de acordo com o inverso da distância entre as partículas, não permitindo
o efeito de aglomeração (Figura 5).
A força de viscosidade surge devido à fricção entre as moléculas do fluido
causando a transformação de energia cinética do fluido em calor. Se a Equação
PUC-Rio - Certificação Digital Nº 0511019/CA
SPH em Animação de Fluido 35
4.20 ou a Equação 4.21 fossem usadas no cálculo da força de viscosidade, em
certos pontos, o Laplaciano das funções seria negativo (Figura 4 e Figura 5),
fazendo com que a força de viscosidade acrescentasse energia ao fluido, ao invés
de retirar a energia do fluido. Por causa disso, Müller et al.[24] propuseram uma
terceira função de suavização, usada no cálculo da força de viscosidade do fluido.
≤≤
−++
−
=
contráriocaso
hr
r
h
h
r
h
r
h
W
ityvis
0
0
1
2
2
2
15
2
2
3
3
3
cos
π
( 4.22 )
Figura 6 Gráfico da função de suavização W
Viscosity
e suas derivadas usada por Müller et
al [24]. O eixo das ordenadas representa o valor da função de suavização e o eixo das
abscissas representa a distância entre os pontos.O gráfico mais à esquerda representa a
função de suavização.O gráfico do meio representa o gradiente de W
Viscosity
e o mais à
direita representa o Laplaciano de W
Viscosity
.
Esta função de suavização possui o Laplaciano positivo em todo o
domínio (Figura 6), garantindo que a força de viscosidade gerada seja sempre
dissipativa.
PUC-Rio - Certificação Digital Nº 0511019/CA
5
Implementação
Neste capítulo iremos descrever os aspectos de maior relevância na criação
e implementação da biblioteca desenvolvida para simulação de fluidos. A
biblioteca implementa as idéias apresentadas por Müller et al. [24, 26].
5.1.
Estruturação da biblioteca
A biblioteca foi desenvolvida utilizando a linguagem de programação C++
e está dividia em 3 módulos (Figura 7).
Figura 7 Estruturação da biblioteca desenvolvida.
O módulo de Simulação Física é responsável por simular e animar as
entidades de um mundo físico. As entidades do mundo físico podem ser uma
única partícula ou um sistema de partículas (várias partículas ). A biblioteca,
através de um método numérico (Euler, Leap Frog [7]), calcula a nova posição da
entidade no tempo tt
∆
+
0
devido à ação dos agentes (forças aplicadas) sobre as
entidades.
Após o cálculo das novas posições das entidades, verifica-se a validade
dessas novas posições (ausência de colisão com obstáculos). Caso a entidade
esteja numa posição inválida, traz-se a entidade para uma posição válida.
O módulo de Grid Espacial é responsável por determinar as vizinhanças
das partículas. A simulação de fluidos descrita neste trabalho é baseada em
Simulação Grid Espacial
Visualização
PUC-Rio - Certificação Digital Nº 0511019/CA
Implementação 37
sistemas de partículas e essas partículas interagem umas com as outras. A maneira
mais trivial de descobrir quais partículas estão interagindo, é fazer uma busca
quadrática (O(n
2
)) entre todas as partículas dentro do sistema. Essa solução não é
adequada uma vez que estamos buscando simulações a taxas interativas. Com o
método SPH, funções de suavização possuem um domínio de influência finito, ou
seja, apenas uma fração das partículas estarão interagindo com uma dada
partícula.
Uma maneira de evitar a busca quadrática entre as partículas é aplicar um
método de subdivisão espacial capaz de diminuir a complexidade do problema de
busca. Neste trabalho optou-se pelo o uso da estrutura espacial de Grade Regular
(grid) [9]. Outras estruturas espaciais existem e também são usadas para resolver
o problema de busca [9].
O uso do grid para resolver o problema de busca diminui a complexidade
do problema. Com o uso do grid a ordem do problema deixa de ser quadrática
O(n
2
) e passa a ser O(mn) onde m é o número médio de partículas dentro das
células vizinhas à célula que contém a partícula de interesse. Um dos fatores que
determina a eficiência da estrutura de grid é o tamanho das células. Para o
problema de busca na simulação de fluidos, o tamanho para as células é escolhido
de forma que seja igual ao tamanho do raio de influência das funções de
suavização (assumindo que a área de influência de todas as funções de suavização
seja a mesma).
Com a subdivisão em células, dado que a partícula de interesse encontra-se
na célula {i,j,k} (estamos considerando o caso 3D), para saber quais partículas
estão interagindo com a partícula de interesse, basta testar as partículas presentes
nas células adjacentes nas 3 dimensões. No caso de um grid 2D, temos 9 células
para fazer a busca e no caso 3D temos 27 células. A Figura 8 exemplifica a busca
no caso de um grid em 2D.
Figura 8 A célula cinza clara possui a partícula de interesse. A busca para determinar
quais partículas estarão interagindo com a partícula de interesse será feita na própria
célula e nas células em cinza escuro.
PUC-Rio - Certificação Digital Nº 0511019/CA
Implementação 38
O módulo de visualização é responsável pela renderização dos resultados da
simulação. Este módulo utiliza a API gráfica OpenGL e uma estrutura de grafo de
cena [6]. Duas formas para visualização dos fluidos foram desenvolvidas.
Podemos visualizar as partículas do fluido como esferas ou pontos. Outra forma
de visualizar a simulação de fluidos é através do algoritmo Marching Cubes [19],
não implementado neste trabalho. O algoritmo Marching Cubes extrai uma iso
superfície do fluido e a desenha como uma malha de triângulos.
5.2.
Evolução do Sistema
A simulação de fluidos implementada pode ser descrita em 5 passos.
Destacaremos os pontos mais importantes de cada passo. A Tabela 1 enumera os
atributos associados às partículas do fluido.
Posição )(x
r
Velocidade )(v
r
Massa )(m
Massa específica
)(
ρ
Pressão )(P
Cs
Normal
)(n
r
Lista de partículas vizinhas
Coeficiente de Viscosidade )(
µ
Massa específica de Repouso )(
0
ρ
Constante de Gás )(k
Tabela 1 Os principais atributos das partículas de um fluido presentes neste trabalho. Os
últimos 3 atributos são necessários para o caso de interação entre fluidos (Seção 4.1.6).
1º passo: Preenchimento da estrutura de grid.
No início de cada quadro, a estrutura de grid é limpa e preenchida com as
novas posições das partículas. Após o preenchimento do grid, a lista de
vizinhança de cada partícula é atualizada.
PUC-Rio - Certificação Digital Nº 0511019/CA
Implementação 39
2º passo: Cálculo da Massa específica.
A massa específica das partículas é essencial na formulação do método
SPH. O cálculo da massa específica das partículas é feito após o
preenchimento da lista de vizinhança das partículas. Para o cálculo da
massa específica é usada a Equação 4.7. Após o cálculo da massa
específica, é calculada a pressão de cada partícula através da Equação
4.11.
3º passo: Cálculo do valor de Cs e da normal de superfície
Neste passo é calculado e armazenado o valor do campo de cor da
partícula (Cs), assim como a normal da superfície (gradiente de Cs) e o
divergente de Cs ( Laplaciano de Cs). Esses valores são armazenados para
mais adiante serem utilizados no cálculo da tensão superficial do fluido.
Além disso, com o valor da normal de superfície de cada partícula, é
possível identificar as partículas que compõem a superfície do fluido. A
identificação das partículas que compõem a superfície do fluido é útil
tanto para o efeito de visualização (em casos que haja a extração da iso
superfície) quanto para a colisão contra obstáculos.
4º passo: Aplicação das Forças
Aplicam-se todos os agentes (gradiente de pressão, força de viscosidade,
forças externas) presentes na Equação de Navier-Stokes sobre as partículas
do fluido, guardando a força total exercida sobre cada partícula do fluido.
5º passo: Evolução
Com a força resultante em cada partícula calculada no passo anterior,
obtemos as acelerações (Equação 4.5) e conseqüentemente, as novas
posições das partículas. Para fazer a evolução do sistema o método
numérico Leap Frog foi utilizado. Testes de colisão contra obstáculos são
executados e são aplicadas as correções para evitar que as partículas
estejam em colisão contra os obstáculos.
PUC-Rio - Certificação Digital Nº 0511019/CA
Implementação 40
Após o cálculo das novas posições das partículas, as novas posições são
enviadas para o módulo responsável pela visualização.
5.3.
Discussão
Neste trabalho a implementação da tensão superficial apresentada por
Müller et al.[24] apresentou problema. Na teoria, o valor do vetor normal à
superfície do fluido || n
r
é diferente de zero apenas na superfície do fluido. Por
isso, os autores sugerem que somente deve ser calculada e aplicada a tensão
superficial às partículas que tiverem || n
r
maior que um limiar, já que o cálculo de
n
n
r
r
onde || n
r
é próximo de zero, pode causar problemas de precisão numérica. O
resultado obtido na implementação da biblioteca, demonstrou que o valor deste
limiar não pode ser muito próximo de zero.
Foi feito um teste em que um fluido representando água, com 8000
partículas e massa específica de repouso de 1000 kg/m
3
encheria uma caixa em
forma de cubo. As dimensões da caixa vão de -0.12 a 0.12 nos três eixos X,Y e Z.
Foi esperado que o fluido chegasse ao repouso para amostrar alguns valores das
normas dos vetores normais à superfície e alguns valores de Cs. Os valores foram
amostrados ao longo de uma linha reta traçada cruzando o meio da caixa. A
Figura 9 ilustra o teste feito.
Figura 9 Uma caixa enchida com um fluido representado por 8000 partículas. A linha reta
(em branco) foi utilizada para a amostragem dos valores da norma do vetor normal à
superfície e de Cs.
PUC-Rio - Certificação Digital Nº 0511019/CA
Implementação 41
-1
1
3
5
7
9
11
13
15
17
19
21
23
25
27
29
31
33
35
37
39
41
43
45
-0.12 -0.095 -0.07 -0.045 -0.02 0.005 0.03 0.055 0.08 0.105
Figura 10 O gráfico com os valores da norma do vetor normal à superfície amostrados
em diferentes pontos do fluido. O eixo das abscissas do gráfico representa o eixo X.
0
0.5
1
1.5
-0.12 -0.095 -0.07 -0.045 -0.02 0.005 0.03 0.055 0.08 0.105
Figura 11 O gráfico com os valores do campo de cor (Cs) amostrados em diferentes
pontos do fluido. O eixo das abscissas do gráfico representa o eixo X.
A Figura 10 apresenta o gráfico com os valores das normas dos vetores
normais à superfície das amostras feitas em diferentes posições do fluido sobre a
linha traçada. A Figura 11 apresenta os valores do campo de cor (Cs) das amostras
feitas em diferentes posições do fluido sobre a linha traçada. Através do gráfico da
Figura 10, vemos que no interior do fluido a norma dos vetores normal à
superfície é diferente de zero e que através do gráfico da Figura 11 os valores do
PUC-Rio - Certificação Digital Nº 0511019/CA
Implementação 42
campo de cor oscilam próximos de 1. Esse fato ocorre devido a problemas de
amostragem na aplicação do método SPH. Se fosse possível conseguir uma
enorme quantidade de amostras ao se determinar o campo de cor das partículas, os
valores amostrados seriam muito próximos do valor esperado 1 e
conseqüentemente a norma do vetor normal à superfície do fluido dado pelo
gradiente do campo de cor seria bem próxima de zero.
Na prática, uma quantidade muito grande de amostra compromete a
eficiência do método SPH em relação ao seu custo computacional, tornando
inviável em certos tipos de aplicações. Na biblioteca desenvolvida usamos como
limiar de || n
r
para determinar se uma partícula pertence à superfície do fluido o
valor 10.0. Este valor foi obtido através de alguns experimentos mudando os
parâmetros de configuração de amostragem, como quantidade de partículas totais
e raio de suporte das funções de suavização.
PUC-Rio - Certificação Digital Nº 0511019/CA
6
Resultados
Neste capítulo serão descritos alguns resultados obtidos com a biblioteca
desenvolvida. Algumas cenas foram montadas para demonstrar a capacidade da
biblioteca em animar diversas configurações de fluidos. Em todos os testes foi
utilizado um computador Core 2 Duo E6600 com 3GB de RAM e com uma placa
de vídeo GeForce 7800 GTX 256 MB.
6.1.
Propriedades de Fluidos
Nesta seção de teste iremos demonstrar como os parâmetros do fluido como
volume, viscosidade e tensão superficial afetam o comportamento do fluido.
6.1.1.
Volume
O primeiro parâmetro a ser estudado será o volume do fluido. A escolha do
volume do fluido é dado como parâmetro de entrada para simulação. Neste
experimento será demonstrado como a escolha do volume afeta a simulação do
fluido. O experimento consiste em encher uma caixa com volume total de
aproximadamente 55 litros. Foram usados dois fluidos idênticos, porém com
volumes diferentes. O primeiro fluido usado para encher a caixa possui 8000
partículas, massa específica de repouso de 1000 kg/m
3
e volume de 50 litros. O
segundo fluido usado para encher a caixa possui 8000 partículas, massa específica
de repouso de 1000 kg/m
3
e volume de 10 litros. O passo de simulação foi de
1/300 segundo e o raio de influência da função de suavização foi de 0.025 metros
em ambos os casos. A Figura 12 ilustra o resultado da comparação entre os dois
fluidos.
O tempo gasto por passo de simulação com o fluido de maior volume foi em
média 0.041 segundos, enquanto o tempo gasto por passo de simulação com o
fluido de menor volume foi em média 0.101 segundos. Essa diferença de tempo
PUC-Rio - Certificação Digital Nº 0511019/CA
Resultados 44
gasto por passo de simulação é conseqüência do aumento do número de vizinhos
interagindo por partícula. Um volume pequeno com uma grande quantidade de
partículas compromete a eficiência da busca de vizinhos das partículas, uma vez
que teremos muitas partículas por célula do grid. Podemos contornar esse
problema de duas maneiras. A primeira é diminuindo a quantidade de partículas
que representam o fluido. A segunda é diminuindo o raio de influência da função
de suavização e conseqüentemente o tamanho das células do grid, a fim de
diminuir o número de vizinhos por partícula.
Figura 12 Comparação entre o volume dos fluidos. O fluido da imagem à esquerda
possui volume de 50 litros e o fluido da imagem à direita possui volume de 10 litros.
6.1.2.
Viscosidade
A viscosidade do fluido é uma propriedade dissipativa de energia. Quanto
mais viscoso o fluido, maior a resistência à mudança de movimento. Neste
experimento, um fluido, na configuração inicial mostrada na Figura 13, foi
inicialmente posicionado a uma certa altura e submetido à ação da gravidade, até
cair dentro de uma caixa. Foi comparado o tempo de simulação gasto por dois
fluido idênticos a não ser pelo valor da constante de viscosidade, para chegar ao
estado de repouso. O fluido utilizado é representado por 6000 partículas, massa
específica de repouso de 1000 kg/m
3
. O passo de simulação utilizado foi de 1/300
segundo e a razão da viscosidade entre os fluidos é de 10.
Figura 13 Configuração inicial do fluido utilizado no experimento.
PUC-Rio - Certificação Digital Nº 0511019/CA
Resultados 45
O tempo gasto de simulação pelo fluido menos viscoso foi de
aproximadamente 7 segundos para chegar ao repouso, enquanto o fluido mais
viscoso levou aproximadamente 3,5 segundos.
A viscosidade interfere também no comportamento do fluido. A Figura 14
mostra a diferença de comportamento entre os dois fluidos. A imagem à esquerda
da Figura 14 ilustra o comportamento do fluido menos viscoso instantes após o
fluido se chocar com a caixa. A imagem à direita da Figura 14 ilustra o fluido
mais viscoso instantes após se chocar com a caixa. Comparando as duas imagens
da Figura 14 é possível notar a diferença de comportamento dos dois fluidos.
Figura 14 Diferença de comportamento entre o fluido menos viscoso (imagem à
esquerda) e o fluido mais viscoso (imagem à direita).
6.1.3.
Tensão Superficial
Por último iremos analisar como a tensão superficial afeta o fluido. A tensão
superficial é responsável por minimizar a curvatura da superfície do fluido. A
tensão superficial é responsável pelas criações de gotas e pelo seus formatos quase
esféricos, uma vez que a esfera é a área de menor superfície.
Neste experimento, foram comparados dois fluidos idênticos a menos da
constante de curvatura usada no cálculo da tensão superficial. Foi criado um
fluido com 1000 partículas e massa específica de repouso de 1000 kg/m
3
. Foi
aplicado a equação de Navier–Stokes ao fluido, porém sem aplicar a gravidade e
nenhuma outra força externa com exceção da tensão superficial. A configuração
do fluido é semelhante à configuração inicial ilustrada na Figura 13. O passo de
simulação foi de 1/300 segundo. O resultado é ilustrado na Figura 15. Na imagem
à esquerda da Figura 15, o valor do coeficiente de tensão superficial é zero, ou
seja, não há a contribuição da tensão superficial. As partículas se atraem devido à
força de pressão mas não formam uma superfície esférica. Na imagem à direita o
PUC-Rio - Certificação Digital Nº 0511019/CA
Resultados 46
valor do coeficiente de tensão superficial é 1. A tensão superficial aplicada ao
fluido, faz com que o fluido tenha o formato esférico.
Figura 15 A imagem à esquerda ilustra um fluido com o valor do coeficiente de tensão
superficial 0 e à direita o mesmo fluido com o valor do coeficiente de tensão superficial 1.
Num segundo experimento foi comparado como a tensão superficial afeta
o comportamento do fluido. Foi utilizado a mesma configuração do experimento
anterior, porém com a força de gravidade atuando nos fluidos. O fluido cai devido
à ação da gravidade até se chocar com uma caixa (não desenhada). A Figura 16
mostra a diferença de comportamento entre os dois fluidos.
Figura 16 A imagem de cima, ilustra um fluido com o valor do coeficiente de tensão
superficial 0 e a imagem de baixo o mesmo fluido com o valor do coeficiente de tensão
superficial 1.
PUC-Rio - Certificação Digital Nº 0511019/CA
Resultados 47
Através da Figura 16 é possível perceber a diferença de comportamento entre os
fluidos devido à tensão superficial. Na Figura 16 a imagem de baixo é possível
notar uma superfície arredondada e suave no fluido que sofre a atuação da tensão
superficial. O mesmo não acontece com o fluido que não sofre a atuação da tensão
superficial.
6.2.
Copo de Água
Neste teste serão reproduzidos os testes apresentados por Müller et al.[24].
O teste consiste em encher um copo de água representado por um obstáculo
cilíndrico. Além disso, o teste permite o usuário interagir com o fluido aplicando
uma força externa controlada via dispositivo de entrada (no caso, o mouse).
A Figura 17 mostra alguns instantes da simulação do fluido enchendo o copo.
Nesta simulação a quantidade de partículas usada foi de 3000. O fluido possui
massa específica de 1000 kg/m
3
. O passo de simulação foi de 1/300 segundo e
tempo gasto por passo foi de 0.012 segundos.
Neste teste foi possível a reprodução dos testes propostos em Müller et
al.[24]. No trabalho de Müller et al.[24], é realizada uma simulação com 3000
partículas a uma taxa de 5 quadros por segundo visualizada com o algoritmo de
Marching Cubes[19]. Neste trabalho a simulação de fluido com 3000 partículas
visualizadas por esferas ocorreu a uma taxa de 44 quadros por segundo.
PUC-Rio - Certificação Digital Nº 0511019/CA
Resultados 48
Figura 17 Três instantes da simulação. As imagens da esquerda mostram os resultados
obtidos no trabalho desenvolvido. As imagens da direita foram retiradas do trabalho de
Müller et al.[24].
PUC-Rio - Certificação Digital Nº 0511019/CA
Resultados 49
6.3.
Interação entre Fluidos
Neste teste serão reproduzidos os resultados obtidos em Müller et al.[26]. O
teste consiste em misturar dois fluidos com massas específicas de repouso
diferentes. Primeiramente enche-se um compartimento com um fluido de massa
específica de repouso de 1000 kg/m
3
, de cor azul. Ao acabar de encher o
recipiente com o fluido azul, começamos a derramar o segundo fluido de cor
vermelha, cuja massa específica de repouso é de 500 kg/m
3
. A Figura 18
demonstra os resultados obtidos, comparando com os resultados obtidos em [26].
Figura 18 Comparação entre os resultados obtidos e os resultados obtidos em [26]. As
imagens da esquerda foram obtidas através da biblioteca desenvolvida, e as imagens da
direita foram retiradas de [26].
PUC-Rio - Certificação Digital Nº 0511019/CA
Resultados 50
Para este teste foi utilizado um passo de simulação de 1/300 segundo e foi
utilizado a força de empuxo artificial (Equação 4.19) com o coeficiente igual a 2.
No trabalho de [26] não são dados muitos detalhes sobre os parâmetros dos
fluidos utilizados nos testes. Por isso, neste teste, foram experimentados diferentes
parâmetros até alcançar um comportamento parecido com o obtido em [26].
6.4.
Quebra de Barragem
Neste experimento foi feita uma comparação visual entre os resultados
obtidos na implementação deste trabalho e os resultados apresentados em [17].
Este experimento foi primeiramente introduzido por Monaghan et al. [22]. O
objetivo deste experimento é comparar o modelo físico desenvolvido neste
trabalho com um modelo físico mais preciso.
O cenário de teste é a simulação da quebra de uma barragem também
conhecido como Dam collapse problem. O cenário consiste em encher uma
barragem e subitamente abrir as comportas, liberando toda a água acumulada. As
figuras 19, 20, 21, ilustram a comparação visual entre os resultados obtidos través
o uso da biblioteca desenvolvida e os resultados apresentados em [17]. Devemos
ressaltar que a simulação presente no experimento apresentado em [17] é em 2D
enquanto neste experimento a simulação obtida é em 3D.
Através das figuras 19, 20 e 21, podemos perceber que o comportamento
do fluido simulado pela biblioteca desenvolvida e do fluido apresentado em [17]
são parecidos com exceção da Figura 21. A diferença de comportamento ilustrada
pela Figura 21 pode ser causada pela forma distinta que as duas implementações
tratam o problema de colisão contra obstáculos. Além disso, em [17] não são
divulgados os parâmetros utilizados (passo de simulação, quantidade de
partículas, etc...) no experimento conduzido.
Neste experimento o fluido modelado pela biblioteca desenvolvida utilizou
8000 partículas, com passo de simulação de 1/300 segundo. O fluido simulado
possuía massa específica de repouso de 1000 kg/m
3
o mesmo da água e o tempo
médio por passo de simulação foi de 0.0194 milisegundos.
Mesmo usando um modelo físico mais simplificado, a biblioteca
desenvolvida foi capaz de animar o experimento de quebra de barragem próximo
PUC-Rio - Certificação Digital Nº 0511019/CA
Resultados 51
da animação feitas por modelos mais sofisticados focados para o uso em
aplicações de engenharia.
Figura 19 Instante inicial da simulação. Um fluido está confinado por uma barragem. A
imagem mais acima foi retirada de [17]. A imagem do meio mostra a visualização
ortográfica da simulação e a imagem mais abaixo mostra a vista perspectiva da
simulação.
PUC-Rio - Certificação Digital Nº 0511019/CA
Resultados 52
Figura 20 Instante após a abertura da barragem em que o fluido encontra-se com o
obstáculo (rampa). A imagem mais acima foi retirada de [17]. A imagem do meio mostra
a visualização ortográfica da simulação e a imagem mais abaixo mostra a vista
perspectiva da simulação.
PUC-Rio - Certificação Digital Nº 0511019/CA
Resultados 53
Figura 21 O fluido interagindo com a rampa (obstáculo) formando uma onda. A imagem
mais acima foi retirada de [17]. A imagem mais acima foi retirada de [17]. A imagem do
meio mostra a visualização ortográfica da simulação e a imagem mais abaixo mostra a
vista perspectiva da simulação.
PUC-Rio - Certificação Digital Nº 0511019/CA
Resultados 54
6.5.
Descarga de Água
Neste experimento será feita uma comparação visual entre os resultados
obtidos pela biblioteca desenvolvida e os resultados obtidos em [17]. O objetivo
deste experimento é comparar o modelo físico apresentado em [17] com o modelo
físico desenvolvido neste trabalho. Note que os resultados apresentados em [17]
são em 2D, enquanto os resultados apresentados neste experimento são em 3D. O
experimento consiste em submeter o fluido a uma situação de alta pressão. Para
isto, uma barragem será enchida assim como o experimento anterior, porém ao
invés de abrir a comporta por inteiro, apenas uma pequena parte será aberta (por
volta de 12 %). Nesta configuração, as partículas de água próximas à abertura são
expelidas para fora da barragem devido à forte força de pressão exercida pelo
resto do fluido. As figuras abaixo ilustram a comparação dos resultados obtidos
neste trabalho e os resultados obtidos em [17].
Figura 22 Instante que a comporta é aberta para vazão da água. A imagem da esquerda
foi retirada de [17], a imagem do meio e a imagem da direita são a vista ortográfica e
perspectiva, respectivamente, da simulação obtida neste trabalho.
PUC-Rio - Certificação Digital Nº 0511019/CA
Resultados 55
Figura 23 Neste instante devido à forte força de pressão as partículas de água são
expelidas pelo compartimento aberto formando uma cavidade. A imagem da esquerda foi
retirada de [17], a imagem do meio e a imagem da direita são a vista ortográfica e
perspectiva, respectivamente, da simulação obtida neste trabalho.
Figura 24 Diferença entre os resultados obtidos em [17] e os resultados obtidos neste
trabalho. A imagem da esquerda foi retirada de [17], a imagem do meio e a imagem da
direita são a vista ortográfica e perspectiva, respectivamente, da simulação obtida neste
trabalho.
Analisando as Figuras 22, 23 e 24 vemos que a simulação obtida pela
biblioteca se aproximou da simulação obtida em [17]. A simulação física obtida
pela biblioteca desenvolvida conseguiu simular o efeito de alta pressão sobre as
partículas expelidas. Analisando a imagem da esquerda da Figura 24 as partículas
expelidas pela comporta sobem e depois caem devido à ação da gravidade,
formando uma cavidade abaixo do fluxo expelido. Através das vistas ortogonal e
perspectiva da simulação (imagem central e mais à direita da Figura 24
respectivamente) obtidas neste experimento, podemos perceber tal cavidade.
Para este experimento o fluido foi modelado por 10648 partículas, com
massa específica de 1000 kg/m
3
. O passo usado na simulação foi de 1/500
segundo e o tempo médio gasto por passo de simulação foi de 0.1186 segundo.
PUC-Rio - Certificação Digital Nº 0511019/CA
Resultados 56
6.6.
Custo Computacional
Nesta seção foram levantados alguns pontos relacionados ao custo
computacional da simulação de fluidos e do modelo físico adotado. Será abordada
a relação entre quantidade de partículas e tempo de simulação por passo. Outro
aspecto importante é a determinação do gargalo da simulação no modelo físico
desenvolvido neste trabalho.
O primeiro aspecto a ser discutido é a relação entre quantidade de partículas
e o tempo de simulação por passo. Para isto, foi enchida repetidamente uma caixa
com um fluido representado por diferentes quantidades de partículas, deixando a
simulação correr por 2000 passos tomando o tempo gasto por cada passo. A
Figura 25 e a Tabela 2 mostram a relação entre quantidade de partículas e tempo
médio de simulação por passo. No experimento todas as propriedades do fluido
permaneceram constantes, exceto o volume do fluido que varia para manter
aproximadamente uma quantidade de 50 vizinhos em média por partícula.
Número de
Partículas
Tempo de
Simulação (seg)
1000 0.005936
2000 0.013446
3000 0.020412
4000 0.028696
5000 0.037857
6000 0.042766
7000 0.04922
8000 0.058672
9000 0.064943
10000 0.072568
Tabela 2 Valores obtidos na comparação entre o número de partículas e o tempo médio
de simulação por passo.
PUC-Rio - Certificação Digital Nº 0511019/CA
Resultados 57
Numero de particulas x Tempo de Simulacao
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0 2000 4000 6000 8000 10000 12000
Numero de Particulas
Tempo de Simulacao
Series1
Figura 25 Gráfico dos dados apresentados na Tabela 2.
Através do gráfico ilustrado na Figura 25, podemos notar um
comportamento linear entre o número de partículas e o tempo de simulação gasto
por passo. Esse comportamento já era esperado, devido ao uso do grid para busca
da vizinhança das partículas, uma vez que com o uso do grid, a complexidade da
busca passa a ser ordem O(nm), onde m é a quantidade média de partículas nas
células vizinhas à partícula de interesse.
O segundo aspecto a ser discutido é determinar o gargalo da simulação de
fluidos. Para determinar o gargalo da simulação foi utilizado o programa VTune
da Intel. O VTune é um programa utilizado para analisar o desempenho de
aplicativos.
O teste feito se baseia em encontrar o gargalo da simulação com diferentes
quantidades de partículas. Para se determinar o gargalo da simulação foi utilizado
o mesmo cenário do teste anterior.
PUC-Rio - Certificação Digital Nº 0511019/CA
Resultados 58
0 5 10 15 20 25 30 35
Calculate_Particles_Surface
Apply_Agent_In_Particle
Extract_Neighbors_Threshold
Calculate_Particles_Density
Função
Tempo total gasto (%)
Figura 26 Gargalo da simulação com 1000 partículas.
Para a simulação com a 1000 partículas, o gargalo encontra-se na
identificação das partículas de superfície do fluido (Figura 26). O tempo gasto
pelo método responsável pela identificação da superfície do fluido consome
aproximadamente 33% do tempo total da simulação. O segundo método que mais
consome tempo de simulação (cerca de 30%) é o método responsável por aplicar
os agentes às partículas.
Aumentando a quantidade de partículas para 4000, o gargalo da simulação
muda. Analisando a Figura 27, temos que o gargalo passa a ser o método
responsável por determinar os vizinhos de todas as partículas com
aproximadamente 28% do tempo total da simulação. O método de identificação da
superfície do fluido passa a ocupar aproximadamente 27% do tempo total.
PUC-Rio - Certificação Digital Nº 0511019/CA
Resultados 59
0 5 10 15 20 25 30 35
Extract_Neighbors_Threshold
Calculate_Particles_Surface
Apply_Agent_In_Particle
Calculate_Particles_Density
Funçaõ
Tempo total gasto (%)
Figura 27 Identificação do gargalo com 4000 partículas.
Por último foi analisado o gargalo com um fluido representado por 8000
partículas. Com 8000 partículas o gargalo novamente fica com a determinação dos
vizinhos das partículas (Figura 28). O método passou a gastar aproximadamente
35% do tempo total da simulação, enquanto o segundo método que mais gasta
tempo de simulação (o método de identificação da superfície do fluido) gasta
aproximadamente 26% do tempo total da simulação.
0 10 20 30 40
Extract_Neighbors_Threshold
Calculate_Particles_Surface
Apply_Agent_In_Particle
Calculate_Particles_Density
Função
Tempo total gasto (%)
Figura 28 Identificação do gargalo com 8000 partículas.
PUC-Rio - Certificação Digital Nº 0511019/CA
Resultados 60
Através do gráfico ilustrado na Figura 25, fica evidente o comportamento
linear entre a quantidade de partículas e o tempo gasto por passo de simulação. O
gargalo da simulação depende do número de partículas usado para representar o
fluido. Ao se usar poucas partículas (menos que 4000) o gargalo da simulação fica
sendo o método responsável por identificar as partículas que compõem a
superfície do fluido seguido pelo método responsável por aplicar os agentes às
partículas. Ao se aumentar o número de partículas que representam o fluido, o
gargalo passa a ser o método responsável por determinar os vizinhos de cada
partícula seguido pelo método responsável por identificar as partículas que
compõem a superfície do fluido.
PUC-Rio - Certificação Digital Nº 0511019/CA
7
Conclusão
O objetivo dessa dissertação foi fazer um estudo compreensivo sobre
simulação de fluidos baseada em partículas utilizando o método conhecido como
Smoothed Particle Hydrodynamics (SPH). Além disso, buscamos a validação do
método SPH para simulação de fluidos através da biblioteca desenvolvida. Para a
simulação de fluidos para efeito de animação em tempo interativo, seguimos as
idéias propostas em Müller et al.[24, 26].
Através dos testes desenvolvidos, fomos capazes de criar cenários para a
validação do nosso trabalho. No Capítulo 6, Seção 2, reproduzimos os resultados
obtidos em Müller et al.[24]. Apesar de não haver a possibilidade de uma
comparação de eficiência entre o nosso trabalho e o trabalho apresentado em [24],
a comparação visual da animação se demonstrou parecida.
Nos testes apresentados no Capítulo 6, Seções 4 e 5, comparamos o modelo
físico mais simplificado desenvolvido neste trabalho, contra um modelo físico
mais sofisticado apresentado em [17]. Em ambos os modelos é usado o método
SPH. As diferenças do modelo físico mais sofisticado em relação ao mais
simplificado utilizado nesta dissertação são: cálculo da massa específica, funções
de suavização, cálculo da força de viscosidade e tratamento de fronteira.
Apesar de todas essas diferenças entre os modelos físicos, fomos capazes de
reproduzir alguns experimentos descritos em [17]. No experimento “Quebra de
Barragem”(Dam Collapse) o comportamento obtido com o fluido na animação foi
bem próximo ao comportamento apresentado em [17], e no experimento
“Descarga De Água” (Water Discharge) o modelo físico simplificado também
conseguiu um resultado satisfatório, conseguindo simular a cavidade que surge
devido à força de alta pressão.
Incorporamos também neste trabalho a interação entre fluidos baseada nas
idéias apresentadas em Müller et al.[26]. Apesar de não terem sido implementadas
todas as propostas apresentadas em [26], nossa biblioteca é capaz de interagir dois
PUC-Rio - Certificação Digital Nº 0511019/CA
Conclusão 62
fluidos com massas específicas diferentes, de forma que o fluido menos denso
permaneça acima do fluido mais denso.
O estudo feito para determinar o custo computacional e o gargalo da
simulação demonstrou que o maior problema do método SPH para simulação de
fluidos (representados por mais de 4000 partículas) neste trabalho está na busca da
vizinhança das partículas. Não foi implementado nenhum outro algoritmo para
busca de vizinhança para efeito de comparação.
Em relação a trabalhos futuros, podemos destacar a implementação de todas
as idéias apresentadas em [26] para interação entre fluidos, permitindo a criação
de partículas de ar, interação entre ar e água, e interação entre fluidos permitindo
mudanças de fases.
Para conseguir simulações estáveis, o passo de simulação usado ainda é
muito baixo. Para alcançarmos simulações em tempo real, é necessária a
investigação de métodos capazes de permitir passos de simulação maiores.
Algumas propostas, como raio de suavização variável por partícula [17],
extensões sobre o método SPH [29, 16] e tratamento de fronteira [17], são
apresentadas na literatura, permitindo simulações mais estáveis.
A implementação do método SPH para simulação de fluidos em GPU é um
tópico a ser explorado. Com o aumento do poder de processamento dos atuais
processadores gráficos e com a possibilidade de utilizá-los para processamentos
genéricos, acreditamos obter um aumento de performance e também um aumento
na quantidade de partículas que representam o fluido, se utilizarmos a GPU para a
simulação de fluidos. Alguns trabalhos [13, 33] já utilizam a GPU para simulação
de fluidos utilizando modelos Eulerianos, obtendo resultados bastantes
satisfatórios.
Além disso, a interação entre fluido e corpos rígidos também é um tópico
importante e muito explorado na literatura. Alguns trabalhos [25, 35] estendem o
trabalho de Müller et al.[24] para permitir a interação entre fluidos e corpos
rígidos. Por último, o uso de outras estruturas de subdivisão espacial para tentar
melhorar a performance e diminuir o custo computacional da simulação e a
investigação de métodos de visualização de fluidos para a extração de uma malha
triangular [19], podem ser investigados.
PUC-Rio - Certificação Digital Nº 0511019/CA
8
Bibliografia
[1] Paiva A., Petronetto F., Lewiner T., Tavares G.. Particle-based non-Newtonian
fluid animation for melting objects, 2006.
[2] Carlson M., Mucha P. J., III R. Brooks Van Horn, and Greg Turk. Melting and
flowing. In Proceedings of the ACM SIGGRAPH symposium on Computer
animation, p. 167–174, 2002.
[3] Carlson M., Mucha P. J., Turk G.: Rigid fluid: animating the interplay
between rigid bodies and fluid. ACM Trans. Graph. 23, 3, p. 377–384, 2004.
[4] Clavet S., Beaudoin P., Poulin P. Particle based Viscoelastic Fluid Simulation.
Eurographics ACM SIGGRAPH Symposium on Computer Animation, 2005.
[5] Desbrun M., Gascuel M. P. . Smoothed particles: A new paradigm for
animating highly deformable bodies. In Computer Animation and Simulation´96
(Proceedings of EG Workshop on Animation and Simulation), p. 61-76, 1996.
[6] Eberly H. D. . 3D Game Engine Architecture: Engineering Real Time
Applications with Wild Magic. Morgan Kaufmann, 2005.
[7] Eberly H. D. . Game Physics. Morgan Kaufmann, 2004.
[8] Enright D., Marschner ,S., Fedkiw R.. Animation and rendering of complex
water surfaces. In Proceedings of the 29th annual conference on Computer
graphics and interactive techniques, ACM Press, p. 736–744, 2002.
[9] Ericson C. Real Time Collision Detection. Morgan Kaufmann, 2005.
[10] Foster N., Metaxas D. . Realistic animation of liquids. Graphical Models and
Image Processing 58,5, p. 471-483, 1996.
[11] Foster N., Metaxas D. .Modeling the motion of hot, turbulent gas. In
SIGGRAPH, p.181-188, 1997.
[12] Gingold R. A., Monaghan J. J.. Smoothed particle hydrodynamics: theory
and application to non-spherical starts. Monthly Notices of the Royal
Astronomical Society, p. 182:375-398, 1977.
[13] Harris J. M. . Fast Fluid Dynamics Simulation on the GPU. GPU Gems.
Addison Wesley, p. 637-665, 2004.
[14] Kass M., Miller G. . Rapid stable fluid dynamics for computer graphics.
Computer Graphics Proceedings of SIGGRAPH, p. 49-57, 1990.
[15] Keiser R., Adams B., Gasser D., Bazzi P., Dutré P., Gross M. . A Unified
Lagrangian Approach to Solid-Fluid Animation, 2005.
[16] Koshizuka S., Oka Y. . Moving-particle semi implicit method for
fragmentation of incompressible fluid. Nuclear Science Engineering 123, p. 421-
434, 1996.
PUC-Rio - Certificação Digital Nº 0511019/CA
Bibliografia 64
[17] Liu G.R., Liu M.B, Smoothed Particle Hydrodynamics, a meshfree particle
method. World Scientific Publishing , 2003.
[18] Lombardo J. C., Puech C. .. Oriented particles: A tool for shape memory
objects modelling. In Graphics Interface, p. 255–262, 1995.
[19] Lorensen W. E., Harvey E. C. . Marching cubes: A high resolution 3d surface
construction algorithm. In Proceedings of the 14th annual conference on
Computer graphics and interactive techniques, p. 163–169, 1987.
[20] Lucy L. B. . A numerical approach to the testing of the fission hypothesis.
The Astronomical Journal, 82: p. 1013-1024, 1977.
[21] Miller G. , Pearce A. . Globular dynamics: A connected particle system for
animating viscous fluids. Computer and Graphics 13, 3, p. 305-309, 1989.
[22] Monaghan J. J. . Simulating free surface flow with SPH. Journal of
Computational Physics, 110: p. 399-406, 1994.
[23] Morris J. P. Simulating surface tension with smoothed particle
hydrodynamics. International Journal for Numerical Methods in Fluids, 33(3):
p. 333–353, 2000.
[24] Müller M, Charypar D., Gross M, Particle based fluid simulation for
interactive applications. Proceedings of 2003 ACM SIGGRAPH Symposium on
Computer Animation, p. 154-159, 2003.
[25] Müller M., Schirm S., Teschner M., Heidelberger B., Gross M. . Interaction
of fluids with deformable solids. Journal of Computer Animation and Virtual
Worlds (CAVW) 15, 3-4, p. 159-171, 2004.
[26] Muller M, Solenthaler B, Keiser R, Gross M, Particle based Fluid-Fluid
Interaction. Eurographics ACM Synposium on Computer Animation, 2005.
[27] Nixon D., Lobb R.. A fluid-based soft-object model. IEEE Computer
Graphics and Applications, p. 68–75, 2002.
[28] Neto A., Rodrigues P. Sergio, Giraldi G., Apolinário A. . Animação
Computacional de Fluidos via Smoothed Particle Hydrodynamics, 2005.
[29] Premože S., Tasdizen T., Bigler J., Lefohn A., Whitaker R. T. . Particle-
based simulation of fluids. Computer Graphics Forum 22, 3, p. 401-410, 2003.
[30] Pozrikidis C. . Numerical Computation in Science and Engineering. Oxford
Univ. Press, NY, 1998.
[31] Reeves W. T. .Particle systems — a technique for modeling a class of fuzzy
objects. ACM Transactions on Graphics 2(2), p. 91–108, 1983.
[32] Stam J. . Stable fluids In SIGGRAPH, p. 121-128, 1999.
[33] Scheidegger C. E. , Comba J. L. , Da Cunha R. D. . Navier Stokes on
programmable graphics hardware using smac. In Proceedings of Sibgrapi, 2004.
[34] Takahashi T., Heihachi U., Kunimatsu A., Fujii H. . The simulation of fluid-
rigid body interaction. ACM Siggraph Sketches & Applications, 2002.
[35] Takashi A. Real-Time Particle-Based Fluid Simulation with Rigid Body
Interaction. Game Programming Gems 6, Charles River Media, p. 189-205, 2006.
PUC-Rio - Certificação Digital Nº 0511019/CA
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
