Universidade Federal de Santa Maria

Ci. e Nat., Santa Maria, v. 41, e9, 2019.

DOI: http://dx.doi.org/10.5902/2179460X33534

Received: 09/07/2018 Accepted: 22/01/2019

Section Mathematics

Análise do comportamento oscilatório de gotas de fluidos de van der Waals sob condições de microgravidade usando o método Smoothed Particle Hydrodynamics

Analysis of the oscillatory behavior of van der Waals fluid droplets under microgravity conditions using the Smoothed Particle Hydrodynamics method

Géssica Ramos da Silva^I

Maicon de Paiva Torres^II

Marciana Lima Goés^III

Helio Pedro Amaral Souto^IV

^IUniversidade do Estado do Rio de Janeiro – UERJ, Rio de Janeiro, RJ, Brasil

grsilva@iprj.uerj.br

^IIUniversidade do Estado do Rio de Janeiro – UERJ, Rio de Janeiro, RJ, Brasil

mptorres.g1@gmail.com

^IIIUniversidade do Estado do Rio de Janeiro – UERJ, Rio de Janeiro, RJ, Brasil

marcianaamorimlima@gmail.com

^IVUniversidade do Estado do Rio de Janeiro – UERJ, Rio de Janeiro, RJ, Brasil

helio@iprj.uerj.br

Resumo

No presente trabalho, o método Smoothed Particle Hydrodynamics (SPH) é utilizado na modelagem do problema bidimensional da formação de uma gota líquida, sem a presença de uma atmosfera externa, considerando um fluido de van der Waals. Além disso, também estuda-se a influência do número de Reynolds, da razão de aspecto inicial e do número de Péclet no comportamento oscilatório da gota, em condições que simulam experimentos no espaço com gravidade zero. Tendo em vista as pesquisas mais recentes na área, utiliza-se um kernel hiperbólico em todas as simulações, garantindo, assim, a formação de gotas de líquido mais uniformes, ou seja, sem aglomeração de partículas. O uso dessa função de suavização reduz a instabilidade de tensão, de modo que não é necessário empregar nenhum tratamento especial quando da aplicação do método numérico.

Palavras-chave: Oscilação. Número de Péclet. Número de Reynolds. Razão de aspecto. Smoothed Particle Hydrodynamics.

Abstract

In the present work, the Smoothed Particle Hydrodynamics method (SPH) is employed in the modeling of a two-dimensional drop formation problem, without an external atmosphere, using a van der Waals fluid. In addition, influences of Reynolds number, initial aspect ratio and Péclet number on drops oscillatory behavior are studied in open space conditions with zero gravity assumption. In view of the most recent researches in drops formation field using the SPH method, a hyperbolic kernel is applied to all simulations in order to ensure uniform drops liquid formation, i.e., without particle agglomeration. Through the usage of this smoothing function the tensile instability is reduced, avoiding unnecessary numerical treatments.

Keywords: Oscillation. Péclet number. Reynolds number. Aspect ratio. Smoothed Particle Hydrodynamics.

1 Introdução

A dinâmica da formação de gotas de líquido vêm sendo investigada há anos devido à sua vasta ocorrência na natureza e na produção industrial. Como exemplo, pode-se citar: a formação de gotas de chuva, a impressão a jato de tinta, a atomização da injeção de combustível, dentre outras (Yang et al., 2014).

A formação e a deformação de gotas de líquido estão associadas à morfologia evolutiva de interfaces livres variáveis e interfaces móveis, que representam um grande desafio às simulações numéricas. No decorrer das últimas décadas, diferentes métodos foram propostos com o intuito de modelar esse tipo de problema. Dentre eles, pode-se citar o Smoothed Particle Hydrodynamics (SPH). O SPH é um método numérico puramente lagrangiano do tipo meshfree, ou seja, que prescinde do uso de uma malha computacional, e que pode ser considerado como o método de partículas mais utilizado na atualidade. Foi inicialmente proposto por Gingold e Monaghan (1977) e Lucy (1977) para a resolução de problemas astrofísicos no espaço tridimensional e, desde então, vem sendo aplicado nas mais diversas áreas, como, por exemplo, na Dinâmica dos Fluidos Computacional.

O processo da formação de gotas, empregando o método SPH e a equação de estado de van der Waals (vdW), foi inicialmente investigado por Nugent e Posch (2000). Os autores estudaram a formação de uma gota em duas dimensões, cercada pelo seu vapor, para um fluido de van der Waals. Dentre os diversos testes realizados, foi possível concluir que, para que houvesse a formação de gotas de líquido estáveis empregando esse modelo, o termo coesivo de pressão na equação de estado de van der Waals deveria, por hipótese, exceder a todas as outras forças presentes nas equações da conservação da quantidade de movimento (momentum) e da energia.

