Solução computacional do problema da cavidade cúbica através das equações de Navier-Stokes tridimensionais

Computational solution of the cubic cavity problem using the tridimensional Navier-Stokes equations

ARTIGOS GERAIS

Solução computacional do problema da cavidade cúbica através das equações de Navier-Stokes tridimensionais

Computational solution of the cubic cavity problem using the tridimensional Navier-Stokes equations

Darlon Vasata; Guilherme Galante; Rogério Luis Rizzi; Reginaldo A. Zara¹

Centro de Ciências Exatas e Tecnológicas, Universidade Estadual do Oeste do Paraná, Cascavel, PR, Brasil

RESUMO

Este trabalho apresenta uma abordagem para a solução de problemas de escoamento tridimensional de fluidos com aplicação no problema da cavidade cúbica com tampa deslizante. A quantidade de movimento é modelada pelas equações de Navier-Stokes, as quais são discretizadas utilizando o método de diferenças finitas. A solução para o problema foi obtida utilizando o método SOR (Successive Over Relaxation) para resolver a equação de Poisson para pressão para cada passo de tempo. Apresentam-se os resultados de simulações numéricas na forma de gráficos da velocidade, pressão e superfícies de corrente. É observada a formação vórtices a partir do cantos da cavidade próximo à tampa deslizante. Os resultados generalizam a solução do problema da caixa bidimensional e mostram que as estratégias numérico-computacionais empregadas são válidas e possuem boa qualidade numérica.

Palavras-chave: Navier-Stokes, cavidade cúbica, método SOR.

ABSTRACT

This paper presents an approach to solving problems of three-dimensional fluid flow with application to the problem cubic cavity with sliding covering. The fluid flow is modeled by the Navier-Stokes equations which are discretized using the finite difference method. The solution to the problem was obtained using the SOR method (Successive Over Relaxation) to solve the Poisson equation for pressure for each time step. The results of numerical simulations are presented as graphs of velocity, pressure and surface current. It is observed the formation of tridimensional vortices rising from the corners of the cavity near the covering. The results generalize the solution of the two-dimensional case and show the numerical-computational strategies employed are appropriated and exhibit good numerical quality.

Keywords: Navier-Stokes, cubic cavity, SOR method.

1. Introdução

Escoamentos de fluidos podem ser encontrados em diferentes ramos da ciência e da tecnologia, envolvendo a indústria, a agricultura, as diversas engenharias, ciências ambientais entre outros [1, 2]. A área da ciência que investiga os fluidos e suas propriedades, chamada de mecânica dos fluidos, desenvolveu-se inicialmente de forma experimental, com o estudo do movimento dos fluidos em canais e tubos. O início da investigação de caráter matemático pode ser atribuído a Euler e sua descrição, através de equações, das relações entre o movimento dos fluidos e as forças envolvidas neste movimento. A hidrodinâmica introduzida pelas equações de Euler ganhou impulso no século XIX através dos trabalhos de Claude Navier, Simeon Poisson e George Stokes, com a proposição de um conjunto de equações conhecidas como equações de Navier-Stokes [1].

O modelo hidrodinâmico de Navier-Stokes é um conjunto de equações diferenciais que descrevem o transporte de massa, momento e energia. Quando resolvidas com as condições de contorno apropriadas, as soluções destas equações fornecem uma descrição detalhada do processo de escoamento, sendo que as propriedades macroscópicas do sistema tais como a velocidade de escoamento, as quedas de pressão e os coeficientes de transporte são calculadas como médias das quantidades microscópicas correspondentes. Estas médias devem ser tomadas em um volume pequeno quando comparado ao tamanho do sistema mas grande o bastante para que as equações de movimento sejam válidas. Assim, realizando médias sobre todos os pontos do espaço pode-se obter as grandezas macroscópicas do sistema as quais, em geral, podem ser relaciondas através da lei empírica de Darcy. No entanto, resolver as equações de movimento é um problema difícil, mesmo para os sistemas mais simples. Além disso, sistemas demasiadamente simplificados não são suficientemente precisos para a determinação das quantidades físicas relevantes da maioria das situações reais de interesse.

