Uma nova metodologia para a simulação de escoamentos de fluidos viscoelásticos

A new approach for simulation of viscoelastic fluid flows

ARTIGO TÉCNICO CIENTÍFICO

Uma nova metodologia para a simulação de escoamentos de fluidos viscoelásticos

A new approach for simulation of viscoelastic fluid flows

André R. Muniz; Argimiro R. Secchi; Nilo S. M.Cardozo

Departamento de Engenharia Química, UFRGS

^{Endereço para correspondência}

RESUMO

É proposta neste trabalho uma nova metodologia para resolução das equações governantes de fluidos viscoelásticos, baseada no método dos volumes finitos, usando o arranjo co-localizado para as variáveis e malhas estruturadas. São utilizadas aproximações de alta ordem para os fluxos lineares e não-lineares médios nas interfaces dos volumes, e para os termos não-lineares que surgem da discretização das equações constitutivas. Nesta metodologia, os valores médios das variáveis nos volumes são usados durante todo o procedimento de resolução, e os valores pontuais são obtidos ao final, através da deconvolução dos valores médios. O sistema de equações discretizadas é resolvido de forma simultânea, pelo método de Newton. A metodologia é exemplificada para um problema clássico em mecânica de fluidos computacional, o escoamento stick-slip, usando como equação constitutiva o modelo de Oldroyd-B. As soluções obtidas apresentaram boa precisão, sendo livres de oscilações mesmo em regiões de grandes gradientes das variáveis.

Palavras-chave: Fluidos viscoelásticos, simulação, método dos volumes finitos.

ABSTRACT

In this work, a new methodology to solve the governing equations of viscoelastic fluid flows is proposed. This methodology is based on the finite-volume method with co-located arrangement of the variables, using high-order approximations for the linear and nonlinear average fluxes in the interfaces and for the nonlinear terms resulting from the discretization of the constitutive equations. In this methodology, the average values of the variable in the volumes are used during the resolution, and the point values are recovered in the post-processing step by deconvolution of the average values. The nonlinear equations, resulting from the discretization technique, are solved simultaneously, using the Newton's method. The solutions obtained are oscillation-free and accurate, as can be seen in the solution of the stick-slip flow, used as an illustrative example.

Keywords: Viscoelastic fluids, simulation, finite-volume method.

Introdução

Nos dias de hoje, a dinâmica de fluidos computacional (CFD – computational fluid dynamics) vem sendo cada vez mais utilizada em diversos segmentos na indústria (automotiva, aeroespacial, processos químicos, geração de energia, metalurgia, etc.), para os mais diversos fins. Também nas indústrias de transformação de polímeros verifica-se um interesse crescente no uso de softwares de CFD. Um exemplo importante de aplicação nesta área é o projeto de equipamentos, onde estes softwares reduzem a necessidade de execução de experimentos e criação de protótipos em escala de bancada, tarefas que consomem muito tempo e envolvem altos custos. Os mesmos podem ainda ser empregados para a otimização de um processo já existente, visando o aumento da produção e/ou melhoria da qualidade do produto, entre muitas outras possibilidades.

Na concepção de um software deste tipo, dois pontos fundamentais devem ser focados: a modelagem do problema físico e a resolução das equações do modelo.

O procedimento de modelagem consiste na descrição matemática do problema físico a ser analisado. No caso de escoamentos de fluidos, o modelo matemático consiste nas equações de conservação (massa, energia e quantidade de movimento), nas condições iniciais e de contorno, e em uma equação constitutiva mecânica que estabeleça a relação entre o campo de tensões e o campo de velocidades no escoamento.

Para escoamentos de fluidos poliméricos, o ponto chave é a seleção de uma equação constitutiva que represente adequadamente o comportamento reológico de fluidos poliméricos, o qual é muito complexo^[1,2]. Em simulações numéricas tem-se dado preferência ao uso de modelos de fluido Newtoniano generalizado e modelos diferenciais e integrais não-lineares para fluidos viscoelásticos. Os primeiros levam em conta apenas a dependência da viscosidade com a taxa de deformação, deixando de lado os efeitos decorrentes da elasticidade do fluido. Já os modelos diferenciais e integrais não-lineares para fluidos viscoelásticos permitem contemplar uma grande variedade das características reológicas apresentadas pelos fluidos poliméricos. Em relação aos modelos integrais, os modelos diferenciais são os mais usados em simulações, pois são mais fáceis de serem implementados e resolvidos, comparados com os integrais, e levam a predições consistentes, no mínimo, qualitativamente^[1].