A simulação da formação de gotas a partir do modelo proposto por Nugent e Posch (2000) apresenta um problema inerente ao método SPH, conhecido como a instabilidade de tensão, que resulta na aglomeração de partículas em algumas situações específicas. Posteriormente, Meleán et al. (2004) realizaram algumas simulações numéricas na tentativa de minimizar/solucionar esse problema. Nele, os autores propuseram a adição de uma força viscosa artificial e de um termo fonte de energia nas equações padrão que governam o movimento e a transferência de energia, respectivamente, no método SPH. Esse modelo pôde ser utilizado com sucesso na simulação de sólidos elásticos ou quebradiços e os autores mostraram que, a partir da escolha adequada dos parâmetros essenciais para o termo de tensão viscosa artificial, o problema da aglomeração de partículas seria amenizado.

Lopez e Sigalotti (2006) estudaram as oscilações não-lineares de gotas líquidas viscosas e condutoras de calor, no vácuo e com gravidade zero. As gotas também foram modeladas como sendo um fluido de van der Waals. O trabalho restringiu-se ao estudo das oscilações, em duas dimensões espaciais, de pequenas e grandes amplitudes de gotas elípticas, inicialmente em uma condição estática. Nesse trabalho, os autores propuseram a utilização de um comprimento de suavização variável (que influencia diretamente na eficiência computacional e na acurácia dos resultados) na estimativa da massa específica, visando melhorar a acurácia e a estabilidade do método de partículas para essa classe de problemas.

Uma outra alternativa, para solucionar ou amenizar o problema da aglomeração de partículas, foi recentemente proposta por Yang et al. (2014). Nele, os autores utilizaram uma função kernel hiperbólica, que possui as suas derivadas segundas não-negativas, em conjunto com uma aproximação específica para a discretização do termo de condução de calor, presente na equação do balanço de energia. Tais modificações levam à obtenção de um campo de temperatura mais uniforme na gota e também possui a vantagem de acarretar num menor esforço computacional do que os modelos propostos anteriormente.

No presente trabalho, o método Smoothed Particle Hydrodynamics é utilizado na simulação numérica do problema da oscilação de uma gota de fluido de van der Waals em condições de microgravidade. Como diferencial, com relação aos trabalhos realizados previamente, incorporou-se a função kernel hiperbólica, introduzida por Yang et al. (2014), e realizou-se um estudo sobre a influência de alguns dos principais parâmetros físicos que afetam o comportamento oscilatório das gotas. Com relação ao termo de condução de calor, a sua discretização é feita de maneira análoga à sugerida por Yang et al. (2014). Diversos testes foram realizados variando-se a razão de aspecto, que é a relação entre os eixos principal e secundário da gota elíptica, o número de Reynolds e o coeficiente de condutividade térmica, que está diretamente relacionado ao número de Péclet, a fim de monitorar-se o movimento das gotas considerando diferentes conjuntos desses parâmetros.

Além disso, algumas simulações foram realizadas a fim de mostrar-se que com o modelo aqui proposto, usando a função kernel hiperbólica, foi possível obter gotas de líquido estáveis e com uma distribuição uniforme de partículas em todas as simulações, sem a necessidade de que nenhuma outra correção ou termo adicional fosse introduzido no método SPH.

2 Modelagem matemática

No problema da dinâmica da formação e da oscilação de uma gota de líquido, em um ambiente de vácuo, as equações que governam o deslocamento das partículas de fluido são a equação de conservação da massa (Slattery, 1999)

(1)

a equação que expressa a conservação da quantidade de movimento

(2)

e a equação do balanço da energia

(3)

onde é o tempo, é a massa específica, é a velocidade, é o tensor de tensão, representa as forças de corpo por unidade de massa, é a energia interna específica e é o vetor do fluxo de calor.

O tensor de tensão, para um fluido newtoniano, é escrito como (Slattery, 1999)

(4)

onde é a pressão interna, é o tensor identidade e é o tensor das tensões viscosas, dado por (Slattery, 1999)

(5)

onde e são os coeficientes de cisalhamento e de viscosidade em massa, respectivamente.

O vetor fluxo de calor pode ser escrito, da lei de Fourier, como (Slattery, 1999)

(6)

onde é o coeficiente de condutividade térmica e é a temperatura do fluido.

No caso de um fluido compressível, ou quase-incompressível, além do sistema de Eqs. (1)-(3) que governam o escoamento, necessita-se também de equações adicionais relacionando a pressão e a energia interna com a temperatura. Para um fluido de van der Waals têm-se que (Nugent e Posch, 2000)

(7)

(8)

onde , e . Sendo a constante de Boltzmann, m a massa da partícula, a uma medida da ação coesiva responsável pelas forças atrativas () de curto alcance, entre as moléculas vizinhas, e um parâmetro constante devido ao tamanho finito das moléculas, responsável pelas forças repulsivas (). Portanto, a Eq. (7) pode ser reescrita, em termos dos termos atrativo e repulsivo, na seguinte forma:

(9)

3 Smoothed Particle Hydrodynamics

De acordo com Liu e Liu (2003), alguns procedimentos básicos fazem parte da metodologia que deve ser empregada ao se aplicar o método SPH na resolução numérica de um determinado fenômeno físico (como, por exemplo, a dinâmica dos fluidos):

• Discretização do domínio físico, representando-o por um conjunto de partículas arbitrariamente distribuídas;

• Aproximação de uma função empregando a Representação Integral (ou Aproximação de Kernel) para a determinação das grandezas de campo do problema, tais como: a massa específica, a pressão e a velocidade;

• Substituição da representação integral por um somatório, ou seja, a Aproximação por Partículas, sobre todos os valores correspondentes das partículas vizinhas dentro do chamado suporte compacto da função de suavização. Essa etapa tem por finalidade gerar um sistema de equações diferenciais ordinárias (EDO’s) com dependência somente no tempo;

• Resolução do sistema de EDO’s utilizando um método numérico adequado para a obtenção dos campos de interesse do problema.

A aproximação de Kernel e a Aproximação por Partículas são os principais passos na obtenção de uma formulação numérica usando o método SPH e, portanto, ambas serão detalhadas na sequência.

3.1 Aproximação de kernel de uma função

A aproximação de kernel de uma função , no método SPH, tem como ponto de partida a identidade

(10)

onde é uma função da posição e δ é a função delta de Dirac, definida como

(11)

Na Eq. (10), é o volume do domínio de integração. Substituindo-se a função delta de Dirac pela função de suavização , com suporte compacto, a igualdade é trocada pela aproximação de kernel (Liu e Liu, 2003)

(12)

onde é o comprimento de suavização que define a influência ou o suporte compacto da função de suavização W.

Já a aproximação da divergência de uma função, , é obtida simplesmente pela substituição de por na Eq. (12), obtendo-se

(13)

onde as derivadas na integral são com relação à coordenada principal. Considerando-se a relação algébrica expressa como

(14)

a seguinte igualdade pode ser obtida

(15)

A primeira integral do lado direito da Eq. (15) pode ser convertida, usando o Teorema da Divergência, em uma integral sobre a superfície S do domínio de integração , de modo que

(16)

onde é o vetor normal unitário à superfície . Como a função de suavização possui um suporte compacto, o valor de na superfície do domínio de integração na Eq. (16) é nulo, caso o domínio de influência da função de suavização não intercepte a fronteira do domínio. Consequentemente, a aproximação de kernel para a divergência pode ser escrita como

(17)

Uma representação análoga pode ser obtida para o operador gradiente (Liu e Liu, 2003).

3.2 Aproximação por partículas

Após a discretização do domínio computacional, incluindo as suas fronteiras, empregando um número finito de partículas, a forma contínua da aproximação de kernel, expressa pela Eq. (12), pode ser escrita numa forma discretizada como um somatório (Liu e Liu, 2003), considerando-se que (onde , e representa o volume finito da partícula )

(18)

onde é o número total de partículas dentro da área de influência da partícula posicionada em . Esse procedimento do somatório sobre as partículas vizinhas é conhecido como a aproximação por partículas. Ele determina o valor da função, em uma dada partícula, como sendo dado pela média dos valores da função em todas as partículas que estão dentro do suporte compacto da função de suavização. Seguindo o mesmo procedimento da aproximação de kernel, a aproximação por partículas do divergente de uma função vetorial pode ser obtida como

(19)

onde o gradiente é avaliado na partícula . O mesmo procedimento pode ser empregado no que diz respeito ao gradiente de uma função (Liu e Liu, 2003).

3.3 Equações governantes na forma discretizada

Utilizando-se a aproximação de kernel, a aproximação por partículas e algumas identidades algébricas (Liu e Liu, 2003), o sistema de Eqs. (1)-(3) pode ser expresso, no formalismo do método SPH, pelas seguintes equações (Góes, 2016):

(20)

(21)

(22)

com os subscritos e denotando as partículas .

Neste trabalho, o método Leap Frog (Liu e Liu, 2003), utilizando um passo de tempo constante e igual a (escolhido de forma a garantir a condição de estabilidade numérica de Courant–Friedrichs–Lewy (CFL)), é utilizado na resolução numérica do sistema formado pelas Eqs. (20)-(22), acrescido da equação

(23)

para a determinação das partículas.

De acordo com Nugent e Posch (2000) o termo coesivo de pressão na Eq. (7), , deve ser considerado separadamente na forma discretizada das equações de balanço do momentum e da energia. Portanto, a Eq. (21) é reescrita como

(24)

onde . De forma análoga, os termos coesivos aparecerão explicitamente na equação de energia (Eq. (22)),