As equações de Navier-Stokes são equações diferenciais parciais não lineares para as quais soluções analíticas são conhecidas apenas para alguns casos específicos. Como a solução de problemas hidrodinâmicos importantes não podem aguardar pela solução analítica das equações matemáticas correspondentes, métodos alternativos de solução devem ser empregados. Muitas aplicações tecnológicas hoje essenciais foram desenvolvidas através de ensaios experimentais realizados em tanques de água ou túneis de vento. Porém, além de dispendiosos e de tomar tempo considerável para execução, os ensaios nem sempre são suficientes para dirimir dúvidas sobre a dinâmica do escoamento ou sobre as propriedades físicas do sistema.

A dinâmica dos fluidos computacional (DFC) é uma abordagem contemporânea para pesquisas em dinâmica dos fluidos [1]. Ela surgiu com o advento dos computadores e integra a física e a matemática com a computação de forma a proporcionar solução de vários problemas complexos. O emprego da DFC permite que problemas para os quais a realização de experimentos é economicamente inviável possam ser investigados numericamente sob as mais diferentes condições. Assim, a partir das informações fornecidas pela DFC, pesquisas teóricas podem ser validadas ou novas abordagens podem ser sugeridas e os testes experimentais necessários podem ser melhor planejados e executados [2].

Neste trabalho são construídos e implementados modelos matemáticos aplicados à DFC baseados nas equações de Navier-Stokes tridimensionais visando a solução numérica para o problema da cavidade cúbica com tampa deslizante. A complexidade matemática das equações resultantes impossibilita encontrar soluções analíticas a partir da análise teórica, havendo a necessidade da utilização de métodos numéricos de solução. Na seção 2 são apresentadas as condições físicas do fluido investigado, o conjunto de equações de Navier-Stokes resultante para tais condições. Estratégias de solução das equações são discutidas descrevendo o processo de discretização das equações, as aproximações e as condições de fronteira utilizadas. A seção 3 detalha o funcionamento do procedimento numérico utilizado para a solução das equações de escoamento. Na seção 4 é apresentado o problema da cavidade cúbica, os parâmetros da simulação são apresentados e os resultados obtidos são discutidos.

2. Equações de Navier-Stokes

As equações de Navier-Stokes são equações diferenciais parciais que permitem determinar os campos de velocidade e pressão de um fluido [1], descrevendo a física de um grande número de fenômenos de interesse acadêmico, tecnologicos e economicamente importantes, podendo ser usadas para modelar desde o escoamento de fluidos em dutos até correntes oceânicas e atmosféricas.

Para formular uma solução das equações de Navier-Stokes princípios básicos de conservação da massa, momento e energia devem ser levados em conta. Considerando o fluido Newtoniano, contínuo e incompressível as propriedades do escoamento podem ser descritas em termos de grandezas macroscópicas como velocidade, pressão, temperatura e densidade [3].

Para obter uma representação da conservação das propriedades físicas do sistema é necessário resolver a equação da continuidade, a qual é uma expressão matemática da conservação da massa. Essa equação pode ser escrita matematicamente como

onde ρ é a densidade do fluido e u, υ e w são as componentes da velocidade do fluido nas direções coordenadas x, y e z respectivamente.

As equações para a conservação da quantidade de movimento nas direções x, y e z podem ser escritas, respectivamente, como

onde p representa a pressão, e µ_H e µ_υ representam os coeficientes de viscosidade cinemática horizontal e vertical, respectivamente.

Uma estratégia de solução para as equações de Navier-Stokes é reescrevê-las, após certa manipulação algébrica, na forma de uma equação de Poisson para a pressão a cada passo de tempo. A representação da equação de Poisson em sua forma analítica é mostrada na Eq. (5), a seguir

onde g agrega as contribuições dos operadores explícitos para a velocidade e demais termos, sendo que p designa pressão [1].

Resolver as equações de Navier-Stokes significa encontrar, para um dado conjunto de condições iniciais e de contorno, as funções soluções u(x, y, z, t), v(x, y, z, t), w(x, y, z, t), p(x, y, z, t) e ρ(x, y, z, t) nas direções x, y e z respectivamente. Como, até o momento, soluções analíticas para as equações de Navier-Stokes ainda não são disponíveis a investigação numérica tem sido utilizada para solucionar as situações de interesse. Para isto, métodos numérico-computacionais têm sido desenvolvidos e testados em diferentes condições para análise de acurácia e condições de aplicação.