O outro ponto, e não de menor importância, é a etapa de simulação, que consiste na resolução numérica das equações do modelo matemático. Para tanto, deve-se ter um código computacional que leve a resultados confiáveis de forma eficiente. Existem, na literatura, diversos métodos numéricos para a simulação de escoamentos viscoelásticos^[3-7], cada qual com uma determinada combinação de características dentre as muitas desejadas. Não existe um método de aplicação global, sendo que há uma busca constante no desenvolvimento e aperfeiçoamento dos diferentes métodos, no que se refere a precisão das soluções obtidas, facilidade de implementação e de entendimento, versatilidade e eficiência.

Este trabalho apresenta uma nova metodologia, baseada no método dos volumes finitos, para a simulação de escoamentos de fluidos viscoelásticos, usando equações constitutivas diferenciais não-lineares. Esta metodologia é descrita em maiores detalhes a seguir, mostrando quais as vantagens que possui em relação à metodologia convencional usada em volumes finitos. Ao final, são mostrados resultados obtidos para o escoamento de um fluido viscoelástico na geometria "stick-slip", usando o modelo de Oldroyd-B.

Modelo Matemático

As equações governantes do escoamento estacionário, incompressível e isotérmico de um fluido viscoelástico, caso considerado neste trabalho, são as equações de conservação de massa e de quantidade de movimento. Quando adimensionalizadas, estas equações assumem, respectivamente, a seguinte forma:

onde n é o vetor velocidade, p a pressão, Re = rUL/h₀ o número de Reynolds, r a massa específica e U e L, a velocidade e comprimento característicos do problema, respectivamente, h_V = h_v / (h_v +h_e) é a contribuição viscosa para a viscosidade total h₀ = h_v + h_e, que é dada pela soma das "contribuições" viscosa e elástica. O tensor das tensões é dado pela soma de uma contribuição newtoniana, dada pela lei da viscosidade de Newton, , com, sendo a taxa de deformação, e uma polimérica , para a qual se deve utilizar uma equação constitutiva adequada.

A metodologia apresentada pode, a princípio, ser empregada para qualquer modelo diferencial não-linear. Pelo fato de ser muito utilizado na literatura, é utilizado para descrição da metodologia, o modelo diferencial de Oldroyd-B:

onde We = lU/L é o número de Weissenberg (adimensional), l um parâmetro da equação constitutiva, chamado de tempo de relaxação, h_E = h_e / (h_v +h_e) é a "contribuição elástica" na viscosidade, tal que h_E + h_V = 1 e é a derivada convectiva do tensor , dada por:

Metodologia Numérica

A metodologia de resolução do sistema de equações do modelo, composto por equações diferenciais parciais, é formada pelas etapas de discretização das equações, aproximação dos fluxos médios e resolução das equações discretizadas.

Discretização das equações

Para a discretização das equações utiliza-se o método dos volumes finitos^[8], que consiste em integrar as equações diferenciais parciais em volumes elementares do domínio. Chega-se, assim, a um sistema de equações algébricas, escritas em termos dos fluxos médios nas interfaces dos volumes.

Neste trabalho serão considerados problemas bidimensionais, usando malhas estruturadas. Neste caso, o valor médio do fluxo f de uma quantidade f, para uma interface i de um dado volume, paralela à direção y, é definido como:

O valor médio da quantidade f no volume é definido por:

sendo que nestas expressões o par (i, j) representa um vértice do volume considerado e o par (i+1/2,j+1/2) representa o centróide deste volume.

Geralmente, tomam-se aproximações de 2ª ordem para os fluxos médios nas interfaces e para os valores médios nos volumes. Desta forma, o fluxo médio é igual ao fluxo no ponto central da interface e o valor médio da variável no volume é igual ao valor no ponto central, de forma a se trabalhar apenas com valores pontuais (nós e centro das interfaces). Os fluxos nos pontos centrais das interfaces são obtidos a partir de valores das variáveis nos centros dos volumes, através de esquemas de interpolação. Com o uso desta aproximação de 2ª ordem para os fluxos médios, mesmo que se utilizem esquemas de interpolação de ordem mais elevada, a ordem global permanece 2, não havendo ganho significativo com o uso destas^[9].

Na metodologia proposta, trabalha-se com valores médios no volume das variáveis^[10,11], sendo usados esquemas de interpolação que relacionam os fluxos médios nas interfaces com os valores médios das variáveis no volume. Com isto, a utilização de aproximações de alta ordem para os fluxos médios nas interfaces se torna mais simples, evitando a necessidade de integrações dos fluxos nas interfaces para obtenção de aproximações de alta ordem^[10].