(25)

A escolha do comprimento de suavização H, nas Eqs. (24) e (25), é determinada levando-se em consideração alguns critérios de estabilidade. Em particular, gotas circulares e estáveis são formadas quando (Nugent e Posch, 2000), ou seja, quando assume-se que o alcance da interação das forças coesivas (atrativas) excede o de todas as outras forças que aparecem nas equações de movimento.

No método SPH, a determinação do padrão de interpolação é de extrema importância e isso deve ser feito com a utilização de uma função de suavização adequada para cada problema. Neste trabalho, conforme já mencionado, emprega-se uma função de suavização conhecida como kernel hiperbólico (Yang et al., 2014):

(26)

com , sendo a distância entre os pares de partículas e o fator de normalização que, para o caso bidimensional, assume o valor . A escolha dessa função de suavização será responsável pela obtenção de gotas estáveis com uma distribuição mais uniforme de partículas. Segundo Yang et al. (2014), o seu uso contribui de forma significativa para a atenuação dos efeitos da instabilidade de tensão.

Ainda em consonância com o trabalho de Yang et al. (2014), o termo de condução de calor foi discretizado na forma (Góes, 2016)

(27)

ou ainda, para uma condutividade térmica constante,

(28)

onde e a fim de evitar-se uma singularidade.

Tabela 1 – Parâmetros adotados nas simulações

Parâmetro	Valor
Massa das partículas ()
Temperatura crítica do fluido ()
Densidade crítica do fluido ()
Pressão crítica do fluido ()
Constante de Boltzmann ()
Comprimento de suavização ()

4 Resultados e discussões

Na Tabela 1 encontram-se listados os principais parâmetros adotados nas simulações numéricas, aqui realizadas, e que correspon- dem aos mesmos valores adimensionais que os utilizados por Nugent e Posch (2000) no estudo da formação de uma gota líquida cercada pelo seu vapor.

Os modelos de gota utilizados nas simulações do fenômeno de oscilação possuem um formato elíptico, obtido a partir da deformação de uma gota circular estável. Portanto, inicialmente será feita uma breve apresentação do problema que é utilizado como referência para o caso da oscilação, ou seja, o problema a partir do qual se obtém uma gota circular estável. A formação de uma gota, desconsiderando a existência de uma atmosfera externa, tem como condição inicial um quadrado de lado , relativamente menor que o domínio do problema , discretizado mediante o uso de partículas igualmente espaçadas entre si ) e que se encontram em repouso, conforme pode ser visto na Figura 1.

Figura 1 – Configuração inicial para o problema da formação de uma gota estável.

As constantes e presentes na Eq. (7) são obtidas considerando-se as seguintes relações (Bonilla e Herrera, 2006)

(29)

que são verificadas caso

(30)

(31)

De acordo com Meleán e Sigalotti (2005), a temperatura e a massa específica iniciais devem satisfazer às seguintes desigualdades

(32)

(33)

conhecidas como condições de estabilidade termodinâmica para que ocorra a formação de uma gota estável. A primeira desigualdade garante que o calor específico, o módulo de compressibilidade adiabática e o produto entre o coeficiente de expansão térmica e o coeficiente de Grüneisen sejam todas quantidades positivas. Já a segunda desigualdade garante que a pressão cinética, dada pelo primeiro termo do lado direito da Eq. (7), seja sempre positiva. A combinação dessas restrições com a positividade da massa específica, da temperatura absoluta e do coeficiente de Grüneisen definem o domínio de fase para esse modelo de fluido de van der Waals (Meleán et al., 2004).

Na Figura 2 são apresentadas as configurações, no estado de equilíbrio, para gotas com diferentes temperaturas iniciais (). Nessas simulações, foram considerados os seguintes valores: , e . A fim de respeitar as condições de estabilidade termodinâmica (Eqs. (32) e (33)), o seguinte valor para a massa específica inicial foi adotado: .

(a) . (b) .

Figura 2 – Configuração de equilíbrio para o problema da formação de uma gota, para diferentes valores da temperatura inicial.

Conforme pode ser observado nas Figuras 2(a)–2(d), todas as gotas formadas apresentam um formato circular, que corresponde ao estado de mais baixa energia (equilíbrio termomecânico). Entretanto, no que diz respeito à distribuição uniforme de partículas, vê-se que apesar da utilização do kernel hiperbólico, a uniformidade da distribuição de partículas torna-se mais difícil de ser obtida quanto mais próxima a temperatura inicial for da temperatura crítica do fluido em questão, . Portanto, todos os resultados que serão apresentados na sequência foram obtidos levando-se em consideração a gota formada cuja temperatura inicial é igual a 0,2 (Figura 2(a)).

