XVII Encontro de Modelagem Computacional
V Encontro de Ciência e Tecnologia de Materiais
Universidade Católica de Petrópolis (UCP), Petrópolis/RJ, Brasil. 15-17 out. 2014
UMA FORMULAÇÃO IMPLÍCITA NO TEMPO PARA O MÉTODO SMOOTHED PARTICLE
HYDRODYNAMICS
Ricardo Dias dos [email protected]
Programa de Pós-Graduação em Engenharia Mecânica
Universidade Federal Fluminense
Rua Passo da Pátria, 156, 24210-240 Niterói, RJ, Brasil
Helio Pedro Amaral [email protected]
Departamento de Modelagem Computacional
Instituto Politécnico, Universidade do Estado do Rio de Janeiro
Caixa Postal 97282, 28610-974 Nova Friburgo, RJ, Brasil
Resumo. Algumas formulações implı́citas foram aplicadas para a integração temporal no método Smoothed Particle Hydrodynamics (SPH) de modo a possibilitar o uso de maiores incrementos de tempo e
uma forte estabilidade no processo de marcha temporal. Devido ao alto custo computacional exigido
pela busca das partı́culas vizinhas a cada passo no tempo, esta implementação só é viável se forem
aplicados algoritmos eficientes para o tipo de estrutura matricial considerada, tais como os métodos
do subespaço de Krylov. Portanto, fez-se um estudo para a escolha apropriada dos métodos que mais
se adequavam a este problema, sendo escolhidos os métodos: Bi-Conjugate Gradient (BiCG), o BiConjugate Gradient Stabilized (BiCGSTAB) e o Quasi-Minimal Residual (QMR). Testes foram utilizados
a fim de validar as soluções numéricas obtidas com a versão implı́cita do método SPH.
Palavras-chave: Dinâmica dos fluidos computacional, Métodos implı́citos, Métodos de Runge-Kutta,
Métodos do subespaço Krylov, Smoothed Particle Hydrodynamics (SPH)
1.
INTRODUÇÃO
Em boa parte dos problemas da dinâmica dos fluidos é de interesse obter-se soluções em regime
transiente e, portanto, deve-se empregar técnicas de integração temporal. Uma primeira possibilidade,
muito usual, seria a de utilizar-se métodos explı́citos com esta finalidade, devido à sua simplicidade
e eficiência computacional. Entretanto, estes métodos frequentemente são apenas condicionalmente
estáveis e estão sujeitos a severas restrições na escolha do passo no tempo. Para problemas advectivos governados por equações hiperbólicas, esta restrição é comumente conhecida como a condição de
Courant-Friedrichs-Lewy (CFL). Quando da necessidade de buscar-se soluções numéricas para grandes
perı́odos de tempo, essa condição torna-se um empecilho uma vez que o valor máximo do incremento de
tempo não deve violá-la. A fim de contornar-se essa barreira, métodos implı́citos (geralmente incondicionalmente estáveis) de integração temporal são utilizados, apesar do seu maior custo computacional.
Atualmente, na literatura, há bem poucos relatos de trabalhos que utilizaram o método Smoothed
Particle Hydrodynamics conjuntamente com métodos implı́citos de integração temporal. Inicialmente,
XVII Encontro de Modelagem Computacional
V Encontro de Ciência e Tecnologia de Materiais
Universidade Católica de Petrópolis (UCP), Petrópolis/RJ, Brasil. 15-17 out. 2014
Monaghan (1987) propôs uma formulação implı́cita para a resolução de problemas de dinâmica dos gases. Por outro lado, Knapp (2000) apresentou uma formulação implı́cita, empregando métodos de Krylov e o método de Newton-Raphson na resolução de um problema não-linear. A transferência de calor
também é resolvida implicitamente via o método de Crank-Nicolson em Rook et al. (2007). Recentemente, Laibe e Price (2012) também utilizaram uma formulação implı́cita e o método SPH na resolução
de problemas astrofı́sicos. Por último, Lanzafame (2013) introduziu uma metodologia implı́cita na qual
não é necessária a inversão da matriz Jacobiana dos coeficientes para a obtenção da solução numérica.
Entretanto, não é do conhecimento dos autores desse trabalho a existência de aplicações do método
SPH que utilizem métodos de Runge-Kutta implı́citos na resolução do sistema de equações diferenciais
ordinárias com relação ao tempo, oriundo da discretização espacial das equações de balanço.
Sabe-se que a utilização dos métodos implı́citos implica na resolução de um sistema algébrico,
cujo número de equações depende do tamanho da malha computacional utilizada, o que acarreta num
custo computacional extra. No método SPH, que usa uma abordagem Lagrangiana e não utiliza malhas espaciais, o tamanho do sistema linear a ser resolvido, para a obtenção das grandezas fı́sicas de
interesse, será da ordem do número de partı́culas utilizadas na simulação. Este valor, que pode chegar a centenas de milhares ou mesmo milhões, acaba sendo um fator desestimulador para o uso dos
métodos implı́citos (Domı́nguez et al., 2011). Em contrapartida, como o maior custo computacional do
método SPH reside na determinação das partı́culas próximas vizinhas (Góes, 2011), quando é empregada
a técnica da Busca Direta, entende-se que o uso de passos de tempo maiores (violando a condição de
CFL) podem resultar em um ganho de eficiência computacional apesar do custo adicional da resolução
de um sistema algébrico de equações, em função da diminuição da necessidade do emprego dos algoritmos de busca pelas partı́culas vizinhas. Além disso, o método SPH é altamente paralelizável, conforme
pode ser visto, por exemplo, em Góes (2011).
2.
CONCEITOS BÁSICOS DO MÉTODO SPH
Quando o método da Representação Integral e a técnica da Aproximação por Partı́culas (Liu e Liu,
2003) do método SPH são aplicadas às equações diferenciais parciais que governam o fenômeno fı́sico,
obtém-se um sistema semi-discreto de equações diferenciais ordinárias (EDOs) dependentes da variável
temporal, cuja resolução numérica será o foco deste trabalho. Algumas ideias fundamentais inerentes ao
método SPH são (Liu e Liu, 2003):
• O domı́nio do problema é representado por uma distribuição arbitrária de partı́culas, não necessitando de uma conectividade entre elas. O que acarreta na natureza do método livre de malhas;
• O método da Representação Integral é utilizado para aproximar as funções de campo, conhecido
também como aproximação de núcleo. Este método é fundamental para que se possa garantir a
estabilidade necessária para o método SPH;
• A Aproximação por Partı́culas se constitui na substituição da aproximação de núcleo por uma
aproximação usando partı́culas, que consiste em substituir a Representação Integral das funções
e suas derivadas por somatórios, ponderados pela função núcleo, dos valores correspondentes de
todas as partı́culas dentro do domı́nio de suporte da função núcleo;
• A aproximação das partı́culas é executada em cada passo no tempo e depende da distribuição
espacial atual das partı́culas;
XVII Encontro de Modelagem Computacional
V Encontro de Ciência e Tecnologia de Materiais
Universidade Católica de Petrópolis (UCP), Petrópolis/RJ, Brasil. 15-17 out. 2014
• A Aproximação por Partı́culas é aplicada a todos os termos nas EDPs de forma a produzir um
sistema semi-discreto de EDOs em relação ao tempo;
• O sistema de EDOs é resolvido utilizando-se métodos de integração temporal, de forma a obter-se
os valores para as variáveis dependentes em todas as partı́culas ao longo da evolução no tempo.
2.1 Equações de Balanço
O Postulado da Conservação da Massa diz que a massa de um corpo material (composto sempre pelas mesmas partı́culas) é constante, ou seja, independente do tempo. Partindo desse princı́pio e aplicando
o Teorema do Transporte de Reynolds (Slattery, 1999) para uma grandeza especı́fica φ = ρ, obtém-se a
equação de conservação da massa ou equação da continuidade
Dρ
+ ρ∇ · v = 0
Dt
(1)
com ρ sendo a massa especı́fica do corpo e v o vetor velocidade.
Segundo o Princı́pio da Conservação da Quantidade de Movimento Linear, a taxa de variação da
quantidade de movimento de um corpo é igual à soma das forças externas agindo sobre o corpo. Assim
sendo, a Primeira Lei de Cauchy do Movimento na sua forma diferencial, para o escoamento incompressı́vel de um fluido Newtoniano, resulta na equação de Navier-Stokes (Slattery, 1999)
ρ
Dv
= −∇p + µ∇ · (∇v) + fext
Dt
(2)
com p sendo a pressão, µ a viscosidade e fext uma força externa.
2.2 Discretização das Equações de Balanço
Na aplicação do método SPH a implementação de duas etapas fundamentais são necessárias, a
primeira é a Representação Integral, também conhecida como Aproximação de Núcleo e a segunda é a
Aproximação por Partı́culas (Liu e Liu, 2003).
A seguir são apresentadas as equações de balanço discretizadas usando-se as formulações do método
SPH. Com esta finalidade são utilizados alguns operadores diferenciais, escolhidos em função da acurácia
e eficiência computacional (dos Santos, 2014). A equação da continuidade, apresentada na Eq. (1), tem
a seguinte forma discretizada (dos Santos, 2014)
N
Dρi X
=
mj vij · ∇i Wij
Dt
j=1
(3)
com vij = vi − vj , Wij = W (xi − xj , h) sendo W a função de suavização (ou núcleo) e h representa
o comprimento de suavização que define o domı́nio de influência da função núcleo (Liu e Liu, 2003). A
forma final discretizada da equação de Navier-Stokes para escoamentos incompressı́veis é dada por (dos
Santos, 2014)
N
N
X
Dvi
1 X
2mj (µi + µj )
xij
1
=− 2
mj pij ∇i Wij +
vij 2
·
∇
W
+
fext
i
ij
2)
Dt
ρi j=1
ρ
ρ
(r
+
η
ρ
i
j
ij
j=1
(4)
XVII Encontro de Modelagem Computacional
V Encontro de Ciência e Tecnologia de Materiais
Universidade Católica de Petrópolis (UCP), Petrópolis/RJ, Brasil. 15-17 out. 2014
sendo η um parâmetro de regularização, normalmente dado por η 2 = 10−5 h2 (Morris et al., 1997).
Como alternativa à formulação WCSPH (Weak Compressible Smoothed Particle Hydrodynamics)
(Liu e Liu, 2003), que considera o fluido como sendo quase-incompressı́vel, considera-se neste estudo
a versão incompressı́vel do método SPH conhecida como PCISPH (Predictive-Corrective Incompressible SPH), onde uma equação para a determinação da pressão é introduzida e solucionada numericamente. Nesta linha, de Freitas e Souto (2014) propuseram um algoritmo iterativo baseado no Método da
Projeção, que emprega a decomposição de Helmholtz-Hodge, no qual o campo de pressões é obtido ao
final de duas etapas. Primeiro resolve-se a equação de Navier-Stokes para escoamentos incompressı́veis
(Eq. (2)) sem as forças de pressão, obtendo-se um campo intermediário de velocidades

µ
fext
Dv∗


=
∇ · (∇v∗ ) +
em Ω

Dt
ρ
ρ
(5)


 v∗ = 0
em δΩ
sendo δΩ a fronteira do domı́nio.
Calculado o campo de velocidades v∗ , a próxima etapa consiste em obter-se o campo de velocidades
n+1
v
que satisfaça à condição de incompressibilidade, através do seguinte problema

Dv
1


= − ∇pn+1 em Ω


Dt
ρ




(6)
∇·v = 0
em Ω







 ∂p = 0
em δΩ
∂n
com pn+1 sendo o campo de pressões no instante de tempo tn + ∆t. Assim, para que o novo campo de
velocidades tenha divergência nula (∇ · vn+1 = 0) deve-se satisfazer à equação de Poisson
1 n+1
1
∇·
∇p
=
∇ · v∗
(7)
ρ
∆t
sendo a condição ∂p/∂n = 0 em δΩ imposta mediante uso das partı́culas fantasmas na fronteira. Então,
o processo iterativo faz uso da equação (de Freitas e Souto, 2014)
N
N
X
mj
4
xij
1 X
∗
pnj 2
·
∇
W
−
∆t
mj vij
· ∇i Wij
i ij
2
ρ
ρ
+
ρ
(r
+
η
)
ρ
j
i
j
i
ij
j=1
j=1
p∗i =
(8)
N
X mj
4
xij
· ∇i Wij
2
ρj ρi + ρj (rij
+ η2)
j=1
e de pn+1
= pni − α∆pi , onde ∆pi = p∗i − pni e α é um parâmetro de sub-relaxação. Calculado o campo
i
de pressões pn+1 , a velocidade para o próximo instante de tempo pode ser obtida por (Eq. (6))
vn+1 = v∗ −
∆t n+1
∇p
ρ
(9)
XVII Encontro de Modelagem Computacional
V Encontro de Ciência e Tecnologia de Materiais
Universidade Católica de Petrópolis (UCP), Petrópolis/RJ, Brasil. 15-17 out. 2014
Por último, obtida a velocidade de cada partı́cula no domı́nio de interesse, faz-se necessário a
determinação da posição futura das respectivas partı́culas, que pode ser calculada resolvendo-se numericamente Dxi /Dt = vi .
2.3 Condições de Contorno
Quando se deseja obter soluções para fenômenos fı́sicos modelados por equações diferenciais é
de extrema importância impor corretamente, para a geometria do domı́nio, as condições de contorno,
de forma a assegurar que os problemas sejam bem-postos. Representar corretamente estas condições
no método SPH, de forma a impor-se o comportamento desejado, pode não ser trivial, visto que as
partı́culas podem se mover indefinidamente desde que dentro das limitações impostas pelas equações
de balanço. Portanto, para um domı́nio finito de simulação, neste método são acrescentadas forças
externas com a finalidade de impedir que as partı́culas atravessem as fronteiras do domı́nio imposto. Ao
se fazer isso, o truncamento do domı́nio de suporte perto da fronteira pode acarretar na não verificação
da condição de normalização da função núcleo (Liu e Liu, 2003), além do fato de que haverá uma
deficiência de partı́culas dentro do domı́nio de suporte, fazendo com que as aproximações resultem em
soluções incorretas (Neto, 2007). O tratamento adequado das condições de contorno, incluindo o uso de
partı́culas fantasmas dos tipos I e II, bem como condições de contorno periódicas, pode ser encontrado
em de Freitas e Souto (2014) e dos Santos (2014).
2.4 Estrutura da Matriz dos Coeficientes
Na formulação implı́cita os valores das variáveis dependentes avaliados para as partı́culas vizinhas
também devem ser determinados no próximo instante de tempo n + 1. Este procedimento é similar ao
utilizado nos métodos que usam malhas computacionais e implica na resolução de um sistema algébrico
de equações. Como no método SPH as partı́culas dentro do domı́nio de simulação não possuem qualquer
conectividade entre elas, ao contrário da conectividade existente entre os nós de uma malha computacional, não é possı́vel saber antecipadamente quais posições serão ocupadas por estas partı́culas em um
determinado instante de tempo futuro. Portanto, a vizinhança de uma dada partı́cula poderá conter, em
princı́pio, qualquer uma das demais partı́culas. Assim, para cada passo de tempo, o sistema linear a ser
resolvido, do tipo Abn+1 = cn , possui uma matriz dos coeficientes A esparsa (Fig. 1) e é da ordem do
número de partı́culas utilizadas na simulação. Por esse motivo, decidiu-se usar os métodos do subespaço
Krylov, visto a comprovada eficiência computacional quando da resolução de sistemas esparsos.
3.
MÉTODOS DO SUBESPAÇO DE KRYLOV
Algumas questões devem ser levadas em consideração quando se trata de resolver um sistema linear
algébrico. Dentre as principais destaca-se a arquitetura fı́sica (hardware) disponı́vel, que irá definir
o balanço entre a memória necessária a ser alocada e o número de operações primordiais exigidas pelo
algoritmo. Para os métodos aqui estudados são elas: o produto interno, a multiplicação de uma matriz por
um vetor, a soma entre vetores e a resolução do sistema obtido com o pré-condicionamento. Claramente,
esta última vai depender do pré-condicionador que melhor irá se adequar à estrutura matricial oriunda
do sistema a ser resolvido. A questão envolvendo o espaço de memória exigido é de suma importância,
pois pretende-se resolver sistemas lineares esparsos da ordem do número de partı́culas empregadas na
discretização do domı́nio de simulação. Uma caracterı́stica importante dos métodos do subespaço de
XVII Encontro de Modelagem Computacional
V Encontro de Ciência e Tecnologia de Materiais
Universidade Católica de Petrópolis (UCP), Petrópolis/RJ, Brasil. 15-17 out. 2014
Figura 1: Estrutura da matriz dos coeficientes em função do posicionamento das partı́culas.
Krylov é que eles podem ser aplicados em ambientes de execução computacional em paralelo.
Na Tabela 1 mostra-se, para os métodos de Krylov que podem ser aplicados a matrizes não-simétricas,
a quantidade de memória necessária por iteração i. Nessa tabela matriz denota o espaço necessário para
a alocação da matriz e n a ordem da matriz dos coeficientes. Os métodos de Krylov não-estacionários que
podem ser aplicados a matrizes não-simétricas são: o Generalized Minimal Residual (GMRES), o BiConjugate Gradient (BiCG), o Quasi-Minimal Residual (QMR), o Conjugate Gradient Squared (CGS)
e o BiConjugate Gradient Stabilized (Bi-CGSTAB). Levando-se em consideração todas as questões relativas aos prós e contras desses métodos (dos Santos, 2014), optou-se por utilizar os métodos BiCG,
Bi-CGSTAB e QMR.
Tabela 1: Espaço de memória exigido por iteração i.
Método
Espaço em memória
GMRES
matriz + (i + 5)n
BiCG
matriz + 10n
QMR
matriz + 16n
CGS
matriz + 11n
Bi-CGSTAB
matriz + 10n
XVII Encontro de Modelagem Computacional
V Encontro de Ciência e Tecnologia de Materiais
Universidade Católica de Petrópolis (UCP), Petrópolis/RJ, Brasil. 15-17 out. 2014
4.
MÉTODOS DIAGONALMENTE IMPLÍCITOS DE RUNGE-KUTTA
A famı́lia dos métodos de Runge-Kutta pode ser representada seguindo o princı́pio descrito por Butcher (1987), no qual os coeficientes são estrategicamente colocados em uma tabela para facilitar o cálculo
de suas propriedades. A forma geral para se determinar a solução numérica é dada por
yn+1 = yn + ∆x
m
X
bi k i
i=1
para 1 ≤ i ≤ s, sendo que nesta equação ki = f (xn + ci ∆x, Yi ), i = 1, 2, . . . , s e Yi é dado por
Yi = yn + ∆x
s
X
aij kj
j=1
sendo s o número de estágios do método.
Os coeficientes caracterı́sticos (aij , bi e ci ), que vão determinar a acurácia e a estabilidade desses
métodos, são determinados convenientemente através da tabela de Butcher (1987)
c1 a11 a12 · · · a1s
c2 a21 a22 · · · a2s
..
..
..
..
...
.
.
.
.
cs as1 as2 . . . ass
c A
bT
b1
b2
···
bs
Os coeficientes do vetor c estão relacionados com os elementos da matriz A
ci =
s
X
aij
i = 1, 2, ..., s.
j=1
e o vetor b deve ser de tal que
s
X
bi = 1
i=1
As formulações para as quais aij = 0 se j > i são chamadas de formulações Diagonalmente
Implı́citas ou DIRK(e,o) (Diagonally Implicit Runge-Kutta), onde e indica o número de estágios e o a
ordem do método
c A
bT
c1 a11 0 · · · 0
c2 a21 a22 · · · 0
..
..
..
.
..
. ..
.
.
.
cs as1 as2 . . . ass
b1
b2
···
bs
XVII Encontro de Modelagem Computacional
V Encontro de Ciência e Tecnologia de Materiais
Universidade Católica de Petrópolis (UCP), Petrópolis/RJ, Brasil. 15-17 out. 2014
Nesse trabalho formulações estáveis (dos Santos, 2014) do tipo DIRK são empregadas na resolução do
sistema de equações diferenciais ordinárias oriundo da discretização das equações de balanço. Especificamente, utiliza-se os métodos DIRK(1,2), DIRK(2,3), DIRK(3,4) e DIRK(3,3).
Na sequência apresenta-se os coeficientes que determinam esses diferentes métodos implı́citos.
Inicia-se considerando a formulação DIRK(1,2), cujos valores de c1 , a11 e b1 são
1
2
1
2
(10)
1
enquanto que no casos dos métodos DIRK(2,3) e DIRK(3,4) os coeficientes são dados por
1
1
+ √
2 2 3
1
1
− √
2 2 3
1
1
+ √
2 2 3
1
−√
3
1/2
1
1
+ √
2 2 3
1/2
(11)
e
(1 + α)/2 (1 + α)/2
1/2
−α/2
(1 + α)/2
(12)
(1 − α)/2
1+α
−(1 + 2α) (1 + α)/2
1/(6α2 ) 1 − 1/(3α2 ) 1/(6α2 )
√
sendo α = 2 cos(π/18)/ 3. Em se tratando da última formulação empregada, a DIRK(3,3), tem-se que
α
α
c2 c2 − α α
b1
b2 α
1
b1
b2 α
onde α é a raiz de x3 − 3x2 + 32 x − 16 = 0, no intervalo
−(6α2 − 16α + 1)/4 e b2 = (6α2 − 20α + 5)/4.
5.
(13)
1 2
,
6 2
, sendo que c2 = (1 + α)/2, b1 =
O PROBLEMA DA CAVIDADE
O problema tratado nesse trabalho consiste no escoamento incompressı́vel de um fluido Newtoniano viscoso no interior de uma cavidade quadrada. O escoamento ocorre graças ao fato de que a
parede superior da cavidade é posta em movimento com uma velocidade constante (Vtop ) enquanto as
demais paredes permanecem fixas. Na discretização espacial do domı́nio de resolução constatou-se
que a utilização de 1.600 partı́culas igualmente espaçadas gerariam resultados suficientemente acurados,
Fig. 2. Na representação das fronteiras fı́sicas desse domı́nio foram utilizadas 160 partı́culas virtuais do
tipo I (dos Santos, 2014).
A massa especı́fica é determinada a partir da equação do somatório normalizado e uma função
quı́ntica por partes é utilizada como função de suavização e não considera-se os efeitos da força gravitacional (dos Santos, 2014). Os principais parâmetros utilizados na simulação numérica foram altura e
XVII Encontro de Modelagem Computacional
V Encontro de Ciência e Tecnologia de Materiais
Universidade Católica de Petrópolis (UCP), Petrópolis/RJ, Brasil. 15-17 out. 2014
Figura 2: Distribuição inicial das partı́culas: cavidade
comprimento L = 10−3 m; viscosidade ν = 10−3 m2 /s; massa especı́fica ρ = 103 kg/m3 ; velocidade
inicial v = (0, 0) m/s; velocidade da parede superior da cavidade Vtop = 10−3 m/s; espaçamento inicial entre partı́culas ∆x = 2, 5 × 10−5 m; comprimento de suavização h = 1, 2 ∆x; parâmetro para a
correção XSPH ε = 0, 3 (dos Santos, 2014) e número de Reynolds Re = Vtop L/ν = 1.
A fim de obter-se valores para o campo de velocidades que poderiam ser confrontados com os
resultados numéricos obtidos com a versão PCISPH, utilizou-se um simulador desenvolvido a partir do
método dos Volumes Finitos (MVF) que emprega esquemas do tipo Total Variation Diminishing (TVD)
na determinação dos fluxos convectivos nas faces dos volumes de controle (Schulz e Amaral Souto,
2010).
Como o objetivo principal desse trabalho é a proposição de uma versão implı́cita para o método SPH,
os resultados calculados são provenientes da implementação das formulações implı́citas DIRK(1,2),
DIRK(2,3), DIRK(2,2), DIRK(3,4) e DIRK(3,3) e optou-se por manter o incremento de tempo ∆t constante e igual 5, 0 × 10−3 s. Primeiramente são mostrados, na Fig 3, os valores correspondentes aos
perfis da componente horizontal adimensionalizada do vetor velocidade considerando os métodos MVF
e PCISPH. Dessa figura nota-se que embora os resultados reproduzam adequadamente o perfil de velocidade, nem todos são suficientemente acurados e resultam em aproximações que são mesmo inferiores
aos resultados fornecidos pela implementação explı́cita (dos Santos, 2014), principalmente próximo à
região de fronteira superior do domı́nio de resolução.
Com relação ao perfil da componente vertical adimensionalizada do vetor velocidade, os resultados
são comparados entre si e os seus valores, tanto os do método PCISPH quanto o do MVF, mostrados na
Fig. 4. Também para esse perfil de velocidade pode-se averiguar que nem todas as formulações implı́citas
foram capazes de fornecer valores que aproximam corretamente o perfil fornecido pelo método dos
Volumes Finitos, principalmente próximo das fronteiras laterais aonde as partı́culas de fluido deveriam
verificar a condição de não-deslizamento.
XVII Encontro de Modelagem Computacional
V Encontro de Ciência e Tecnologia de Materiais
Universidade Católica de Petrópolis (UCP), Petrópolis/RJ, Brasil. 15-17 out. 2014
1.0
0.9
0.8
y(adim ensional)
0.7
0.6
0.5
0.4
0.3
MVF
DIRK(1,2)
DIRK(2,3)
DIRK(2,2)
DIRK(3,4)
DIRK(3,3)
0.2
0.1
0.0
-0.4
-0.2
0.0
0.2
0.4
0.6
0.8
1.0
Vx/Vt op
Figura 3: Comparação entre as soluções obtidas com o Método de Volumes Finitos e o PCISPH para as
formulações implı́citas em relação a componente horizontal adimensionalizada da velocidade
0.25
MVF
DIRK(1,2)
DIRK(2,3)
DIRK(2,2)
DIRK(3,4)
DIRK(3,3)
0.20
0.15
0.10
Vy/Vt op
0.05
0.00
-0.05
-0.10
-0.15
-0.20
-0.25
0.0
0.1
0.2
0.3
0.4
0.5
0.6
x(adim ensional)
0.7
0.8
0.9
1.0
Figura 4: Comparação entre as soluções obtidas com o Método de Volumes Finitos e o PCISPH para as
formulações implı́citas em relação a componente vertical adimensionalizada da velocidade
XVII Encontro de Modelagem Computacional
V Encontro de Ciência e Tecnologia de Materiais
Universidade Católica de Petrópolis (UCP), Petrópolis/RJ, Brasil. 15-17 out. 2014
De forma a poder comparar os diferentes resultados numéricos obtidos com o MVF e as formulações
explı́cita ERK(4,4) (∆t = 1, 0 × 10−4 s) e implı́citas (DIRK), a Tabela 2 mostra o valor calculado do erro
considerando a norma euclidiana das respectivas aproximações. Nota-se que o melhor resultado para o
perfil horizontal foi obtido com o DIRK(2,3) enquanto que o DIRK(3,3) forneceu o resultado com o
menor erro para o perfil vertical.
Tabela 2: Erro comparativo entre as soluções dos métodos MVF e Runge-Kutta.
Formulações
ERK(4,4)
DIRK(1,2)
DIRK(2,3)
DIRK(2,2)
DIRK(3,4)
DIRK(3,3)
Perfil horizontal
7, 34 × 10−1
7, 34 × 10−1
7, 21 × 10−1
7, 39 × 10−1
7, 43 × 10−1
7, 45 × 10−1
Perfil vertical
5, 26 × 10−2
5, 34 × 10−2
5, 70 × 10−2
6, 41 × 10−2
1, 28 × 10−1
5, 24 × 10−2
Em função dos valores obtidos, pôde-se constatar que os métodos BiCG e Bi-CGSTAB tiveram o
seu critério de convergência satisfeito, para uma tolerância tol = 1, 0 × 10−9 , com somente uma única
iteração para todas as formulações. O método Bi-CGSTAB apresentou uma maior redução do resı́duo
do que o método Bi-CG (dos Santos, 2014). Por último, o método QMR necessitou de 7 a 14 iterações,
em função da formulação implı́cita, para que a convergência fosse alcançada (dos Santos, 2014).
6.
CONCLUSÃO
Nesse trabalho foi desenvolvida uma versão implı́cita no tempo para o método Smoothed Particle Hydrodynamics empregando uma famı́lia de métodos do tipo Diagonally Implicit Runge-Kutta. Na
validação dos resultados obtidos, foi considerado o problema do escoamento bidimensional no interior
de uma cavidade quadrada, sendo que neste caso as soluções numéricas foram comparadas com àquelas
fornecidas por um simulador construı́do baseado no método dos Volumes Finitos.
As formulações mostraram-se comparáveis entre si, com pouca variação nos valores dos erros para
as duas componentes adimensionalizadas do vetor velocidade. Deve-se observar que estes resultados
foram alcançados com as formulações explı́cita e implı́citas empregando incrementos de tempo iguais a
∆t = 1, 0 × 10−4 e 5, 0 × 10−3 , respectivamente. Portanto, um passo de tempo cinquenta vezes maior
que o da formulação explı́cita foi utilizado com êxito na formulação implı́cita.
Na resolução dos sistemas lineares, todos os métodos de Krylov empregados convergiram. No
entanto, o método QMR precisou de um número mais alto de iterações, quando comparado aos métodos
BiCG e Bi-CGSTAB, para que a convergência fosse alcançada.
REFERÊNCIAS
Butcher, J. C., 1987. The Numerical Analysis of Ordinary Differential Equations: Runge-Kutta and General Linear
Methods. New York-USA: Wiley-Interscience.
Domı́nguez, J. M.; Crespo, M.; Dost, S., 2011. Optimization strategies for parallel CPU and GPU implementations
of a meshfree particle method. Computing Research Repository, abs/1110.3711.
XVII Encontro de Modelagem Computacional
V Encontro de Ciência e Tecnologia de Materiais
Universidade Católica de Petrópolis (UCP), Petrópolis/RJ, Brasil. 15-17 out. 2014
Gingold, R. A. & Monaghan, J. J., 1977. Smoothed Particle Hydrodynamics: Theory and application to nonspherical stars. Physical Review, v. 181, p. 375-389.
de Freitas, M. M. & Amaral Souto, H. P., 2014. Simulação de Escoamentos Incompressı́veis Empregando o Método
SPH Utilizando Um Algoritmo Iterativo na Determinação do Campo de Pressões. Revista de Engenharia da
Universidade Católica de Petrópolis, v. 8, 2, p. 1-20.
Góes, J. F., 2011. Resolução numérica de escoamentos compressı́veis empregando um método de partı́culas livre
de malhas e o processamento em paralelo (CUDA). Dissertação de Mestrado, Universidade do Estado do Rio
de Janeiro, Instituto Politécnico, Nova Friburgo, Brasil.
Knapp, C. E., 2000. An Implicit Smooth Particle Hydrodynamic Code. Tese de Doutorado, University of New
Mexico, Albuquerque, New Mexico.
Laibe, G.; Price, D. J., 2012. Dusty gas with Smoothed Particle Hydrodynamics - II. Implicit timestepping and
astrophysical drag regimes. Monthly Notices of the Royal Astronomical Society, v. 420, p. 2365-2376.
Lanzafame, G., 2013. Implicit integrations for sph in semi-lagrangian approach: Application to the accretion disc
modeling in a microquasar. Computer Physics Communications, v. 184, p. 671-688.
Liu, G. R.; Liu, M. B., 2003. Smoothed Particle Hydrodynamics: A Meshfree Particle Method. Singapore: World
Scientific.
Lucy, L. B., 1977. Numerical approach to testing the fission hypothesis. Astronomical Journal, v. 82, p. 1013-1024.
Monaghan, J. J., 1987. Implicit SPH drag and dusty gas dynamics. Journal of Computational Physics, v. 138,
p. 801-820.
Morris, J. P.; Fox, P.; Zhu, Y., 1997. Modeling low reynolds number incompressible flows using SPH. Journal of
Computation Physics, v. 136, p. 214-226.
Neto, A. P., 2007. Uma Abordagem Lagrangeana para Simulação de Escoamentos de Fluidos Viscoplásticos e
Multifásicos. Tese de Doutorado, Pontifı́cia Universidade Católica, Rio de Janeiro, Brasil.
dos Santos, R. D., 2014. Uma Formulação Implı́cita para o Método Smoothed Particle Hydrodynamics.
Dissertação de Mestrado, Universidade do Estado do Rio de Janeiro, Instituto Politécnico, Nova Friburgo,
Brasil.
Rook, R.; Yildiz, A. J. C. & Gesteira, M. G., 2007. Modeling transient heat transfer using SPH and implicit time
integration. Numerical Heat Transfer, Part B: Fundamentals, v. 51, p. 1-23.
Slattery, J. C., 1999. Advanced Transpot Phenoma. Cambridge University Press.
Schulz, J. A. T. & Amaral Souto, H. P., 2010. Numerical Simulation of a Two-Phase Flow in a Porous Medium Employing the Finite Volume Method and a Semi-Implicit Approach. International Review of Chemical
Engineering, v. 2, p. 728-735.
AN IMPLICIT-TIME FORMULATION FOR THE SMOOTHED PARTICLE HYDRODYNAMICS METHOD
Abstract. In this study, some implicit formulations for time integration are used in the Smoothed Particle Hydrodynamics (SPH) method to enable the use of larger time increments and obtain a strong stability in the time
evolution process. Due to the high computational cost required by the particles tracking at each time step, the
implementation will be feasible only if efficient algorithms were applied for this type of matrix structure such as
Krylov subspace methods. Therefore, we carried out a study for the appropriate choice of methods best suited to
this problem, and the methods chosen were the Bi-Conjugate Gradient (BiCG), the Bi-Conjugate Gradient Stabilized (BiCGSTAB) and the Quasi-Minimal Residual(QMR). Tests were used to validate the numerical solutions
obtained with the implicit version of the SPH method.
Keywords: Computational fluid dynamics, Smoothed Particle Hydrodynamics, Implicit methods, Runge-Kutta
methods, Krylov subspace methods