Para exemplificar o procedimento de discretização utilizado, são mostradas a seguir, a equação de conservação da quantidade de movimento, na direção x, e a equação constitutiva para o tensor das tensões, na forma discretizada.

Aproximação dos fluxos nas interfaces e dos termos não-lineares

Os valores dos fluxos médios nas interfaces e os valores das médias de produtos no volume devem ser expressos em função do valor médio das variáveis no volume que, como já dito, são as variáveis do problema a ser resolvido.

Para a determinação dos fluxos médios nas interfaces, são usados esquemas de interpolação de alta ordem, que levam a soluções com melhor precisão, possibilitando a utilização de malhas mais grosseiras, reduzindo o custo computacional da simulação. Detalhes da metodologia de obtenção destes esquemas de interpolação podem ser encontrados nos trabalhos de Pereira et al.^[10] e de Kobayashi ^[11]. A grande redução do tempo computacional que se tem com o uso desta metodologia, comparada com a baseada em valores pontuais, foi analisada por Muniz^[12], para alguns exemplos de aplicação.

A notação para os índices mostrada na figura 1 é usada para descrever os esquemas de interpolação que serão mostrados na seqüência.

Para a determinação dos fluxos advectivos , utiliza-se uma aproximação de Lagrange de 3ª ordem, que leva em conta a direção do escoamento, tomando dois volumes a montante e um a jusante da interface:

Para os fluxos difusivos , utiliza-se uma aproximação de Lagrange centrada de 4ª ordem:

Como é bem sabido, o uso de esquemas de alta ordem para os fluxos advectivos pode levar a oscilações na solução em regiões com gradientes elevados^[9,13]. Estas oscilações, além de levar a soluções irreais, podem levar a instabilidade numérica, causando problemas de convergência.

Para contornar este problema, esquemas de 3ª ordem, para diferentes estênceis , foram usados para formar um esquema WENO (Weighted-Essentially-Non-Oscillatory)^[14,15] . Um esquema WENO consiste em uma combinação ponderada de esquemas de interpolação de mesma ordem de aproximação, mas com diferentes estênceis, levando a um esquema com ordem superior. Sua característica principal é evitar a presença de oscilações irreais na solução em regiões de elevados gradientes ou descontinuidades. Este tipo de esquema é amplamente usado em problemas de diferenças finitas, mas pouco explorado em volumes finitos.

Para aproximações de Lagrange de 3ª ordem para os fluxos advectivos em uma interface genérica i paralela a y, considerando v_x > 0, três possíveis estênceis levam aos seguintes esquemas f_k:

Desta forma, o fluxo na interface i é dado por uma combinação convexa destes esquemas:

onde r é o número de estênceis e w _k o peso de cada esquema no resultante, calculado de acordo com a suavidade na solução junto à região de grande variação ou próximo a pontos de singularidade^[14,15] e para resultar em uma aproximação de 5 ª ordem.

Para finalizar a etapa de aproximações, devem ser obtidos os fluxos não-lineares nas interfaces, como os termos , por exemplo, e as médias no volume de produtos entre variáveis ou entre variáveis e derivadas, tais como . Estas aproximações merecem uma atenção especial. Na metodologia convencional, a aproximação de um produto de variáveis é dado pelo produto das aproximações de cada termo, resultando em uma aproximação de 2ª ordem:

Para obter uma alta ordem global, uma aproximação de mais alta ordem deve ser usada também para estes termos não-lineares. Para isto, utiliza-se a idéia desenvolvida por Pereira et al.^[10], que consiste em calcular os termos de ordem inferior da aproximação dada pelas equações (13) e (14). Isto é feito ao comparar as expansões em série de Taylor para uma ordem desejada, de ambos os lados destas equações. Assim, para uma aproximação de 4ª ordem dos fluxos não lineares nas interfaces, os termos resultantes são:

E para a média dos produtos de variáveis no volume, tem-se:

Torna-se necessário, então, obter as aproximações dos termos do lado direito da equação, que são expressos em termos das médias dos volumes na vizinhança.

Condições de contorno

As condições de contorno são impostas naturalmente nas equações discretizadas. Quando se conhece o valor do fluxo ou das variáveis no contorno, estes são diretamente substituídos na equação, como por exemplo, velocidade nula nas paredes, tensão de cisalhamento nula em linhas de simetria, entre outros casos. Em situações, onde o valor no contorno deve ser calculado a partir de valores internos, são usados esquemas de interpolação, com a mesma ordem de aproximação dos fluxos nas interfaces internas, de modo a manter a ordem global. Isto se aplica para todos os tipos de fluxos citados.