O modelo empregado no presente trabalho permite que gotas estáveis, com uma distribuição uniforme de partículas, sejam obtidas para uma relação entre os valores de h e H igual a 1:2. Tal fato implica que o custo computacional para a determinação das partículas vizinhas na Busca Direta (Liu e Liu, 2003), na aproximação por partículas, não será elevado.

A evolução temporal da temperatura (), da pressão (), da energia cinética () e da energia interna () do sistema são apresentadas nas Figura 3(a)–3(d). Os resultados foram obtidos tomando-se uma média dos valores individuais de cada partícula. Conforme pode ser observado na Figura 3(a), a temperatura oscila bastante nos instantes de tempo até e, a partir de então, tende a assumir um valor constante de aproximadamente . Isso se deve ao fato do alto valor adotado para o coeficiente de condutividade térmica, , que garante um rápido ajuste de temperatura até que a condição de equilíbrio seja alcançada. A pressão, Figura 3(b), apresenta um comportamento semelhante ao da temperatura e os valores negativos obtidos são devidos à tensão superficial.

Já a energia cinética, Figura 3(c), e a energia interna, Figura 3(d), estão diretamente relacionadas. Conforme pode ser visto, a energia interna aumenta de acordo com a diminuição da energia cinética (em função da conservação de energia). A simulação é interrompida quando a energia cinética tende a zero e a energia interna tende a um valor constante negativo, indicando que a gota encontra-se em equilíbrio termomecânico ().

O valor teórico da tensão superficial de uma gota pode ser calculado através da equação de Laplace (Nugent e Posch, 2000).

Para uma gota de líquido estável, a equação de Laplace é dada por

(34)

onde é a tensão superficial, é a pressão no centro da gota, é a pressão de vapor, bem distante da gota, e é o raio da gota. Na Figura 4 é apresentado o gráfico da pressão no centro da gota variando com o inverso do raio para gotas de líquido de tamanhos diferentes. Os raios das gotas foram determinados calculando-se o centro de massa das mesmas e comparando-os com a posição da partícula mais externa. A curva de ajuste linear que melhor corresponde aos resultados obtidos é dada pela seguinte equação:

(35)

obtida com o uso da ferramenta Basic Fitting do Matlab®. Das Eqs. (34) e (35) obtém-se, mediante uma simples comparação, os valores e . Conforme a previsão teórica, a tangente da inclinação da reta no gráfico deve fornecer o valor da tensão superficial, que deveria ser constante para as gotas com diferentes valores do seu raio. Portanto, considera-se que, apesar das pequenas variações, os resultados numéricos apresentam uma boa concordância com a previsão teórica.

4.1 Oscilação das gotas de líquido viscosas

Devido à importância da tecnologia de processamento sem recipiente usada em aplicações espaciais, diversos estudos voltados para as oscilações não-lineares de gotas têm se concentrado ou em cálculos numéricos detalhados ou em experimentos com gravidade zero (Lopez e Sigalotti, 2006). Tal constatação reforça a importância do presente estudo acerca da oscilação de gotas viscosas de fluidos de van der Waals.

A gota circular estável obtida previamente, mostrada na Figura 2(a), é tirada do seu estado de equilíbrio e convertida para o formato elíptico por meio de um aplanamento homogêneo sob tensão de cisalhamento puro (Twiss e Gefell, 1990). Isso envolve uma transformação de coordenadas, implicando na preservação de área, dada por

(36)

onde os parâmetros e correspondem, respectivamente, à deformação e ao alongamento sofridos pela gota circular de modo a torná-la elíptica. A transformação expressa pela Eq. (36) é similar àquela que foi utilizada por Lopez e Sigalotti (2006) e irá converter a gota circular em uma elipse cujo eixo principal está alinhado com o eixo x, conforme as Figuras. 5(a) e 5(b). O eixo

(a) temperatura (b) pressão

Figura 3 – Variação temporal da temperatura, pressão, energia cinética e energia interna da gota.

principal e o eixo secundário das gotas elípticas possuem os comprimentos e , respectivamente, onde é o raio da gota circular estável.

Os modelos oscilatórios empregados são parametrizados pelos números adimensionais de Reynolds

(37)

pela razão de aspecto inicial

(38)

e pelo número de Péclet, que é uma medida da importância da condução de calor,

(39)

onde é o calor específico a pressão constante e é a tensão superficial dada pela equação de Laplace (Eq. (34)).

Com o intuito de estudar-se a influência desses parâmetros no comportamento oscilatório da gota de fluido, os seguintes testes foram realizados:

Figura 4 – Pressão no centro da gota em função do inverso do raio das gotas. “ ” representa a pressão no centro da gota e “–” representa a curva de ajuste obtida a partir dos resultados numéricos.

(a) razão de aspecto = 2 (b) razão de aspecto = 4

Figura 5 – Configuração inicial para o problema da oscilação de gotas de líquido viscosas, com diferentes razões de aspecto.