2.1. Discretização do domínio e das equações

Para a solução do problema da cavidade cúbica, utilizando-se as equações de Navier-Stokes tridimensionais, empregou-se o método de diferenças finitas. Na discretização das derivadas espaciais utilizou-se a malha de Akarawa tipo C. Nesse tipo de malha o posicionamento das variáveis tem um arranjo deslocado, que inibe as oscilações de pressão decorrentes da independência entre a pressão e a velocidade no mesmo ponto quando da solução da equação de conservação da massa. Para a abordagem de solução empregada esta escolha tem se mostrado adequada, já que alguns métodos que resolver as equações de Navier-Stokes para a pressão em malha co-localizada nem sempre garantem conservação da massa nas células da malha. Na malha descentrada as velocidades são armazenadas nas faces enquanto a pressão e as outras propriedades do fluido são armazenadas no centro do elemento da malha, como mostrado na Fig. 1.

Em diversas abordagens encontradas na literatura para soluções de escoamentos, as velocidades u, υ e w e a pressão p são representadas no mesmo ponto, ao centro de cada célula da malha. Porém, essa representação pode levar a instabilidades na solução numérica do problema, existindo assim a necessidade de uma nova representação física das velocidades e pressão na malha [4]. Assim, de modo a evitar este tipo de inconveniente, neste trabalho utiliza-se a abordagem de malha tridimensional deslocada.

A malha tridimensional deslocada consiste em representar a velocidade u nas paredes inferior e superior, a velocidade υ nas paredes esquerda e direita, e velocidade w nas paredes de frente e fundo e a pressão p no centro da célula, como mostrado na Fig. 1.

Dessa forma, a Eq. (2) é discretizada no ponto (i + 1/2, j, k) utilizando as seguintes aproximações para as derivadas espaciais e

para a derivada temporal

Porém nas Eqs. (6), (7) e (8) os termos , e não estão definidos na malha. Estes valores devem ser calculados através da interpolação dos valores de u, v e w das faces das células, como e

Nesse caso, ainda sem substituir os termos que não estão definidos na malha, a versão discreta da Eq. (2) pode ser escrita como

Reorganizando os termos da Eq. (18) de forma a isolar obtém-se

onde representa-se como os termos convectivos no passo de tempo n, já discretizados

Sendo ν_H = µ_h/ρ e ν_υ = µ_υ/ρ os coeficientes de viscosidade cinemática, para as componentes nas direções x, y e z, respectivamente, representa os termos viscosos, também já discretizados para o passo de tempo n

Substituindo na Eq. (20) as aproximações para os termos não definidos na malha

Torna-se conveniente, para efeitos da discretização, isolar os termos convectivos, viscosos e a velocidade para separar os operadores explícitos a serem utilizados em cálculos posteriores, sendo sua soma representada como . Tem-se então

Agora, a Eq. (19) pode ser expressa de forma mais suscinta por

Analogamente, a Eq. (3) é discretizada no ponto (i, j + 1/2, k) utilizando aproximações análogas, obtendo

onde

sendo e expressões análogas às Eqs. (20) e (21), porém para o ponto (i, j + 1/2, k). Da mesma forma, a Eq. (4) é discretizada no ponto (i, j, k +1/2) utilizando aproximações análogas, obtendo

onde

e sendo e expressões análogas às Eqs. (20) e (21), porém estas para o ponto (i, j, k +1/2).

2.2. Aproximação para pressão

Diferentemente das equações das velocidades a equação para a pressão é discretizada no ponto (i, j, k). Aplicando a Eq. 1 para o ponto (i, j, k) e utilizando diferenças centrais de segunda ordem, obtém-se

ou seja,

A Eq. (24) implica, por um deslocamento no índice i, como

tem-se também que a Eq. (25) implica

que e a Eq. (27) implica que

Substituindo as Eqs. (31), (32) e (33) na Eq. (30), tem-se a equação de Poisson para a pressão tridimensional discretizada, como

2.3. Condições de fronteira

As condições de fronteira são artifícios utilizados para representar a interação do sistema com o meio externo. Consiste em adicionar valores nas células de borda do domínio computacional, sendo assim denominadas células de fronteira.