Resolução do sistema de equações discretizadas

As equações discretizadas formam um sistema de equações não lineares, que é resolvido de forma simultânea, usando o método de Newton, devido ao forte acoplamento entre as variáveis. A matriz Jacobiana é calculada via perturbação em diferenças finitas, e mantida constante durante um certo número de iterações. O sistema linear que surge a cada iteração é resolvido por um método iterativo, o GMRES^[16], com pré-condicionamento ILU (fatorização LU incompleta). Durante a montagem da matriz Jacobiana, a esparsidade é levada em conta, de modo a minimizar o custo computacional envolvido na construção e resolução dos sistemas lineares, diminuindo o tempo de resolução.

Deconvolução dos valores médios

Após a resolução das equações discretizadas, chega-se aos valores médios das variáveis nos volumes. Os valores das variáveis em pontos de interesse (centro dos volumes, vértices, etc.) podem ser obtidos via deconvolução dos valores médios. Isto pode ser feitos de diferentes maneiras. Neste trabalho utiliza-se, para este fim, um esquema de interpolação de alta ordem, para obter valores nos vértices dos volumes, a partir dos valores médios nas interfaces do volume^[10]. Deve-se ressaltar que a aproximação usada deve manter a ordem de aproximação usada no procedimento de resolução.

Resultados e Discussão

No presente trabalho são mostrados, para a exemplificação do uso da metodologia, alguns resultados para o escoamento na geometria "stick-slip", um problema muito usado para testes de métodos numéricos e esquemas de interpolação devido à dificuldade de resolução, sendo este o motivo da escolha deste exemplo. Este problema consiste no escoamento de um fluido entre placas paralelas, seguido de um escoamento no qual a superfície do fluido encontra-se livre de cisalhamento (Figura 2). Nestas condições, existe uma singularidade dada pela mudança súbita na condição de contorno (stick: não-deslizamento na parede ® slip: deslizamento), que traz dificuldades na resolução numérica do problema. Este é uma primeira aproximação para simulação do escoamento na saída de uma matriz de extrusão, negligenciando o efeito de inchamento do extrudado.

Por questões de simetria, será usada somente a metade superior da geometria. Foram resolvidos escoamentos de um fluido newtoniano e do fluido de Oldroyd-B. No exemplo considerado, utilizou-se H = 1, L₁ = 3, L₂ = 7. Os valores dos parâmetros adimensionais foram: Re = 1, We = 0,2 e h_E = 0,5.

Para demonstrar o uso e as vantagens da metodologia apresentada, foram implementados três esquemas distintos: QUICK, LAG34 e WENO.

• QUICK: uso do esquema QUICK^[9] para os termos advectivos e diferenças centrais de 2ª ordem para os termos difusivos^[9], um conjunto muito usado na literatura, baseado em valores pontuais (aproximações de 2ª ordem).

• LAG34: uso do esquema de Lagrange de 3ª ordem (equação 8) para os termos advectivos e de 4ª ordem (equação 9) para os termos difusivos.

• WENO: uso do esquema WENO (equações 10 e 11), para os termos advectivos e Lagrange de 4ª ordem (equação 9) para os termos difusivos.

A aproximação para os termos não-lineares é de 4ª ordem (equações 14 e 15), exceto para o QUICK, que é de 2ª ordem (equações 12 e 13).

Comparando-se os esquemas QUICK e LAG34, pode-se observar as vantagens do uso da metodologia baseada em valores pontuais e a baseada em valores médios, respectivamente, no que diz respeito à solução obtida em diferentes tamanhos de malha. O comportamento da solução obtida próximo a singularidades pode ser comparado pelos esquemas LAG34 e WENO.

Na Figura 3 podem ser vistas as linhas de corrente para o escoamento do fluido de Oldroyd-B, obtidas em uma malha 80'20, uniformemente espaçada.

Na Figura 4, faz-se uma comparação entre o perfil de velocidade n_y na direção vertical, para escoamento newtoniano, em c = 0,6667, obtidos pelos conjuntos de aproximações LAG34 e QUICK, em diferentes tamanhos de malhas, variando o número de volumes na direção vertical.

Observa-se que as soluções para o esquema LAG34 são praticamente independentes do tamanho da malha, e que as soluções trazidas pelo QUICK, tendem a solução do LAG34 com o refinamento da malha. Porém, mesmo com 40 volumes na direção y, não se consegue uma solução com mesma qualidade obtida pelo LAG34. Além disto, a metodologia proposta tem um custo computacional da mesma ordem de grandeza que a metodologia convencional, para um mesmo tamanho de malha. Portanto, pode-se obter soluções mais precisas com menor tempo computacional usando o esquema LAG34.