1. Variação da razão de aspecto inicial de a , para números de Reynolds na faixa: ;

2. Variação do número de Reynolds de a , para uma mesma razão de aspecto inicial; e

3. Variação do coeficiente de condutividade térmica de a , para dois valores distintos da razão de aspecto inicial: e , e para um mesmo valor do número de Reynolds.

Os valores para a massa específica e a pressão no centro da gota, no estado de equilíbrio (), são obtidos por interpolação a partir da Eq. (18), obtendo-se e , respectivamente. Novamente, o raio da gota é determinado calculando-se o centro de massa da mesma e comparando-o com a posição da partícula mais externa, de modo que . Logo, da equação de Laplace obtém-se , que apresenta uma variação menor do que 3% com relação ao valor determinado pelo ajuste de curva da Eq (35).

Na Figura 6 é apresentada a distribuição de partículas, para diferentes instantes de tempo, do problema de oscilação de uma gota de fluido de van der Waals, com (Figura 5(b)) e , durante o primeiro período de oscilação. Conforme pode ser observado, a gota oscila ao longo dos eixos e (Figuras 6(a)-6(h)), até que, no instante de tempo (Figura 6(i)), ela retorna à sua posição inicial. Porém, como ocorre uma dissipação de energia, devido às tensões viscosas (Silva, 2018), a razão de aspecto da gota vai diminuindo com o passar do tempo, até que ela atinja a condição de equilíbrio estável (formato circular).

(a) t = 10 (b) t = 20 (c) t = 30

(d) t = 40 (e) t = 50 (f) t = 60

(g) t = 70 (h) t = 80 (i) t = 90

Figura 6 – Distribuição de partículas de uma gota de fluido de vdW durante o primeiro período de oscilação, com e .

Na Figura 7 encontra-se o perfil de velocidade do sistema de partículas correspondente à gota mostrada na Figura 6, para diferentes instantes de tempo. A gota, nos modelos oscilatórios, parte de uma condição inicial na qual as partículas encontram-se todas em repouso. Como a gota não se encontra no estado de mais baixa energia, à medida que o tempo avança, as partículas começam a se mover, promovendo o movimento oscilatório da mesma ao longo dos eixos x e y.

O perfil de velocidade para t = 10 pode ser visto na Figura 7(a), ou seja, bem no início da simulação. Conforme pode ser visto, as partículas mais externas ao longo do eixo x são as que apresentam as maiores velocidades (em módulo) e os vetores estão direcionados para o centro da gota, o que causa o posterior alongamento da mesma ao longo do eixo y. Já na Figura 7(b), tem-se o momento de máxima deformação sofrida pela gota durante o primeiro período de oscilação. Para esse instante de tempo (t = 50), tem-se que a maioria das partículas do sistema retorna à condição inicial de repouso, já que ela não mais irá se alongar nessa direção.

Em um instante posterior (t = 70), mostrado na Figura 7(c), observa-se novamente que as partículas mais externas apresentam as maiores velocidades e os vetores apontam para a direção em que ocorrerá o alongamento. Por último, na Figura 7(d) é apresentado o momento em que a gota completa o primeiro período de oscilação (t = 90). Constata-se, para esse instante de tempo, que praticamente todas as partículas atingem a condição de repouso, antes de iniciar um novo ciclo oscilatório.

(a) t = 10 (b) t = 50

Figura 7 – Perfil de velocidade das partículas da gota para diferentes instantes de tempo.

Nas Figuras 8(a) e 8(b) são apresentadas as variações da razão de aspecto da gota em função do tempo para diferentes formatos elípticos iniciais, considerando um número de Reynolds igual a 9,09 (Figura 8(a)) e 90,88 (Figura 8(b)). Da observação dos gráficos das figuras constata-se que à medida que as razões de aspecto aumentam também aumentam a amplitude das oscilações, e o equilíbrio é alcançado após transcorrido um intervalo de tempo cada vez maior, embora ele não possa ser observado na Figura 8(b), devido à simulação ter sido interrompida em t = 500. A dissipação de energia ocorre mais rapidamente para os menores valores da razão de aspecto.

Levando-se em consideração a influência do número de Reynolds, é possível observar que, para uma mesma razão de aspecto inicial, quanto maior for o número de Reynolds mais pronunciado será o comportamento oscilatório da gota. Isso já era de se esperar, uma vez que o número de Reynolds está diretamente relacionado com o coeficiente de cisalhamento (η), um dos principais parâmetros responsáveis pela dissipação da energia do sistema (Silva, 2018).