Neste trabalho, as células de fronteira aplicam a técnica da reflexão, podendo comportar-se das formas escorregadia, ou fronteira de simetria e não escorregadia ou fronteira sólida. Assim, as condições de fronteira, pela técnica de reflexão são:

• Para células de fronteira escorregadias nas bordas laterais da malha: O valor referente a u da célula interna vizinha recebe o valor de u da célula de fronteira, com seu sinal invertido.

• Para células de fronteira não-escorregadias nas bordas laterais da malha: O valor referente a u da célula interna vizinha recebe o valor de u da célula de fronteira.

• Para células de fronteira escorregadias nas bordas superior ou inferior da malha: O valor referente a υ da célula interna vizinha recebe o valor de υ da célula de fronteira, com seu sinal invertido.

• Para células de fronteira escorregadias nas bordas superior ou inferior da malha: O valor de υ da célula interna vizinha recebe o valor de υ da célula de fronteira.

3. Algoritmo de solução

O procedimento para o cálculo das variáveis e da pressão no nível de tempo n + 1 consiste em resolver os passos descritos na Fig. 2 e descrito a seguir.

Inicialmente é feito todo pré-processamento, que inclui a leitura dos arquivos de entrada são e a inicialização das variáveis e parâmetros da simulação.

Na sequência define-se as condições de fronteira, de acordo com o descrito na seção 2.3. Nos dois passos subsequentes, efetua-se as operações relacionadas ao cálculo das variáveis F, G e H no interior e nas bordas do domínio.

No quarto passo resolve-se o lado direito da Eq. (34), que envolve o cálculo das pressões, para todos os pontos internos da malha, gerando um sistema de equações. Para resolver tal sistema utiliza-se o método numérico SOR (Successive Over Relaxation). Assim, a Eq. (35) deve a ser aplicada para cada ponto interno da malha

onde (k) representa a iteração do SOR, β = , γ = , w é um fator de relaxação, que deve estar entre 0 < w < 2, e

A abordagem de solução utilizada, que é obter as componentes da velocidade a partir dos valores atualizados da pressão, obriga a que em cada iteração do SOR seja efetuada uma correção da pressão com o intuito de estimar o campo de velocidades, de forma a satisfazer a equação da continuidade e obter uma representação física consistente [1].

Considerando que não há no escoamento dentro da cavidade, a pressão deve ser normalizada após cada iteração do método SOR de modo a se evitar que os valores da pressão cresçam ou diminuam arbitrariamente. A normalização consiste em manter constante o valor da pressão em uma determinada célula da malha. Como manter um valor constante pode ocasionar instabilidades na solução, uma alternativa é tomar uma determinada célula como referência e aplicar a correção (Eq. (37)) a todos os pontos da malha a cada iteração do método SOR, indicada pelo índice k. Repare que utilizando a técnica de correção, o valor da pressão na célula de referência é mantido constante em 0 [1].

O método apresentado pode ser aplicado a diversos níveis de tempo, incrementando seu passo de tempo t até que um determinado critério de parada seja satisfeito. Neste artigo o critério de parada escolhido é

onde NI, NJ e NK representam o número de células nas direções x, y e z, respectivamente. A escolha deste critério deve-se a que, quanto menor for o valor da soma do módulo das velocidades internas, mais perto estamos do estado estacionário. A escolha do valor ε deve ser um valor tão pequeno que faça com que a diferença entre as duas soluções seja desprezível. É razoável adotar valores ε 10^-7.

Uma vez finalizada a resolução dos sistemas de equações, é verificado se a simulação atingiu um estado estacionário. Em caso negativo, inicia-se uma nova iteração no tempo t. Em caso positivo encerra-se a aplicação e gera-se os gráficos para a visualização de resultados, como pode ser visto na seção 4.1.

4. O problema da cavidade cúbica

O problema da cavidade cúbica é uma extensão para a dimensão Z do problema de uma cavidade quadrada resolvida no plano XY, o qual é um problema de escoamento bidimensional amplamente estudado e sobre o qual é possível encontrar muitos trabalhos na literatura [5-8]. Este problema é frequentemente utilizado para validações de precisão, condições de contorno e eficiência de soluções numéricas.

