Uma técnica de correção de interface para o método "Incompressible Smoothed Particle Hydrodynamics"

Uma técnica de correção de interface para o método "Incompressible Smoothed Particle Hydrodynamics"

D.F. Cordeiro^I; F.S. Sousa^I,*; A. Castelo Filho^I; J.M. Nóbrega^II

^IInstituto de Ciências Matemáticas e de Computação (ICMC), Universidade de São Paulo (USP), 13566-590 São Carlos, SP, Brasil. E-mails: [e-mail]; [e-mail]; [e-mail]

^IIDepartamento de Engenharia de Polímeros, Universidade do Minho, 4800-058 Guimarães, Portugal. E-mail: [e-mail]

RESUMO

A manutenabilidade de uma boa distribuição de partículas é algo fundamental em simulações através do SPH, sendo um dos fatores responsáveis pela estabilidade do sistema e a acurácia do método. Em soluções através do método ISPH (acrônimo em inglês para Incompressible Smoothed Particle Hydrodynamcs) este resultado pode ser alcançado através de uma técnica baseada na introdução de um deslocamento de partículas, de forma a retirá-las de suas trajetórias lagrangeanas. Entretanto, quando diretamente aplicada a escoamentos multifásicos, esta técnica acaba introduzindo um comportamento não físico na interface entre fluidos, resultando em um incorreto posicionamento da interface, e prejudicando a simulação. Neste trabalho, propomos uma correção baseada em uma função suave, definida inicialmente como a distância normal à interface. Testes numéricos mostram que esta correção consegue recuperar a correta posição da interface. Finalmente, um estudo sobre a conservação de massa desta técnica é apresentado.

Palavras-chave: ISPH, correção de interface, escoamentos multifásicos.

ABSTRACT

The maintainability of a good particle distribution is of prior importance toSPH (Smoothed Particle Hydrodynamics) simulations, since it affects both stability and accuracy of the numerical results. When using ISPH (Incompressible SPH), the particle distribution can be improved by dislocating particles away from their Lagrangian trajectories. However, when applying this technique to multiphase flows, the correct position of the interface between fluids is disturbed, resulting in unphysical flow behavior. In this work, we propose a correction, based upon a constructed smooth distance function. Numerical tests show that this correction can recover the expected position of the interface. Finally, issues regarding mass conservation are assessed and reported.

Keywords: ISPH, interface correction, multiphase flows, particle distribution.

1 INTRODUÇÃO

O método SPH (do inglês, Smoothed Particle Hydrodynamics) é um método de aproximação livre de malha que, através de um conjunto finito de partículas e uma formulação completamente Lagrangeana, permite a solução de diversos tipos de escoamentos. Proposto por Gingold &Monaghan [5 ] e Lucy [8 ], o método foi inicialmente desenvolvido para o tratamento de escoamentos compressíveis. De um modo geral, o SPH aproxima o valor de uma função ou derivada através de conceitos de interpolação aliados a uma função núcleo. As vantagens características do método, tais como a eficiente aproximação entre contínuo e discreto [7 ], a semelhança com problemas de dinâmica molecular, além dos ganhos computacionais, motivaram o seu aperfeiçoamento para o tratamento de escoamentos incompressíveis. Este aperfeiçoamento foi alcançado através de um esquema de aproximação, no qual os fluidos incompressíveis são tratados como fluidos quasi-compressíveis, sendo denominado de WCSPH (do inglês, Weakly Compressible SPH) [9 ], onde a pressão é obtida por uma equação de estado, baseada essencialmente nos valores da velocidade do som e na densidade de partículas. Entretanto, o WCSPH apresenta alguns problemas, como a necessidade de um menor passo de tempo devido aos altos valores de velocidade do som [ ], e a possibilidade de um comportamento não físico no campo de pressão, podendo tornar o sistema instável.

Para resolver os problemas apresentados pelo WCSPH, Cummins e Rudman [4 ] propuseram calcular a pressão através de uma equação de Poisson. Este novo esquema, baseado em conceitos do método da projeção [2 ] e denominado de ISPH (Incompressible SPH), é mais estável que o WCSPH, e permite passos de tempo maiores, além da obtenção de resultados para a pressão com uma maior coerência física [6 ].