Como foi dito, o uso de esquemas de alta ordem em geral tem o inconveniente de levar a solução oscilatórias em regiões de grande variação das variáveis, tais como ocorre junto à singularidade, neste problema. Isto pode ser visualizado pelo perfil de tensão normal t_xx em uma linha horizontal próxima ao contorno "stick-slip". Com o uso do esquema LAG34 este problema ocorre, como pode ser visto na figura 5. Para escoamentos newtonianos, as oscilações não chegam a prejudicar a convergência, estando presentes somente junto à singularidade. Porém, em maiores números de Weissenberg, estas oscilações tendem a se espalhar, desestabilizando o procedimento numérico, levando à problemas de convergência. Com o uso do esquema WENO, foi possível resolver este problema, como também pode ser visto na Figura 5, chegando-se a uma solução totalmente livre de oscilações.

Conclusões

Foi apresentada neste trabalho uma nova metodologia para a resolução das equações governantes de escoamentos de fluidos viscoelásticos, baseada no método dos volumes finitos, na utilização de esquemas de interpolação de alta ordem, e no uso dos valores médios das variáveis durante o procedimento de resolução. Esta metodologia, comparada com a convencional, baseada em valores pontuais, se mostra mais vantajosa em relação ao custo computacional necessário para a obtenção de soluções com boa precisão. A ocorrência de oscilações e instabilidade numérica em escoamentos com gradientes severos é evitada usando um esquema de interpolação WENO. Foram mostrados alguns resultados que comprovam estas características do método, sendo que resultados mais detalhados, para diferentes tipos de problemas, serão apresentados em trabalhos futuros. A metodologia desenvolvida também será adaptada para coordenadas generalizadas, para resolver problemas em geometrias mais complexas.

Enviado: 16/01/04

Reenviado: 17/12/04

Aprovado: 21/12/04

• 1. Bird, R. B.; Armstrong R. C. & Hassager. 0. - "Dynamics of Polymeric Liquids Vol. 1, Fluid Mechanics", John Wiley (1987).
• 2. Macosko, C. - "Rheology: Principles, Measurements and Applications", VCH Publishers (1994).
• 3. Marchal, J.M. & Crochet, M.J.- J. Non-Newt. Fluid Mech., 26, p.77 (1987).
• 4. Huang, X.; Phan-Thien, N. & Tanner, R.I. - J. Non-Newt Fluid Mech, 64, p.71 (1996)
• 5. Oliveira, P.J.; Pinho, F.T. & Pinto, G.A. - J. Non-Newt. Fluid Mech., 79, p.1 (1998).
• 6. Baaijens, F. P. T. - J. Non-Newt. Fluid Mech., 79, p.361 (1998).
• 7. Aboubacar, M. & Webster, M.F. - J. Non-Newt. Fluid Mech., 98, p.83 (2001).
• 8. Patankar, S. V. - "Numerical Heat Transfer and Fluid Flow", McGraw-Hill, New York (1980).
• 9. Ferziger, J. H. & Peric M. - "Computational Methods for Fluid Dynamics", Springer, Berlin (1999).
• 10. Pereira, J. M. C.; Kobayashi, M. H. & Pereira, J. C. F. - J. Comp. Phys., 167, p.217 (2001).
• 11. Kobayashi, M. H. - J. Comp. Phys., 156, p.137 (1999).
• 12. Muniz, A. R. - "Desenvolvimento de um Método de Volumes Finitos de Alta Ordem para a Simulação de Escoamentos de Fluidos Viscoelásticos", Dissertação de Mestrado, Universidade Federal do Rio Grande do Sul, Brasil (2003).
• 13. Gaskell, P. H. & Lau, A. K. C. - Int. J. Num. Methods Fluids, 8, p.617 (1988).
• 14. Liu, X.; Osher, S. & Chan, T. - J. Comp. Phys., 115, p.200 (1994).
• 15. Jiang, G. S. & Shu, C. W. - J. Comp. Phys., 126, p. 202 (1996).
• 16. Saad, Y. & Schultz, M. H. - SIAM J. Sci. Stat. Comp., 7, p. 856 (1986).

Datas de Publicação

• Publicação nesta coleção
  14 Jul 2005
• Data do Fascículo
  Mar 2005

Histórico

• Aceito
  21 Dez 2004
• Revisado
  17 Dez 2004
• Recebido
  16 Jan 2004