No problema da cavidade cúbica um fluido encontra-se inicialmente em repouso, contido em uma caixa de formato cúbico de dimensão linear L, preenchendo a caixa por completo. As paredes da caixa são rígidas e impermeáveis, as paredes laterais e inferior são não-escorregadias e a superior é móvel e escorregadia. Instantaneamente, a tampa superior do cubo desloca-se com uma velocidade u₀ ao longo da coordenada x. O resultado esperado é a movimentação do fluido no interior da caixa, devido ao contato do mesmo com a parede superior em movimento. As equações de Navier-Stokes são resolvidas numericamente determinando-se para cada instante de tempo o campo de velocidade e pressão, permitindo obter uma descrição completa da movimentação do fluido no interior da caixa. Na próxima seção os parâmetros da simulação são apresentados e os resultados obtidos são discutidos.

4.1. Resultados

O algoritmo para o cálculo tridimensional da pressão e das velocidades foi implementado utilizando a linguagem de programação C com o compilador GCC (Gnu Compiler Collection) sob sistema operacional Linux, kernel 2.6.5. O objetivo é simular o escoamento de fluidos em um determinado domínio computacional tridimensional regular através da solução do conjunto des equações de Navier-Stokes associadas à caixa cúbica. Para que a solução possa ser visualizada, o software disponibiliza arquivos de saída com dados sobre os valores da pressão e velocidades, que podem ser analisados através da ferramenta VisIT [9].

Neste trabalho são apresentados resultados da simulação para problema da caixa cúbica de lado L = 40 efetuada com o seguinte conjunto de como parâmetros de entrada: Δt = 1.0 ×10^-1, Δx = 1.0 ×10^-1, Δ y = 1.0 ×10^-1, Δz = 1.0 ×10^-1, u₀ = 1.541667 ×10^-4, ν_H = 1.850 ×10^-3, ν_υ = 1.850 ×10^-4. O coeficiente de relaxação empregado no SOR foi w = 0.7 e, como critério de convergência, admitiu-se que a diferença entre grandezas físicas não fosse superior a = 1.0 ×10^-5.

Os resultados obtidos são apresentados na forma de gráficos tridimensionais para o deslocamento do fluido em dois instantes (iterações do método) de tempo o que permite compreender a dinâmica do fluido no interior da caixa nas três direções de interesse. Para isto são utilizados gráficos tridimensionais para os campos de velocidade u, v e w (respectivamente nas direções x y e z), pressão no interior do fluido e para as superfícies de corrente. Estas últimas mostram a trajetória do fluido no domínio.

A Fig. 3 mostra as componentes da velocidade nas direções X, Y e Z, denotadas por u, v e w respectivamente, em uma escala de tons na qual os tons mais claros correspondem a velocidades mais baixas, conforme indicado nas legendas. Esta figura permite descrever o movimento do fluido pois os valores de u, v e w retratam os valores da velocidade do fluido ao longo de cada um dos eixos, enquanto seus sinais indicam o sentido o movimento ao longo dos eixos. Ao longo do eixo X, nas regiões mais próximas às paredes a velocidade do fluido é menor devido à reflexão nas paredes da caixa. De fato, próximo às paredes observa-se valores negativos para a componente u da velocidade indicando a colisão do fluido com as paredes localizadas no eixo Z. Ao analisar, em conjunto com a componente u, o valores e a direção da componente da velocidade ao longo do eixo Y (υ), obtém-se uma imagem bidimensional do movimento: o fluido produz vórtices no plano XY ao colidir com as paredes da caixa, sendo que, em uma das paredes a componente υ da velocidade assume valores negativos. Os resultados obtidos para o movimento do fluido no plano XY são consistents com outros trabalhos disponíveis na literatura para o caso bidimensional. Resultados para o caso tridimensional são escassos sendo que, muitas vezes, o caso tridimensional é analisado através de extrapolações do caso em duas dimensões. Como na simulação aqui apresentada as três componentes da velocidade são numericamente calculadas, é possível também analisar a componente da velocidade ao longo do eixo Z. Estas componentes podem ser visualizadas na Fig. 3. Os resultados mostram que no plano XZ também ocorre a formação de vórtices, porém ao contrário do plano XY onde dois os vórtices possuem o mesmo sentido, no plano XZ são formados quatro vórtices a partir dos cantos superiores da caixa.

