Embed Size (px)
Citation preview
Tiago de Holanda Cunha Nobrega
Sistema de Partículas para Modelagem e SimulaçãoInterativa de Fluidos
Florianópolis, Santa Catarina
Fevereiro de 2007
Tiago de Holanda Cunha Nobrega
Sistema de Partículas para Modelagem e SimulaçãoInterativa de Fluidos
Trabalho de conclusão de curso apresentadocomo parte dos requisitos para obtenção do graude Bacharel em Ciências da Computação
Orientador:
Prof. Dr. rer.nat. Aldo von Wangenheim
Co-orientador:
Thiago Ramos dos Santos
BACHARELADO EM CIÊNCIAS DA COMPUTAÇÃO
DEPARTAMENTO DE INFORMÁTICA E ESTATÍSTICA
CENTRO TECNOLÓGICO
UNIVERSIDADE FEDERAL DE SANTA CATARINA
Florianópolis, Santa Catarina
Fevereiro de 2007
Dedico esta monografia à
minha família, por seu
amor incondicional, e
aos meus amigos,
tanto do Rio quanto de Floripa,
tanto antigos quanto novos,
pelo apoio na jornada.
Resumo
Este trabalho discute e implementa um método para a simulação de fluidos através de sis-temas de partículas. Vários métodos existentes na literatura são apresentados. Tópicos alheiosao método, como detecções de colisão e busca rápida de vizinhança, são abordados e imple-mentados. Por fim, os resultados são apresentados e diferentes implementações comparadas ediscutidas.
Sumário
Lista de Figuras
Lista de Tabelas
1 Introdução p. 12
1.1 Objetivo Geral . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 13
1.2 Objetivos Específicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 14
2 Notação e Conceitos Matemáticos p. 15
3 Dinâmica dos Fluidos Computacional p. 17
3.1 Métodos baseados em Grids para Representação de Fluidos . . . . . . . . . . p. 17
3.1.1 Método Semi-Lagrangeano . . . . . . . . . . . . . . . . . . . . . . . p. 18
3.1.2 Marker-and-Cell (MAC) Grids . . . . . . . . . . . . . . . . . . . . . p. 18
3.1.3 Aproximações HeightField . . . . . . . . . . . . . . . . . . . . . . . p. 20
3.2 Sistemas de Partículas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 21
3.2.1 Smoothed Particle Hydrodynamics . . . . . . . . . . . . . . . . . p. 22
3.2.2 Método Moving Particle Semi-Implicit (MPS) . . . . . . . . . . . . . p. 24
3.3 Outras Metodologias de Simulações de Fluidos . . . . . . . . . . . . . . . . p. 24
4 Simulação de Fluidos utilizando Smoothed Particle Hydrodynamics p. 26
4.1 Núcleos de Suavização . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 27
4.1.1 Núcleo de Suavização Poly6 . . . . . . . . . . . . . . . . . . . . . . p. 29
4.1.2 Núcleo de Suavização Spiky . . . . . . . . . . . . . . . . . . . . . . p. 30
4.1.3 Núcleo de Suavização Viscosity . . . . . . . . . . . . . . . . . . . . p. 30
4.2 As grandezas que governam os fluidos em SPH . . . . . . . . . . . . . . . . p. 31
4.2.1 Densidade das massas . . . . . . . . . . . . . . . . . . . . . . . . . p. 32
4.2.2 Pressão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 32
4.2.3 Força de Pressão . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 35
4.2.4 Força de Viscosidade . . . . . . . . . . . . . . . . . . . . . . . . . . p. 36
4.2.5 Forças externas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 37
4.2.5.1 Força da Gravidade . . . . . . . . . . . . . . . . . . . . . p. 37
4.2.5.2 Força de Superfície . . . . . . . . . . . . . . . . . . . . . p. 37
4.3 Integração Temporal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 40
4.4 Implementação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 42
5 Detecção e Tratamento de Colisões p. 45
5.1 Definições de Planos e suas Propriedades . . . . . . . . . . . . . . . . . . . p. 46
5.1.1 Distância Perpendicular Ponto-Plano . . . . . . . . . . . . . . . . . p. 47
5.1.2 Intersecção Reta-Plano . . . . . . . . . . . . . . . . . . . . . . . . . p. 48
5.2 Intersecções Partículas - Planos . . . . . . . . . . . . . . . . . . . . . . . . . p. 49
5.2.1 Obtenção do Raio de um Partícula . . . . . . . . . . . . . . . . . . . p. 49
5.2.2 Intersecção sem trajetória . . . . . . . . . . . . . . . . . . . . . . . p. 50
5.2.3 Intersecção com trajetória . . . . . . . . . . . . . . . . . . . . . . . p. 51
5.2.4 Interações entre dois ou mais planos . . . . . . . . . . . . . . . . . . p. 59
5.3 Balanço do recipiente . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 63
5.4 Implementação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 65
6 Busca por Vizinhos p. 66
6.1 Força Bruta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 66
6.1.1 Otimização . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 67
6.1.2 Vantagens . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 68
6.1.3 Desvantagens . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 68
6.2 Busca em Grid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 68
6.2.1 Inserção e Remoção no Grid . . . . . . . . . . . . . . . . . . . . . . p. 70
6.2.2 Otimização . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 71
6.2.3 Grids Fora da Origem . . . . . . . . . . . . . . . . . . . . . . . . . p. 72
6.2.4 Vantagens . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 73
6.2.5 Desvantagens . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 73
6.3 Implementação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 74
7 Resultados e Conclusões p. 76
7.1 Busca de Vizinhança . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 78
7.2 Conclusões . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 79
7.3 Trabalhos Futuros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 79
Referências p. 82
Lista de Figuras
1 Exemplo da simulação de um líquido interagindo com uma esfera (FOSTER;
FEDKIW, 2001). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 13
2 Um exemplo do grid estático de células que contém o fluido (STAM, 1999). . . p. 17
3 Retrocesso no tempo de execução para obtenção de ”x” (STAM, 1999) . . . . p. 18
4 Uma célula com a densidade p avaliada em seu centro, e a velocidade ~v de-
compostas em u, v e w (BRIDSON; MULLER-FISCHER; GUENDELMAN, 2006). . p. 19
5 Heightfield em duas dimensões: as colunas variam de altura h ao longo da
simulação (KASS; MILLER, 1990). . . . . . . . . . . . . . . . . . . . . . . . . p. 20
6 Uma Octree funcionando como um grid de fluido (LOSASSO; GIBOU; FEDKIW,
2004). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 25
7 O núcleo de suavização pode ser visto como uma esfera ao redor de uma
partícula. (KELAGER, 2006) . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 27
8 Duas partículas em proximidade. Por causa da propriedade de reflexão, W (~r, h) =
W (~q, h) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 28
9 O núcleo gaussiano com h = 1(KELAGER, 2006) . . . . . . . . . . . . . . . p. 28
10 Os três núcleos utilizados neste trabalho: Poly6, Spiky e Viscosity, respecti-
vamente. As linhas grossas representam o valor dos núcleos, as finas o valor
dos gradientes, e as tracejadas os laplacianos (MULLER; CHARYPAR; GROSS,
2003). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 31
11 A relação entre a densidade das massas e a pressão. No primeiro caso, as
partículas destacadas estão em equilíbrio. No segundo, a concentração de
partículas aumenta a densidade das massas e consequentemente, a pressão na
região. Por fim, no terceiro caso a baixa densidade das massas acarreta em
uma pressão baixa, porém existente (KELAGER, 2006). . . . . . . . . . . . . p. 34
12 A força exercida na superfície do fluido aponta para seu interior.(KELAGER,
2006) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 38
13 O erro introduzido pelo processo de integração temporal. A curva cheia re-
presenta a trajetória ideal, contínua da partícula A. Integrando a aceleração ~a
durante um intervalo de tempo t, a nova posição da partícula é B. Com 2t, a
partícula pára em C, visivelmente distante da trajetória correta. . . . . . . . . p. 41
14 No método Leap-Frog, a velocidade (u) ”pula” por cima da posição (r), e
vice-versa (KELAGER, 2006). . . . . . . . . . . . . . . . . . . . . . . . . . . p. 42
15 O diagrama das classes envolvidas no passo do método SPH. . . . . . . . . . p. 44
16 As situações do recipiente aberto e fechado. Em (A), a partícula não tem
como alcançar as porções dos planos infinitos que não compõem as paredes.
Em (B), uma partícula pode escapar pelo topo do recipiente e ainda assim
colidir com o plano que forma a parede da esquerda. . . . . . . . . . . . . . p. 46
17 O ponto A está do lado positivo de plano, enquanto B se encontra no lado
negativo. Consequentemente, ”a” é positivo e ”b” é negativo. . . . . . . . . . p. 47
18 No caso mais simples, se a distância |d| do centro C ao plano for menor do
que o raio r, então a partícula está colidindo com a superfície. . . . . . . . . p. 50
19 Uma partícula cruzando um plano ao longo de sua trajetória em um intervalo
de tempo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 50
20 Os três casos possíveis para o valor da distância d de C1 ao plano. Em (A), o
fato de d ser maior que o raio da partícula evidencia a não colisão. Em (B), d
é positiva e menor que r . Em (C), d é negativa pois C1 está do lado negativo
do plano. Nos últimos dois casos houve colisão. . . . . . . . . . . . . . . . . p. 53
21 O vetor ~u liga a posição final atual e a correta da partícula. Na posição inicial
C0, a esfera já toca o plano. . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 54
22 A trajetória de três dos pontos da partícula. Como todos viajam à mesma
velocidade e o ponto P0 é o que tem a menor distância (x) a percorrer até
chegar ao plano, ele é o que entra em contato primeiro. . . . . . . . . . . . . p. 54
23 A distância a de P1 ao plano vale r + |d|. Como d é garantidamente negativa
no caso da figura, o valor desta distância pode ser escrito como r − d. . . . . p. 55
24 Mesmo quando C1 se encontra do lado positivo do plano, a distância de P1 à
ele ainda vale r − d. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 55
25 Cf = C1 + ~u, onde o tamanho de ~u vale duas vezes o valor da distância do
ponto mais distante ao plano. . . . . . . . . . . . . . . . . . . . . . . . . . . p. 56
26 O eixo E a partir do qual começa a reflexão do centro passa pelo ponto CI . . p. 57
27 Utilizando o eixo E como eixo de reflexão, o comprimento do vetor ~u ainda
é r − d. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 57
28 ~v2 é o resultado da reflexão de ~v1 no plano com normal ~n. . . . . . . . . . . . p. 58
29 O que pode acontecer quando dois ou mais planos existem no ambiente. Por
causa da escolha do plano B para o teste, a posição final Cf não representa o
que aconteceria na realidade. . . . . . . . . . . . . . . . . . . . . . . . . . . p. 59
30 O resultado esperado da situação da figura 29. A partícula deve primeiro
colidir com o plano A, para então colidir com o plano B. . . . . . . . . . . . p. 59
31 A distância d de C1 ao plano indica se PI deve ser obtido, pela presença ou
não de colisão. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 61
32 A distância entre a posição inicial da partícula e sua posição no instante da
colisão vale |I − P0|. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 62
33 C0′ é o novo início da trajetória da partícula, e pode ser obtido por C0+(PI−P0). p. 63
34 A classe Box é a responsável pela detecção e o tratamento das colisões. . . . p. 65
35 Na esquerda, as estruturas que representam as partículas ocupam regiões con-
secutivas em um vetor. Na direita, cada estrutura sabe o endereço da próxima
em uma lista. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 67
36 Em um determinado momento (A), está sendo efetuada a busca por vizinhos
da partícula A. Testa-se esta relação entre A e B. Quando, em (B), for feita
a busca por vizinhos de B, não é necessário testar com A, pois a relação de
vizinhança é simétrica e este teste já foi efetuado, sendo elas vizinhas ou não. p. 68
37 Com um Grid, todas as partículas potencialmente vizinhas à partícula A se
encontram na célula de A e nas células diretamente adjacentes. . . . . . . . . p. 69
38 Uma partícula com centro C = (9, 3) fica na posição (4, 1) do Grid da imagem. p. 71
39 A região hachurada em azul contém todas as partículas candidatas à vizi-
nhança de A e B. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 71
40 Uma partícula com centro C = (12, 5) fica na célula (4, 1) em um Grid
deslocado (3, 2) da origem. . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 73
41 A classe Grid possibilita a busca rápida por vizinhança. . . . . . . . . . . . . p. 75
42 Cenas de seis pontos diferentes da simulação. . . . . . . . . . . . . . . . . . p. 77
43 A gravidade atua sempre nas partículas, inclusive nas respostas de colisões. . p. 81
44 Uma artéria modelada como uma malha de triângulos pode possuir dezenas
de milhares de polígonos. . . . . . . . . . . . . . . . . . . . . . . . . . . . p. 81
Lista de Tabelas
1 Valores utilizados na simulação. . . . . . . . . . . . . . . . . . . . . . . . . p. 76
2 Comparação entre os diferentes métodos de busca de vizinhança. Os valores
na coluna indicam o número de ciclos que cada método leva em média no
processador, dividido por mil. . . . . . . . . . . . . . . . . . . . . . . . . . p. 78
3 Porcentagem de processamento com e sem busca por Grid. . . . . . . . . . . p. 78
12
1 Introdução
A simulação realista de fluidos não é um problema novo. A maneira como líquidos se
comportam é intuitivamente natural para nós - É simples prever com um certo grau de precisão o
que acontece quando uma torneira é aberta e água é despejada numa bacia. Métodos e equações
para modelar fluidos já existem há várias décadas, como as equações elaboradas por Claude
Navier em 1822 e, posteriormente, George Stokes, em 1845 (WRENNINGE, 2003). Stam (2003)
afirma que há um consenso entre cientistas da área de que estas equações, conhecidas como
equações Navier-Stokes, são um modelo bastante apropriado para fluidos. As equações de
Navier-Stokes cumprem um papel no objetivo deste trabalho, de forma que serão brevemente
apresentadas no capítulo 2.
Estas equações são geralmente difíceis (ou impossíveis) de resolver analiticamente (STAM,
2003), de forma que desde então vários métodos foram propostos no campo da Dinâmica dos
Fluidos Computacional (Computational Fluid Dynamics - CDF) para simular fluidos resol-
vendo estas equações numericamente, tais como os métodos de Kass e Miller (1990), Foster e
Metaxas (1996), Stam (1999) e Losasso, Gibou e Fedkiw (2004).
O propósito deste trabalho é discutir e implementar um destes métodos, conhecido como
Smoothed Particle Hydrodynamics (SPH), inicialmente concebido por Lucy (1977) e Gingold
e Monaghan (1977), para simular eficientemente um recipiente retangular, em três dimensões e
com uma determinada quantidade de líquido. Este método foi escolhido por apresentar resul-
tados promissores (MULLER; CHARYPAR; GROSS, 2003), particularmente na interatividade com
o usuário, que é um dos objetivos principais deste trabalho. Muller, Charypar e Gross (2003)
comentam que o método é flexível o bastante para simular qualquer tipo de fluido, mas este
13
trabalho foca na simulação de líquidos.
O método SPH em si possui alguns pontos ”abstratos” em seu algoritmo: Não há definição
precisa da evolução do tempo de simulação, ou da estrutura de dados que representa os obstá-
culos do ambiente. Não é mencionada a forma de detecção de colisões entre as várias entidades
integrantes do sistema, e nem o que fazer no momento da colisão. Também fica em aberto a
questão do valor dos parâmetros iniciais ideais (como velocidade, massa, pressão, etc). Este
trabalho se propõe a solucionar estas questões, tendo em vista a eficiência da simulação e da
interatividade com o usuário.
A motivação do trabalho vem da necessidade de simulação realista de fluidos em diversas
áreas: Muller, Schirm e Teschner (2004) afirmam que a modelagem interativa de fluidos como a
água e o sangue humano é um componente essencial na cirurgia computacional, enquanto Stam
(2003) observa que é desejável incluir modelagens de fluidos em engines de jogos pois um dos
maiores objetivos dos mesmos é a imersão do jogador num mundo virtual plausível.
Figura 1: Exemplo da simulação de um líquido interagindo com uma esfera (FOSTER; FEDKIW,2001).
1.1 Objetivo Geral
Implementar uma aplicação consistindo de um recipiente em três dimensões, com posição
variável e definível interativamente pelo usuário da aplicação, e com uma determinada quanti-
dade de líquido em seu interior.
14
1.2 Objetivos Específicos
Para alcançar o objetivo geral, este trabalho assume os seguintes objetivos específicos:
• Analisar brevemente os métodos existentes na literatura para simulação de fluidos - em
particular os métodos baseados em Grid (Eulerianos) e os baseados em partículas (La-
grangeanos).
• Analisar mais especificamente o método conhecido como Smoothed Particle Hydrodyna-
mics (SPH) para a simulação de fluidos.
• Projetar um sistema para a implementação eficiente do SPH, tratando os componentes
de um sistema de simulação física que o método não aborda, como valores iniciais dos
parâmetros do sistema, busca de vizinhança eficiente e detecção eficiente de colisão entre
as partículas e o ambiente externo (neste caso, o recipiente).
15
2 Notação e Conceitos Matemáticos
Em 1822, Claude Louis Marie Henri Navier derivou as equações que regem o movimento
de um fluido viscoso. Ele partiu do trabalho de Euler de 1755, que já havia derivado as equações
para fluidos ideais não viscosos. Após isso, em 1845 George Gabriel Stokes derivou as equações
na maneira como elas são conhecidas hoje.
Fluidos incompressíveis são fluidos cuja densidade não varia. Líquidos são fluidos com-
pressíveis, mas como a compressão só se dá sobre grande pressão eles são geralmente simplifi-
cados como incompressíveis(BRIDSON; MULLER-FISCHER; GUENDELMAN, 2006).
A evolução do estado de um fluido incompressível no tempo é dada pelas seguintes especi-
alizações das equações de Claude Navier e George Stokes (M.S.CRAMER, 2002):
∂ρ
∂t+∇. (ρ~v) = 0 (2.1)
ρ
(∂~v
∂t+ ~v.∇~v
)= −∇p + ~fext + µ∇2~v (2.2)
A equação 2.1 rege a conservação da massa no líquido, e a equação 2.2 dita a conservação do
momento. Em ambas:
ρ é a densidade do líquido no ponto considerado;
p é a pressão do líquido no ponto considerado;
t é o tempo;
16
v é a velocidade do líquido no ponto considerado;
~fext é o vetor das forças externas ao líquido, e.g. a gravidade.
É útil também definir as seguintes convenções de notação:
∂f∂x
é a derivada de f em relação a x ;
∇f é o gradiente de uma função f de n variáveis, e é definido como o vetor(
∂f∂x1
, ∂f∂x2
, . . . , ∂f∂xn
);
∇.f é o operador vetorial de divergência de uma função f de n variáveis, e é definido
como a soma dos componentes de ∇f ;
∇2f é o operador laplaciano escalar, definido como a divergência do gradiente de f , i.e.
∇.∇f . No espaço cartesiano Rnele é definido como a soma das derivadas parciais
de segunda ordem não mistas, ou seja, ∇2f = ∂2f∂x2
1+ ∂2f
∂x22
+ . . . + ∂2f∂x2
n;
∇2~v é o operador laplaciano vetorial.
17
3 Dinâmica dos Fluidos Computacional
Os modelos matemáticos gerados pelas equações descritas no capítulo 2 só possuem solu-
ção analítica nos casos mais simples. Com a introdução da ferramenta computacional come-
çaram a surgir métodos para resolver estas equações numéricamente na tentativa de simular o
comportamento de fluidos. Tais métodos podem ser divididos em dois tipos distintos (BRIDSON;
MULLER-FISCHER; GUENDELMAN, 2006): Métodos baseados em Grids (ou Eulerianos), e mé-
todos baseados em Partículas (ou Lagrangeanos). Como o sistema proposto por este trabalho é
baseado em partículas, uma discussão mais elaborada será feita sobre esta segunda categoria de
métodos na seção 3.2.
3.1 Métodos baseados em Grids para Representação de Flui-dos
Os métodos baseados em Grids (Malhas), têm sua base na discretização do espaço onde
ocorre a simulação em subespaços de tamanho fixo (Células), transformando o espaço inteiro
em uma Malha n-dimensional. Nestes métodos, a evolução do fluido é dada pela evolução do
espaço de velocidade v, o espaço de densidade ρ e o espaço de pressão p, avaliados em cada
célula. A figura 2 ilustra o grid.
Figura 2: Um exemplo do grid estático de células que contém o fluido (STAM, 1999).
18
A maneira como as três quantidades são avaliadas na malha varia de acordo com o método.
3.1.1 Método Semi-Lagrangeano
Stam (1999) considera o valor das grandezas de interesse no centro de cada célula, e obtém
este valor usando o valor que a ”porção de líquido” em questão possuía no início da simulação.
Por exemplo, para obter a velocidade do fluido em uma célula qualquer do grid, o método
percorre o campo vetorial de velocidade no sentido inverso (”voltando” no tempo), e obtém a
velocidade que aquela porção de líquido possuía no instante inicial. Segundo Stam, a velocidade
na nova célula é igual à velocidade na célula original, no tempo inicial, conforme ilustra a figura
3.
Figura 3: Retrocesso no tempo de execução para obtenção de ”x” (STAM, 1999)
Este método é chamado de Semi-Lagrangeano porque os métodos originais lagrangeanos
para a simulação de fluidos são baseados em partículas, e não em grids de posições fixas. Neste
método, a ”partícula” é apenas uma abstração para a resolução do sistema euleriano, não apare-
cendo explicitamente na solução em código.
Stam aponta que seu método é isento de instabilidade numérica, e o considera extrema-
mente eficiente na simulação realista em tempo real de gases, enquanto Foster e Metaxas (1996)
indicam que a dissipação de massa resultante o torna impróprio para a simulação de líquidos.
3.1.2 Marker-and-Cell (MAC) Grids
O método Marker-and-Cell, originalmente proposto por Harlow e Welch (1965), tem como
característica principal o fato de avaliar diferentes quantidades em diferentes regiões das células,
19
ao contrário do método anterior que sempre as avalia nos centros.
A pressão p continua sendo avaliada no centro de cada célula, porém o vetor de velocidade
~v é decomposto em suas três coordenadas, u, v, e w, e cada coordenada é avaliada em um dos
três pares de faces opostos da célula, conforme ilustra a figura 4.
Figura 4: Uma célula com a densidade p avaliada em seu centro, e a velocidade ~v decompostasem u, v e w (BRIDSON; MULLER-FISCHER; GUENDELMAN, 2006).
O nome do método deriva da necessidade de saber aonde se encontra a superfície do fluido
(por exemplo, a superfície do líquido). Para resolver este problema, Harlow e Welch (1965)
incluem no conjunto de células um grupo de partículas, chamadas de Marker Particles (partícu-
las marcadoras), inicialmente posicionadas em células que sabidamente contém o fluido. Estas
partículas se movem com o fluido, utilizando uma interpolação linear das velocidades nas várias
faces da célula que contém a partícula para determinar a velocidade da mesma.
Desta forma, a qualquer momento da simulação, células que contém partículas são células
com fluido, e células sem partículas estão vazias. Células que possuem vizinhas vazias são
células de superfície.
Este método simplifica a obtenção da superfície do líquido, sob custo da perda da garantia
da estabilidade numérica de Stam (1999).
20
3.1.3 Aproximações HeightField
As aproximações baseadas em Heightfields para grids surgiram do fato de que, para algumas
aplicações, é desnecessário e custoso levar em conta todo o volume do fluido considerado para
se obter uma simulação realista. Bridson, Muller-Fischer e Guendelman (2006) apontam como
exemplo a simulação de um lago. Do ponto de vista de um observador externo, só interessa a
evolução da superfície do lago.
Nestes métodos, o grid Euleriano é visto como um conjunto de ”colunas” de fluido. Ao
longo da simulação, a altura de cada coluna varia e o conjunto de alturas h das colunas ca-
racteriza a superfície do fluido. As desvantagens apontadas por Kass e Miller (1990) incluem
a inabilidade de simular ”respingos” de água e ondas quebrando. A vantagem marcante é a
grande redução da complexidade computacional, conforme indicam Kass e Miller (1990) em
seu trabalho relacionado.
Figura 5: Heightfield em duas dimensões: as colunas variam de altura h ao longo da simulação(KASS; MILLER, 1990).
Os métodos baseados em Grids, por possuírem pontos fixos no espaço onde avaliar o estado
do fluido, possuem como vantagem a simplicidade de implementação (de fato, Stam (2003)
desafia o leitor a implementar um algoritmo mais simples que o seu, que possui em torno de
100 linhas de código escrito na linguagem C). Estes métodos se estabeleceram na literatura
21
por serem úteis para simulações estáticas e sem variações durante a execução. São, porém,
computacionalmente caros e de difícil uso em aplicações que requerem interatividade com o
usuário em tempo real.
3.2 Sistemas de Partículas
Sistemas de partículas foram originalmente propostos por Reeves (1983), como uma ma-
neira de modelar objetos que o autor denominou fuzzy (”nebulosos”), como o fogo, água, fu-
maça, etc. As características principais deste tipo de objeto são as seguintes:
1. Objetos deste tipo não são representados por um conjunto de polígonos que delimitam sua
superfície, mas por uma nuvem de tipos primitivos denomidados Partículas que definem
o seu volume.
2. Tais objetos não são estáticos - as partículas se movimentam, alteram sua forma, saem e
entram do sistema com o passar do tempo.
3. Geralmente há um fator estocástico (aleatório) que influencia o comportamento geral do
objeto, isto é, o comportamento do objeto não é determinístico.
Desde então, sistemas de partículas foram refinados para representar uma grande variedade de
elementos - sólidos deformáveis, como por exemplo órgãos do corpo humano (JAILLET; SHA-
RIAT; VANDORPE, 1997), fogo e fumaça (STAM, 2000), líquidos (MULLER; CHARYPAR; GROSS,
2003; PREMOZE et al., 2003; MUELLER et al., 2004), entre outros. Cada tipo de elemento pos-
sui suas características distintas - objetos sólidos, por exemplo, atribuem menos liberdade de
movimento às partículas, enquanto é comum modelar o fogo como um conjunto de partículas
onde cada partícula ”nasce”, se desloca para cima com algum grau de variação (o elemento
estocástico), e ”morre”.
Para líquidos, que são o foco deste trabalho, o sistema é bem definido: A massa do líquido é
fixa, e como a massa é distribuída entre as partículas, o número de partículas é fixo. A evolução
22
do sistema ao longo do tempo é regida por equações conhecidas da mecânica dos fluidos, de
forma que não existe o elemento aleatório.
Sistemas de partículas possuem vantagens sobre a maneira clássica baseada em malhas para
simulação de fluidos: São geralmente mais computacionalmente eficientes, por avaliarem os pa-
râmetros que regem o estado do fluido apenas em pontos específicos do espaço, ao contrário de
em todo ele. Por serem mais eficientes, eles são mais facilmente implementados em aplicações
que requerem interatividade com o usuário. Em contrapartida, uma desvantagem significativa é
a dificuldade da geração da superfície do líquido representado pelas partículas.
3.2.1 Smoothed Particle Hydrodynamics
O método Smoothed Particle Hydrodynamics (SPH) (LUCY, 1977; GINGOLD; MONAGHAN,
1977; MONAGHAN, 1992) foi originalmente concebido para a solução de modelos astrofísicos
em três dimensões, mas é genérico a ponto de ter sido adaptado para vários outros campos da
engenharia e computação, particularmente a simulação de fluidos compressíveis e não compres-
síveis.
A síntese do método está em seu nome: As grandezas relevantes ao sistema são discre-
tizadas nas partículas através de uma interpolação ponderada do valor destas grandezas nas
partículas vizinhas. Ou seja, o estado de cada partícula depende do estado de suas partículas
vizinhas. Generalizando para uma grandeza A na posição espacial r:
As (r) =n∑
j=1
mjAj
ρj
W (r − rj, h) (3.1)
Onde:
mj é a massa da partícula j;
ρj é a densidade da partícula j;
W é a função de suavização utilizada, arbitrária;
23
h é uma constante que define o raio máximo de atuação da função W.
A função de suavização W (Smoothing Kernel) é o que caracteriza o método. É prático pensar
neste núcleo como uma função que leva todo valor a uma distância r − rj do centro de uma
esfera de raio h a um valor de interpolação.
Um exemplo da aplicação da equação 3.1 é o cálculo da densidade em uma partícula que se
encontra na posição r:
ρs (r) =n∑
j=1
mjρj
ρj
W (r − rj, h) =n∑
j=1
mjW (r − rj, h) (3.2)
Como esta equação pode auxiliar na simulação realista de fluidos? Por praticidade, apresenta-
se novamente a equação da conservação do momento de Navier-Stokes ( 2.2 na página 15):
ρ
(∂v
∂t+ v.∇v
)= −∇p + ~fext + µ∇2~v
Muller, Charypar e Gross (2003) afirmam que a expressão ∂v∂t
+ v.∇v pode ser simplificada
a ∂v∂t
pois, ao contrário dos sistemas estáticos baseados em malhas, em sistemas de partículas
as mesmas se deslocam junto com o fluido, significando que a variação do espaço vetorial de
velocidade do fluido é representado simplesmente pela variação no tempo da velocidade das
partículas que o compõem. Desta forma, dividindo esta simplificação de 2.2 por ρ, obtém-se:
∂v
∂t=−∇p + ρ~fext + µ∇2~v
ρ(3.3)
Que lembra intuitivamente a segunda lei de Newton, aplicada para líquidos:
~a =~F
ρ
O método então consiste em determinar as grandezas que compõem F através dos núcleos de
suavização, utilizando a equação 3.1. O método SPH foi o escolhido para a simulação do fluido
neste trabalho e será discutido com mais detalhes no capítulo 4.
24
3.2.2 Método Moving Particle Semi-Implicit (MPS)
Este método, originalmente proposto por Koshizuka, Nobe e Oka (1998) também utiliza
a idéia de que o estado de uma partícula é influenciado por suas partículas vizinhas. Desta
forma, existe neste método uma função de peso w que, semelhantemente ao método SPH, in-
terpola um valor a partir de uma distância r − rj. A diferença básica entre os dois métodos é o
comportamento da função de peso e o algoritmo de atualização das grandezas do líquido.
Premoze et al. (2003) citam diversas aplicações do método MPS, como transições de fase,
fluxo multifásico, fluidos com sedimentação e estruturas elásticas.
3.3 Outras Metodologias de Simulações de Fluidos
Além dos métodos apresentados nas seções anteriores, vários outros foram propostos ao
longo dos anos e com o avanço da capacidade computacional, geralmente com a intenção de
solucionar eficientemente um problema específico.
Métodos híbridos, que misturam em sua solução partículas e grids eulerianos estáticos, têm
se tornado bastante populares por apresentarem as vantagens de cada abordagem. Um exemplo
de uma abordagem híbrida é o uso de grids para cálculo das grandezas do fluido, e partículas
para derivação da superfície implícita (FOSTER; METAXAS, 1996; FOSTER; FEDKIW, 2001).
Losasso, Gibou e Fedkiw (2004) ilustram um método que implementa a abordagem de grid
utilizando uma Octree. Enquanto Foster e Metaxas (1996) consideravam as faces da célula
como ponto de interesse do cálculo de certas grandezas do líquido (figura 4), a abordagem com
Octree decompõe a célula em nodos, tomando assim as grandezas nas ”quinas” da célula para
interpolação, conforme ilustra a figura 6.
Finalmente, avanços na tecnologia de hardware gráfico possibilitam otimizações de perfor-
mance que exploram os recursos dedicados da placa aceleradora de vídeo, conforme ilustram
Amada et al. (2004), que indicam conseguirem o dobro de velocidade de processamento em sua
versão utilizando a Unidade de Processamento Gráfico, em comparação ao processamento em
25
Figura 6: Uma Octree funcionando como um grid de fluido (LOSASSO; GIBOU; FEDKIW, 2004).
CPU.
26
4 Simulação de Fluidos utilizandoSmoothed Particle Hydrodynamics
O método Smoothed Particle Hydrodynamics (SPH) foi apresentado brevemente na seção
3.2.1. Por ter sido o método escolhido para este trabalho, é válido discutí-lo com maior profun-
didade.
Conforme mencionado, o coração do método SPH está em seus núcleos de suavização. São
eles que descrevem a forma como as partículas interagem entre si e o comportamento geral do
fluido simulado. A equação 3.1 na página 22 dita o valor de uma grandeza A qualquer na posição
representada pelo ponto r. Muller, Charypar e Gross (2003) afirmam que as derivadas destas
grandezas, comuns em equações de fluidos, afetam apenas o núcleo de suavização envolvido.
Desta forma, podemos simplificar o gradiente de uma grandeza A, que é o vetor formado pelas
derivadas parciais de A, conforme mencionado no capítulo 2:
∇As (r) =n∑
j=1
mjAj
ρj
∇W (r − rj, h) (4.1)
E seu laplaciano, que é a soma das derivadas parciais de segunda ordem não mistas:
∇2As (r) =n∑
j=1
mjAj
ρj
∇2W (r − rj, h) (4.2)
Por que esta propriedade é interessante, ou até mesmo relevante?
A equação 3.3 na página 23 apresenta a variação da velocidade em um fluido, representado
por um sistema de partículas. O método apresentado por Muller, Charypar e Gross (2003)
27
consiste em resolver a expressão−∇p+ ~fext+µ∇2~v através de SPH, de forma a obter a variação
da velocidade e simular o fluido. O fato da derivação das grandezas poder ser realizada em
função dos núcleos de suavização simplifica bastante o processo, conforme será demonstrado
na seção 4.2. Mas antes, é preciso abordar os núcleos de suavização.
4.1 Núcleos de Suavização
A maneira mais simples e intuitiva de abordar os núcleos de suavização é pensar nas partí-
culas no espaço. Se isolarmos uma destas partículas, a qual chamaremos de p, pensarmos em
uma esfera hipotética de raio h ao redor de p e considerarmos que cada partícula dentro desta
esfera hipotética contribui para o comportamento de p, então um núcleo de suavização nada
mais faz do que explicitar a maneira pela qual esta contribuição é feita. A região compreendida
dentro desta esfera hipotética é chamada de vizinhança de p e todas as partículas dentro desta
região são ditas na vizinhança de p.
Figura 7: O núcleo de suavização pode ser visto como uma esfera ao redor de uma partícula.(KELAGER, 2006)
Um núcleo de suavização é uma função. Segundo Muller, Charypar e Gross (2003), se
esta função assumir as propriedades de reflexão e normalização, ela é de segunda ordem de
precisão, o que significa dizer que o erro obtido ao interpolar A(r) utilizando o núcleo é da
ordem de O (h2)ou menor (KELAGER, 2006), onde h é o raio da esfera.
A primeira propriedade diz que o núcleo deve ser uma função reflexiva, ou seja, W (~r, h)
deve ser igual a W (−~r, h). Isto permite que as partículas sejam relacionadas em qualquer
ordem, como ilustra a figura 8.
28
Figura 8: Duas partículas em proximidade. Por causa da propriedade de reflexão, W (~r, h) =W (~q, h)
A segunda propriedade diz que o núcleo deve ser normalizado, para que as constantes pre-
sentes na função da grandeza A sejam interpoladas corretamente (GINGOLD; MONAGHAN, 1977),
ou seja:
∫W (~r, h) d~r = 1
Adicionalmente, existe um limite ao módulo de ~r, acima do qual a função sempre vale zero.
Na analogia de que o núcleo de suavização representa uma esfera, este limite é o raio da mesma,
h. Este valor é também chamado de raio de suporte do núcleo, e serve para limitar o número de
partículas em uma vizinhança. Isto é computacionalmente interessante, visto que uma partícula
não dependerá de todas as outras do sistema. A figura 9 exemplifica o gráfico do núcleo de
suavização gaussiano com h = 1.
Figura 9: O núcleo gaussiano com h = 1(KELAGER, 2006)
Kelager (2006) afirma que a escolha de um valor para o raio de suporte h apropriado é
vital para a simulação robusta e estável do fluido. O autor afirma que um valor grande demais
diminuiria a contribuição de partículas próximas ao centro no núcleo, e propõe uma maneira
29
simples de obter um valor aceitável. Dado o volume total V do líquido, o número n de partículas
que o representa, e o número médio x de partículas dentro de uma região de vizinhança qualquer,
o valor do raio de suporte pode ser obtido por:
h =3
√3V x
4πn
Esta fórmula é obtida do fato que nV
é a densidade das partículas no líquido, ou seja, o
número de partículas por unidade de volume. Multiplicando este valor pelo volume de uma
esfera de raio h, obtém-se o número de partículas que ocupam aquele volume:
x =(
n
V
)4
3πh3
Ambas equações possuem duas incógnitas não disponíveis imediatamente, x e h. Kela-
ger (2006) sugere que x seja testado até o menor número possível que simule o fluido e suas
propriedades convincentemente.
Nas seguintes subseções serão apresentados os núcleos de suavização utilizados neste tra-
balho, conforme o modelo de Muller, Charypar e Gross (2003). Em todos os casos, se |~r| > h,
W (~r, h) = 0:
4.1.1 Núcleo de Suavização Poly6
O primeiro núcleo proposto por Muller, Charypar e Gross (2003) é atraente por não requerer
o cálculo de raízes, computacionalmente caros:
Wpoly6 (~r, h) =315
64πh9
(h2 − |~r|2
)3(4.3)
O gradiente de núcleo Poly6 vale:
∇Wpoly6 (~r, h) = − 945
32πh9~r(h2 − |~r|2
)2(4.4)
30
E seu laplaciano:
∇2Wpoly6 (~r, h) = − 945
32πh9
(h2 − |~r|2
) (3h2 − 7|~r|2
)(4.5)
4.1.2 Núcleo de Suavização Spiky
Este núcleo é utilizado no cálculo da força de pressão sobre as partículas, e sua utilidade
será explicada na subseção 4.2.2:
Wspiky (~r, h) =15
πh6(h− |~r|)3 (4.6)
Com gradiente:
∇Wspiky (~r, h) = − 45
πh6
~r
|~r|(h− |~r|)2 (4.7)
E laplaciano:
∇2Wspiky (~r, h) = − 90
πh6
1
|~r|(h− |~r|) (h− 2|~r|) (4.8)
4.1.3 Núcleo de Suavização Viscosity
O terceiro núcleo proposto por Muller, Charypar e Gross (2003) é utilizado no cálculo da
força de viscosidade atuando sobre as partículas, por uma razão que será discutida na subseção
4.2.4:
Wviscosity (~r, h) =15
2πh3− |~r|3
2h3+|~r|2
h2+
h
2|~r|− 1 (4.9)
Com gradiente:
31
∇Wviscosity (~r, h) =15
2πh3~r
(−3|~r|2h3
+2
h2− h
2|~r|3
)(4.10)
E laplaciano:
∇2Wviscosity (~r, h) =45
πh6(h− |~r|) (4.11)
Figura 10: Os três núcleos utilizados neste trabalho: Poly6, Spiky e Viscosity, respectivamente.As linhas grossas representam o valor dos núcleos, as finas o valor dos gradientes, e as tracejadasos laplacianos (MULLER; CHARYPAR; GROSS, 2003).
4.2 As grandezas que governam os fluidos em SPH
Finalmente é possível aplicar o método SPH sobre a equação da variação da velocidade no
fluido (equação 3.3). Para isso é preciso utilizar os núcleos de suavização em cada grandeza
envolvida na expressão −∇p + ~fext + µ∇2~v.
Observando novamente a equação geral do método (equação 3.1 na página 22),
As (r) =n∑
j=1
mjAj
ρj
W (r − rj, h)
nota-se que, além do valor da própria grandeza A, a equação também depende da massa
das partículas (m) e da densidade das massas (ρ) em uma vizinhança. A massa é conhecida
e constante, definida no início da simulação. A densidade das massas varia de acordo com a
32
vizinhança de uma partícula, de forma que esta grandeza sempre deve ser calculada primeiro,
antes de ser possível obter as outras. Esta é a única restrição de ordem no método, todas as
outras grandezas podem ser calculadas na ordem que for mais conveniente.
4.2.1 Densidade das massas
O valor da densidade das massas em partícula pode ser entendido como a quantidade de
massa em sua vizinhança. Quanto maior o número de partículas em uma determinada vizi-
nhança, maior a densidade das massas.
A equação para o cálculo da densidade das massas é obtida simplesmente substituindo-se
a variável A por ρ na equação 3.1. O resultado desta substituição foi obtido na seção 3.2.1 na
página 22:
ρs (r) =n∑
j=1
mjρj
ρj
W (r − rj, h) =n∑
j=1
mjW (r − rj, h) (4.12)
Vale notar que a densidade das massas em uma partícula depende unicamente do valor da
massa das partículas em sua vizinhança e de sua própria massa.
4.2.2 Pressão
De posse da densidade das massas, é possível calcular a pressão de cada partícula direta-
mente, sem a necessidade da interpolação do método SPH. Este valor pode ser obtido através
da fórmula do gás ideal (KELAGER, 2006):
pV = nRT (4.13)
onde V é o volume por unidade de massa, n é o número de partículas de gás em 1 mol, R
é a constante universal dos gases, e T é a temperatura. Se a massa do fluido é constante (ou
seja, não existe introdução ou perda de fluido durante a simulação) e a temperatura também, a
expressão nRT pode ser abreviada para uma constante k. O resultado se torna:
33
pV = k (4.14)
V é o volume em uma unidade de massa. A densidade é definida como a quantidade de
massa por volume. Em uma unidade de massa, obtém-se:
ρ =m
V
ρ =1
V
V =1
ρ
Substituindo na equação 4.14, obtém-se:
p1
ρ= k
p = kρ (4.15)
O valor p obtido na equação 4.15 é utilizado posteriormente no cálculo da força de pressão
exercida sobre uma partícula (subseção 4.2.3). Ele pode ser visto como o ”potencial” que uma
partícula apresenta para influenciar suas vizinhas. Quanto maior a densidade das massas de uma
partícula, maior o valor da pressão p, e consequentemente maior a força que ela exercerá para
repelir suas vizinhas. A figura 11 ilustra esta relação.
Desbrun e Gascuel (1996) comentam que o uso da equação 4.15 resulta em forças pura-
mente repulsivas. Em aplicações em astrofísica, a força de pressão é comumente combinada
com forças gravitacionais que balanceiam a repulsão. Forças de pressão puramente repulsivas
fazem sentido em gases, que tendem a se expandir e ocupar todo o espaço, mas líquidos ten-
dem a manter uma coesão interna (KELAGER, 2006). Desbrun e Gascuel (1996) propõem uma
34
Figura 11: A relação entre a densidade das massas e a pressão. No primeiro caso, as partículasdestacadas estão em equilíbrio. No segundo, a concentração de partículas aumenta a densidadedas massas e consequentemente, a pressão na região. Por fim, no terceiro caso a baixa densidadedas massas acarreta em uma pressão baixa, porém existente (KELAGER, 2006).
alternativa à equação 4.15:
p = k (ρ− ρ0) (4.16)
Onde ρ0 é chamado de densidade de repouso e depende do tipo de fluido. O resultado é
que partículas distantes tenderão a se atrair, enquanto que partículas muito próximas tenderão
a se afastar. Desbrun e Gascuel (1996) identificam duas vantagens significativas na inclusão de
ρ0 à equação 4.15. Primeiro, se as partículas possuírem a mesma massa (ou seja, se a massa
do líquido for proporcionalmente distribuída entre as partículas), elas tenderão a se distribuir
uniformemente no espaço. Segundo, como as partículas tenderão a encontrar um estado de
equilíbrio entre a atração e a repulsão, a densidade das massas tenderá a ser constante, o que re-
sulta em um volume total aproximadamente constante. Desta forma, o fluido tenderá a assumir
este estado de equilíbrio após uma deformação.
Kelager (2006) chama a atenção para um cuidado em especial com a constante k. Este
valor, como mencionado, representa as três constantes nRT da equação 4.13, mas os valores
físicos reais dessas grandezas não podem ser utilizados na obtenção de k por gerarem valores
absurdamente grandes, fato que também foi observado durante a implementação deste sistema.
Desta forma, k deve ser obtido empiricamente, idealmente sendo o maior possível que ainda
permita uma simulação realista. Quanto maior este valor, menor deve ser o intervalo de tempo
entre iterações da simulação, para evitar instabilidade no fluido.
Uma hipótese para os valores exageradamente grandes de k obtidos por Kelager (2006) é
35
que a lei do gás ideal, além de ser um modelo estritamente teórico, diz respeito apenas a substân-
cias no estado gasoso (ATKINS, 1999). Como fluidos no estado líquido possuem pressão interna
e temperatura intuitivamente mais baixas do que gases, faz sentido que k deva ser amenizada
quando simulando líquidos.
De posse da densidade e da pressão nas partículas, pode-se extrair a força de pressão resul-
tante em cada uma delas.
4.2.3 Força de Pressão
Retornando à equação 3.3, a parcela relevante à força de pressão(−∇p) é obtida pelo gra-
diente do núcleo de suavização Spiky, da seção 4.1.2 na página 30:
−∇p (r) = −n∑
j=1
mjpj
ρj
∇W (r − rj, h) (4.17)
Por que este núcleo deve ser utilizado, ao invés do núcleo padrão Poly6 (subseção 4.1.1 na
página 29)?
Observando novamente a equação 4.17, nota-se que o gradiente do núcleo de suavização
deve ser utilizado. O gráfico do gradiente do núcleo Poly6 (ilustrado na figura 10 na página 31)
mostra que conforme |~r| tende a 0, ∇W (~r, h) também o faz. Espacialmente, isto significa que
quanto mais próximas duas partículas estiverem de si, menor a força de pressão gerada entre
elas. Isto é exatamente o oposto do comportamento intuitivo para uma força de pressão: partícu-
las muito próximas deveriam se repelir com uma intensidade maior do que partículas distantes.
Desbrun e Gascuel (1996) propuseram então o núcleo Spiky, cujo gráfico da figura 10 ilustra
que |∇W (~r, h) | tende a aumentar conforme |~r| tende a 0. Tal estratégia foi posteriormente
adotada por Muller, Charypar e Gross (2003) e Kelager (2006) .
Desbrun e Gascuel (1996) atentam ao fato de a força resultante da equação 4.17 não ser
simétrica. Em outras palavras, a força que uma partícula exerce sobre sua vizinha não será igual
à força que ela sofre em retorno, a não ser que as densidades de ambas sejam iguais. Como
resultado, a lei de ação-reação não é garantida e o resultado não é realista.
36
Estes mesmos autores propõem a seguinte variação à equação, que garante a simetria:
−∇p (~ri) = −ρi
n∑j=1
mj
(pi
ρ2i
+pj
ρ2j
)∇W (~r − ~rj, h)
Muller, Charypar e Gross (2003) propõem uma outra alternativa, computacionalmente mais
eficiente:
−∇p (~ri) = −n∑
j=1
(pi + pj
2
)mj
ρj
∇W (~r − ~rj, h)
4.2.4 Força de Viscosidade
A viscosidade em um fluido representa a capacidade do fluido de resistir à deformação.
Conforme o fluido se deforma, suas moléculas atritam entre si, gerando calor e reduzindo a
energia cinética (KELAGER, 2006).
A parcela na equação 2.2 referente à viscosidade é µ∇2~v. Interpolando com o método SPH
obtém-se:
µ∇2~v = µn∑
j=1
~vj∇2W (r − rj, h) (4.18)
Esta equação gera forças assimétricas, tal qual a equação 4.17, mas pelo motivo de que a
velocidade varia entre as partículas. Uma alternativa é utilizar no lugar da velocidade absoluta
a velocidade relativa entre as partículas (MULLER; CHARYPAR; GROSS, 2003):
µ∇2~v = µn∑
j=1
(~vj − ~vi)∇2W (r − rj, h)
Muller, Charypar e Gross (2003) afirmam que, se ∇2W (~r − ~rj, h) não for positivo para
todo |~r| < h, a força de viscosidade ao invés de agir como um atenuante irá introduzir energia
no sistema e aumentar a velocidade das partículas. Desta forma, tanto o núcleo Poly6 quanto
o Spiky não podem ser utilizados no cálculo da força de viscosidade, por possuírem valores
negativos em determinados pontos (figura 10 na página 31). Os autores então propuseram o
37
núcleo Viscosity, cujo laplaciano é positivo em todo o domínio, como a própria figura 10 ilustra.
4.2.5 Forças externas
A última parcela restante na expressão −∇p + ~fext + µ∇2~v diz respeito às forças externas
que influenciam o fluido. São elas a força gerada pela aceleração da gravidade e a força gerada
na superfície do fluido, ou seja:
~fext = ~fgrav + ~fsup
onde:
• ~fgrav equivale a força exercida pelo campo gravitacional da Terra (subseção 4.2.5.1),
• ~fsup é a força exercida na superfície do líquido (subseção 4.2.5.2).
4.2.5.1 Força da Gravidade
A força que a aceleração gravitacional exerce sobre o fluido é igual em todas as partículas,
e vale (MULLER; CHARYPAR; GROSS, 2003):
~fi = ρi~g
onde ~g é a aceleração da gravidade, geralmente denotada como [ 0.0 −9.82 0.0 ]T .
4.2.5.2 Força de Superfície
A segunda força externa atuando sobre o fluido é aquela agindo sobre a superfície do
mesmo. Esta força aponta sempre para o interior do fluido e serve para suavizar curvaturas
acentuadas (figura 12) (MORRIS, 2000).
Para simular a força de superfície Muller, Charypar e Gross (2003) utilizam as idéias de
Morris (2000), que diz:
38
Figura 12: A força exercida na superfície do fluido aponta para seu interior.(KELAGER, 2006)
fsup(r) = −σ∇2c(r)~n(r)
|~n(r)|(4.19)
onde:
• σ é um coeficente de tensão que depende das duas substâncias que interagem na superfície
(como o ar agindo na superfície da água);
• ~ni é o vetor normal que aponta para dentro do fluido na posição da partícula i;
• ci é o valor do campo de cor na posição da partícula i.
O campo de cor é um quantidade adicional que sempre vale exatamente 1 nas posições do
espaço que contém fluido e 0 em todos os outros pontos. Assim, as partículas representam os
pontos não-nulos deste campo. Com a interpolação do método SPH, este campo vale
c(r) =∑j
cjmj
ρj
W (r − rj, h)
=∑j
mj
ρj
W (r − rj, h)
E o gradiente deste valor, que vale
~n(r) = ∇c(~r)
39
=∑j
mj
ρj
∇W (r − rj, h)
É o vetor normal que aponta para dentro da superfície do fluido na posição ri (MULLER;
CHARYPAR; GROSS, 2003). Este vetor só possui módulo diferente de zero na região próxima
da superfície, de forma que |~n| tende a zero conforme as partículas se distanciam das bordas
do fluido. E conforme |~n| tende a zero, ~n|~n| tende a um vetor com módulo infinitamente grande,
introduzindo instabilidade no fluido (KELAGER, 2006). Uma maneira de evitar esta instabilidade
é calcular a força de superfície apenas para partículas onde |~n| > ε , onde ε é algum limite
observado empiricamente. Kelager (2006) recomenda o uso de
ε =
√ρ
x
onde ρ é a densidade do material e x é o número médio de partículas em uma vizinhança.
Outra possibilidade é o limite de Morris (2000) , que vale
ε =0.01
h
onde h é o raio de suporte dos núcleos de suavização. Finalmente, o laplaciano do campo
de cor c vale
∇2c(r) =∑j
mj
ρj
∇2W (r − rj, h)
Agora é possível obter o valor da força de superfície atuando sobre as partículas:
fsup(r) = −σ∇2c(r)~n(r)
|~n(r)|
fsup(r) = −σ(∑
j
mj
ρj
∇2W (r − rj, h))(∑
jmj
ρj∇W (r − rj, h))
|(∑jmj
ρj∇W (r − rj, h))|
(4.20)
40
4.3 Integração Temporal
O resultado da aplicação do método SPH sobre a equação de conservação de momento ( 3.3
na página 23) é um vetor de aceleração ~a para cada partícula. O passo seguinte é aplicar este
vetor sobre as partículas para atualizar suas velocidades e posições e assim simular a dinâmica
do fluido.
Idealmente, todo este cálculo seria feito instantâneamente e a posição de cada partícula seria
atualizada continuamente, em um intervalo de tempo absolutamente pequeno. Infelizmente isto
não é possível, o processamento das várias grandezas do fluido leva tempo. Assim, é preciso
aplicar o vetor de aceleração obtido em um instante qualquer durante um intervalo de tempo
consideravelmente grande, para ”anular” pelo menos uma parcela do custo do processamento
do estado do fluido.
A escolha do tamanho deste intervalo de tempo (ou tamanho do frame) deve ser cuidadosa:
Valores pequenos demais significarão uma simulação lenta, e valores grandes demais introdu-
zirão instabilidade na simulação. A razão para esta instabilidade é que o vetor de aceleração é
simplesmente isto, um vetor com direção, módulo e sentido. Por melhor que seja o método de
integração utilizado, sempre é introduzido um erro na posição obtida da partícula com relação
à sua posição ”correta”, da trajetória contínua (figura 13).
O método mais simples de aplicar a aceleração na atualização da velocidade e posição das
partículas é o o esquema de Euler (KELAGER, 2006). Nele:
Rt+4t = Rt +4t~vt
~vt+4t = ~vt +4t~at
onde R é a posição da partícula, ~v a sua velocidade e ~a sua aceleração. Kelager (2006)
comenta que, se a atualização da posição já levar em conta a atualização na velocidade, este
método se torna o método Euleriano Semi-Implícito, mais estável:
41
Figura 13: O erro introduzido pelo processo de integração temporal. A curva cheia representaa trajetória ideal, contínua da partícula A. Integrando a aceleração ~a durante um intervalo detempo t, a nova posição da partícula é B. Com 2t, a partícula pára em C, visivelmente distanteda trajetória correta.
Rt+4t = Rt +4t~vt+4t
Ambos os métodos são extremamente simples de implementar e eficientes, mas dependem
de intervalos de tempo muito pequenos para manterem a estabilidade (HUT; MAKINO, ).
Vários outros métodos foram propostos na literatura, e neste trabalho utiliza-se um método
chamado Leap-Frog. A característica principal do método Leap-Frog é que a posição de uma
partícula é definida em intervalos de tempo ti, ti+1, ti+2, como no método de Euler, mas a ve-
locidade é definida na metade destes intervalos, i.e. ti− 12, ti+ 1
2, ti+ 3
2(HUT; MAKINO, ). O nome
deste método vem do fato que a posição e a velocidade ”pulam” uma sobre a outra (figura 14).
Assim:
Ri = Ri−1 +4t~vi− 12
42
Figura 14: No método Leap-Frog, a velocidade (u) ”pula” por cima da posição (r), e vice-versa(KELAGER, 2006).
~vi+ 12
= ~vi− 12
+4t~ai
Esta maneira de utilizar a aceleração para atualizar a velocidade e posição da partícula
é mais estável, mas requer definir a velocidade em ”meio-intervalos”, que não são diretamente
intuitivos. Hut e Makino () apresentam uma outra maneira de reescrever as equações do método,
agora considerando todas as grandezas em intervalos inteiros:
Ri+1 = Ri + ~vi 4+1
2~ai(4t)2 (4.21)
~vi+1 = ~vi +1
2(~ai + ~ai+1)4 t (4.22)
4.4 Implementação
A classe Particle é a base de todo o sistema. Ela representa uma partícula, e para tal possui
atributos como massa, densidade, pressão, e uma lista de vizinhos. Ela também possui como
atributos sua posição e velocidade e sabe aplicar o método Leap-Frog para a integração tempo-
ral.
A classe System representa o sistema inteiro, com as partículas e o ambiente externo. Ele é
responsável pela inicialização da simulação e pelo seu fluxo. Seus atributos incluem a lista de
partículas, os núcleos de suavização e o recipiente.
43
Finalmente, a classes Poly6Kernel, SpikyKernel e ViscosityKernel representam os núcleos
de suavização. Elas possuem métodos para o cálculo do núcleo, seu gradiente e seu laplaciano
em função de um vetor de distância entre duas partículas.
O algoritmo 1 explicita a principal função do sistema. Nele, deltaT é o intervalo de tempo
entre cada instante da simulação, e P é a lista de partículas do objeto da classe System. A função
busqueVizinhos (linha 3) é discutida no capítulo 6, e a função colidaComObjetos (linha 11) no
capítulo 5.
Algorithm 1 Pseudocódigo da principal função do sistema.1 função atualize(deltaT: Real)2 para cada partícula p em P3 busqueVizinhos(p)4 atualizeDensidadeEPressao(p)5 fim-para6 para cada partícula p em P7 calculeForça(p)8 fim-para9 para cada partícula p em P10 integreComLeapFrog(p, deltaT)11 colidaComObjetos(p)12 fim-para13 para cada partícula p em P14 renderize(p)15 fim-para16 renderizeRecipiente17 fim
As funções atualizeDensidadeEPressao (linha 4) e calculeForça (linha 7) calculam as gran-
dezas do fluido utilizando os núcleos de suavização do método SPH, conforme foi discutido ao
longo deste capítulo. Isto é, a densidade das massas, pressão, forças de pressão, viscosidade,
gravidade e superfície em cada partícula são calculadas utilizando o método SPH sobre a lista de
vizinhos de cada partícula, obtida na função busqueVizinhos. O resultado é um vetor de acele-
ração, ~a. A função integreComLeapFrog (linha 10) utiliza o vetor ~a para atualizar a velocidade
e posição das partículas, implementando as equações 4.21 e 4.22.
Finalmente, a renderização das partículas é feita utilizando chamadas à biblioteca OpenGL
para desenhar uma esfera na posição de cada partícula, e com seu raio (seção 5.2.1 na pá-
gina 49). O recipiente é renderizado como seis quadriláteros com transparência. A figura 15
44
ilustra o diagrama das classes envolvidas nestes estágios da simulação.
Figura 15: O diagrama das classes envolvidas no passo do método SPH.
45
5 Detecção e Tratamento de Colisões
O método SPH do capítulo 4 fornece a maneira para representar e simular fluidos utilizando
partículas, mas não trata dos outros componentes que fazem parte de um sistema de simulação
física, como por exemplo a interação do fluido com o seu ambiente.
Esta interação é essencial: na ausência de objetos no ambiente, o máximo que é possível
simular é um líquido em queda livre, ou um gás subindo na atmosfera. Na maioria das aplica-
ções o interesse é tanto pelo comportamento realista do fluido quanto pela interação do mesmo
com outros objetos. Exemplos incluem frascos contendo líquidos, cubos de gelo boiando em
água, artérias se expandindo devido à pressão do sangue que corre por dentro delas, etc. E isto
significa que o contato entre o fluido e os objetos deve ser devidamente detectado (detecção de
colisão), e a reação correspondente deve ser então tratada.
A proposta deste trabalho é implementar um recipiente retangular que contém uma deter-
minada quantidade de líquido, de forma que as interações entre estes dois tipos de objetos (o
recipiente e o fluido) em particular serão discutidas aqui.
O recipiente pode ser visto como sendo composto por seis planos, se ele for fechado, ou
cinco, caso seja aberto. Se o recipiente for fechado de forma que não haja maneira do líquido
escapar, pode-se utilizar as equações do plano clássico da geometria, que é infinito, para tra-
tar as colisões. Caso contrário, é preciso tomar o cuidado adicional de considerar apenas o
quadrilátero que compõe cada parede. A figura 16 ilustra ambos os casos.
A seção 5.2 discutirá o método para detectar e tratar as colisões entre partículas e planos,
mas para isso é preciso primeiro apresentar definições e propriedades relevantes da Geometria.
46
Figura 16: As situações do recipiente aberto e fechado. Em (A), a partícula não tem comoalcançar as porções dos planos infinitos que não compõem as paredes. Em (B), uma partículapode escapar pelo topo do recipiente e ainda assim colidir com o plano que forma a parede daesquerda.
5.1 Definições de Planos e suas Propriedades
Um plano é totalmente definido por um ponto X0, pertencente ao plano, e um vetor ~n não-
nulo e perpendicular à ele. O conjunto de pontos pertencentes ao plano é aquele que satisfaz a
equação
~n. (X −X0) = 0 (5.1)
Por convenção:
• ~n = [a, b, c]t
• X0 = [x0, y0, z0]t
• X = [x, y, z]t
Aplicando estas convenções à equação 5.1 obtém-se:
ax + by + cz − ax0 − by0 − cz0 = 0
47
A expressão −ax0 − by0 − cz0 é comumente abreviada para uma constante d, resultando
em
ax + by + cz + d = 0 (5.2)
que é a equação geral do plano (WEISSTEIN, 2002).
5.1.1 Distância Perpendicular Ponto-Plano
A medida da distância entre um ponto e um plano é vital para a interação entre as partículas
e o recipiente, por razões que ficarão explícitas na seção 5.2. A distância D entre um ponto
P = [xp, yp, zp] e um plano qualquer é definida pela seguinte equação:
D =axp + byp + czp + d√
a2 + b2 + c2(5.3)
Caso o vetor normal ~n esteja normalizado (ou seja, seu módulo é 1), o denominador da
equação 5.3 pode ser omitido, já que é a equação do módulo do vetor.
Um aspecto interessante é que o valor D obtido na equação 5.3 não é absoluto: ele é
positivo caso o ponto P esteja do lado do plano para qual o vetor normal aponta, e negativo
caso contrário (WEISSTEIN, 2003). A figura 17 ilustra ambos casos.
Figura 17: O ponto A está do lado positivo de plano, enquanto B se encontra no lado negativo.Consequentemente, ”a” é positivo e ”b” é negativo.
48
5.1.2 Intersecção Reta-Plano
Assim como a distância ponto-plano, o cálculo da intersecção entre uma reta e um plano é
essencial para o método de detecção e tratamento de colisão deste trabalho. As razões para tal
serão discutidas na seção 5.2.
Uma reta L pode ser definida por um ponto P0 e um vetor ~u que indica sua direção. Desta
forma, a equação paramétrica
P (s) = P0 + s~u (5.4)
representa todos os pontos pertencentes à reta. Ela indica que todo ponto da reta pode ser
obtido somando-se um vetor paralelo à ~u escalado por um fator s ao ponto conhecido P0. A
seguinte equação dá o valor de s que resulta no ponto PI da intersecção entre a reta e um plano:
sI =−(ax0 + by0 + cz0)
~n.~u(5.5)
PI = P0 + sI~u (5.6)
Vale notar que é preciso verificar antes se a reta L não é paralela ao plano. Se este for o
caso, a reta nunca toca o plano, ou nele está contida completamente. O teste de paralelismo
pode ser efetuado verificando-se se os vetores ~n e ~u são perpendiculares entre si:
~n.~u = 0
Em caso positivo, a reta é paralela ao plano. Para verificar se a reta está então contida no
plano, basta obter a distância entre o ponto L0 e o plano, com a equação 5.3. Se esta for zero, o
ponto está no plano, e consequentemente também a reta.
49
5.2 Intersecções Partículas - Planos
Para finalmente ser possível a discussão da colisão entre uma partícula e um plano, é preciso
definir como uma partícula é representada no espaço. Uma partícula pode ser abstraída como
uma esfera com centro C e raio r. O centro é o ponto utilizado ao longo do método SPH, de
forma que falta apenas obter o valor do raio.
5.2.1 Obtenção do Raio de um Partícula
O raio de uma partícula pode ser obtido através da densidade das massas e a massa da
partícula. O volume V de uma esfera vale
V =4
3πr3
e também
V =m
ρ
sendo m a massa da partícula e ρ a densidade das massas ao redor do centro C. Substituindo
uma equação na outra, obtém-se
4
3πr3 =
m
ρ
r3 =3m
4πρ
r = 3
√3m
4πρ(5.7)
que é o raio da partícula.
50
5.2.2 Intersecção sem trajetória
A maneira mais simples e intuitiva de verificar se uma partícula está em colisão com um
plano é obtendo a distância d entre o centro C da partícula e o plano utilizando a equação 5.3.
Sendo |d| menor do que o raio, a partícula está atravessando o plano, como ilustra a figura 18.
Figura 18: No caso mais simples, se a distância |d| do centro C ao plano for menor do que oraio r, então a partícula está colidindo com a superfície.
Esta maneira de verificar a colisão é simples e rápida, mas insuficiente. Dados dois instantes
de tempo consecutivos t0 e t1, é completamente possível que a partícula atravesse totalmente
o plano entre os instantes, mas nas posições inicial e final não o intersecte. A figura 19 ilustra
esta situação.
Figura 19: Uma partícula cruzando um plano ao longo de sua trajetória em um intervalo detempo.
51
Uma solução possível neste caso é diminuir o tamanho do intervalo entre os instantes de
simulação, forçando a partícula a colidir com o plano em um instante observável. Como o
tempo de processamento não será proporcionalmente diminuído, o resultado é uma simulação
mais lenta.
Outra solução é considerar a trajetória da partícula entre os instantes, e usar também esta
informação no teste de colisão. Esta é a solução adotada neste trabalho, e será discutida a seguir.
5.2.3 Intersecção com trajetória
Os pontos inicial e final de uma partícula durante um intervalo de tempo podem ser utiliza-
dos para um teste (e tratamento) de colisões mais realista e robusto. Para isto, convém adotar
algumas convenções:
• r é o raio da partícula.
• C0 é o centro da partícula no início de sua trajetória, ou seja, no instante t0.
• C1 é o centro da partícula no final da sua trajetória, no instante t1. Esta é a posição obtida
após todo o cálculo do método SPH e subsequente integração no tempo, sem o teste de
colisão ainda.
• Cf é o centro da partícula no final da sua trajetória, mas após o teste e tratamento de
colisão. Esta é a posição final correta e desejada. Se não houve colisão, Cf = C1.
Assumem-se também duas premissas:
1. No instante inicial t0, a partícula sempre está totalmente dentro do recipiente que a con-
tém, ou seja, a distância entre C0 e qualquer plano do recipiente é sempre maior ou igual
a r.
2. O lado de dentro do recipiente é o lado positivo de todos os planos que o formam (segundo
a definição da subseção 5.1.1 na página 47), e os vetores normais destes planos estão
sempre normalizados.
52
O primeiro passo é verificar o valor da distância d entre o centro C1 da partícula e um dado
plano, usando a equação 5.3 na página 47. O resultado obtido cai em um de três possíveis
casos, ilustrados na figura 20:
• d é positiva, e maior ou igual a r. Neste caso a partícula, em sua posição final, não
intercepta o plano (ou no máximo o toca, no caso em que d = r ). Como na origem
da trajetória a partícula está totalmente do lado positivo do plano (premissa 2) e no final
também, e como a trajetória é uma linha reta, conclui-se que a partícula não colidiu com
o plano durante este intervalo de tempo. Nada mais precisa ser feito para este plano e esta
partícula.
• d é positiva, mas menor do que r. Neste caso, a partícula, em sua posição final, está
intersectando o plano, mas sua maior parte ainda está dentro da caixa, do lado positivo do
plano. Houve colisão.
• d é negativa. Agora a partícula, em sua posição final, está toda ou parcialmente fora da
caixa, e houve colisão.
Na figura 20, no primeiro caso não houve colisão, e Cf = C1. Nos dois últimos houve
colisão, e o tratamento em ambos os casos é o mesmo.
O próximo passo é colocar a partícula em sua posição final correta, simulando a reflexão da
mesma na parede da caixa. Para isso é preciso obter um vetor ~u que desloca a partícula de C1 a
Cf . Intuitivamente, este vetor tem a mesma direção e sentido que o vetor normal do plano, ~n,
como mostra a figura 21. Basta então encontrar o seu módulo. Esta perspectiva do tratamento
da colisão é análoga ao método híbrido de Impulso-Projeção de Kelager (2006).
O primeiro ponto da partícula que, durante a sua trajetória, tocou o plano é aquele que está
mais próximo do mesmo na posição inicial. Este é o ponto P0 da figura 22.
P1 é o ponto final da trajetória que parte de P0. A distância a de P1 ao plano é r − d, onde
d é a distância do centro da partícula ao plano (equação 5.3) e r é o raio da mesma. A figura 23
ilustra estas grandezas.
53
Figura 20: Os três casos possíveis para o valor da distância d de C1 ao plano. Em (A), o fato ded ser maior que o raio da partícula evidencia a não colisão. Em (B), d é positiva e menor quer . Em (C), d é negativa pois C1 está do lado negativo do plano. Nos últimos dois casos houvecolisão.
54
Figura 21: O vetor ~u liga a posição final atual e a correta da partícula. Na posição inicial C0, aesfera já toca o plano.
Figura 22: A trajetória de três dos pontos da partícula. Como todos viajam à mesma velocidadee o ponto P0 é o que tem a menor distância (x) a percorrer até chegar ao plano, ele é o que entraem contato primeiro.
55
Figura 23: A distância a de P1 ao plano vale r + |d|. Como d é garantidamente negativa no casoda figura, o valor desta distância pode ser escrito como r − d.
O valor de a sempre pode ser obtido por r − d, independentemente da posição final atual
da partícula ser dentro ou fora da caixa. Para provar, vale lembrar que o valor da distância de
um ponto ao plano da seção 5.1.1 é sinalizado.
Se a partícula estiver no lado positivo do plano, d será positivo e a expressão r − d repre-
sentará exatamente a porção da partícula que está do lado de fora dele (figura 24). Estando a
partícula no lado negativo, d será negativo e r − d será r + |d|, ou seja, o raio da partícula mais
o valor absoluto da distância do plano ao centro (figura 23).
Figura 24: Mesmo quando C1 se encontra do lado positivo do plano, a distância de P1 à eleainda vale r − d.
Em seu comportamento correto e no caso de uma colisão totalmente elástica (ou seja, uma
colisão aonde não haja perda de energia e consequente desaceleração), o ponto P1 ainda estaria
a uma distância a do plano, mas em seu outro lado. Assim, o valor procurado do módulo do
vetor ~u é 2a:
56
a = r − d
~u = ~n.2a (5.8)
Somando-se ~u à C1, obtém-se Cf . A figura 25 ilustra estes valores.
Figura 25: Cf = C1 + ~u, onde o tamanho de ~u vale duas vezes o valor da distância do pontomais distante ao plano.
Outra maneira de entender porquê ~u = ~n.2a é pensar neste deslocamento por ~u como uma
reflexão ao redor do plano, como se este fosse um espelho e Cf fosse o reflexo de C1. A questão
é que o eixo ao redor do qual é feita a reflexão não está contido no plano, mas sim acima dele.
Mais especificamente, ele passa pelo ponto CI , que é o centro da partícula no exato instante da
colisão (figura 26). É daquele momento que começa a reflexão, de modo que toda a distância
que a partícula percorreria sem a presença do plano deve ainda ser percorrida, na direção da
reflexão.
57
Figura 26: O eixo E a partir do qual começa a reflexão do centro passa pelo ponto CI .
Assim, tal qual em um reflexo de espelho, a distância a entre C1 e o eixo E da figura 26
deve ser a mesma que a distância deste mesmo eixo à Cf . Este valor é a distância de C1 ao plano
mais a distância do plano ao eixo. Como o eixo é paralelo ao plano e passa por CI , esta segunda
distância é igual ao raio da partícula! Assim, a distância de C1 ao eixo vale r − d sendo d a
distância negativa de C1 ao plano. Este valor é idêntico ao obtido utilizando P1 como referência
(figura 23), de modo que |~u| realmente vale 2(r − d) (equação 5.8). A figura 27 ilustra todos
estes valores.
Figura 27: Utilizando o eixo E como eixo de reflexão, o comprimento do vetor ~u ainda é r− d.
Vale notar que é possível alterar o tamanho de ~u para simular perdas de energia na colisão.
Basta alterar o coeficiente que multiplica a para qualquer valor entre 1 e 2. 1 simula uma
colisão totalmente inelástica (como um punhado de lama arremessado contra uma parede), e
58
2 uma colisão totalmente elástica. Na prática valores como 1,5 e 1,6 foram observados como
úteis.
Finalmente, quando a partícula colidiu com o plano, a direção de sua velocidade se alterou,
de forma que o último passo é atualizar a velocidade da partícula.
O vetor ~v de velocidade da partícula pode ser decomposto em outros dois: Um vetor ~vn
paralelo ao vetor normal do plano (normalizado), e outro vetor ~vp perpendicular ao mesmo
(MOLLER ERIC HAINES, 2002):
~vn = ~n(~v.~n)
~vp = ~v − ~vn
Assim, o método para refletir a velocidade ao redor do plano consiste em inverter ~vn e
somá-lo novamente a ~vp:
~vf = ~vp − ~vn
A figura 28 ilustra a idéia.
Figura 28: ~v2 é o resultado da reflexão de ~v1 no plano com normal ~n.
59
5.2.4 Interações entre dois ou mais planos
O método apresentado na subseção anterior é suficiente na presença de apenas um plano.
Quando existem vários planos (como as seis paredes da caixa deste trabalho), a escolha arbitrá-
ria da ordem dos planos para o teste pode produzir resultados absolutamente errados, conforme
a figura 29 exemplifica. A figura 30 mostra o resultado correto esperado.
Figura 29: O que pode acontecer quando dois ou mais planos existem no ambiente. Por causada escolha do plano B para o teste, a posição final Cf não representa o que aconteceria narealidade.
Figura 30: O resultado esperado da situação da figura 29. A partícula deve primeiro colidir como plano A, para então colidir com o plano B.
Assim, a escolha cuidadosa de qual plano será testado primeiro é essencial. O critério para
seleção do primeiro plano a ser verificado é simples: O plano que possuir ponto de contato com
60
a partícula (em qualquer ponto de sua trajetória) mais próximo da posição inicial dela é o ideal.
Traduzindo, a parede na qual a partícula ”bateria” primeiro é a candidata correta.
Como este plano pode ser obtido?
Na seção 5.2.3 foi convencionado que o ponto P0 é o primeiro ponto a entrar em contato
com o plano. A trajetória deste ponto ao longo do intervalo de tempo é um segmento de reta
paralelo à trajetória do centro da partícula, e de mesmo tamanho.
Assim, os dois pontos que formam esta nova trajetória são obtidos diretamente somando-se
um vetor ~r aos dois pontos da trajetória do centro da partícula. Intuitivamente, ~r é paralelo ao
vetor ~n do plano, com sentido contrário e módulo igual ao raio da partícula:
~r = −r~n
Com este vetor e os dois pontos da trajetória do centro, obtém-se agora a trajetória do ponto
P0 ao longo do intervalo de tempo. As equações 5.5 e 5.6 na página 48 dão o ponto de contato
PI entre esta trajetória (um segmento de reta) e o plano. Este é o ponto onde a partícula tocou
pela primeira vez o plano.
Essas duas equações, a 5.5 e a 5.6 na página 48, são para retas infinitas mas podem ser
utilizadas com as trajetórias (segmentos de retas), se for tomado um cuidado em potencial: O
ponto PI de intersecção da trajetória só deve ser calculado se, e somente se, a distância d entre
C1 e o plano indicar uma colisão (vide seção 5.2.3). Se d indicar que não houve colisão, não
é necessário calcular a trajetória de P0 e obter PI , já que a partícula não colidiu com o plano
em momento algum. O resultado disso é que PI só deve ser obtido para aqueles planos que
potencialmente colidem com a partícula.
Tomando este cuidado em especial, tem-se a garantia de que o ponto PI resultante das duas
equações de intersecção estará contido na trajetória de P0. Se não estivesse, a partícula não
colidiria com o plano, e logo este nem deveria ter sido considerado. A figura 31 ilustra esta
justificativa.
61
Figura 31: A distância d de C1 ao plano indica se PI deve ser obtido, pela presença ou não decolisão.
62
O módulo do vetor PI−P0 é a distância entre o ponto de contato e o ponto mais próximo do
plano no instante inicial. Como as duas trajetórias (a do centro e a do ponto P0) são paralelas,
este valor é também a distância entre o centro das partículas nos dois instantes. A figura 32
ilustra esta idéia.
Figura 32: A distância entre a posição inicial da partícula e sua posição no instante da colisãovale |I − P0|.
O método consiste em obter o valor desta distância para todos os planos, tomar o plano com
a menor distância, e aplicar o procedimento da seção 5.2.3 na página 51. Em seguida, repetir
até que não existam mais colisões. Mas antes de repetir é preciso tomar uma última medida:
atualizar os pontos da trajetória da partícula.
Após a colisão e consequente reflexão a trajetória da partícula ao longo do intervalo não é
mais a mesma do início do processo. O ponto inicial da nova trajetória é o ponto do centro da
partícula no momento de primeiro contato, e o ponto final é o ponto obtido na seção anterior,
após o tratamento da colisão.
Obter o novo ponto inicial é simples: Ele está contido na trajetória antiga, a uma distância
|PI − P0| do centro inicial (figura 33):
C0novo = C0antigo + (PI − P0)
O fato deste procedimento ter que ser repetido até que não existam mais colisões implica
em um custo computacional razoável. Para poucos planos (como as 6 paredes do recipiente),
este custo é pequeno em relação ao custo maior de busca por vizinhança e o método SPH, mas
63
Figura 33: C0′ é o novo início da trajetória da partícula, e pode ser obtido por C0 + (PI − P0).
conforme o número de planos cresce este método tem o potencial de se tornar um dos gargalos
da performance da simulação.
5.3 Balanço do recipiente
O método da seção 5.2 é suficiente para o caso mais comum em uma simulação, que é
um ambiente com objetos estacionários. O tratamento para os casos nos quais os obstáculos do
fluido também podem se deslocar é diferente; é preciso considerar os obstáculos como entidades
físicas também, com massa, velocidade e etc.
Quando dois objetos se chocam ocorre uma troca instantânea de energia entre eles. Esta
troca se traduz na mudança na velocidade dos objetos. O modelo de colisão para esta situação
utilizado neste trabalho é chamado ”Lei de Newton de Restituição para Colisões Instantâneas
sem Fricção” (HECKER, 1997).
O primeiro passo é avaliar se a colisão entre dois objetos realmente ocorreu. Neste trabalho,
isto significa verificar a colisão entre os planos do recipiente e as partículas do fluido. Como
neste caso não apenas as partículas, mas também os planos, estão se movendo, não é válido
utilizar o modelo de colisão com trajetória da seção 5.2.3 puro. Se o deslocamento de cada plano
do recipiente for pequeno o suficiente para que nenhuma partícula chegue a sair totalmente do
recipiente, é suficiente utilizar o modelo de colisão sem trajetória da seção 5.2.2 na página 50.
Na prática foi observado que este modelo realmente é suficiente.
64
Havendo uma colisão, o segundo passo é determinar a velocidade relativa dos objetos em
questão (no caso, a velocidade relativa entre uma partícula e um dos planos) e o vetor normal
da colisão (o vetor normal do plano, ~n). A velocidade relativa entre dois objetos ”a” e ”b” pode
ser obtida com a seguinte fórmula:
~vab = ~va − ~vb
Ou seja, a velocidade relativa é a velocidade de um objeto em relação ao outro. É a veloci-
dade na qual ”b” vê ”a” se afastar ou se aproximar dele. Outras grandezas novas são o impulso
(j) e o coeficiente de restituição (e) , que aparecem na seguinte equação:
j =−(1 + e)~vab.~n
~n.~n( 1ma
+ 1mb
)
Onde m é a massa de ”a” ou ”b”. O coeficiente de restituição é um escalar que representa
a quantidade de energia dissipada na colisão. Semelhantemente à colisão totalmente elástica da
seção 5.2.3 na página 51, ele pode variar de 0 (um punhado de lama que atinge o solo e nele
fica) a 1 (uma esfera perfeitamente reflexiva).
O valor j é o impulso que cada objeto sofre como consequência da colisão. Este impulso
é diretamente na direção do vetor normal da colisão (~n), com sentidos opostos para ”a” e ”b”,
representando a nova velocidade dos objetos:
~v2a = ~v1
a +j
ma
~n
~v2b = ~v1
b −j
mb
~n
O processo para obtenção destas equações é demonstrado por Hecker (1997). Assim, cada
colisão detectada entre as partículas e os planos pode ser tratata atualizando as velocidades de
ambos e expulsando a partícula da superfície do plano imediatamente ( 5.2.2 na página 50).
65
5.4 Implementação
A detecção e o tratamento de colisões são feitos na função colidaComObjetos, na linha 11
do algoritmo 1 na página 43. Esta função é a implementação direta dos métodos da seções 5.2.3
na página 51 e 5.2.4 na página 59. Para tal, uma classe extra chamada Box foi acrescentada ao
sistema. Esta classe sabe detectar e tratar as colisões com partículas corretamente, levando em
conta pontos como a ordem de colisão.
Esta classe é ainda responsável pela implementação do balanço do recipiente, utilizando
diretamente as equações da seção 5.3. O evento ”balanço da caixa” é disparado por uma tecla,
que faz o recipiente transladar em uma direção específica. A classe System utiliza este vetor de
translação e simula o deslocamento em um curto intervalo de tempo. Com o deslocamento e o
tempo, obtém-se a velocidade do recipiente:
~v =~r
t
Onde ~r é o vetor de translação e t é o intervalo de tempo. O objeto Box do sistema utiliza
estas informações e a lista de partículas do sistema para implementar o impulso e consequente
troca de velocidade das partículas que colidem com o recipiente. Esta implementação é sim-
plista mas na prática funciona para deslocamentos pequenos da caixa.
Com o acréscimo da classe Box, o diagrama de classes do sistema é ilustrado na figura 34
Figura 34: A classe Box é a responsável pela detecção e o tratamento das colisões.
66
6 Busca por Vizinhos
Os núcleos de suavização (seção 4.1 na página 27) são o coração do método Smoothed
Particle Hydrodynamics. Como mencionado, existe um valor chamado raio de suporte que
limita a distância em que uma partícula se encontra de outra para que existam contribuições
entre elas. A seção 4.1 explicou um núcleo de suavização como uma esfera ao redor de uma
partícula, e o raio de suporte como o raio desta esfera hipotética. Assim, apenas partículas
que se encontram dentro desta esfera hipotética possuem contribuições não nulas. Encontrar as
partículas na vizinhança de uma partícula p é o primeiro passo do método SPH.
A seguir serão apresentados e discutidos dois métodos para a busca de vizinhos, salientando-
se as vantagens e desvantagens de cada um.
6.1 Força Bruta
As partículas existem na memória de um computador como estruturas de dados que contém
a posição da partícula, sua densidade, e todas as suas outras propriedades. As maneiras mais
simples de organizar estas estruturas na memória são em um vetor, onde cada partícula conse-
cutiva ocupa uma região também consecutiva na memória, ou em uma lista, onde cada estrutura
representativa da partícula possui o endereço da região que a estrutura seguinte ocupa. A figura
35 mostra estas duas estruturas de dados.
Nenhuma destas estruturas mantém a relação espacial entre as partículas explícita na me-
mória do computador. Duas partículas que ocupam posições consecutivas em um vetor podem
estar em extremidades opostas do ambiente da simulação.
67
Figura 35: Na esquerda, as estruturas que representam as partículas ocupam regiões consecuti-vas em um vetor. Na direita, cada estrutura sabe o endereço da próxima em uma lista.
A maneira mais simples de encontrar as partículas (na memória) que estão em uma vizi-
nhança (no ambiente simulado) é verificar o valor da posição de cada uma delas em relação às
outras. Para cada partícula p , percorre-se inteiramente a estrutura de dados que armazena as
partículas procurando partículas vizinhas a p. Esta maneira de encontrar a vizinhança é chamada
de Força Bruta por percorrer o vetor ou lista cegamente várias vezes, sem aproveitar nenhuma
propriedade que o problema possa possuir.
6.1.1 Otimização
O método Força Bruta original pode ser facilmente acelerado se for levado em consideração
o fato de que, se a partícula A é vizinha de B , então B também é vizinha de A. Assim, quando
se identifica que uma partícula A é vizinha de B , a informação de que B é uma das vizinhas
de A pode também ser registrada garantidamente.
No nível das estruturas armazenadoras das partículas, isto significa que o vetor ou lista não
mais precisa ser percorrido inteiramente para cada partícula, apenas da região de memória que
ela habita em diante, como explica a figura 36.
68
Figura 36: Em um determinado momento (A), está sendo efetuada a busca por vizinhos dapartícula A. Testa-se esta relação entre A e B. Quando, em (B), for feita a busca por vizinhosde B, não é necessário testar com A, pois a relação de vizinhança é simétrica e este teste já foiefetuado, sendo elas vizinhas ou não.
6.1.2 Vantagens
• O método Força Bruta é o mais fácil de todos de compreender e programar. Mesmo com
a otimização a implementação é trivial e requer poucos recursos do computador.
• O tempo total de busca por vizinhos é longo (em relação a outros métodos), mas sempre
constante.
6.1.3 Desvantagens
• Tempo total de busca imenso, em comparação à outros métodos. O tempo de processa-
mento do método Força Bruta inibe a interatividade da simulação com poucos milhares de
partículas. A otimização reduz este tempo consideravelmente, mas ainda assim o tempo
obtido é considerado insatisfatório.
6.2 Busca em Grid
Um Grid é uma matriz, tridimensional no caso de uma simulação em três dimensões. Ele
é uma maneira de se dividir o espaço da simulação em espaços menores, de mesmo tamanho.
Se um recipiente contendo fluido for visto como um Grid, cada partícula que compõe o fluido
estará contída em uma das células da matriz.
Indo um passo além, é possível aproveitar as propriedades do método SPH e dividir um
69
recipiente em um Grid onde cada célula tem o lado exatamente igual ao raio de suporte dos
núcleos de suavização. Assim, as partículas potencialmente vizinhas à uma partícula em uma
célula X estarão necessariamente dentro da célula X ou nas células diretamente adjacentes a
ela (MULLER; CHARYPAR; GROSS, 2003). A figura 37 ilustra esta situação.
Figura 37: Com um Grid, todas as partículas potencialmente vizinhas à partícula A se encon-tram na célula de A e nas células diretamente adjacentes.
Um detalhe importante é que apenas a posição do centro de uma partícula importa para
determinar a qual célula ela pertence. A justificativa é que todos os núcleos de suavização do
método SPH utilizam apenas as posições dos centros nos cálculos, de forma que o raio (se-
ção 5.2.1 na página 49) não influi. Ele está implicitamente contabilizado no valor da densidade
das massas das partículas e suas forças de pressão resultantes.
O método por Grid foi utilizado neste trabalho para acelerar a busca por vizinhos, de forma
que é válido discutí-lo com mais detalhes. As subseções a seguir assumem que o ponto de
origem do Grid está na origem do universo do sistema, e que o Grid é alinhado aos eixos x ,
y e z . O caso mais genérico aonde o Grid se encontra em qualquer posição será tratado na
seção 6.2.3 na página 72.
70
6.2.1 Inserção e Remoção no Grid
Para inserir uma partícula com centro C = (xp, yp, zp) em um Grid é necessário encontrar
os índices da célula que a conterá. Estes índices podem ser obtidos diretamente através de C e
do tamanho dos lados das células.
Seja G um Grid tridimensional onde cada célula possui hx de base, hy de altura e hz de
profundidade. A célula onde uma partícula com centro C = (xp, yp, zp) se encontra possui
índices i, j e k que podem ser obtidos com as seguintes equações:
i =⌊xp
hx
⌋(6.1)
j =
⌊yp
hy
⌋
k =⌊zp
hz
⌋
Onde⌊
ab
⌋é o valor da divisão de a por b arredondado para baixo.
Como exemplo, seja o Grid bidimensional da figura 38 na página seguinte. A partícula da
figura possui centro C = (9, 3), de forma que a célula que ela ocupará possui os índices
i =⌊9
2
⌋= 4
j =⌊3
2
⌋= 1
O método para busca de vizinhos em Grid consiste então em obter o índice da célula que
contém a partícula de interesse para então obter o conjunto de todas as partículas que habitam
esta célula e as adjacentes.
71
Figura 38: Uma partícula com centro C = (9, 3) fica na posição (4, 1) do Grid da imagem.
6.2.2 Otimização
O método de busca em Grid pode ser acelerado à uma maneira semelhante do método de
busca por Força Bruta. As partículas candidatas à vizinhança de uma partícula qualquer P ,
obtidas pelas células de P e adjacentes a ela, são também as candidatas à vizinhança de todas
as partículas na célula de P , como exemplifica a figura 39.
Figura 39: A região hachurada em azul contém todas as partículas candidatas à vizinhança deA e B.
Assim, quando se testa a vizinhança de P pode-se imediatamente testar também a vizi-
72
nhança de todas as outras partículas da célula de P . A abordagem passa de considerar cada
partícula individualmente para considerar cada célula que contém partículas individualmente.
Este ”truque” economiza novas buscas no Grid e melhora a performance.
6.2.3 Grids Fora da Origem
É absolutamente possível que a origem do Grid não se encontre na origem do universo
(como o ponto (0, 0, 0) do plano cartesiano). Neste trabalho, isto acontece frequentemente com
o balanço do recipiente ( 5.3 na página 63), que se desloca livremente pelo espaço. O Grid deve
se deslocar junto com o recipiente para continuar representando todo seu interior.
Por sorte, o tratamento para deslocar o Grid é simples. Basta manter controle do ponto
de origem do Grid (por exemplo o ponto (0, 0, 0) no início da simulação). Seja este ponto o
ponto A. Quando o recipiente é deslocado por um vetor de translação ~x, aplica-se também esta
translação a A. No momento de inserção e remoção leva-se em conta então também o ponto de
origem do Grid nas equações 6.1 na página 70:
i =⌊xp − xA
hx
⌋
j =
⌊yp − yA
hy
⌋
k =⌊zp − zA
hz
⌋
Como exemplo, seja novamente o Grid da figura 38 na página anterior, mas deslocado da
origem de forma que a origem do recipiente está no ponto A = (3, 2) (figura 40 na próxima
página). A partícula da figura possui centro C = (12, 5), de forma que a célula que ela ocupará
possui os índices
i =⌊12− 3
2
⌋= 4
73
j =⌊5− 2
2
⌋= 1
Figura 40: Uma partícula com centro C = (12, 5) fica na célula (4, 1) em um Grid deslocado(3, 2) da origem.
6.2.4 Vantagens
• Busca por Grids é muito mais rápida do que a busca por Força Bruta. Em testes foram
observadas melhorias de até 7 vezes no tempo de busca total.
• Apesar de mais difícil de implementar do que a Força Bruta, o Grid ainda é relativamente
simples de modelar e implementar.
6.2.5 Desvantagens
• A grande desvantagem do Grid é seu consumo de memória. O espaço inteiro da simu-
lação é dividido em células de tamanho igual (e geralmente pequeno). Mesmo com uma
quantidade razoável de fluido a maioria das células sempre vai estar vazia, ocupando
espaço desnecessário.
• O tempo total de busca de vizinhos geralmente não é constante ao longo da simulação. Ele
depende da quantidade de partículas nas células. Células com mais partículas representam
74
um conjunto maior de candidatas a vizinhança.
6.3 Implementação
A classe Grid é fundamental para a busca eficiente de vizinhos na simulação. No algo-
ritmo do ciclo principal de simulação ( 1 na página 43), esta classe é responsável pela função
busqueVizinhos (linha 3).
No início da simulação, o objeto da classe System cria um objeto Grid com as dimensões
do recipiente, e de forma que cada célula tenha como lado o tamanho do raio de suporte dos
núcleos de suavização. Em seguida, ele insere todas as partículas da lista no Grid.
Na função busqueVizinhos o Grid primeiramente obtém o índice de uma partícula p em
questão e retorna todas as partículas contidas naquela célula e em suas vizinhas. O System
então verifica a relação de vizinhança correta naquele grupo, isto é, apenas partículas a uma
distância menor ou igual ao raio de suporte dos núcleos de suavização são realmente vizinhas
de p. Como otimização, o System aproveita o grupo de partículas obtido pelo Grid e testa todas
as partículas da célula de p.
75
O diagrama de classes final, acrescido da classe Grid é ilustrado na figura 41.
Figura 41: A classe Grid possibilita a busca rápida por vizinhança.
76
7 Resultados e Conclusões
Neste capítulo serão apresentados e discutidos os resultados da implementação do sistema
de partículas proposto. Os valores utilizados na simulação da água são demonstrados na tabela
1. Os valores físicos de referência são os utilizados por Kelager (2006). O tempo de cada passo
da simulação (um frame) é de 0,01 segundos.
Parâmetro Unidade ValorNo de Partículas - 1000
Volume de Líquido m3 0,02Massa kg 0,02
Raio de Suporte (h) m 0,0415Partículas na Vizinhança - 20
Densidade de Repouso (ρ0)kgm3 998,29
Viscosidade (µ) Pa.s 3,5Tensão de Superfície (σ) N
m0,0728
Limite (l) - 7,065Constante k J 3
Tabela 1: Valores utilizados na simulação.
A figura 42 ilustra algumas cenas durante a simulação. Utilizando o método de busca de
vizinhança mais eficiente (Grid com otimização, seção 6.2 na página 68), a simulação ocorre
a uma média de 25 quadros por segundo. Como o intervalo de tempo utilizado na simulação é
de 0,01 segundos (acima disto, o método LeapFrog perde sua estabilidade, como foi observado
experimentalmente), e cada quadro leva 125
= 0, 04 segundos para ser processado, o resultado é
que o fluido parece estar se movendo a uma velocidade quatro vezes menor que a normal. Isto é
longe de ser uma simulação em tempo real, de forma que otimizações consideráveis de métodos
e código teriam que ser feitas antes de se obter uma simulação totalmente realista.
77
Figura 42: Cenas de seis pontos diferentes da simulação.
78
7.1 Busca de Vizinhança
A tabela 2 compara todos os métodos de busca de vizinhança discutidos no capítulo 6.
No deIterações
Com Grid,Otim.
Com Grid,Não Otim.
Força Bruta,Otim.
Força Bruta,Não Otim.
10 20 31 84 16720 25 39 88 16630 23 36 83 16040 19 23 80 15950 20 32 85 16060 24 29 80 16170 26 31 82 15780 22 32 79 16190 20 31 83 162100 22 30 82 164
Média 22,1 31,4 82,6 161,7
Tabela 2: Comparação entre os diferentes métodos de busca de vizinhança. Os valores na colunaindicam o número de ciclos que cada método leva em média no processador, dividido por mil.
Como a tabela demonstra, o método com Grid e otimização é em média quase oito vezes
mais eficiente que o pior de todos, o método Força Bruta sem otimização. De fato, a utiliza-
ção do método Força Bruta na prática inibiu completamente qualquer sensação de realismo na
simulação do fluido. Contudo, ele foi utilizado como base para teste e validação da busca por
Grid, por ser de fácil implementação.
A tabela 3 ilustra a diferença entre os dois extremos pela porcentagem que cada um ocupa
do ciclo total de simulação.
Passo Força Bruta GridBusca de Vizinhos 80% 20%
SPH 15% 70%Integração no Tempo < 1% < 1%
Tratamento de Colisões 2% 5%Renderização 2% 5%
Tabela 3: Porcentagem de processamento com e sem busca por Grid.
79
7.2 Conclusões
O método Smoothed Particle Hydrodynamics é um dos métodos possíveis para a simulação
de fluidos, mas não é o único componente de uma simulação física. Este trabalho tratou a
implementação do método SPH para a água e também apresentou métodos para a interação do
fluido com seu ambiente (capítulo 5 na página 45) e para a busca rápida de vizinhos (capítulo 6
na página 66).
O método SPH, em teoria, simula líquidos corretamente. A observação feita é que o método
é pesado e possui um custo computacional alto, de forma que não foi possível simular líquidos
em tempo real.
A interação do fluido com seu ambiente é essencial, e o método para detecção e resposta à
colisões com o recipiente deste trabalho foi avaliado como robusto e suficiente, não implicando
em perdas de performance proibitivamente grandes. Esta observação leva em conta o fato de
que o número de objetos no ambiente é pequeno (seis paredes).
A busca acelerada de vizinhos é necessária para a simulação interativa, ainda que não sufi-
ciente. O método por Grid apresentado neste trabalho é de simples implementação e apresenta
uma melhoria no tempo de busca por vizinhos de quase uma ordem de grandeza, se comparado
ao método mais ineficiente.
7.3 Trabalhos Futuros
Como mencionado, não foi possível a simulação de fluido com 1000 partículas em tempo
real. Medidas que poderiam ser tomadas para acelerar a simulação incluem:
• Utilização de hardware mais poderoso.
• Busca de métodos mais eficientes para ”gargalos” da simulação. Por exemplo, um mé-
todo mais eficiente que o Grid poderia ser concebido para a busca de vizinhança, ou o
procedimento para detecção de colisão poderia ser mais inteligente. A propriedade da
80
simetria da vizinhança entre partículas poderia ser explorada para otimizar o cálculo das
contribuições do método SPH.
• Otimizações em nível de código e hardware. Características da plataforma utilizada
(compilador, sistema operacional, processador) poderiam ser exploradas para a geração
de código mais eficiente, sacrificando portabilidade.
Além de otimizações, vários capacidades adicionais poderiam ser adicionadas à simulação,
como por exemplo:
• Algoritmo de colisões: Generalizar o algoritmo para outros tipos de objetos, diferentes
das seis paredes do recipientes. Permitir que o recipiente seja aberto.
• Algoritmo de colisões: Incluir propriedades físicas do domínio da aplicação. Por exem-
plo, levar em conta a força da gravidade na reflexão de uma partícula (figura 43)
• Extrair a superfície do fluido antes de renderizá-lo, para um efeito visual mais realista.
Possibilidades para isto incluem métodos como Marching Cubes (LORENSEN; CLINE,
1987) e Surface Splatting (MULLER; CHARYPAR; GROSS, 2003).
• Interações do fluido com sólidos deformáveis. Isto permitiria por exemplo a simulação
de sangue correndo através de uma artéria. Este problema seria interessante mesmo que
a artéria não se expandisse com o fluxo sanguíneo, pelo simples acréscimo do número de
polígonos para colisão (figura 44).
81
Figura 43: A gravidade atua sempre nas partículas, inclusive nas respostas de colisões.
Figura 44: Uma artéria modelada como uma malha de triângulos pode possuir dezenas demilhares de polígonos.
82
Referências
AMADA, T.; IMURA, M.; YASUMURO, Y.; MANABE, Y.; CHIHARA, K. Particle-basedfluid simulation on gpu. ACM Workshop on General-Purpose Computing on GraphicsProcessors and SIGGRAPH 2004 Poster Session, 2004.
ATKINS, P. W. Fisico-Quimica. 6. ed. [S.l.]: LTC, 1999.
BRIDSON, Robert; MULLER-FISCHER, Matthias; GUENDELMAN, Eran. Fluid simulation.In: SIGGRAPH 2006 Course Notes. [S.l.: s.n.], 2006.
DESBRUN, Mathieu; GASCUEL, Marie-Paule. Smoothed particles: a new paradigm foranimating highly deformable bodies. In: Proceedings of the Eurographics workshop onComputer animation and simulation ’96. New York, NY, USA: Springer-Verlag New York,Inc., 1996. p. 61–76. ISBN 3-211-82885-0.
FOSTER, Nick; FEDKIW, Ronald. Practical animation of liquids. In: SIGGRAPH ’01:Proceedings of the 28th annual conference on Computer graphics and interactive techniques.New York, NY, USA: ACM Press, 2001. p. 23–30. ISBN 1-58113-374-X.
FOSTER, Nick; METAXAS, Dimitri. Realistic animation of liquids. Graphical mo-dels and image processing: GMIP, v. 58, n. 5, p. 471–483, 1996. Disponível em:<citeseer.ist.psu.edu/foster95realistic.html>.
GINGOLD, R.A.; MONAGHAN, J.J. Smoothed particle hydrodynamics - theory andapplication to non-spherical stars. mnras, v. 181, p. 375–389, nov. 1977.
HARLOW, F. H.; WELCH, J. E. Numerical calculation of time-dependent viscousincompressible flow of fluid with free surface. In: The Physics of Fluids. [S.l.: s.n.], 1965. p.2182–2189.
HECKER, Chris. Collision response. Game Developer Magazine, June 1997.
HUT, Piet; MAKINO, Jun. The art of computer science.
JAILLET, F.; SHARIAT, B.; VANDORPE, D. Deformable volume object mode-ling with a particle-based system for medical applications. 1997. Disponível em:<citeseer.ist.psu.edu/jaillet97deformable.html>.
KASS, Michael; MILLER, Gavin. Rapid, stable fluid dynamics for computer graphics.In: SIGGRAPH ’90: Proceedings of the 17th annual conference on Computer graphicsand interactive techniques. New York, NY, USA: ACM Press, 1990. p. 49–57. ISBN0-201-50933-4.
KELAGER, Micky. Lagrangian fluid dynamics using smoothed particle hydrodynamics. What.January 2006. Disponível em: <www.opentissue.org>.
83
KOSHIZUKA, Seiichi; NOBE, Atsushi; OKA, Yoshiaki. Numerical analysis of breakingwaves using the moving particle semi-implicit method. International Journal for NumericalMethods in Fluids, v. 26, p. 751–769, 1998.
LORENSEN, William E.; CLINE, Harvey E. Marching cubes: A high resolution3d surface construction algorithm. In: SIGGRAPH ’87: Proceedings of the 14thannual conference on Computer graphics and interactive techniques. New York, NY,USA: ACM Press, 1987. v. 21, n. 4, p. 163–169. ISSN 0097-8930. Disponível em:<http://portal.acm.org/citation.cfm?id=37422>.
LOSASSO, Frank; GIBOU, Frederic; FEDKIW, Ron. Simulating water and smoke with anoctree data structure. ACM Trans. Graph., ACM Press, New York, NY, USA, v. 23, n. 3, p.457–462, 2004. ISSN 0730-0301.
LUCY, Leon B. A numerical approach to testing the fission hypothesis. v. 82, n. 12, p.1013–1924, dez. 1977.
MOLLER ERIC HAINES, Tomas Akenine-Moller Tomas. Real-Time Rendering. 2Âa. ed.[S.l.]: AK Peters, Ltd., 2002.
MONAGHAN, J.J. Smoothed particle hydrodynamics. Annual Review of Astronomy andAstrophysics, p. 543–574, 1992.
MORRIS, J. P. Simulating surface tension with smoothed particle hydrodynamics. InternationalJournal for Numerical Methods in Fluids, v. 33, p. 333–353, jun. 2000.
M.S.CRAMER. Navier-Stokes Equations - Incompressible Flows. 2002. Disponível em:<http://www.navier-stokes.net/nsinc.htm>.
MUELLER, Matthias; SCHIRM, Simon; TESCHNER, Matthias; HEIDELBERGER, Bruno;GROSS, Markus. Interaction of Fluids with Deformable Solids. 2004. Disponível em:<citeseer.ist.psu.edu/mueller04interaction.html>.
MULLER, Matthias; CHARYPAR, David; GROSS, Markus. Particle-based fluid simulation forinteractive applications. In: SCA ’03: Proceedings of the 2003 ACM SIGGRAPH/Eurographicssymposium on Computer animation. Aire-la-Ville, Switzerland, Switzerland: EurographicsAssociation, 2003. p. 154–159. ISBN 1-58113-659-5.
MULLER, Matthias; SCHIRM, Simon; TESCHNER, Matthias. Interactive blood simulationfor virtual surgery based on smoothed particle hydrodynamics. Technol. Health Care, IOSPress, Amsterdam, The Netherlands, The Netherlands, v. 12, n. 1, p. 25–31, 2004. ISSN0928-7329.
PREMOZE, S.; TASDIZEN, T.; BIGLER, J.; LEFOHN, A.; WHITAKER, R. Particle-based si-mulation of fluids. 2003. Disponível em: <citeseer.ist.psu.edu/premoze03particlebased.html>.
REEVES, W. T. Particle systems: a technique for modeling a class of fuzzy objects. ACMTrans. Graph., ACM Press, New York, NY, USA, v. 2, n. 2, p. 91–108, 1983. ISSN 0730-0301.
STAM, Jos. Stable fluids. In: SIGGRAPH ’99: Proceedings of the 26th annual conference onComputer graphics and interactive techniques. New York, NY, USA: ACM Press/Addison-Wesley Publishing Co., 1999. p. 121–128. ISBN 0-201-48560-5.
84
STAM, Jos. Interacting with smoke and fire in real time. Commun. ACM, ACM Press, NewYork, NY, USA, v. 43, n. 7, p. 76–83, 2000. ISSN 0001-0782.
STAM, J. Real-time fluid dynamics for games. 2003. Disponível em:<citeseer.ist.psu.edu/stam03realtime.html>.
WEISSTEIN, Eric W. Plane. MathWorld–A Wolfram Web Resource., 2002. Disponível em:<http://mathworld.wolfram.com/Plane.html>.
WEISSTEIN, Eric W. Point-plane distance. MathWorld–A Wolfram Web Resource., 2003.Disponível em: <http://mathworld.wolfram.com/Point-PlaneDistance.html>.
WRENNINGE, Magnus. Fluid Simulation for Visual Effects. Dissertação (Mestrado) —Linköping University, Linköpings universitet 581 83 LINKÖPING, Sweden, 2003.