Na Figura 9 tem-se o gráfico da variação da razão de aspecto em função do tempo transcorrido para diferentes valores do número de Reynolds. Conforme é do conhecimento geral, o número de Reynolds expressa a razão entre as forças de inércia e as forças de dissipação viscosas e, portanto, quanto menor for o seu valor, maior será a importância da dissipação viscosa, que é responsável pelo decaimento assimptótico das oscilações. À medida que o número de Reynolds aumenta, as forças viscosas já não são mais suficientes para impedir o aparecimento do movimento oscilatório e, quanto maior for esse número adimensional, mais lento será o decaimento da amplitude das oscilações.

Os resultados apresentados na Figura 9 também podem ser analisados do ponto de vista da mudança no regime de oscilações:

(a) Re = 9,72

(b) Re = 97,25

Figura 8 - Evolução temporal da razão de aspecto da gota para diferentes deformações iniciais.

para as oscilações apresentam um comportamento periódico (sub-amortecido); para , as oscilações apresentam um comportamento aperiódico, ou seja, a gota retorna rapidamente ao seu estado de equilíbrio circular; e, para , tem-se um comportamento típico de um sistema superamortecido, que corresponde a um modo de decaimento aperiódico mais lento, que domina os movimentos oscilatórios em tempos mais avançados. Para oscilações de amplitudes infinitesimais, os mesmos efeitos foram previstos por Prosperetti (1980).

Finalmente, a Tab. 2 contém os valores da razão de aspecto para diferentes valores do período de oscilação e dos coeficientes de condutividade térmica (), que correspondem a números de Péclet na faixa . Percebe-se que aumentando o valor do coeficiente de condutividade térmica, mais alongada é a forma da gota no final do período de cada oscilação para uma razão de aspecto inicial igual a . Já para as oscilações de grandes amplitudes (), o coeficiente de condutividade térmica só passa a exercer uma influência apreciável no comportamento oscilatório à medida que o tempo avança, ou seja, nos instantes de tempo iniciais o efeito é quase imperceptível.

É importante destacar que o coeficiente de condutividade térmica não tem um efeito tão significativo quanto a razão de aspecto inicial da gota e o número de Reynolds, que está diretamente relacionado com o coeficiente de cisalhamento. Isso se deve ao fato da gota encontrar-se em um ambiente de vácuo, ou seja, não existe troca de calor com o ambiente externo.

Fazendo-se uma análise conjunta de todos os parâmetros que influenciam a dinâmica do comportamento oscilatório de gotas sob condições de microgravidade, pode-se concluir que, para oscilações de baixas amplitudes, a dissipação de energia ocorre pelo efeito combinado da viscosidade e da condução de calor, enquanto que para oscilações de grandes amplitudes, a dinâmica oscilatória é predominantemente governada pela interação das forças inerciais e viscosas.

Por último, é feita uma breve análise da variação da tensão superficial de gotas submetidas a oscilações de baixas amplitudes, a

Figura 9 - Efeito do número de Reynolds nas oscilações de grandes amplitudes e para .

Tabela 2 – Evolução temporal da razão de aspecto da gota para diferentes coeficientes de condutividade térmica

		χ₀ = 1,5			χ₀ = 3,0
Período	κ = 1,0	κ = 5,0	κ = 10,0	κ = 1,0	κ = 5,0	κ = 10,0
τ₀	1,5000	1,5000	1,5000	3,0000	3,0000	3,0000
τ₁	1,3585	1,3955	1,3981	2,5528	2,5868	2,5868
τ₂	1,2198	1,2429	1,2594	2,2019	2,2514	2,2755
τ₃	1,0925	1,1299	1,1519	1,9302	1,9792	1,9981
τ₄	0,9934	1,0282	1,0557	1,7264	1,7637	1,8000
τ₅	0,9057	0,9463	0,9792	1,5566	1,6049	1,6415
τ₆	0,8462	0,8842	0,9085	1,4038	1,4915	1,5283

fim de comparar-se os resultados com aqueles obtidos pela equação de Laplace (Eq. (34)). Para oscilações de baixas amplitudes de gotas em duas dimensões, o período de oscilação ( ) pode ser obtido a partir da seguinte relação (Yang et al., 2014)

(40)

de modo que

(41)

onde denota a tensão superficial que é obtida no trabalho realizado por Rayleigh (1879).

Na Figura 10 é apresentado o gráfico da tensão superficial, obtida empregando-se a Eq. (41), para gotas com diferentes raios. As gotas simuladas partem de uma condição inicial na qual , que corresponde a uma oscilação de baixa amplitude. A linha tracejada corresponde ao valor médio da tensão superficial () obtida para as diferentes gotas. Comparando-se esse resultado com o que foi obtido pela equação de Laplace, obtém-se um erro relativo, dado por , de aproximadamente 0,001068, ou seja, menor do que 1%, mostrando uma boa concordância entre os valores obtidos com essas duas equações distintas, para o caso específico aqui tratado. Deve-se levar em consideração que erros são introduzidos quando do cálculo do centro de massa, do raio da gota e da determinação do período de oscilação.