A partir das componentes da velocidade é possível construir o gráfico da velocidade resultante do fluido no interior da caixa, conforme mostrado na Fig. 4. Pode-se perceber que a interação entre os vórtices formados nos planos XY e XZ produz um vórtice central o qual forma-se a partir dos cantos superiores da caixa. Este efeito também pode ser visualizado a partir das superfícies de correntes na Fig. 5. As superfícies circulares mostradas na figura indicam alguns pontos os valores das velocidades u e υ possuem o mesmo módulo. Através destas, é possível acompanhar o movimento no interior da caixa e visualizar o vórtice interior. A figura também ilustra a influência da componente vertical da velocidade sobre o sistema. Nota-se ainda que os valores da componente w da velocidade possui valores consideravelmente menores do que as demais velocidades.

Os resultados mostrados na Fig. 5 também ilustram como as componentes da velocidade são dependentes dos valores das pressões no sistema, evidenciado matematicamente pelas Eqs. (24), (25) e (27). Comparando com as Figs. 3 e 5 é possível visualizar que os movimentos circulares internos que ocorrem a partir dos cantos superiores da caixa, onde estão os pontos que contêm os valores extremos das pressões. Ressalta-se ainda que estes pontos são justamente onde a força exterior que rege o movimento na tampa da caixa é aplicada. Pode-se observar ainda que nos pontos onde acontece a colisão do fluido com a parede a pressão é maior, como mostra a Fig 6. Percebe-se também que no momento de estabilidade do sistema a pressão é distribuída em todo o sistema, fazendo com que os valores de tais variáveis sejam menores nos cantos superiores.

5. Conclusões

Neste trabalho foi apresentada uma abordagem para a solução de problemas de escoamento do fluido com atenção ao problema da cavidade cúbica com tampa deslizante. As equações de Navier-Stokes que governam o movimento do fluido foram discretizadas utilizando o método de diferenças finitas e resolvidas utilizando o método de relaxação SOR (Successive Over Relaxation) para a equação de Poisson para pressão. Os resultados das simulações numéricas foram apresentados na forma de gráficos para a velocidade e pressão no fluido.Além disso, foram apresentadas as superfícies de corrente que permitem construir a representação do movimento do fluido na cavidade, na qual percebe-se a formação de vórtices tridimensionais a partir dos cantos da cavidade próximos à tampa.

[9] VisIt é uma ferramente livre e interativa para visualização e análise gráfica de dados científicos. Permite ao usuário gerar a visualição de seus dados, compor animações temporais e manipulá-los de forma gráfica.

Recebido em 9/4/2010; Aceito em 27/5/2011; Publicado em 8/7/2011

• [1] A. Fortuna, Técnicas Computacionais para Dinâmica dos Fluidos: Conceitos Básicos e Aplicações (Editora EDUSP, São Paulo, 2000).
• [2] R.L. Rizzi, Modelo Computacional Paralelo para a Hidrodinâmica e para o Transporte de Substâncias Bidimensional e Tridimensional Tese de Doutorado, UFGRS, 2002.
• [3] J.M. McDonough, Lectures on Computational Numerical Analysis of Partial Differencial Equations, Disponível em http://www.engr.uky.edu/~acfd/, acessado em 7/2/2007.
• [4] J.D. Anderson, Computational Fluid Dynamics (McGraw-Hill, New York, 1995).
• [5] A.C.H. Marques e J.L. Doricio, in: XXVIII Congresso Brasileiro de Matemática Aplicada e Computacional, São Paulo, 2005.
• [6] A.C.H. Marques J.L. Doricio, in: XXVI Iberian Latin American Congress on Computational Methods in Engineering, Guarapari, 2005, p. 411.
• [7] C.H. Bruneau and M. Saad, Computers and Fluids, 35, 326 (2006).
• [8] T.P. Chiang and W.H. Sheu, Computational Mechanics 20, 379 (1997).

Datas de Publicação

• Publicação nesta coleção
  02 Set 2011
• Data do Fascículo
  Jun 2011

Histórico

• Recebido
  09 Abr 2010
• Aceito
  27 Maio 2011