Apesar destas vantagens, o ISPH não é capaz de resolver um problema característico do SPH: a possibilidade de formação de regiões de alta densidade de partículas em escoamentos com um elevado movimento inercial, relacionado ao número de Reynolds. Para tratar este problema em soluções através do ISPH, torna-se necessária a aplicação de uma técnica que garanta uma boa distribuição de partículas ao longo de toda a simulação. Uma das primeiras propostas para resolver este problema foi apresentada por Chaniotis et al. [1 ], onde as partículas tem suas posições periodicamente reiniciadas com base em uma malha uniforme, e os valores de suas variáveis interpolados para as novas posições. Uma desvantagem desta técnica, é o fato de o sistema perder a característica de totalmente livre de malhas. Outra técnica foi apresentada por Xu et al. [10 ], onde a é aplicado um deslocamento nas partículas e realizada uma correção das variáveis hidrodinâmicas por série de Taylor. Embora ambas estratégias alcancem bons resultados, quando aplicadas a escoamentos multifásicos, acabam gerando um comportamento não físico na interface entre as fases, prejudicando a simulação. Desta forma, propomos aqui uma correção da interface em escoamentos multifásicos a ser aplicada ao método de deslocamento de Xu et al. [10 ], que se baseia na definição e advecção de uma função distância suave, bem conhecida de métodos de Level-set.

Este trabalho é organizado da seguinte forma: na Seção 2 apresentamos a formulação SPH para escoamentos incompressíveis, seguida da descrição do método ISPH na Seção 3. Na Seção 4 descrevemos o método de correção das trajetórias proposto por Xu et al., e em seguida, na Seção 5, definimos o problema e sua correção através da resolução de um problema de transporte de partículas em um campo de velocidades conhecido. Resultados são apresentados ainda na Seção 5, e conclusões na Seção 6.

2 FORMULAÇÃO SPH

O SPH se propõe a resolver as equações de Navier-Stokes em sua forma lagrangeana, dadas por

onde ρ é a densidade, p a pressão, µ a viscosidade dinâmica, v o vetor velocidade e g as forças externas.

Neste sentido, o SPH segue a seguinte formulação: seja f uma função qualquer, sabe-se que esta pode ser representada na seguinte forma integral, dita rigorosa ou exata,