Esses últimos testes foram realizados apenas com o intuito de demonstrar que as gotas de líquido de van der Waals e a tensão superficial podem ser eficientemente modeladas e calculadas através do emprego do método SPH na resolução do problema da dinâmica do movimento das gotas.

Figura 10 – Tensão superficial da gota, obtida pela Eq. 41, para gotas de tamanhos diferentes. “ ” representa a tensão superficial da gota e “– –” representa o valor médio dentre os resultados obtidos.

5 Conclusões

No presente trabalho, o método numérico livre de malhas Smoothed Particle Hydrodynamics (SPH) é utilizado na simulação do problema da oscilação de uma gota de fluido de van der Waals (vdW) sob condições de microgravidade. Os parâmetros adimensionais que governam o comportamento oscilatório da gota, cuja influência no movimento é estudada, são o número de Reynolds (), a razão de aspecto inicial () e o número de Péclet ().

Os resultados numéricos mostram que para oscilações de baixas amplitudes () a dissipação de energia ocorre de forma mais rápida, mesmo para valores elevados do número de Reynolds. Já para oscilações de grandes amplitudes (), é possível determinar para quais faixas do número de Reynolds ocorre a transição do comportamento oscilatório periódico () para o aperiódico (). Por último, conclui-se que o coeficiente de condutividade térmica, associado ao número de Péclet, afeta principalmente o comportamento dos modelos oscilatórios de baixas amplitudes.

Em relação a alguns aspectos numéricos e do modelo empregado no método SPH, vale a pena destacar que, com a implementação de modificações relativamente simples, como o uso do kernel hiperbólico e da discretização do termo de condução de calor baseado em uma expansão em série de Taylor, é possível obter gotas de líquido estáveis e com uma distribuição uniforme de partículas para todas as simulações contempladas. Além disso, o método mostra-se ser perfeitamente capaz de modelar satisfatoriamente as gotas de líquido de van der Waals e os seus respectivos efeitos de tensão superficial.

Agradecimentos

Os autores agradecem o apoio financeiro concedido pelas agências de fomento CAPES, CNPq e FAPERJ.

Referências

Anderson, D., Tannehill, J. C., Pletcher, R. H. (2016). Computational fluid mechanics and heat transfer. CRC Press.

Bonilla, B., Herrera, J. (2006). Revisando la ecuación de van der Waals. Revista Mexicana de Física E, 52(1), 65–77.

Gingold, R. A., Monaghan, J. J. (1977). Smoothed Particle Hydrodynamics: theory and application to non-spherical stars. Monthly Notices of the Royal Astronomic Society, 181, 375–389.

Góes, M. L. (2016). Estudo de interfaces livres empregando o método Smoothed Particle Hydrodynamics (SPH) e as equações de estado de van der Waals e de Martin. Tese de Doutorado, Universidade do Estado do Rio de Janeiro.

Liu, G. R., Liu, M. B. (2003). Smoothed Particle Hydrodynamics: a meshfree particle method. World Scientific.

Lopez, H., Sigalotti, L. D. G. (2006). Oscillation of viscous drops with Smoothed Particle Hydrodynamics. Physical Review E, 73(5), 051,201.

Lucy, L. B. (1977). Numerical approach to testing the fission hypothesis. The Astronomical Journal, 82, 1013–1024.

Meleán, Y., Sigalotti, L. D. G. (2005). Coalescence of colliding van der Waals liquid drops. International Journal of Heat and Mass Transfer, 48(19), 4041–4061.

Meleán, Y., Sigalotti, L. D. G., Hasmy, A. (2004). On the SPH tensile instability in forming viscous liquid drops. Computer Physics Communications, 157(3), 191–200.

Nugent, S., Posch, H. (2000). Liquid drops and surface tension with smoothed particle applied mechanics. Physical Review E, 62(4), 4968.

Prosperetti, A. (1980). Free oscillations of drops and bubbles: the initial-value problem. Journal of Fluid Mechanics, 100(2), 333–347.

Rayleigh, L. (1879). On the capillary phenomena of jets. Proc R Soc London, 29(196-199), 71–97.

Silva, G. R. (2018). Estudo da dinâmica da formação de gotas empregando o método Smoothed Particle Hydrodynamics e a equação de estado de van der Waals. Dissertação de Mestrado, Universidade do Estado do Rio de Janeiro.

Slattery, J. C. (1999). Advance Transport Phenomena. Cambridge University Press.

Twiss, R. J., Gefell, M. J. (1990). Curved slickenfibers: a new brittle shear sense indicator with application to a sheared serpentinite. Journal of Structural Geology, 12(4), 471–481.

Yang, X., Liu, M., Peng, S. (2014). Smoothed Particle Hydrodynamics modeling of viscous liquid drop without tensile instability. Computers & Fluids, 92, 199–208.