onde δ(x- x') é o Delta de Dirac. Como não é possível definir uma função que satisfaça as propriedades da distribuição delta de Dirac, usa-se uma aproximação por uma função núcleo W:

A integral de uma determinada função aplicada a um ponto do domínio do problema, denominado Θ, pode ser discretizada através de um somatório sobre as partículas do sub-espaço Ω_i de Θ definido pelo suporte compacto, onde Ω_i = x_j, |x_i - x_j| < κh, sendo h o comprimento do núcleo e κ uma constante que define a área não nula da função núcleo, e o volume infinitesimal dx aproximado por δV, o qual refere-se ao volume ocupado por uma partícula, tal que

onde m_j é a massa da partícula e W_ij = W(x_i -x_j,h). Neste artigo, em todos os testes, foi utilizada uma função núcleo de quinta ordem.

Os operadores para cálculo de derivadas podem ser obtidos de maneira similar [7 ]. Neste sentido, o operador gradiente SPH é dado pela seguinte expressão

sendo similar ao operador divergente. Estes operadores são utilizados no ISPH durante o cálculo do gradiente de pressão e do divergente de velocidade. O termo viscoso é obtido como

Por outro lado, o laplaciano da pressão utiliza um operador proposto em [3 ], mais adequado para escoamentos multifásicos, onde é observado descontinuidade em relação a densidade, dado por

3 O MÉTODO ISPH

Diferentemente do WCSPH, no ISPH o valor de densidade é constante, sendo associada às partículas na fase inicial da simulação, o que respeita diretamente a condição de incompressibilidade (2.1). Assim, a dinâmica do sistema é obtida através da equação de conservação da quantidade de movimento (2.2). Neste sentido, a primeira tarefa do ISPH é obter a posição intermediária x* advectando-se a posição para um passo de tempo intermediário

e com isso obter um campo de velocidade livre de divergência, denominado velocidade intermediária v*, o qual é calculado no espaço intermediário Ω* através da equação de quantidade de movimento sem a contribuição do gradiente de pressão

A partir disso, a pressão no tempo n+1 é obtida resolvendo-se uma equação de Poisson

Com os valores de pressão calculados, o campo de velocidade livre de divergência, no tempo n+1, pode ser obtido da correção

e com isso, a nova posição das partículas é finalmente obtida de

Durante a aplicação do ISPH, as condições de contorno de Dirichlet para a velocidade e de Neumann para a pressão são definidas utilizando-se partículas fantasmas [7 ].

4 A TÉCNICA DE DESLOCAMENTO DE PARTÍCULAS

A formação de aglomerados de partículas é um problema crítico em soluções de escoamentos através do ISPH, tanto pelo comportamento não físico apresentando quanto por questões relacionadas ao próprio método ISPH, tais como o mal condicionado da matriz referente à equação de Poisson (3.3), podendo ser observado em experimentos clássicos como, por exemplo, na resolução de vórtices de Taylor-Green. Na Figura 1, pode ser visto o campo de velocidade que define os vórtices de Taylor-Green (a), e o resultado do transporte de partículas neste campo, sem qualquer tratamento (b). Para corrigir este problema, introduziu-se em [10 ] uma técnica que consiste em realizar um deslocamento nas partículas, de forma a retirá-las de suas linhas de corrente originais, e posteriormente corrigir as variáveis hidrodinâmicas, fazendo com que a distribuição de partículas seja mais uniforme (como ilustrado na Fig. 1 (c)). Para isso, após o cálculo da posição no tempo n+1 através do ISPH (Equação 3.5), este deslocamento é calculado como

sendo C uma constante entre 0.01 e 0.1, α a magnitude do deslocamento, onde α = V_maxΔt, com V_max sendo a velocidade máxima do sistema e X_i é o vetor deslocamento calculado por

onde M_i é o número de vizinhos da partícula i. A partir disso, as variáveis hidrodinâmicas, no caso monofásico, velocidade e pressão, são então corrigidas através de uma expansão de primeira ordem em série de Taylor:

onde ϕ é uma váriável a ser corrigida e δx_ii' é o vetor distância entre a antiga posição da partícula x_i e nova posição da partícula x_i'. A correção possui um erro da ordem de O(). Experimentos realizados com a técnica demonstram sua robustez em relação à manutenabilidade da disposição de partículas e, além disso, a introdução de um ganho de acurácia no método ISPH [10 ].

5 UM MÉTODO DE CORREÇÃO DE INTERFACE NO ISPH

Para introduzir o problema causado pelo deslocamento das partículas lagrangeanas no ISPH, consideremos um teste numérico realizado em um domínio bidimensional [0,1 ]× [0,1 ], com um campo de velocidades prescrito e fixo, que dá origem aos chamados vórtices de Taylor-Green,

representado na Figura 1 (a).

Considerando partículas marcadoras representando dois fluidos diferentes, realiza-se a integração temporal das partículas imersas neste campo, até um determinado tempo. Em seguida, inverte-se o campo de velocidades para, integrando a mesma quantidade de passos, recuperar a configuração inicial.

Este teste foi realizado para a distribuição de partículas apresentada na Figura 2(a-d). Nesta figura, partículas de cor mais escura estão posicionadas dentro de um círculo de centro em (0.5,0.5) e raio r = 0.4. Quando as trajetórias são integradas aplicando um método Runge-Kutta de ordem 4, sem qualquer tipo de perturbação adicional, a posição inicial e a interface entre as partículas são recuperadas quase exatamente, com o erro causado apenas pelo método RK4, quase imperceptível.

Em um contexto do método ISPH, as distribuições intermediárias das partículas causariam problemas de estabilidade no método. Desta forma, foi realizado o mesmo teste, agora com a correção nas trajetórias proposta por Xu et al. [10 ]. Após cada passo de integração das trajetórias, as partículas são perturbadas de acordo com a Equação (4.1). A Figura 2(e-h) mostra a evolução no tempo deste problema. Apesar da boa distribuição de partículas no domínio em todos os passos de tempo, a interface não foi bem recuperada (Fig. 3), com a presença de partículas que cruzaram a interface, destruindo a configuração original esperada.

Assim, uma correção da interface é necessária. Para tanto, foi definida uma função suave, que em cada partícula é calculada como a distância normal à interface, ou seja, para a interface descrita acima, tem-se

Tal campo é negativo no interior do círculo, positivo no exterior, e portanto, a interface pode ser definida, a exemplo dos métodos de Level-set, como o conjunto de nível zero deste campo, mais especificamente como Γ = {(x,y): Φ(x,y) = 0}. Como este campo é suave, pode-se aplicar a mesma correção utilizada para corrigir as variáveis hidrodinâmicas por uma expansão em série de Taylor de primeira ordem, a saber,

e a nova configuração da interface é então computada, baseando-se no sinal do campo resultante. O resultado ao final de 800 passos de integração no teste de Taylor-Green é uma interface bem definida entre os "fluidos" (Fig. 4).

Outro caso, onde o problema da interface pode ser facilmente observado são os escoamentos com vórtice simples. Para isso, seja um domínio bi-dimensional [0,1 ]× [1,0 ] e o campo de velocidades prescrito e fixo, dado em termos da função de corrente

Considerando a distribuição de partículas apresentada na Figura 5(a-d), sendo as partículas marcadoras de cor mais escura localizadas em um círculo de raio 1.5 e centro em (0.5,0.75), onde as trajetórias de todas as partículas do sistema são integradas com RK4 durante 1000 passos, quando o campo de velocidades é invertido, e novamente integrada mais 1000 passos. Neste sentido, nota-se na Figura 5(a-d) a ocorrência de regiões com baixa densidade de partículas, assim como aglomerados de partículas. Aplicando-se a ténica de correção de deslocamento, obtém-se uma boa distribuição de partículas, porém a interface entre os fluidos é prejudicada (Fig. (5(e-h))). Por outro lado, quando a técnica de correção de interface é utilizada, a disposição final da interface entre os fluidos encontra-se mais próxima do esperado (Fig. 5(k-l)).

Na Figura 6 é possível ver com maior clareza a diferença que a técnica de correção de interface proporciona nos resultados finais, o que destaca a importância na resolução de problemas com interfaces.

5.1 Análise da conservação de volume

Apesar das vantagens obtidas com a correção de interface, a alteração de partículas de uma fase para outra provoca efeitos indesejados na conservação de volume (área em 2D). Em escoamentos de fluidos incompressíveis, por exemplo, a alteração do volume das fases introduz erros na simulação numérica. Utilizando o problema modelo de Taylor-Green, estudou-se este tipo de erro.

Como pode ser observado nas figuras anteriores (veja por exemplo a Fig. 6 (c)), a correção proposta diminui a quantidade de partículas posicionadas erroneamente, mas não elimina totalmente o problema. Ainda, as trocas de fase impostas pela função Level-set corrigida não conservam o volume inicial, devido ao erro de ordem O() das correções sucessivas. Uma análise mais detalhada pode ser vista na Figura 7 (esquerda), onde é apresentada a evolução do volume normalizado (da "bolha"), definido como

sendo V₀ o volume inicial e V(t) o volume no instante de tempo t.

Da Figura 7 (esquerda), percebe-se que o volume pode ser alterado para mais ou para menos ao longo da evolução da interface, sendo que uma medida pontual pode mascarar o erro cometido ao longo da integração numérica. Uma análise mais justa pode ser obtida calculando-se o acúmulo do erro em função do tempo, dado por

ou seja, a função E_V(t) mede o erro de volume acumulado até o instante de tempo t. Esta expressão pode ser aproximada através da regra trapezoidal, no passo t_n, de maneira que

Esta expressão fornecerá o erro acumulado devido às sucessivas trocas de partículas durante a evolução da interface. Na Figura 7 (direita) é apresentado o erro acumulado para o problema dos vórtices de Taylor-Green em função dos passos de integração das posições. Observa-se um crescimento linear de E_V(t) ao longo do tempo. Se a técnica de correção de interface não alterasse o volume em nenhum passo de tempo, E_V(t) seria idealmente constante. Este comportamento linear fornecido pela expressão (5.6), indica que a fração de volume |1-V_n(t)| é uma função constante no tempo, sendo sua integral uma função linear. Isso mostra que, apesar de pequenas variações em V_n(t), ao longo do tempo, na média, a fração de volume parece se manter constante, não indicando um aumento catastrófico do volume de uma das fases. Mesmo assim, espera-se obter no futuro técnicas que consigam reduzir o erro acumulado, em comparação com o apresentado aqui.

6 APLICAÇÃO EM ESCOAMENTOS BIFÁSICOS

Para ilustrar a aplicabilidade do método desenvolvido aqui na simulação de escoamentos bifásicos usando o método ISPH, foi resolvido o problema da cavidade impulsionada bifásica. Tal problema consiste em um domínio quadrado [0,1 ]× [0,1 ], com condições de contorno u = (0,0), exceto na tampa superior, onde é imposto u = (1,0). Esta simulação foi realizada até tempo t = 8 s, com massa específica constante, mas razão de viscosidade µ₁/µ₂ = 200 sendo o fluido mais viscoso posicionado em cima do menos viscoso (Fig. 8 (a)). A Figura 8 ilustra a evolução no tempo deste escoamento, comprovando o bom comportamento e não espalhamento da interface entre os fluidos, devido ao método aqui proposto.

7 CONCLUSÃO

Um método de correção de interface para aplicação do método ISPH com perturbação de trajetórias proposto por Xu et al. em escoamentos incompressíveis multifásicos foi apresentado. O método baseia-se na utilização de uma função distância suave, como feito nos métodos Level-set, e na correção das variáveis hidrodinâmicas do método de Xu et al. Resultados mostram que a correção efetivamente elimina o problema causado quando o método de Xu é aplicado em escoamentos multifásicos, o que foi demonstrado utilizando-se um teste de escoamento prescrito, e posteriormente ilustrado através da simulação de uma cavidade impulsionada multifásica. Apesar de efetivo, o método não conserva volume exatamente, mas análises indicam que não há indícios de acúmulo catastrófico de erro.

Recebido em 29 setembro, 2012

Aceito em 10 outubro, 2013

• [1 ] A.K. Chaniotis, D. Poulikakos & P. Koumoutsakos. Remeshed smoothed particle hydrodynamics for the simulation of viscous and heat conducting flows. Journal of Computational Physics, 182 (2002), 67-90.
• [2 ] A.J. Chorin. Numerical Solution of the Navier-Stokes Equations. Mathematics of Computation, 22 (1968), 745-762.
• [3 ] P. Cleary & J.J. Monaghan. Conduction modelling using smoothed particle hydrodynamics. Journal of Computational Physiscs, 148 (1999), 227-264.
• [4 ] S.J. Cummins & M. Rudman. An SPH Projection Method. Journal of Computational Physics, 152 (1999), 584-607.
• [5 ] R.A. Gingold & J.J. Monaghan. Smoothed particle hydrodynamics. Theory and application to non-spherical stars. Royal Astron. Society, 181 (1977), 375-389.
• [6 ] E.-S Lee, C. Moulinec, R. Xu, D. Violeau, D. Laurence & P. Stansby. Comparisons of weakly compressible and truly incompressible algorithms for the SPH mesh free particles method. IEEE Trans. on Visual. and C. Graphics, 227 (2008), 8417-8436.
• [7 ] G.R. Liu, M.B & M.B. Liu. "Smoothed Particle Hydrodynamics - a meshfree particle method",World Scientific, Singapure, (2003).
• [8 ] L.B. Lucy. A numerical approach to the testing of the fission hypothesis. Astron. Journal, 82 (1977), 1013-1024.
• [9 ] J.J. Monaghan. Simulating free surface flow with SPH. Journal of Computational Physics, 110 (1994), 399-406.
• [10 ] R. Xu, P. Stansby & D. Laurence. Accuracy and stability in incompressible SPH (ISPH) based on the projection method and a new approach. Journal of Computational Physics, 228 (2009), 6703-6725.

Datas de Publicação

• Publicação nesta coleção
  07 Mar 2014
• Data do Fascículo
  Dez 2013

Histórico

• Recebido
  29 Set 2012
• Aceito
  10 Out 2013