Estudos numéricos da formulação natural do tensor: aplicação em escoamentos viscoelásticos com escorregamento Fabiano Ruano Neto Orientador: Prof. Dr. Cassio Machiaveli Oishi Programa: Matemática Aplicada e Computacional Presidente Prudente, Março de 2021 UNIVERSIDADE ESTADUAL PAULISTA Faculdade de Ciências e Tecnologia de Presidente Prudente Programa de Pós-Graduação em Matemática Aplicada e Computacional Estudos numéricos da formulação natural do tensor: aplicação em escoamentos viscoelásticos com escorregamento Fabiano Ruano Neto Orientador: Prof. Dr. Cassio Machiaveli Oishi Dissertação apresentada ao Programa de Pós-Graduação em Matemática Aplicada e Computacional da Faculdade de Ciências e Tecnologia da UNESP para obtenção do tí- tulo de Mestre em Matemática Aplicada e Computacional. Presidente Prudente, Março de 2021 R894e Ruano Neto, Fabiano Estudos numéricos da formulação natural do tensor : aplicação em escoamentos viscoelásticos com escorregamento / Fabiano Ruano Neto. -- Presidente Prudente, 2021 114 f. : il., tabs. Dissertação (mestrado) - Universidade Estadual Paulista (Unesp), Faculdade de Ciências e Tecnologia, Presidente Prudente Orientador: Cassio Machiaveli Oishi 1. Formulação natural do tensor. 2. Condição de contorno com escorregamento. 3. Solução numérica. 4. Diferenças finitas. I. Título. Sistema de geração automática de fichas catalográficas da Unesp. Biblioteca da Faculdade de Ciências e Tecnologia, Presidente Prudente. Dados fornecidos pelo autor(a). Essa ficha não pode ser modificada. CERTIFICADO DE APROVAÇÃO UNIVERSIDADE ESTADUAL PAULISTA Estudos numéricos da formulação natural do tensor: aplicação em escoamentos viscoelásticos com escorregamento TÍTULO DA DISSERTAÇÃO: AUTOR: FABIANO RUANO NETO ORIENTADOR: CÁSSIO MACHIAVELI OISHI Aprovado como parte das exigências para obtenção do Título de Mestre em MATEMÁTICA APLICADA E COMPUTACIONAL, pela Comissão Examinadora: Prof. Dr. CÁSSIO MACHIAVELI OISHI (Participaçao Virtual) Departamento de Matemática e Computação / Faculdade de Ciencias e Tecnologia de Presidente Prudente Prof. Dr. ANTONIO CASTELO FILHO (Participaçao Virtual) Departamento de Ciências de Computação e Estatística / Universidade de São Paulo Prof. Dr. LUIS JORGE LIMA FERRÁS (Participaçao Virtual) Universidade do Minho Presidente Prudente, 15 de março de 2021 Faculdade de Ciências e Tecnologia - Câmpus de Presidente Prudente - Rua Roberto Simonsen, 305, 19060900, Presidente Prudente - São Paulo http://www.fct.unesp.br/pos-graduacao/--matematica-aplicada-e-computacional/CNPJ: 48.031.918/0009-81. Dedico este trabalho à minha família. Agradecimentos A Deus pela vida e por todas as oportunidades e experiências até aqui. À minha mãe e ao meu pai por todo suporte e incentivo desde sempre. Aos professores da graduação e pós-graduação pela contribuição em minha formação acadêmica. Em especial ao professor José Roberto, meu primeiro orientador, e à professora Gilcilene por me inserir na área de matemática computacional e por toda motivação inspirada. Ao professor e orientador Cassio por toda dedicação, paciência, disponibilidade e aten- ção no desenvolvimento, detalhes e discussões do trabalho. Ao Hugo por toda ajuda e disponibilidade com as implementações no código. Às amigas Karina e Marcela pela amizade desde a graduação, por todas as risadas e descontrações. À Cinthia por toda ajuda e atenção com a parte burocrática do programa. A todos que estiveram comigo durante o curso e contribuíram de forma direta ou indireta. O presente trabalho foi realizado com apoio da Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Código de Financiamento 001. Resumo Neste projeto investiga-se numericamente escoamentos com escorregamento utilizando uma formulação alternativa para descrever o tensor que representa os efeitos viscoelásticos do fluido. São analisados o caso Newtoniano e, no caso viscoelástico, os modelos cons- titutivos Oldroyd-B e PTT (Phan-Thien–Tanner) linear, os últimos tanto na formulação cartesiana quanto na alternativa. Em particular, nesta formulação, conhecida como for- mulação natural do tensor, as equações constitutivas que definem o modelo matemático do fluido viscoelástico são construídas utilizando uma decomposição diádica do tensor polimérico. Estudos recentes indicaram que a aplicação dessa formulação pode melhorar a precisão numérica da solução em escoamentos com singularidades de tensão, como é o caso da contração abrupta e do stick-slip. Entretanto, até o presente momento, nenhum estudo foi realizado com a formulação natural do tensor aplicada na solução de escoa- mentos viscoelásticos com escorregamento. Dessa forma, neste projeto são analisados os efeitos da condição de contorno Navier linear combinada à formulação natural do tensor na solução numérica de escoamentos com singularidades de tensão, como a contração 4:1 e a expansão 1:4. Para este estudo é utilizado um método de projeção e as equações são discretizadas através de fórmulas de diferenças finitas para malhas não-uniformes. Os re- sultados obtidos foram positivos no contexto do estudo numérico fornecendo dados ainda não presentes na literatura. Verificou-se que a variação no coeficiente de escorregamento, com ambas as formulações, provoca alterações nas propriedades do escoamento bem como, em alguns casos, impede a convergência. Em contrapartida, é necessário um referencial teórico para uma conclusão mais concreta dos resultados como, por exemplo, considerar a condição slip no estudo assintótico, o que ainda não foi publicado na literatura. Palavras-Chave: Formulação natural do tensor, Condição de contorno com escorrega- mento, Solução numérica, Diferenças finitas. Abstract In this project we numerically investigate slip flows using an alternative formulation to describe the tensor which represents the fluid viscoelastic effects. We analyze both New- tonian and non-Newtonian cases, considering Oldroyd-B and linear PTT (Phan-Thien– Taner) models with cartesian and alternative formulations. In particular, in this formu- lation which is known as Natural Stress Formulation (NSF), the constitutive equations that define the mathematical modelling of viscoelastic materials are constructed by using a dyadic decomposition of the polymeric tensor. According to recent studies, this formu- lation can improve the accuracy of the numerical solution for solving stress singularity flows, as for instance, the contraction flow and the stick-slip problem. However, until now, there is no study using the NSF for solving viscoelastic fluid slip flows. Therefore, in this project, we analyze the effects of linear Navier slip law combined with NSF in the nume- rical solution of viscoelastic flows with stress singularities, as in the 4:1 contraction and 1:4 expansion. This study is conducted using a projection method and the equations are discretized with a finite difference scheme for non-uniform meshes. The obtained results are good in the numerical study context as they provide new data to the literature. It has been verified that varying the slip coefficient in both formulations causes changes in the flow properties and, in some cases, it is not possible to achieve convergence. In contrast, a theoretical reference is still needed in order to confirm and concretize the results as, for instance, considering the slip boundary condition in asymptotic analysis, not yet avaiable in the literature. Keywords: Natural stress formulation, slip boundary condition, numerical solution, finite difference. Lista de Figuras 2.1 Tipos de parede em que a condição slip é empregada. . . . . . . . . . . . . 30 3.1 Pontos de interesse de uma célula. . . . . . . . . . . . . . . . . . . . . . . . 35 3.2 Células do canal de placas paralelas (a) e da contração (b). . . . . . . . . . 35 3.3 Discretização da condição slip em uma parede horizontal. . . . . . . . . . . 39 3.4 Discretização da condição slip em uma parede vertical. . . . . . . . . . . . 41 3.5 Algumas regiões de interesse da contração. . . . . . . . . . . . . . . . . . . 43 3.6 Pontos da discretização da equação em u na Região 1. . . . . . . . . . . . . 43 3.7 Pontos da discretização da equação em v na Região 1. . . . . . . . . . . . . 44 3.8 Pontos da discretização da equação em u na Região 2. . . . . . . . . . . . . 44 3.9 Pontos da discretização da equação em v na Região 2. . . . . . . . . . . . . 45 3.10 Pontos da discretização da equação em u na Região 3. . . . . . . . . . . . . 46 3.11 Pontos da discretização da equação em v na Região 3. . . . . . . . . . . . . 46 4.1 Geometria da expansão 1:4 de comprimento 40L. . . . . . . . . . . . . . . 49 4.2 Regiões em que a malha não-uniforme é mais refinada. . . . . . . . . . . . 50 4.3 Comprimentos de vórtices obtidos por Ferrás [15] comparados aos deste trabalho. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.4 Ângulo θ e distância r para o estudo assintótico. . . . . . . . . . . . . . . . 51 4.5 Comportamento assintótico da velocidade e das componentes do tensor Ts próximo à quina ao longo de θ = π 2 - Re = 0.01. . . . . . . . . . . . . . . . 52 4.6 Comportamento assintótico das propriedades do escoamento para diferen- tes valores de k̄ ao longo de θ = π 2 - Re = 0.01. . . . . . . . . . . . . . . . . 53 4.7 Variação no comprimento dos vórtices com o aumento de k̄ - NSF, M2. . . 55 4.8 Dimensão dos vórtices obtidos por [15] comparados aos deste trabalho - CSF e NSF. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.9 Ângulos tomados para estudo do comportamento assintótico . . . . . . . . 56 4.10 Cortes verticais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.11 Comportamento assintótico das soluções numéricas próximo à quina ao longo de θ = π 2 - Re = 0.01, ε = 0.25, Wi = 1, β = 1 9 . . . . . . . . . . . . . 58 4.12 Comportamento assintótico das soluções numéricas próximo à quina ao longo de θ = 3 4 π - Re = 0.01, ε = 0.25, Wi = 1, β = 1 9 . . . . . . . . . . . . 59 4.13 Comportamento assintótico das soluções numéricas próximo à quina ao longo de θ = π - Re = 0.01, ε = 0.25, Wi = 1, β = 1 9 . . . . . . . . . . . . . 60 4.14 Variação assintótica das variáveis naturais ao longo de diferentes ângulos θ - Re = 0.01, ε = 0.25, Wi = 1, β = 1 9 . . . . . . . . . . . . . . . . . . . . . . 61 4.15 Comportamento assintótico próximo à quina ao longo de θ = π 2 - CSF, Re = 0.01, ε = 0.25, Wi = 1, β = 1 9 . . . . . . . . . . . . . . . . . . . . . . . 62 4.16 Comportamento assintótico próximo à quina ao longo de θ = π 2 - NSF, Re = 0.01, ε = 0.25, Wi = 1, β = 1 9 . . . . . . . . . . . . . . . . . . . . . . . 63 11 LISTA DE FIGURAS 12 4.17 Comportamento assintótico da pressão ao longo de θ = π 2 , 3π 4 e π - Re = 0.01, ε = 0.25, Wi = 1, β = 1 9 . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.18 Perfis das propriedades do escoamento com a variação de k̄ em x ≈ 10 - NSF, Re = 0.01, ε = 0.25, Wi = 1, β = 1 9 . . . . . . . . . . . . . . . . . . . 65 4.19 Perfis das propriedades do escoamento com a variação de k̄ em x = 20 - NSF, Re = 0.01, ε = 0.25, Wi = 1, β = 1 9 . . . . . . . . . . . . . . . . . . . 66 4.20 Perfis das propriedades do escoamento com a variação de k̄ em x ≈ 30 - NSF, Re = 0.01, ε = 0.25, Wi = 1, β = 1 9 . . . . . . . . . . . . . . . . . . . 67 4.21 Representação da contração 4:1, regiões de formação de vórtices e singula- ridade onde é feita a análise assintótica. . . . . . . . . . . . . . . . . . . . . 68 4.22 Regiões mais refinadas da malha não-uniforme. . . . . . . . . . . . . . . . . 68 4.23 Resultados do caso Newtoniano obtidos por [16] comparados aos deste tra- balho. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.24 Cortes horizontal e vertical na geometria (L = 1). . . . . . . . . . . . . . . 70 4.25 Comportamento assintótico da velocidade, das componentes do tensor Ts e da pressão próximo à quina - Re = 0.01. . . . . . . . . . . . . . . . . . . 71 4.26 Comportamento assintótico da velocidade, das componentes do tensor Ts e da pressão próximo à quina - Re = 0.01. . . . . . . . . . . . . . . . . . . 72 4.27 Perfis das propriedades do escoamento com a variação de k̄ em x ≈ 10 - Re = 0.01. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.28 Perfis das propriedades do escoamento com a variação de k̄ em x = 20 - Re = 0.01. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.29 Perfis das propriedades do escoamento com a variação de k̄ em x ≈ 30 - Re = 0.01. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.30 Perfil horizontal de (a) u e (b) p em y = 4 - Re = 0.01. . . . . . . . . . . . 76 4.31 Comparação do comportamento assintótico das propriedades do escoa- mento utilizando as formulações CSF e NSF na malha M1 - Re = 0.01, Wi = 1 e β = 1 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 4.32 Propriedades do escoamento ao se aproximarem da quina da contração - CSF, Re = 0.01, Wi = 1 e β = 1 2 . . . . . . . . . . . . . . . . . . . . . . . . 78 4.33 Propriedades do escoamento ao se aproximarem da quina da contração - NSF, Re = 0.01, Wi = 1 e β = 1 2 . . . . . . . . . . . . . . . . . . . . . . . . 79 4.34 Redução da dimensão dos vórtices com o aumento de k̄ - Oldroyd-B, NSF, M1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.35 Comportamento assintótico das variáveis com β = 1 9 - CSF, M1, Re = 0.01, Wi = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 4.36 Comportamento assintótico das variáveis com β = 1 9 - CSF, M2, Re = 0.01, Wi = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.37 Redução da dimensão dos vórtices com o aumento de k̄ - PTT, NSF, M1. . 84 4.38 Comportamento assintótico das variáveis com β = 1 2 - CSF, M2, Re = 0.01, Wi = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.39 Comportamento assintótico das variáveis com β = 1 2 - NSF, M2, Re = 0.01, Wi = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 Lista de Tabelas 4.1 Detalhes das malhas utilizadas nas simulações. . . . . . . . . . . . . . . . . 50 4.2 Tamanhos de vórtice, XR - Newtoniano, expansão. . . . . . . . . . . . . . . 51 4.3 Tamanhos de vórtice, XR - CSF e NSF - PTT, expansão. . . . . . . . . . . 54 4.4 Tamanhos de vórtice, XR - Newtoniano, contração. . . . . . . . . . . . . . 69 4.5 Tamanhos de vórtice, XR - PTT, contração. . . . . . . . . . . . . . . . . . 81 4.6 Tamanhos de vórtice, XR - CSF e NSF - PTT, contração. . . . . . . . . . . 84 13 Lista de Siglas CSF: Cartesian Stress Formulation NSF: Natural Stress Formulation PTT: Phan-Thien–Tanner 15 Sumário Resumo 7 Abstract 9 Lista de Figuras 10 Lista de Tabelas 12 Lista de Siglas 15 Capítulos 1 Introdução 19 2 Formulação Matemática 23 2.1 Equações governantes para fluidos Newtonianos e viscoelásticos . . . . . . 23 2.2 Adimensionalização . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.3 Natural Stress Formulation (Formulação natural do tensor) . . . . . . . . . 26 2.4 Condições auxiliares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.4.1 Condições iniciais . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.4.2 Paredes rígidas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.4.3 Entrada de fluido (inflow) . . . . . . . . . . . . . . . . . . . . . . . 28 2.4.4 Saída de fluido (outflow) . . . . . . . . . . . . . . . . . . . . . . . . 29 2.5 Leis de escorregamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.5.1 Adimensionalização das leis de escorregamento . . . . . . . . . . . . 32 3 Método numérico 33 3.1 Método de projeção . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.2 Malha computacional e discretização do domínio . . . . . . . . . . . . . . . 34 3.3 Classificação das células . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.4 Discretização temporal da equação de momento . . . . . . . . . . . . . . . 35 3.5 Atualização do tensor tensão - CSF . . . . . . . . . . . . . . . . . . . . . . 36 3.6 Atualização do tensor tensão - NSF . . . . . . . . . . . . . . . . . . . . . . 37 3.7 Discretização da lei Navier linear . . . . . . . . . . . . . . . . . . . . . . . 39 3.8 Equações de movimento com influência slip . . . . . . . . . . . . . . . . . . 42 3.9 Equações constitutivas - Oldroyd-B e PTT . . . . . . . . . . . . . . . . . . 46 3.10 Etapas do método numérico . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4 Resultados Numéricos 49 4.1 Expansão 1:4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.1.1 Fluido Newtoniano . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.1.2 PTT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 4.2 Contração 4:1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.2.1 Fluido Newtoniano . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.2.2 Oldroyd-B . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.2.3 PTT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 5 Conclusão e trabalhos futuros 87 Referências 88 Apêndices A Diferenças finitas para malhas não-uniformes 93 B Discretização espacial das equações de quantidade de movimento e cons- titutiva 95 B.1 Equação de quantidade de movimento - direção x . . . . . . . . . . . . . . 95 B.2 Equação de quantidade de movimento - direção y . . . . . . . . . . . . . . 96 B.3 Equação de Poisson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 B.4 Equações de atualização da velocidade . . . . . . . . . . . . . . . . . . . . 98 B.5 Equações constitutivas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 B.5.1 Formulação Cartesiana - CSF . . . . . . . . . . . . . . . . . . . . . 99 B.5.2 Formulação tensão natural - NSF . . . . . . . . . . . . . . . . . . . 99 C Equações de movimento em termos das variáveis naturais - detalhes 101 D Leis de escorregamento em termos das variáveis naturais 103 E Análise assintótica da contração: solução Newtoniana 105 F Equações em coordenadas polares 111 Capítulo 1 Introdução A análise matemática sobre a aplicação de condições de contorno em superfícies sólidas nos escoamentos de fluidos é extremamente complexa, combinando as dificuldades da teoria clássica da mecânica do contínuo com efeitos moleculares que surgem na escala microscópica. Dessa forma, um passo inicial nos estudos das equações de Navier-Stokes é assumir a condição de contorno sem escorregamento ou sem deslizamento (no-slip), ou seja, assumir que a velocidade relativa entre a superfície rígida e o fluido é zero na direção tangencial do escoamento. Essa imposição já foi causa de debates tanto do ponto de vista teórico como prático, ainda no contexto de fluidos Newtonianos [31], surgindo a necessidade de um olhar mais cuidadoso para a chamada “condição de deslizamento na parede” (wall slip). O artigo de revisão como em [25] abordou questões importantes sobre esse assunto para fluidos Newtonianos. Os estudos sobre escoamentos com escorregamento ganharam ainda mais destaque no contexto de fluidos não-Newtonianos, já que na falta de uma imposição correta da condição de contorno, pode ocorrer a formação de uma camada limite próxima a parede, resultando na amplificação de efeitos de tensão cisalhante. Os artigos de revisão de [4] e [32] detalham os mecanimos por trás da imposição de condições de parede com efeitos de deslizamento/escorregamento para materiais não-Newtonianos. Os dois parágrafos iniciais dessa Introdução ilustram de forma superficial a complexi- dade do assunto, tanto para fluidos Newtonianos, como para materiais não-Newtonianos. Portanto, a fim de simplificar os estudos sobre escoamentos viscoelásticos com escorre- gamento, neste trabalho foca-se apenas na análise numérica da imposição da condição de deslizamento na parede, analisando os efeitos da aplicação de diferentes tipos de con- dições de escorregamento já consolidadas na literatura, como as leis de Navier-linear e Navier-não-linear. Alguns estudos numéricos sobre a combinação de equações constitutivas de modelos viscoelásticos usando condições de contorno de escorregamento na parede foram apresen- tados na literatura. Por exemplo, no trabalho de [33], efeitos da condição de deslizamento foram analisados para o modelo de Carreau considerando um escoamento compressível e isotérmico enquanto que nos trabalhos de Karapetsas e Mitsoulis [20] e Ferrás et al. [16], os autores investigaram, respectivamente, a relação da lei de Navier-linear em escoamen- tos viscoelásticos modelados pelas equações constitutivas do fluido K-BKZ e sPTT. Outro trabalho interessante foi apresentado por Ferrás et al. [17], onde modelos não-lineares para imposição da condição de contorno de escorregamento na parede foram numericamente analisados no contexto do método de volumes finitos. Já em [19], o modelo de Navier- não-linear foi imposto na solução do escoamento viscoelástico do modelo Oldroyd-B em um tubo. 19 1. Introdução 20 Mais recentemente, pesquisadores têm estendido os estudos de escoamentos com es- corregamentos para solução de problemas de interesse prático. Em [21], os autores utili- zaram um método espectral para estudar escoamentos viscoelásticos com escorregamento em meios porosos, enquanto que Drapaca [5] aplicou o conceito de derivadas fracioná- rias de Caputo combinado com condições de contorno de escorregamento para modelar numericamente microaneurismas cerebrais. Com respeito ao problema de extrusão, am- plamente investigado no contexto industrial, Kwon [22] investigou a relação entre modelos de deslizamento na parede com o surgimento de instabilidade de fraturas, como é o caso do efeito sharkskin. Já no contexto de processamento de polímeros, Fernandes et al. [13] implementaram no pacote computacional OpenFOAM as condições de contorno de deslizamento parcial acopladas as equações constitutivas do modelo PTT. Os trabalhos numéricos citados acima envolvendo escoamentos viscoelásticos com con- dições de contorno de escorregamento na parede utilizaram, para definir as equações constitutivas, a clássica formulação cartesiana do tensor polimérico. Nesta formulação, o tensor que representa a contribuição não-Newtoniana é decomposto em componentes combinadas com os vetores da base canônica do sistema de coordenadas cartesiano. Em particular, essa abordagem no tratamento do tensor polimérico é amplamente adotada em Reologia Computacional [26]. Entretanto, uma dificuldade reconhecida na literatura ao utilizar a formulação cartesiana é o tratamento numérico da singularidade do tensor. Estas singularidades surgem como consequência de uma mudança abrupta nas condições de contorno, como no caso do stick-slip, ou devido a presença de quinas na geometria do problema, como no escoamento da contração. Uma forma alternativa de representar o ten- sor polimérico foi apresentada por Renardy [29], com o objetivo de representar melhor o comportamento da singularidade da tensão. Essa formulação recebeu o nome de formula- ção natural do tensor (Natural Stress Formulation). Porém, essa formulação foi proposta originalmente no contexto de escoamentos estacionários, e como é de conhecimento, os problemas mais complexos são intrinsecamente transientes. A primeira versão das equações constitutivas e transientes da formulação natural do tensor foi apresentada pelo orientador desse projeto e o colaborador Jonathan D. Evans em [8]. Logo após essa publicação, a formulação foi analisada em mais detalhes nos recentes trabalhos [10, 12], que evidenciaram o potencial dessa metodologia na obtenção de soluções mais acuradas em escoamentos com singularidades. Vale ressaltar que o pesquisador J.D. Evans já havia investigado a formulação natural do tensor do ponto de vista da análise assintótica nos trabalhos [7, 9], porém nenhum resultado ainda havia sido publicado no contexto de solução das equações transientes. Desta forma, o orientador, em colaboração com J.D. Evans e outros pesquisadores, vem trabalhando no desenvolvimento e aperfeiçoamento numéricos dessa formulação alternativa. Nos trabalhos citados acima do orientador desse projeto e colaboradores, a formulação natural do tensor foi aplicada apenas em escoamentos viscoelásticos sem escorregamen- tos, ou seja, apenas adotando a condição de contorno do tipo sem deslizamento (no-slip). Neste trabalho de mestrado é analisada e implementada a lei de Navier linear para a mo- delagem de escoamentos com escorregamentos utilizando a formulação natural do tensor para a definição do sistema de equações constitutivas. A lei Navier não-linear é estu- dada brevemente porém não implementada. Vale ressaltar que, em grande maioria dos trabalhos que discutiram efeitos das condições de contorno de escorregamento na parede, utilizaram escoamentos que apresentam singularidade da tensão, como no caso da con- tração 4:1 [16] ou no problema de extrusão (stick-slip e slip-stick) [17, 22]. Portanto, a combinação da formulação natural do tensor com modelos de leis de escorregamento na parede auxiliará na obtenção de uma solução numérica mais precisa nestes escoamentos complexos. 1. Introdução 21 Alguns resultados alcançados no trabalho são: • Implementação da lei de escorregamento no código do grupo de pesquisa; • Testes de validação na contração e expansão com as referências [15, 16] considerando os casos Newtoniano e viscoelástico utilizando condições de contorno com escorrega- mento. Para o caso viscoelástico utilizou-se a formulação tradicional para o tensor polimérico; • Simulações de fluidos viscoelásticos nas geometrias com contração e expansão com- binando a formulação natural do tensor com a condição de escorregamento. O trabalho está organizado da seguinte maneira: • Capítulo 2: são apresentadas as equações governantes, o processo de adimensiona- lização, as equações constitutivas escritas com a formulação natural do tensor, as condições de contorno consideradas, as leis de escorregamento na forma dimensional e adimensional; • Capítulo 3: é apresentado o método numérico empregado, a malha computacional utilizada e a discretização do domínio bem como a discretização das equações (equa- ção de momento, atualização dos tensores em ambas as formulações, lei de Navier linear e consequentes mudanças nas equações); • Capítulo 4: resultados numéricos que tratam de escoamentos na contração e expan- são para fluidos Newtonianos e viscoelásticos; • Capítulo 5: considerações finais sobre o trabalho e etapas futuras. Além disso, encontram-se nos apêndices alguns complementos dos estudos desenvolvi- dos como, por exemplo, as equações de quantidade de movimento escritas em termos das variáveis da formulação NSF, as leis de escorregamento Navier linear e Navier não-linear também escritas com as variáveis naturais e estudos iniciais sobre análise assintótica, nos quais considerou-se o caso Newtoniano. Capítulo 2 Formulação Matemática 2.1 Equações governantes para fluidos Newtonianos e viscoelásticos O escoamento de fluidos incompressíveis e compressíveis, laminares e turbulentos é modelado pelas equações de Navier-Stokes. Essas equações representam, matematica- mente, os conceitos físicos de conservação de massa, balanço do momento e conservação de energia. Neste trabalho são considerados escoamentos incompressíveis e isotérmicos e, desta forma, as equações de quantidade de movimento e continuidade se resumem, respectivamente, a ρ ( ∂u ∂t +∇ · (uuᵀ) ) = −∇p+∇ · τ , (2.1) ∇ · u = 0, (2.2) onde • x = (x, y)ᵀ é o vetor posição de um sistema de coordenadas cartesianas bidimensi- onais; • u = (u, v)ᵀ é o vetor velocidade com componentes u = u(x, y, t) e v = v(x, y, t); • t ≥ 0 é o tempo; • ρ é a densidade do fluido; • p é a pressão; • τ = [ τxx τxy τxy τ yy ] é o tensor tensão extra, com τ = τ (x, y, t). O tensor τ é decomposto em uma parcela que representa a contribuição do solvente (Newtoniana) τ s e uma polimérica τ p e é dado por τ = τ s + τ p. (2.3) A contribuição newtoniana é escrita como τ s = 2ηsD, (2.4) 23 2. Formulação Matemática 24 onde ηs é a viscosidade do solvente e D = 1 2 (∇u + (∇u)ᵀ) é o tensor taxa de deformação. A parcela polimérica depende do modelo viscoelástico utilizado. Os modelos utili- zados nesse trabalho são o Oldroyd-B e o PTT (Phan-Thien–Tanner), cujas equações constitutivas são dadas, respectivamente, por τ p + λp O τ p = 2ηpD, (2.5) e τ p + λp ( O τ p + ε ηp tr(τ p)τ p ) = 2ηpD. (2.6) As constantes λp e ηp são o tempo de relaxação do polímero e viscosidade polimérica, respectivamente, e o parâmetro ε, para o modelo PTT, regula a influência do termo quadrático em (2.6), isto é, controla o comportamento elongacional do modelo. Ainda, a derivada convectiva contravariante O τ p é expressa por O τ p = ∂τ p ∂t + (u · ∇)τ p − (∇u)τ p − τ p(∇u)ᵀ. (2.7) 2.2 Adimensionalização A adimensionalização das equações é feita utilizando os seguintes termos x = Lx̄, t = L U t̄, u = U ū, p = ρU2p̄, τ s = η0U L T̄s, τ p = η0U L T̄, (2.8) em que L é o comprimento característico, U é a velocidade característica e η0 = ηs + ηp é a viscosidade total a uma taxa nula de cisalhamento. Substituindo as expressões de (2.8) nas equações de Navier-Stokes e, a fim de simpli- ficar a notação, omitindo as barras, obtém-se as equações na forma adimensional ∂u ∂t +∇ · (uuᵀ) = −∇p+ β Re ∇2u + 1 Re ∇ ·T, (2.9) ∇ · u = 0. (2.10) Também substituindo (2.8) nas equações dos modelos constitutivos, obtém-se, para o fluido Oldroyd-B, a equação T +Wi ( ∂T ∂t + (u · ∇)T− (∇u)T−T(∇u)ᵀ ) = 2(1− β)D, (2.11) e, para o fluido PTT, a equação na forma adimensional é dada por T +Wi ( ∂T ∂t + (u · ∇)T− (∇u)T−T(∇u)ᵀ + ε 1− β tr(T)T ) = 2(1− β)D. (2.12) Os parâmetros adimensionais presentes nas equações governantes são: • Número de Reynolds (Re): razão entre as forças inerciais e as forças viscosas. É dado por Re = ρUL η0 . (2.13) 2. Formulação Matemática 25 • Número de Weissenberg (Wi): razão entre a escala de tempo característica de um fluido e a escala de tempo característica de um escoamento. Este parâmetro é dado por Wi = λpU L . (2.14) • Razão das viscosidades (β ∈ (0, 1]): essa constante regula a contribuição do solvente Newtoniano. Quando β = 1 e T = 0, o escoamento é Newtoniano. O parâmetro em questão é definido como β = ηs η0 . (2.15) Reescrevendo as equações de Navier-Stokes e as equações constitutivas adimensionais em coordenadas cartesianas, tem-se: • Navier-Stokes ∂u ∂t + ∂(uu) ∂x + ∂(uv) ∂y = −∂p ∂x + β Re ( ∂2u ∂x2 + ∂2u ∂y2 ) + 1 Re ( ∂T xx ∂x + ∂T xy ∂y ) , (2.16) ∂v ∂t + ∂(uv) ∂x + ∂(vv) ∂y = −∂p ∂y + β Re ( ∂2v ∂x2 + ∂2v ∂y2 ) + 1 Re ( ∂T xy ∂x + ∂T yy ∂y ) , (2.17) ∂u ∂x + ∂v ∂y = 0, (2.18) • Oldroyd-B T xx +Wi ( ∂T xx ∂t + ∂(uT xx) ∂x + ∂(vT xx) ∂y − 2 ∂u ∂x T xx − 2 ∂u ∂y T xy ) = 2(1− β) ∂u ∂x , (2.19) T xy +Wi ( ∂T xy ∂t + ∂(uT xy) ∂x + ∂(vT xy) ∂y − ∂v ∂x T xx − ∂u ∂y T yy ) = (1− β) ( ∂u ∂y + ∂v ∂x ) , (2.20) T yy +Wi ( ∂T yy ∂t + ∂(uT yy) ∂x + ∂(vT yy) ∂y − 2 ∂v ∂x T xy − 2 ∂v ∂y T yy ) = 2(1− β) ∂v ∂y . (2.21) • PTT T xx+Wi ( ∂T xx ∂t + ∂(uT xx) ∂x + ∂(vT xx) ∂y − 2 ∂u ∂x T xx − 2 ∂u ∂y T xy + ε 1− β fxx ) = 2(1−β) ∂u ∂x , (2.22) T xy+Wi ( ∂T xy ∂t + ∂(uT xy) ∂x + ∂(vT xy) ∂y − ∂v ∂x T xx − ∂u ∂y T yy + ε 1− β fxy ) = (1−β) ( ∂u ∂y + ∂v ∂x ) , (2.23) 2. Formulação Matemática 26 T yy +Wi ( ∂T yy ∂t + ∂(uT yy) ∂x + ∂(vT yy) ∂y − 2 ∂v ∂x T xy − 2 ∂v ∂y T yy + ε 1− β fyy ) = 2(1− β) ∂v ∂y , (2.24) com fxx = (T xx + T yy)T xx, fxy = (T xx + T yy)T xy, fyy = (T xx + T yy)T yy. (2.25) 2.3 Natural Stress Formulation (Formulação natural do tensor) A Formulação natural do tensor (NSF), proposta por Renardy [29], tem por objetivo evitar instabilidades no cálculo do tensor T próximo às quinas presentes, por exemplo, em escoamentos com contração e no cross-slot. A ideia é representar o tensor numa base natural alinhada às linhas de corrente do escoamento em questão. A seguir são descritos os passos para a obtenção das equações constitutivas em termos das variáveis naturais λ, µ, e ν. O primeiro passo é escrever o tensor polimérico T em função do tensor configuração A, isto é, T = (1− β) Wi (A− I). (2.26) Considerando o modelo Oldroyd-B, substituímos (2.26), com O I = −2D, em (2.11) e, fazendo as simplifações necessárias, obtém-se Wi O A + (A− I) = 0. (2.27) O tensor A é escrito como A = λuuᵀ + µ(uwᵀ + wuᵀ) + νwwᵀ, (2.28) onde w = 1 |u|2 (−v, u)ᵀ um vetor ortogonal a u com |u×w|= |u||w|= |u| |u||u|2 = 1. A próxima etapa consiste em substituir (2.28) em (2.27) e escrever as equações em termos das novas variáveis, cujo resultado é o seguinte Wi [ ∂λ ∂t + 2λ |u| ∂|u| ∂t + 2µ |u|4 ( v ∂u ∂t − u∂v ∂t ) + (u · ∇)λ+ 2µ∇ ·w ] + λ− 1 |u|2 = 0, (2.29) Wi [ ∂µ ∂t + ( λ− ν |u|4 )( u ∂v ∂t − v∂u ∂t ) + (u · ∇)µ+ ν∇ ·w ] + µ = 0, (2.30) Wi [ ∂ν ∂t + 2ν |u| ∂|u| ∂t + 2µ ( u ∂v ∂t − v∂u ∂t ) + (u · ∇)ν ] + ν − |u|2= 0, (2.31) com ∇ ·w = 1 |u|4 [ (v2 − u2) ( ∂u ∂y + ∂v ∂x ) + 4uv ∂u ∂x ] . (2.32) A fim de evitar instabilidades numéricas nos termos divididos por potências de |u|, é feita uma mudança nas variáveis naturais. Para tanto, utiliza-se os termos abaixo 2. Formulação Matemática 27 λ = λ̂ |u|2 , µ = µ̂, ν = ν̂|u|2, (2.33) que, substituídos em (2.29)−(2.31), produzem as equações [ ∂λ̂ ∂t + 2µ̂ |u|2 ( v ∂u ∂t − u∂v ∂t ) + |u|2(u · ∇) ( λ̂ |u|2 ) + 2µ̂|u|2∇ ·w ] + λ̂− 1 Wi = 0, (2.34)[ ∂µ̂ ∂t + ( λ̂− ν̂ |u|2 )( u ∂v ∂t − v∂u ∂t ) + (u · ∇)µ̂+ ν̂|u|2∇ ·w ] + µ̂ Wi = 0, (2.35)[ ∂ν̂ ∂t + 2µ̂ |u|2 ( u ∂v ∂t − v∂u ∂t ) + 1 |u|2 (u · ∇)(ν̂|u|2) ] + ν̂ − 1 Wi = 0, (2.36) com |u|2∇ ·w = 1 |u|2 [ (v2 − u2) ( ∂u ∂y + ∂v ∂x ) + 4uv ∂u ∂x ] . (2.37) No caso do fluido PTT, os passos seguidos são os mesmos realizados para o Oldroyd-B. Substituindo (2.26) com O I = −2D, agora em (2.12), o resultado é Wi O A + (A− I) + ε tr(A− I)(A− I) = 0. (2.38) Escrevendo as equações de (2.38) em termos de λ̂, µ̂ e ν̂, tem-se [ ∂λ̂ ∂t + 2µ̂ |u|2 ( v ∂u ∂t − u∂v ∂t ) + |u|2(u · ∇) ( λ̂ |u|2 ) + 2µ̂|u|2∇ ·w ] + α(λ̂− 1) Wi = 0, (2.39)[ ∂µ̂ ∂t + ( λ̂− ν̂ |u|2 )( u ∂v ∂t − v∂u ∂t ) + (u · ∇)µ̂+ ν̂|u|2∇ ·w ] + αµ̂ Wi = 0, (2.40)[ ∂ν̂ ∂t + 2µ̂ |u|2 ( u ∂v ∂t − v∂u ∂t ) + 1 |u|2 (u · ∇)(ν̂|u|2) ] + α(ν̂ − 1) Wi = 0, (2.41) onde α = 1 + ε(λ̂− 2 + ν̂) e |u|2∇ ·w = 1 |u|2 [ (v2 − u2) ( ∂u ∂y + ∂v ∂x ) + 4uv ∂u ∂x ] . (2.42) Agora, a fim de determinar uma relação entre as variáveis NSF e as CSF, substitui-se (2.28) e (2.33) em (2.26). Desta forma, as variáveis se relacionam pelas expressões abaixo T xx = (1− β) Wi [ 1 |u|2 ( λ̂u2 − 2µ̂uv + ν̂v2 ) − 1 ] , (2.43) T xy = (1− β) Wi|u|2 [ λ̂uv + µ̂(u2 − v2)− ν̂uv ] , (2.44) T yy = (1− β) Wi [ 1 |u|2 ( λ̂v2 + 2µ̂uv + ν̂u2 ) − 1 ] , (2.45) e 2. Formulação Matemática 28 λ̂− 1 = Wi (1− β)|u|2 ( u2T xx + 2uvT xy + v2T yy ) , (2.46) µ̂ = Wi (1− β)|u|2 ( −uvT xx + (u2 − v2)T xy + uvT yy ) , (2.47) ν̂ − 1 = Wi (1− β)|u|2 ( v2T xx + 2uvT xy + u2T yy ) . (2.48) 2.4 Condições auxiliares O uso de condições auxiliares na solução numérica de equações diferenciais parciais é muito importante uma vez que busca-se solução única para um determinado problema. Assim, são necessárias as condições iniciais que são valores das propriedades do escoa- mento conhecidos em todo o domínio no tempo t = 0 e também as condições de contorno que são valores prescritos de propriedades do escoamento em determinadas áreas do do- mínio. A seguir são apresentadas as condições auxiliares utilizadas neste trabalho. 2.4.1 Condições iniciais As condições iniciais adotadas nas simulações para o campo de velocidade u, a pressão p e o tensor polimérico T são nulas. Em decorrência disso, com (2.46)-(2.48), para as variáveis naturais tem-se λ̂ = ν̂ = 1 e µ̂ = 0. 2.4.2 Paredes rígidas Como o fluido não penetra as paredes, a velocidade normal é nula e, representada por un = 0. (2.49) A velocidade tangencial depende do tipo de condição de contorno adotada: sem es- corregamento (no-slip) ou com escorregamento (slip). No caso no-slip o fluido adere à parede e não se movimenta, de modo que a velocidade tangencial é ut = 0. (2.50) Já quando a condição é slip, tem-se ut = uslip. (2.51) Neste trabalho, é adotada a lei de escorregamento Navier linear, cuja expressão é dada por (2.56). Detalhes e a discretização das condições de escorregamento são apresentadas na Seção 2.5. 2.4.3 Entrada de fluido (inflow) Na região por onde o fluido é injetado a velocidade normal é prescrita (uinflow) e consideramos a velocidade tangencial nula,{ un = uinflow ut = 0 . (2.52) 2. Formulação Matemática 29 Para o tensor T é considerada a solução analítica de um escoamento em um canal de placas paralelas totalmente desenvolvido para um fluido Oldroyd-B. As condições são expressas por T xx = 2(1− β)Wi ( ∂u ∂y )2 , T xy = (1− β) ( ∂u ∂y ) , T yy = 0. (2.53) No caso da NSF, as condições no inflow são λ̂ = Wi (1− β) T xx + 1, µ̂ = Wi (1− β) T xy, ν̂ = Wi (1− β) T yy + 1, (2.54) e são obtidas por meio das equações (2.46)−(2.48). 2.4.4 Saída de fluido (outflow) Na região de outflow é aplicada condição de Neumann homogênea tanto para u quanto para T, ou seja,  ∂un ∂n = 0, ∂ut ∂n = 0, ∂T ∂n = 0, (2.55) sendo n o vetor normal ao outflow. 2.5 Leis de escorregamento Por muito tempo a condição de não escorregamento (no-slip) tem sido empregada na modelagem de experimentos macroscópicos e essa condição afirma que, em um es- coamento, o fluido adere à parede e as partículas de fluido próximas da parede não se movimentam [14]. Porém, em algumas situações como o escoamento de gás, de fluidos não-Newtonianos e dinâmica de linha de contato sabe-se que a velocidade do fluido na parede não é nula [23]. Em 1822, Navier [24] apresentou a ideia de que um líquido pode deslizar sobre uma superfície sólida com velocidade de escorregamento uslip oposta a uma força de fricção proporcional ao gradiente de velocidade na direção normal à parede, isto é, a velocidade uslip é proporcional à taxa de cisalhamento na superfície [1]. A Figura 2.1 ilustra os quatro tipos de paredes presentes nas geometrias dos escoamen- tos estudados neste trabalho e as componentes normal e tangencial do vetor velocidade. A lei de escorregamento Navier linear é dada pela expressão uparede = −k{n · τ − [(n · τ ) · nᵀ]n}, (2.56) onde n é o vetor normal ao contorno, k é o coeficiente de fricção que controla a intensidade da velocidade de escorregamento na parede e τ é o tensor tensão extra apresentado em (2.3). A partir dessa expressão, reescrita em coordenadas cartesianas bidimensionais, tem- se as expressões para cada uma das paredes apresentadas na Figura 2.1. Começa-se reescrevendo o tensor tensão extra τ = 2ηsD + τ p, ou seja, 2. Formulação Matemática 30 parede (a) Parede horizontal superior. parede (b) Parede horizontal inferior. pa re de (c) Parede vertical à esquerda. parede (d) Parede vertical à direita. Figura 2.1: Tipos de parede em que a condição slip é empregada. [ τxx τxy τxy τ yy ] = ηs  2 ∂u ∂x ∂v ∂x + ∂u ∂y ∂v ∂x + ∂u ∂y 2 ∂v ∂y + [ τxxp τxyp τxyp τ yyp ] . (2.57) Inicialmente, considera-se uma parede horizontal superior e, neste caso, o vetor normal é dado por n = [0, 1]. Assim, substituindo n e τ em (2.61), a velocidade uparede é expressa por uparede =− k { [0 1] [ τxx τxy τxy τ yy ] − [( [0 1] [ τxx τxy τxy τ yy ])[ 0 1 ]] [0 1] } parede =− k{ [τxy τ yy]− [0 τ yy] }parede =− k [τxy 0]parede =− k [ ηs ( ∂v ∂x + ∂u ∂y ) + τxyp 0 ] parede . (2.58) Como se trata de uma parede paralela ao eixo x, temos que ∂v ∂x = 0. Finalmente, a expressão se reduz a uparede = −k [ ηs ∂u ∂y + τxyp , 0 ] parede . (2.59) Quando a condição de escorregamento é empregada numa parede horizontal inferior o vetor normal é n = [0, −1] e, procedendo de maneira análoga ao que foi feito acima, obtem-se 2. Formulação Matemática 31 uparede = k [ ηs ∂u ∂y + τxyp , 0 ] parede . (2.60) Considerando paredes verticais, a maneira de determinar as expressões para uparede segue o mesmo procedimento utilizado para paredes horizontais, porém, agora a derivada nula é ∂u ∂y . Para uma parede vertical à esquerda tem-se n = [−1, 0] e, então uparede = k [ 0, ηs ∂v ∂x + τxyp ] parede , (2.61) e, para uma parede vertical à direita, como n = [1, 0], uparede = −k [ 0, ηs ∂v ∂x + τxyp ] parede . (2.62) Os mesmos procedimentos são realizados para a lei de Navier não-linear, dada pela expressão ut = −k‖τ t‖m τ t ‖τ t‖ , (2.63) onde τ t é o tensor tensão extra tangente à parede dado por τ t = (I− nᵀ · n) (τ · nᵀ) . (2.64) Desta forma, de acordo com o vetor normal n, ou seja, dependendo da parede que se toma como referência, tem-se as seguintes expressões: • Para n = [0, 1]: uparede = −k [( ηs ∂u ∂y + τxyp )m , 0 ] parede , (2.65) • Para n = [0, −1]: uparede = k [( ηs ∂u ∂y + τxyp )m , 0 ] parede , (2.66) • Para n = [−1, 0]: uparede = k [ 0, ( ηs ∂v ∂x + τxyp )m] parede , (2.67) • Para n = [1, 0]: uparede = −k [ 0, ( ηs ∂v ∂x + τxyp )m] parede . (2.68) 2. Formulação Matemática 32 2.5.1 Adimensionalização das leis de escorregamento Para adimensionalizar as expressões referentes as leis de escorregamento Navier linear (m = 1) e não-linear, são utilizadas as seguintes variáveis adimensionais x = Lx̄, u = U ū, τ p = η0U L T̄, k = U1−m (η0 L )−m k̄. (2.69) Com as substituições e os cálculos necessários são obtidas as seguintes expressões adimensionais para uparede em que as barras são omitidas: • Navier linear (m = 1): uparede =∓ k [ β ∂u ∂y + T xy, 0 ] parede , (2.70) uparede = ∓ k [ 0, β ∂v ∂x + T xy ] parede . (2.71) • Navier não-linear: uparede =∓ k [( β ∂u ∂y + T xy )m , 0 ] parede , (2.72) uparede = ∓ k [ 0, ( β ∂v ∂x + T xy )m] parede . (2.73) Vale lembrar que, para m = 1, tem-se a lei Navier linear e que, no caso de fluidos Newtonianos, basta tomar β = 1 nas expressões acima. Capítulo 3 Método numérico 3.1 Método de projeção Devido à dificuldade de resolver simultaneamente o sistema de equações (2.9)−(2.10), o método de projeção objetiva desacoplar as variáveis u e p desse sistema, resultando, assim, em equações onde u e p são determinadas de maneira isolada. O método de projeção é embasado no Teorema da Decomposição de Helmholtz-Hodge que garante que, para um campo vetorial ũ em uma região Ω com fronteira suave ∂Ω, a decomposição de ũ na forma ũ = u +∇ψ, (3.1) existe e é única, sendo ψ um campo escalar definido em Ω e com ∇ · u = 0. Inicialmente, o método de projeção baseia-se em resolver a equação ∂ũ ∂t +∇ · (ũũᵀ) = −∇p̃+ β Re ∇2ũ + 1 Re ∇ ·T, (3.2) em que tem-se os valores de T e p̃ 6= p, uma pressão tentativa. Como, para calcular u, é necessário obter o valor de ψ, isolamos u em (3.1) e, aplicando o divergente, tem-se ∇ · u = ∇ · ũ−∇2ψ. (3.3) Desconsiderando a equação de incompressibilidade para ũ e aplicando (2.2) em (3.3), chega-se à equação de Poisson ∇2ψ = ∇ · ũ. (3.4) As condições de contorno adotadas para ũ são as mesmas utilizadas para u apresen- tadas na Seção 2.4. Já para ψ as condições de contorno para inflow e fronteiras rígidas, segundo [34], são ∂ψ ∂n = 0, (3.5) e no outflow, ψ = 0. (3.6) 33 3. Método numérico 34 3.2 Malha computacional e discretização do domínio A malha computacional em que as equações são discretizadas é do tipo não-uniforme e as células desta malha são deslocadas. Desta forma, define-se por Nx e Ny a quantidade de intervalos em que o domínio é dividido nas direções x e y, respectivamente. Os pontos da malha distam, na direção x, ∆xi = xi − xi−1, i = 1, 2, · · · , Nx, (3.7) e, na direção y, ∆yj = yj − yj−1, j = 1, 2, · · · , Ny. (3.8) É necessário identificar as células que compõem a malha e os pontos em que as incóg- nitas das equações governantes serão calculadas. Deste modo, as células são identificadas pelo par ordenado (i, j) e seu centro pelo ponto xi,j, matematicamente representado por xi,j = ( xi + 1 2 ∆xi, yj + 1 2 ∆yj ) . (3.9) Nesses pontos são calculadas as componentes do tensor tensão T xx, T xy, T yy, as variá- veis naturais λ, µ e ν, a pressão p e o campo escalar ψ. Agora, é possível classificar pontos que se encontram nas faces das células do domínio. Um ponto localizado na face vertical esquerda de uma célula é indicado por xi− 1 2 ,j e representado por xi− 1 2 ,j = ( xi, yj + 1 2 ∆yj ) , (3.10) e, neste ponto, é calculada a componente u da velocidade. Um ponto que está na face vertical direita de uma célula é da forma xi+ 1 2 ,j = ( xi + ∆xi, yj + 1 2 ∆yj ) . (3.11) Já um ponto posicionado na face horizontal inferior de uma célula é escrito como xi,j− 1 2 e sua representação é xi,j− 1 2 = ( xi + 1 2 ∆xi, yj ) , (3.12) ponto em que a componente v da velocidade é calculada. Por outro lado, um ponto na face horizontal superior é representado por xi,j+ 1 2 = ( xi + 1 2 ∆xi, yj + ∆yj ) . (3.13) A Figura 3.1 ilustra os pontos descritos onde as propriedades u, v, p, ψ, T xx, T xy, T yy, λ, µ e ν são calculadas. 3.3 Classificação das células As células computacionais são classificadas como: • INFLOW (I): células da região por onde o fluido entra; 3. Método numérico 35 x x x x x x x x x Δ Δ Figura 3.1: Pontos de interesse de uma célula. • OUTFLOW (O): células da região por onde o fluido sai; • FULL (F): células da região preenchida por fluido; • BOUNDARY (B): células de regiões de parede rígida (ou contorno). Em escoamentos em um canal de placas paralelas e em uma contração abrupta, por exemplo, as células são classificadas como na Figura 3.2. B B B B B B B B B B B B B B B B B B B B B B B B F F F F F F F F F F F F B B F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F I I I I I O O O O O O O (a) I I I I I I I I I I I I B F B F F F F F F F F F F F B F B F F F F F F F F F F F B F B F F F F F F F F F F F B F B F F F F F F F F F F F B F B F F F F F F F F F F F B F B F F F F F F F F F F F B F B F F F F F F F F F F F B F B F F F F F F F F F F F B F B F F F F F F F F F F F B F B F F F F F F F F F F F B F B F F F F F F F F F F F B F B F F F F F F F F F F F B B B B B B F F F B B F B B B F F F F B B F F F F B B F F F F B B F F F F B B F F F B F B F F F F B B F F F F B B F F F F B B F F F F B O O O O O O (b) Figura 3.2: Células do canal de placas paralelas (a) e da contração (b). 3.4 Discretização temporal da equação de momento Nesta seção é feita a discretização temporal da equação (3.2) de modo que a velocidade u é aproximada de maneira semi-implícita, uma vez que a derivada temporal é aproximada com o método de Euler Implícito e os termos convectivos são aproximados utilizando o método de Euler Explícito. 3. Método numérico 36 Assim tem-se a seguinte equação discretizada ũn+1 − ũn ∆t +∇ · (ũũᵀ)n = −∇p̃n+1 + β Re ∇2ũn+1 + 1 Re ∇ ·Tn. (3.14) Como procura-se determinar ũn+1, é assumido que ũn = un e p̃n+1 = pn, valores conhecidos. Desta forma, a equação anterior torna-se ũn+1 − un ∆t +∇ · (uuᵀ)n = −∇pn + β Re ∇2ũn+1 + 1 Re ∇ ·Tn. (3.15) Também é necessário atualizar a pressão pn+1 e essa atualização é feita inicialmente discretizando a equação (2.9) no tempo un+1 − un ∆t +∇ · (uuᵀ)n = −∇pn+1 + β Re ∇2un+1 + 1 Re ∇ ·Tn, (3.16) e realizando a operação (3.16) − (3.15), que resulta em un+1 − ũn+1 ∆t = −∇pn+1 +∇pn + β Re [ ∇2un+1 −∇2ũn+1 ] . (3.17) Agora, substituindo (3.1) em (3.17) obtem-se un+1 − un+1 −∇ψn+1 ∆t = −∇pn+1 +∇pn + β Re [ ∇2un+1 −∇2 ( un+1 +∇ψn+1 )] . (3.18) Procedendo com as simplificações possíveis e descartando o termo do operador ∇2, tem-se a seguinte equação para atualizar a pressão pn+1 = pn + ψn+1 ∆t . (3.19) 3.5 Atualização do tensor tensão - CSF Agora é necessário calcular as componentes do tensor tensão, isto é, T xx, n+1, T xy, n+1 e T yy, n+1. Isso é feito discretizando, no tempo com o método de Euler Explícito, o conjunto de equações (2.19)-(2.21) se o modelo constitutivo Oldroyd-B for utilizado ou então as equações (2.22)-(2.24) quando o modelo escolhido for o PTT. Para o modelo Oldroyd-B, as equações discretizadas são dadas por T xx, n+1 = T xx, n + ∆t · [ −∂(un+1T xx, n) ∂x − ∂(vn+1T xx, n) ∂y + 2 ∂un+1 ∂x T xx, n +2 ∂un+1 ∂y T xy, n + 2 (1− β) Wi ∂un+1 ∂x − 1 Wi T xx, n ] , (3.20) T xy, n+1 = T xy, n + ∆t · [ −∂(un+1T xy, n) ∂x − ∂(vn+1T xy, n) ∂y + ∂vn+1 ∂x T xx, n + ∂un+1 ∂y T yy, n + (1− β) Wi ( ∂un+1 ∂y + ∂vn+1 ∂x ) − 1 Wi T xy, n ] , (3.21) 3. Método numérico 37 T yy, n+1 = T yy, n + ∆t · [ −∂(un+1T yy, n) ∂x − ∂(vn+1T yy, n) ∂y + 2 ∂vn+1 ∂x T xy, n +2 ∂vn+1 ∂y T yy, n + 2 (1− β) Wi ∂vn+1 ∂y − 1 Wi T yy, n ] . (3.22) Já para o PTT, temos as seguintes equações T xx, n+1 = T xx, n + ∆t · [ −∂(un+1T xx, n) ∂x − ∂(vn+1T xx, n) ∂y + 2 ∂un+1 ∂x T xx, n +2 ∂un+1 ∂y T xy, n − ε 1− β fnxx + 2 (1− β) Wi ∂un+1 ∂x − 1 Wi T xx, n ] , (3.23) onde fnxx = (T xx, n + T yy, n)T xx, n, T xy, n+1 = T xy, n + ∆t · [ −∂(un+1T xy, n) ∂x − ∂(vn+1T xy, n) ∂y + ∂vn+1 ∂x T xx, n + ∂un+1 ∂y T yy, n − ε 1− β fnxy + (1− β) Wi ( ∂un+1 ∂y + ∂vn+1 ∂x ) − 1 Wi T xy, n ] , (3.24) onde fnxy = (T xx, n + T yy, n)T xy, n, T yy, n+1 = T yy, n + ∆t · [ −∂(un+1T yy, n) ∂x − ∂(vn+1T yy, n) ∂y + 2 ∂vn+1 ∂x T xy, n +2 ∂vn+1 ∂y T yy, n − ε 1− β fnyy + 2 (1− β) Wi ∂vn+1 ∂y − 1 Wi T yy, n ] . (3.25) onde fnyy = (T xx, n + T yy, n)T yy, n. 3.6 Atualização do tensor tensão - NSF Assim como feito para a formulação cartesiana, também é necessário atualizar os tensores da formulação NSF, quando esta é a adotada. O método utilizado para discretizar as equações (2.34)-(2.36), do fluido Oldroyd-B, no tempo é o de Euler Explícito e o resultado da discretização é: λ̂n+1 = λ̂n + ∆t · f λ̂NSF (3.26) µ̂n+1 = µ̂n + ∆t · f µ̂NSF (3.27) ν̂n+1 = ν̂n + ∆t · f ν̂NSF (3.28) com (3.29) f λ̂NSF = − 2µ̂n |un+1|2 + tol ( vn+1∂u n+1 ∂t − un+1∂v n+1 ∂t ) −|un+1|2 ( ∂(un+1λn) ∂x + ∂(vn+1λn) ∂y ) −2µ̂n|un+1|2∇·wn+1− 1 Wi ( λ̂n−1 ) , 3. Método numérico 38 onde λn = λ̂n |un+1|2+tol , (3.30) f µ̂NSF = − ( λ̂n − ν̂n |un+1|2 + tol )( un+1∂v n+1 ∂t − vn+1∂u n+1 ∂t ) − ( ∂(un+1µ̂n) ∂x + ∂(vn+1µ̂n) ∂y ) − ν̂n|un+1|2∇ ·wn+1 − µ̂n Wi , (3.31) f ν̂NSF = − 2µ̂n |un+1|2 + tol ( un+1∂v n+1 ∂t − vn+1∂u n+1 ∂t ) − 1 |un+1|2 + tol ( ∂(un+1νn) ∂x + ∂(vn+1νn) ∂y ) − 1 Wi (ν̂n − 1) , onde νn = ν̂n|un+1|2 e |un+1|2∇·wn+1 = 1 |un+1|2 + tol [ [(vn+1)2−(un+1)2] ( ∂un+1 ∂y + ∂vn+1 ∂x ) +4un+1vn+1∂u n+1 ∂x ] . (3.32) Para o modelo PTT, seguindo o mesmo procedimento, também pelo método de Euler Explícito discretiza-se as equações (2.39)-(2.41) e obtem-se λ̂n+1 = λ̂n + ∆t · gλ̂NSF (3.33) µ̂n+1 = µ̂n + ∆t · gµ̂NSF (3.34) ν̂n+1 = ν̂n + ∆t · gν̂NSF (3.35) com (3.36) gλ̂NSF = − 2µ̂n |un+1|2 + tol ( vn+1∂u n+1 ∂t − un+1∂v n+1 ∂t ) −|un+1|2 ( ∂(un+1λn) ∂x + ∂(vn+1λn) ∂y ) −2µ̂n|un+1|2∇·wn+1− αn Wi ( λ̂n−1 ) , onde λn = λ̂n |un+1|2+tol , (3.37) gµ̂NSF = − ( λ̂n − ν̂n |un+1|2 + tol )( un+1∂v n+1 ∂t − vn+1∂u n+1 ∂t ) − ( ∂(un+1µ̂n) ∂x + ∂(vn+1µ̂n) ∂y ) − ν̂n|un+1|2∇ ·wn+1 − αnµ̂n Wi , (3.38) gν̂NSF = − 2µ̂n |un+1|2 + tol ( un+1∂v n+1 ∂t − vn+1∂u n+1 ∂t ) − 1 |un+1|2 + tol ( ∂(un+1νn) ∂x + ∂(vn+1νn) ∂y ) − αn Wi (ν̂n − 1) , 3. Método numérico 39 onde νn = ν̂n|un+1|2, e |un+1|2∇·wn+1 = 1 |un+1|2 + tol [ [(vn+1)2−(un+1)2] ( ∂un+1 ∂y + ∂vn+1 ∂x ) +4un+1vn+1∂u n+1 ∂x ] . (3.39) Note que, para as três últimas equações, temos αn = 1 + ε(λ̂n − 2 + ν̂n). Nessas equações foram acrescidas as parcelas tol em determinados denominadores pois os vetores da base (2.28) se degeneram em regiões próximas ao contorno e pontos onde a velocidade é nula. Assim, para evitar possíveis divisões por zero, a constante tol é adicionada aos denominadores. 3.7 Discretização da lei Navier linear A discretização das condições slip é iniciada pela equação (2.70), referente a uma parede horizontal superior, ilustrada na Figura 3.3 abaixo. B F F F F B B F F Figura 3.3: Discretização da condição slip em uma parede horizontal. Calculando up = ui− 1 2 ,j+ 1 2 , com (2.70), tem-se (3.40) ui− 1 2 ,j+ 1 2 = −k ( β ∂u ∂y + T xy )∣∣∣∣ i− 1 2 ,j+ 1 2 = −k [ β ( aui− 1 2 ,j + bui− 1 2 ,j+ 1 2 + cui− 1 2 ,j+1 ) + T xy i− 1 2 ,j+ 1 2 ] , onde a = − hy2 hy1(hy1 + hy2) , b = hy2 − hy1 hy1hy2 , c = hy1 hy2(hy1 + hy2) , (3.41) e hy1 = 1 2 (∆yj−1 + ∆yj), hy2 = 1 2 (∆yj + ∆yj+1 ). (3.42) Ao aproximar a derivada ∂u ∂y , o ponto ui− 1 2 ,j+1, que está fora do domínio, é inserido na equação. É utilizada extrapolação de ordem 1 para aproximá-lo, isto é, 3. Método numérico 40 ui− 1 2 ,j+1 = 2ui− 1 2 ,j+ 1 2 − ui− 1 2 ,j, (3.43) e, substituindo em (3.40) ui− 1 2 ,j+ 1 2 = −k { β [ aui− 1 2 ,j + bup + c ( 2ui− 1 2 ,j+ 1 2 − ui− 1 2 ,j )] + T xy i− 1 2 ,j+ 1 2 } . (3.44) Com simplificações a expressão acima torna-se ui− 1 2 ,j+ 1 2 = − kβ(a− c) 1 + kβ (b+ 2c) ui− 1 2 ,j − k 1 + kβ (b+ 2c) T xy i− 1 2 ,j+ 1 2 . (3.45) Note que é necessário calcular T xy i− 1 2 ,j+ 1 2 em um ponto onde não está definido. Para isto, é feita a seguinte média T xy i− 1 2 ,j+ 1 2 = T xy i−1,j+ 1 2 − T xy i,j+ 1 2 2 , (3.46) e, em seguida, a extrapolação de grau 0 abaixo T xy i−1,j+ 1 2 = T xyi−1,j, (3.47) T xy i,j+ 1 2 = T xyi,j . (3.48) Finalmente, substituindo em (3.45), a velocidade na parede é ui− 1 2 ,j+ 1 2 = − kβ(a− c) 1 + kβ (b+ 2c) ui− 1 2 ,j − k 2 [1 + kβ (b+ 2c)] ( T xyi−1,j + T xyi,j ) . (3.49) Em relação a uma parede horizontal inferior (Figura 2.1b), a partir da equação (2.65) e, com cálculos semelhantes aos feitos no caso da parede superior, a expressão encontrada para up = ui− 1 2 ,j− 1 2 é ui− 1 2 ,j− 1 2 = kβ(c− a) 1− kβ (b+ 2a) ui− 1 2 ,j + k 2 [1− kβ (b+ 2a)] ( T xyi−1,j + T xyi,j ) , (3.50) em que os termos a, b e c são os mesmos de (3.41) Considera-se, agora, paredes verticais. Em uma parede à direita, por exemplo, re- presentada pela Figura 3.4, a velocidade normal à parede é nula, porém a velocidade tangencial não é. A expressão para calcular vp será determinada de forma semelhante aos casos de paredes horizontais. A equação (2.5) é discretizada entre as células (i, j) e (i, j − 1), (3.51) vi+ 1 2 ,j− 1 2 = −k ( β ∂v ∂x + T xy )∣∣∣∣ i+ 1 2 ,j− 1 2 = −k [ β ( avi,j− 1 2 + bvi+ 1 2 ,j− 1 2 + cvi+1,j− 1 2 ) + T xy i+ 1 2 ,j− 1 2 ] . com a = − hx2 hx1(hx1 + hx2) , b = hx2 − hx1 hx1hx2 , c = hx1 hx2(hx1 + hx2) , (3.52) 3. Método numérico 41 B BF F BF F F F Figura 3.4: Discretização da condição slip em uma parede vertical. em que hx1 = 1 2 (∆xi−1 + ∆xi), hx2 = 1 2 (∆xi + ∆xi+1 ). (3.53) O ponto vi+1,j− 1 2 está fora do domínio e será aproximado por extrapolação linear: vi+1,j− 1 2 = 2vi+ 1 2 ,j− 1 2 − vi,j− 1 2 , (3.54) Substituindo (3.54) em (3.51), vi+ 1 2 ,j− 1 2 = −k { β [ avi,j− 1 2 + bvi+ 1 2 ,j− 1 2 + c ( 2vi+ 1 2 ,j− 1 2 − vi,j− 1 2 )] + T xy i+ 1 2 ,j− 1 2 } . (3.55) A equação (3.55), com as devidas manipulações, é reescrita como vi+ 1 2 ,j− 1 2 = − kβ(a− c) 1 + kβ (b+ 2c) vi,j− 1 2 − k 1 + kβ (b+ 2c) T xy i+ 1 2 ,j− 1 2 . (3.56) Assim como nos cálculos para paredes horizontais, a componente T xy i+ 1 2 ,j− 1 2 de (3.56) é aproximada por uma média aritmética e, em seguida, por extrapolação de grau 0. Desta forma, a expressão final é vi+ 1 2 ,j− 1 2 = − kβ(a− c) 1 + kβ (b+ 2c) vi,j− 1 2 − k 2 [1 + kβ (b+ 2c)] ( T xyi−1,j + T xyi,j ) . (3.57) De modo análogo, a expressão para uma parede à esquerda é vi− 1 2 ,j− 1 2 = kβ(c− a) 1− kβ (b+ 2a) vi,j− 1 2 + k 2 [1− kβ (b+ 2a)] ( T xyi−1,j + T xyi,j ) , (3.58) com os coeficientes a, b e c de (3.52). 3. Método numérico 42 3.8 Equações de movimento com influência slip Um aspecto importante da discretização é notar as mudanças que as condições de contorno geram na equação de movimento discretizada. No Apêndice B são apresentadas várias fórmulas de diferenças utilizadas na discretização das equações de movimento nas direções x e y. Substituindo tais fórmulas em (B.2) e (B.11), obtém-se, respectivamente, equações da forma: a1ũ n+1 i− 3 2 ,j + b1ũ n+1 i− 1 2 ,j−1 + c1ũ n+1 i− 1 2 ,j + d1ũ n+1 i− 1 2 ,j+1 + e1ũ n+1 i+ 1 2 ,j = zn1 , (3.59) a2ṽ n+1 i−1,j− 1 2 + b2ṽ n+1 i,j− 3 2 + c2ṽ n+1 i,j− 1 2 + d2ṽ n+1 i,j+ 1 2 + e2ṽ n+1 i+1,j− 1 2 = zn2 , (3.60) com os coeficientes em (3.59) dados por a1 = −β∆t Re 2 hx1(hx1 + hx2) , b1 = −β∆t Re 2 hy1(hy1 + hy2) , c1 = 1 + β∆t Re ( 2 hx1hx2 + 2 hy1hy2 ) , d1 = −β∆t Re 2 hy2(hy1 + hy2) , e1 = −β∆t Re 2 hx2(hx1 + hx2) , onde hx1 = ∆xi−1 , hx2 = ∆xi , hy1 = 1 2 (∆yj−1 + ∆yj) e hy2 = 1 2 (∆yj + ∆yj+1 ). Já os coeficientes de (3.60) são dados por a2 = −β∆t Re 2 hx1(hx1 + hx2) , b2 = −β∆t Re 2 hy1(hy1 + hy2) , c2 = 1 + β∆t Re ( 2 hx1hx2 + 2 hy1hy2 ) , d2 = −β∆t Re 2 hy2(hy1 + hy2) , e2 = −β∆t Re 2 hx2(hx1 + hx2) , onde hx1 = 1 2 (∆xi−1 + ∆xi), hx2 = 1 2 (∆xi + ∆xi+1 ), hy1 = ∆yj−1 e hy2 = ∆yj . As condições de contorno slip aplicadas às equações (3.59) e (3.60) são do tipo explícita, isto é, a velocidade na parede, up, é conhecida desde o instante de tempo precedente. Assim, sendo, por exemplo, ũn+1 p a velocidade intermediária no tempo atual, em casos de condição de contorno slip, a velocidade unp é utilizada. Considerando a geometria da contração, a discretização das equações de movimento é detalhada em três regiões de interesse: na quina da contração, em uma parede vertical e em uma parede horizontal. Estas regiões são destacadas na figura abaixo. 3. Método numérico 43 BF BF F F BB F F F F F B F F F F FF F B B B B 31 2 Figura 3.5: Algumas regiões de interesse da contração. Inicia-se pela Região 1, próxima à quina da contração. Esta célula é identificada por (i, j) e, assim, a equação de movimento na direção x discretizada na face esquerda desta célula é a1ũ n+1 i− 3 2 ,j + b1ũ n+1 i− 1 2 ,j−1 + c1ũ n+1 i− 1 2 ,j + d1ũ n+1 i− 1 2 ,j+1 + e1ũ n+1 i+ 1 2 ,j = zn1 , (3.61) F B F F F F FF F Figura 3.6: Pontos da discretização da equação em u na Região 1. e a equação na direção y, discretizada na face inferior da célula, é dada por a2ṽ n+1 i−1,j− 1 2 + b2ṽ n+1 i,j− 3 2 + c2ṽ n+1 i,j− 1 2 + d2ṽ n+1 i,j+ 1 2 + e2ṽ n+1 i+1,j− 1 2 = zn2 . (3.62) Observe que, em ambas as equações, pontos do contorno não estão presentes e, por- tanto, as equações não sofrem alterações. Em contrapartida, observa-se que nas Regiões 2 e 3, há modificações nas equações devido à inclusão de pontos das paredes. 3. Método numérico 44 F B F F F F FF F Figura 3.7: Pontos da discretização da equação em v na Região 1. A Região 2 trata-se de uma parede vertical; é preciso atentar-se às diferenças que a condição de contorno implica nas equações. Discretizando a equação de movimento na direção x na face esquerda da célula (i, j+2) tem-se B BF F BF F F F Figura 3.8: Pontos da discretização da equação em u na Região 2. a1ũ n+1 i− 3 2 ,j+2 + b1ũ n+1 i− 1 2 ,j+1 + c1ũ n+1 i− 1 2 ,j+2 + d1ũ n+1 i− 1 2 ,j+3 + e1u n p = zn1 (3.63) Como em paredes verticais up = (0, vp), isto é, pelo fato do fluido não penetrar a parede, a equação anterior reduz-se a a1ũ n+1 i− 3 2 ,j+2 + b1ũ n+1 i− 1 2 ,j+1 + c1ũ n+1 i− 1 2 ,j+2 + d1ũ n+1 i− 1 2 ,j+3 = zn1 . (3.64) Agora, discretizando a equação de movimento na direção y na face inferior da célula (i, j + 2), encontra-se a2ṽ n+1 i−1,j+ 3 2 + b2ṽ n+1 i,j+ 1 2 + c2ṽ n+1 i,j+ 3 2 + d2ṽ n+1 i,j+ 5 2 + e2ṽ n+1 i+1,j+ 3 2 + = zn2 . (3.65) 3. Método numérico 45 B BF F BF F F F Figura 3.9: Pontos da discretização da equação em v na Região 2. Note que o ponto ṽn+1 i+1,j+ 3 2 está fora do domínio e, como feito na Seção 3.7, será apro- ximado por extrapolação linear, ou seja, ṽn+1 i+1,j+ 3 2 = vn i+1,j+ 3 2 = 2vn i+ 1 2 ,j+ 3 2 − vn i,j+ 3 2 = 2vp − vni,j+ 3 2 . (3.66) Substituindo (3.66) em (3.65), a equação é reescrita como a2ṽ n+1 i−1,j+ 3 2 + b2ṽ n+1 i,j+ 1 2 + c2ṽ n+1 i,j+ 3 2 + d2ṽ n+1 i,j+ 5 2 = zn2 − e2 ( 2vp − vni,j+ 3 2 ) , (3.67) onde vp = − kβ(a− c) 1 + kβ (b+ 2c) vi,j+ 3 2 − k 2 [1 + kβ (b+ 2c)] ( T xyi,j+1 + T xyi,j+2 ) (3.68) Na Região 3, quando discretiza-se a equação de movimento na direção x na face esquerda da célula (i+ 2, j), a equação resultante é a1ũ n+1 i+ 1 2 ,j + b1ũ n+1 i+ 3 2 ,j−1 + c1ũ n+1 i+ 3 2 ,j + d1ũ n+1 i+ 3 2 ,j+1 + e1ũ n+1 i+ 5 2 ,j = zn1 , (3.69) e o ponto ũn+1 i+ 3 2 ,j+1 é aproximado por ũn+1 i+ 3 2 ,j+1 = un i+ 3 2 ,j+1 = 2un i+ 3 2 ,j+ 1 2 − un i+ 3 2 ,j = 2up − uni+ 3 2 ,j , (3.70) Ao substituir (3.70) em (3.69), a equação é a1ũ n+1 i+ 1 2 ,j + b1ũ n+1 i+ 3 2 ,j−1 + c1ũ n+1 i+ 3 2 ,j + e1ũ n+1 i+ 5 2 ,j = zn1 − d1 ( 2up − uni+ 3 2 ,j ) (3.71) onde up = − kβ(a− c) 1 + kβ (b+ 2c) ui+ 3 2 ,j − k 2 [1 + kβ (b+ 2c)] ( T xyi+1,j + T xyi+2,j ) (3.72) Finalmente, discretizando a equação de movimento na direção y na face inferior da célula (i+ 2, j) tem-se a2ṽ n+1 i+1,j− 1 2 + b2ṽ n+1 i+2,j− 3 2 + c2ṽ n+1 i+2,j− 1 2 + d2ṽ n+1 i+2,j+ 1 2 + e2ṽ n+1 i+3,j− 1 2 = zn2 , (3.73) 3. Método numérico 46 BBB F F F F F F Figura 3.10: Pontos da discretização da equação em u na Região 3. Neste caso up = (up, 0) e, então, a última equação torna-se a2ṽ n+1 i+1,j− 1 2 + b2ṽ n+1 i+2,j− 3 2 + c2ṽ n+1 i+2,j− 1 2 + e2ṽ n+1 i+3,j− 1 2 = zn2 . (3.74) BBB F F F F F F Figura 3.11: Pontos da discretização da equação em v na Região 3. 3.9 Equações constitutivas - Oldroyd-B e PTT As condições de contorno slip (ou no-slip) do problema afetam também as equações constitutivas. Quando as equações dos modelos Oldroyd-B e PTT são atualizadas, seja pela formulação CSF ou pela NSF, nas derivadas presentes nessas equações encontram-se termos com as incógnitas un+1 e vn+1. Como esta é a última etapa do método numérico, as componentes un+1 e vn+1 trazem consigo informações das condições de contorno levadas em consideração na etapa do cálculo de ũn+1 e, consequentemente, na atualização de un+1. Portanto, no cálculo das derivadas presentes nas equações constitutivas, essas informações são inclusas nas aproximações por diferenças finitas. 3. Método numérico 47 3.10 Etapas do método numérico São descritos nesta seção os passos do método utilizado e consistem em • Passo 1: Calcular ũn+1 com a equação (3.15): ũn+1 − un ∆t +∇ · (uuᵀ)n = −∇pn + β Re ∇2ũn+1 + 1 Re ∇ ·Tn. • Passo 2: Calcular ψn+1 com a equação (3.4): ∇2ψ = ∇ · ũ. • Passo 3: Atualizar un+1 com a equação (3.1): u = ũ−∇ψ. • Passo 4: Atualizar pn+1 com a equação (3.19): pn+1 = pn + ψn+1 ∆t . • Passo 5: Atualizar Tn+1 com as equações: – Oldroyd-B ∗ (3.20)-(3.22) se a formulação ado- tada for a CSF: T xx,n+1 = T xx,n + ∆t · fxxCSF,OB, T xy,n+1 = T xy,n + ∆t · fxyCSF,OB, T yy,n+1 = T yy,n + ∆t · f yyCSF,OB, ∗ (3.26)-(3.28) se a formulação ado- tada for a NSF: λ̂n+1 = λ̂n + ∆t · f λ̂NSF,OB, µ̂n+1 = µ̂n + ∆t · f µ̂NSF,OB, ν̂n+1 = ν̂n + ∆t · f ν̂NSF,OB. – PTT ∗ (3.23)-(3.25) se a formulação ado- tada for a CSF: T xx,n+1 = T xx,n + ∆t · fxxCSF, PTT , T xy,n+1 = T xy,n + ∆t · fxyCSF, PTT , T yy,n+1 = T yy,n + ∆t · f yyCSF, PTT , ∗ (3.33)-(3.35) se a formulação ado- tada for a NSF: λ̂n+1 = λ̂n + ∆t · f λ̂NSF, PTT , µ̂n+1 = µ̂n + ∆t · f µ̂NSF, PTT , ν̂n+1 = ν̂n + ∆t · f ν̂NSF, PTT . Capítulo 4 Resultados Numéricos A seguir são apresentados alguns resultados de validação do código e dos estudos desenvolvidos combinando a formulação NSF com condições de contorno do tipo slip, bem como comparações entre as duas formulações adotadas para o tensor polimérico. As simulações são feitas nas geometrias da expansão e da contração, cujas proporções são 1:4 e 4:1, respectivamente. 4.1 Expansão 1:4 As simulações na geometria com expansão são realizadas em duas malhas não-uniformes, M1 e M2, que são mais refinadas em regiões de interesse da geometria: quinas superior e inferior e quinas que delimitam a transição do canal menor para o canal maior da expan- são. Os detalhes das malhas utilizadas e das regiões onde é mais refinada são apresentados na Tabela 4.1 e Figura 4.2, respectivamente. São apresentados resultados para os casos Newtoniano e viscoelástico, neste último, considerando-se o modelo PTT. 20L20L 1L 3L in flo w ou tfl owr XR θ Figura 4.1: Geometria da expansão 1:4 de comprimento 40L. 49 4. Resultados Numéricos 50 Tabela 4.1: Detalhes das malhas utilizadas nas simulações. Malha Espaçamento mín. Cél. na direção x Cél. na direção y Número de células M1 ∆x = ∆y = 0.008 190 260 32300 M2 ∆x = ∆y = 0.004 220 340 48400 Figura 4.2: Regiões em que a malha não-uniforme é mais refinada. 4.1.1 Fluido Newtoniano O objetivo inicial é validar os resultados seguindo a referência [15] em que os autores estudam a influência da velocidade de escorregamento no tamanho e intensidade dos vórtices, por exemplo. Para as simulações, é considerada a geometria com expansão, Figura 4.1, e os seguintes parâmetros adimensionais: • Comprimento característico: L = 1; • Número de Reynolds: Re = 0.01; • Passo temporal: ∆t = 10−4; • Inflow : perfil de velocidade desenvolvido com u U = 1.5; • Coeficiente de escorregamento: k̄ = 0, 0.001, 0.1, 1, 10, 45, 100, 10000. Conforme Ferrás [15], aumentar o valor do parâmetro k̄ implica na redução da dimen- são do vórtice, XR. No trabalho do autor citado, os resultados são ilustrados por meio das linhas de corrente do escoamento e por uma imagem semelhante à Figura 4.3. O maior valor testado pelos autores foi k̄ = 4500 e, ainda assim, havia presença de recirculação. Os resultados na Figura 4.3 e na Tabela 4.2 indicam que os valores obtidos nas si- mulações deste trabalho estão de acordo com o que foi descrito no artigo de referência, porém, como a malha aqui utilizada é mais grosseira, não foi possível capturar o vórtice nas simulações com k̄ = 45, 100 e 10000. 4. Resultados Numéricos 51 10-2 100 102 104 7k 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 X R Ferr4as [15] Nossos resultados - M2 Figura 4.3: Comprimentos de vórtices obtidos por Ferrás [15] comparados aos deste trabalho. Tabela 4.2: Tamanhos de vórtice, XR - Newtoniano, expansão. k̄ M2 0 1.503234 0.001 1.520757 0.1 1.442187 1.0 1.163254 10 0.531966 45 0.0 100 0.0 10000 0.0 Vários resultados no presente trabalho são a respeito do comportamento assintótico das propriedades próximo à singularidade da geometria adotada. Tal estudo consiste em avaliá-las à medida que se aproximam da quina, ou seja, sendo r a distância de um ponto até a quina, avalia-se como as propriedades se comportam quando r → 0 e verifica-se se estão em conformidade com as previsões teóricas. Além disso, é relevante ressaltar que essa análise pode ser feita considerando um ou mais ângulos θ e, nesta seção, os resultados são construídos tomando θ = π 2 como ilustrado pela Figura 4.4. Os resultados teóricos de [3], para o caso no-slip, indicam que a velocidade e componentes do tensor do solvente próximos à quina são dados por u ∝ r0.5445, Ts ∝ r−0.4555. (4.1) r θ Figura 4.4: Ângulo θ e distância r para o estudo assintótico. Desta forma, a Figura 4.5 traz o resultado da simulação com k̄ = 0 em que observa-se que o comportamento das propriedades está de acordo com a predição teórica enquanto a Figura 4.6, com gráficos para diferentes valores de k̄, exibe um distanciamento da reta para k̄ > 0. Isso era esperado uma vez que, para k̄ 6= 0, os resultados teóricos da variação assintótica são outros pois a condição de contorno envolve uma lei de escorregamento. 4. Resultados Numéricos 52 -2.5 -2 -1.5 -1 -0.5 0 log r -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 lo g ju j 7k = 0 slope= 0:5445 (a) -2.5 -2 -1.5 -1 -0.5 0 log r -1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 lo g jv j 7k = 0 slope= 0:5445 (b) -3 -2.5 -2 -1.5 -1 -0.5 0 log r -2 -1.5 -1 -0.5 0 0.5 1 lo g jT x x j 7k = 0 slope= !0:4555 (c) -3 -2.5 -2 -1.5 -1 -0.5 0 log r -2 -1.5 -1 -0.5 0 0.5 lo g jT x y j 7k = 0 slope= !0:4555 (d) -3 -2.5 -2 -1.5 -1 -0.5 0 log r -2 -1.5 -1 -0.5 0 0.5 1 lo g jT y y j 7k = 0 slope= !0:4555 (e) Figura 4.5: Comportamento assintótico da velocidade e das componentes do tensor Ts próximo à quina ao longo de θ = π 2 - Re = 0.01. 4. Resultados Numéricos 53 -2.5 -2 -1.5 -1 -0.5 0 log r -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 lo g ju j 7k = 0 7k = 0:1 7k = 1 7k = 10 7k = 100 slope= 0:5445 (a) -2.5 -2 -1.5 -1 -0.5 0 log r -3.5 -3 -2.5 -2 -1.5 -1 -0.5 lo g jv j 7k = 0 7k = 0:1 7k = 1 7k = 10 7k = 100 slope= 0:5445 (b) -3 -2.5 -2 -1.5 -1 -0.5 0 log r -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 lo g jT x x j 7k = 0 7k = 0:1 7k = 1 7k = 10 7k = 100 slope= !0:4555 (c) -3 -2.5 -2 -1.5 -1 -0.5 0 log r -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 lo g jT x y j 7k = 0 7k = 0:1 7k = 1 7k = 10 7k = 100 slope= !0:4555 (d) -3 -2.5 -2 -1.5 -1 -0.5 0 log r -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 lo g jT y y j 7k = 0 7k = 0:1 7k = 1 7k = 10 7k = 100 slope= !0:4555 (e) Figura 4.6: Comportamento assintótico das propriedades do escoamento para diferentes valores de k̄ ao longo de θ = π 2 - Re = 0.01. 4. Resultados Numéricos 54 4.1.2 PTT Nesta subseção [15] ainda é utilizada como referência, agora no estudo do fluido PTT. Os parâmetros adimensionais abaixo descritos são os utilizados nas simulações desta sub- seção em que são consideradas as formulações CSF e NSF. • Comprimento característico: L = 1; • Número de Reynolds: Re = 0.01; • Número de Weissenberg: Wi = 1; • ε = 0.25; • Razão das viscosidades: β = 1 9 ; • Passo temporal: ∆tCSF = 10−4, ∆tNSF = 10−5; • Tolerância da formulação NSF: tolNSF = 10−5; • Inflow : perfil de velocidade e de T desenvolvidos com u U = 1.5; • Coeficiente de escorregamento: k̄ = 0, 0.045, 0.18, 0.27, 0.36, 0.39, 0.41, 0.42, 0.45, 0.6, 0.65, 0.7, 45. Os resultados são validados comparando os tamanhos dos vórtices formados no canto superior da expansão. Na Tabela 4.3 e Figura 4.8 vê-se que os valores obtidos pela CSF, principalmente em M2, são bastante próximos dos de [15]. Já quando a formulação NSF é utilizada, os vórtices têm maior dimensão, quando comparados à CSF para um mesmo k̄, e só desaparecem para valores um pouco maiores de k̄. Na Figura 4.7 é ilustrada a redução no comprimento das recirculações por meio das linhas de corrente do escoa- mento considerando-se a formulação NSF, a malha M2 e alguns valores de k̄ selecionados. Observa-se na Figura 4.7e o vórtice de menor dimensão, obtido com k̄ = 0.7, o qual denomina-se neste trabalho “ k̄ de transição”. Tabela 4.3: Tamanhos de vórtice, XR - CSF e NSF - PTT, expansão. k̄ CSF NSF Ferrás [15] - CSFM1 M2 M1 M2 0 1.376331 1.358653 1.449045 1.420133 1.3665 0.045 1.299366 1.295675 1.295188 1.350113 1.2704 0.18 0.911696 0.885275 1.054136 1.061453 0.8991 0.27 0.627049 0.632389 0.911696 0.885275 0.6247 0.36 0.310571 0.349214 0.668600 0.669394 0.2756 0.39 0.203903 0.236603 - - - 0.41 0.057700 0.120444 - - - 0.42 0 0.013288 - - - 0.45 0 0 0.480676 0.472999 0 0.6 - - 0.203903 0.206044 - 0.65 - - 0.121444 0.130898 - 0.7 - - 0.057700 0.060359 - 45 0 0 0 0 0 4. Resultados Numéricos 55 (a) (b) (c) (d) (e) Figura 4.7: Variação no comprimento dos vórtices com o aumento de k̄ - NSF, M2. 4. Resultados Numéricos 56 10-2 10-1 100 101 102 7k 0 0.2 0.4 0.6 0.8 1 1.2 1.4 X R Ferr4as [15] CSF - M2 NSF - M2 Figura 4.8: Dimensão dos vórtices obtidos por [15] comparados aos deste trabalho - CSF e NSF. Os estudos teóricos da variação assintótica do fluido PTT detalhados em [27] indicam como as propriedades do escoamento devem portar-se próximo da quina da expansão considerando a condição de contorno no-slip. No caso da formulação CSF, as propriedades têm como esperado o seguinte comportamento u ∝ r0.5445, T ∝ r−0.3286, p ∝ r−0.4555, (4.2) e, para a NSF, as variáveis naturais têm como previsão λ ∝ r−1.4176, µ ∝ r−0.1268, ν ∝ r1.1638. (4.3) Inicialmente são comparadas as soluções numéricas obtidas com as formulações CSF e NSF, no caso no-slip, que são plotadas ao longo dos ângulos θ = π 2 , 3π 4 e π. As simulações são realizadas na malha M2 e o ângulo θ é medido a partir da parede horizontal do canal estreito no sentido anti-horário, conforme ilustrado pela Figura 4.9. Figura 4.9: Ângulos tomados para estudo do comportamento assintótico Por meio das Figuras 4.11, 4.12 e 4.13 vê-se que a formulação NSF apresenta me- lhor desempenho no cálculo das componentes de T e apresenta melhoria considerável na aproximação da pressão. Analisa-se, agora, a variação assintótica das variáveis naturais λ, µ e ν ao longo dos mesmos ângulos θ. Neste caso, segundo a Figura 4.14, também tem-se soluções adequadas 4. Resultados Numéricos 57 às potências previstas em (4.3) reforçando a validade e coerência dos resultados. No que diz respeito às simulações com diferentes valores de k̄, o comportamento assintótico das propriedades é exibido nas Figuras 4.15 e 4.16, com as formulações CSF e NSF, respecti- vamente. Para ambas as formulações observa-se que as soluções numéricas distanciam-se da reta com o aumento de k̄ mas que o comportamento da pressão no caso NSF, apa- rentemente, acompanha a inclinação da reta. Com atenção especial à pressão, analisada em θ = π 2 , 3π 4 e π, a Figura 4.17 indica que não há discrepâncias significativas entre as soluções numéricas produzidas por ambas as formulações. Vale destacar que as soluções numéricas para k̄ 6= 0 afastam-se da reta pois a inclinação de tal reta é referente ao caso k̄ = 0. Por fim, são apresentados os perfis das propriedades do escoamentos em três regiões da geometria: x ≈ 10, x = 20 e x ≈ 30, como ilustrado na Figura 4.10. 20 20 1 3 x ≈ 10 x = 20 x ≈ 30 Figura 4.10: Cortes verticais O fato de malhas não-uniformes serem utilizadas nas simulações nem sempre possibilita que o corte seja feito exatamente no ponto x desejado. Nesses casos, é tomado o ponto da malha que se encontra mais próximo do ponto desejado. Essa é a razão do símbolo “≈” acima. Como visto nas Figuras 4.18, 4.19 e 4.20 a condição de escorregamento provoca mu- danças nas propriedades nos três cortes estudados. Em particular, em x ≈ 10, as variáveis naturais λ̂ e µ̂ têm seus valores reduzidos com o aumento de k̄ ao passo que ν̂ é prati- camente uma constante de valor 1. Comportamentos bastante semelhantes das variáveis naturais ocorrem em x ≈ 30. Já em x = 20, λ̂ e ν̂ apresentam maiores variações próximo às paredes enquanto os valores de µ̂ variam ao longo do interior do canal. 4. Resultados Numéricos 58 -2.5 -2 -1.5 -1 -0.5 0 log r -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 lo g ju j 7k = 0 - CSF 7k = 0 - NSF slope= 0:5445 (a) -2.5 -2 -1.5 -1 -0.5 0 log r -1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 lo g jv j 7k = 0 - CSF 7k = 0 - NSF slope= 0:5445 (b) -3 -2.5 -2 -1.5 -1 -0.5 0 log r -2 -1.5 -1 -0.5 0 0.5 1 1.5 lo g jT x x j 7k = 0 - CSF 7k = 0 - NSF slope= !0:3286 (c) -3 -2.5 -2 -1.5 -1 -0.5 0 log r -2 -1.5 -1 -0.5 0 0.5 1 1.5 lo g jT x y j 7k = 0 - CSF 7k = 0 - NSF slope= !0:3286 (d) -3 -2.5 -2 -1.5 -1 -0.5 0 log r -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 lo g jT y y j 7k = 0 - CSF 7k = 0 - NSF slope= !0:3286 (e) -3 -2.5 -2 -1.5 -1 -0.5 0 log r 2.5 2.6 2.7 2.8 2.9 3 3.1 3.2 3.3 3.4 3.5 lo g jp j 7k = 0 - CSF 7k = 0 - NSF slope= !0:4555 (f) Figura 4.11: Comportamento assintótico das soluções numéricas próximo à quina ao longo de θ = π 2 - Re = 0.01, ε = 0.25, Wi = 1, β = 1 9 . 4. Resultados Numéricos 59 -2.5 -2 -1.5 -1 -0.5 0 log r -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 lo g ju j 7k = 0 - CSF 7k = 0 - NSF slope= 0:5445 (a) -2.5 -2 -1.5 -1 -0.5 0 log r -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 lo g jv j 7k = 0 - CSF 7k = 0 - NSF slope= 0:5445 (b) -3 -2.5 -2 -1.5 -1 -0.5 0 log r -2 -1.5 -1 -0.5 0 0.5 1 lo g jT x x j 7k = 0 - CSF 7k = 0 - NSF slope= !0:3286 (c) -3 -2.5 -2 -1.5 -1 -0.5 0 log r -2 -1.5 -1 -0.5 0 0.5 1 lo g jT x y j 7k = 0 - CSF 7k = 0 - NSF slope= !0:3286 (d) -3 -2.5 -2 -1.5 -1 -0.5 0 log r -2.5 -2 -1.5 -1 -0.5 0 0.5 1 lo g jT y y j 7k = 0 - CSF 7k = 0 - NSF slope= !0:3286 (e) -3 -2.5 -2 -1.5 -1 -0.5 0 log r 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3 3.1 lo g jp j 7k = 0 - CSF 7k = 0 - NSF slope= !0:4555 (f) Figura 4.12: Comportamento assintótico das soluções numéricas próximo à quina ao longo de θ = 3 4 π - Re = 0.01, ε = 0.25, Wi = 1, β = 1 9 . 4. Resultados Numéricos 60 -2.5 -2 -1.5 -1 -0.5 0 log r -1.3 -1.2 -1.1 -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 lo g ju j 7k = 0 - CSF 7k = 0 - NSF slope= 0:5445 (a) -2.5 -2 -1.5 -1 -0.5 0 log r -1.1 -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 lo g jv j 7k = 0 - CSF 7k = 0 - NSF slope= 0:5445 (b) -3 -2.5 -2 -1.5 -1 -0.5 0 log r -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 lo g jT x x j 7k = 0 - CSF 7k = 0 - NSF slope= !0:3286 (c) -3 -2.5 -2 -1.5 -1 -0.5 0 log r -2 -1.5 -1 -0.5 0 0.5 1 lo g jT x y j 7k = 0 - CSF 7k = 0 - NSF slope= !0:3286 (d) -3 -2.5 -2 -1.5 -1 -0.5 0 log r -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 lo g jT y y j 7k = 0 - CSF 7k = 0 - NSF slope= !0:3286 (e) -3 -2.5 -2 -1.5 -1 -0.5 0 log r 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3 lo g jp j 7k = 0 - CSF 7k = 0 - NSF slope= !0:4555 (f) Figura 4.13: Comportamento assintótico das soluções numéricas próximo à quina ao longo de θ = π - Re = 0.01, ε = 0.25, Wi = 1, β = 1 9 . 4. Resultados Numéricos 61 -3 -2.5 -2 -1.5 -1 -0.5 0 log r -3 -2 -1 0 1 2 3 4 log j6j log j7j log j8j slope= !1:4176 slope= !0:1268 slope= 1:1638 (a) θ = π 2 -3 -2.5 -2 -1.5 -1 -0.5 0 log r -2 -1 0 1 2 3 4 log j6j log j7j log j8j slope= !1:4176 slope= !0:1268 slope= 1:1638 (b) θ = 3π 4 -3 -2.5 -2 -1.5 -1 -0.5 0 log r -3 -2 -1 0 1 2 3 4 log j6j log j7j log j8j slope= !1:4176 slope= !0:1268 slope= 1:1638 (c) θ = π Figura 4.14: Variação assintótica das variáveis naturais ao longo de diferentes ângulos θ - Re = 0.01, ε = 0.25, Wi = 1, β = 1 9 . 4. Resultados Numéricos 62 -2.5 -2 -1.5 -1 -0.5 0 log r -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 lo g ju j 7k = 0 7k = 0:18 7k = 0:27 7k = 0:36 7k = 45 slope= 0:5445 (a) -2.5 -2 -1.5 -1 -0.5 0 log r -2.2 -2 -1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 lo g jv j 7k = 0 7k = 0:18 7k = 0:27 7k = 0:36 7k = 45 slope= 0:5445 (b) -3 -2.5 -2 -1.5 -1 -0.5 0 log r -2 -1.5 -1 -0.5 0 0.5 1 lo g jT x x j 7k = 0 7k = 0:18 7k = 0:27 7k = 0:36 7k = 45 slope= !0:3286 (c) -3 -2.5 -2 -1.5 -1 -0.5 0 log r -2.5 -2 -1.5 -1 -0.5 0 0.5 1 lo g jT x y j 7k = 0 7k = 0:18 7k = 0:27 7k = 0:36 7k = 45 slope= !0:3286 (d) -3 -2.5 -2 -1.5 -1 -0.5 0 log r -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 lo g jT y y j 7k = 0 7k = 0:18 7k = 0:27 7k = 0:36 7k = 45 slope= !0:3286 (e) -3 -2.5 -2 -1.5 -1 -0.5 0 log r 2.2 2.4 2.6 2.8 3 3.2 3.4 lo g jp j 7k = 0 7k = 0:18 7k = 0:27 7k = 0:36 7k = 45 slope= !0:4555 (f) Figura 4.15: Comportamento assintótico próximo à quina ao longo de θ = π 2 - CSF, Re = 0.01, ε = 0.25, Wi = 1, β = 1 9 . 4. Resultados Numéricos 63 -2.5 -2 -1.5 -1 -0.5 0 log r -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 lo g ju j 7k = 0 7k = 0:18 7k = 0:27 7k = 0:36 7k = 45 slope= 0:5445 (a) -2.5 -2 -1.5 -1 -0.5 0 log r -2.2 -2 -1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 lo g jv j 7k = 0 7k = 0:18 7k = 0:27 7k = 0:36 7k = 45 slope= 0:5445 (b) -3 -2.5 -2 -1.5 -1 -0.5 0 log r -0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 lo g j6 j 7k = 0 7k = 0:18 7k = 0:27 7k = 0:36 7k = 45 slope= !1:4176 (c) -3 -2.5 -2 -1.5 -1 -0.5 0 log r -2.5 -2 -1.5 -1 -0.5 0 0.5 lo g j7 j 7k = 0 7k = 0:18 7k = 0:27 7k = 0:36 7k = 45 slope= !0:12689 (d) -3 -2.5 -2 -1.5 -1 -0.5 0 log r -2.5 -2 -1.5 -1 -0.5 0 0.5 lo g j8 j 7k = 0 7k = 0:18 7k = 0:27 7k = 0:36 7k = 45 slope= 1:1638 (e) -3 -2.5 -2 -1.5 -1 -0.5 0 log r 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 lo g jp j 7k = 0 7k = 0:18 7k = 0:27 7k = 0:36 7k = 45 slope= !0:4555 (f) Figura 4.16: Comportamento assintótico próximo à quina ao longo de θ = π 2 - NSF, Re = 0.01, ε = 0.25, Wi = 1, β = 1 9 . 4. Resultados Numéricos 64 -3 -2.5 -2 -1.5 -1 -0.5 0 log r 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 lo g jp j 7k = 0 - CSF 7k = 0 - NSF 7k = 0:27 - CSF 7k = 0:27 - NSF 7k = 45 - CSF 7k = 45 - NSF slope = !0:4555 (a) -3 -2.5 -2 -1.5 -1 -0.5 0 log r 1.8 2 2.2 2.4 2.6 2.8 3 lo g jp j 7k = 0 - CSF 7k = 0 - NSF 7k = 0:27 - CSF 7k = 0:27 - NSF 7k = 45 - CSF 7k = 45 - NSF slope = !0:4555 (b) -3 -2.5 -2 -1.5 -1 -0.5 0 log r 1.6 1.8 2 2.2 2.4 2.6 2.8 3 lo g jp j 7k = 0 - CSF 7k = 0 - NSF 7k = 0:27 - CSF 7k = 0:27 - NSF 7k = 45 - CSF 7k = 45 - NSF slope = !0:4555 (c) Figura 4.17: Comportamento assintótico da pressão ao longo de θ = π 2 , 3π 4 e π - Re = 0.01, ε = 0.25, Wi = 1, β = 1 9 . 4. Resultados Numéricos 65 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 y 0 0.5 1 1.5 u 7k = 0 7k = 0:18 7k = 0:45 7k = 0:7 7k = 45 (a) 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 y -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 v #10-6 7k = 0 7k = 0:18 7k = 0:45 7k = 0:7 7k = 45 (b) 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 y 0 1 2 3 4 5 6 7 6̂ 7k = 0 7k = 0:18 7k = 0:45 7k = 0:7 7k = 45 (c) 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 y -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 7̂ 7k = 0 7k = 0:18 7k = 0:45 7k = 0:7 7k = 45 (d) 3 3.5 4 4.5 5 y 0.9999 0.99995 1 1.00005 1.0001 1.00015 1.0002 1.00025 1.0003 1.00035 8̂ 7k = 0 7k = 0:18 7k = 0:45 7k = 0:7 7k = 45 (e) 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 y 200 400 600 800 1000 1200 1400 1600 1800 2000 2200 p 7k = 0 7k = 0:18 7k = 0:45 7k = 0:7 7k = 45 (f) Figura 4.18: Perfis das propriedades do escoamento com a variação de k̄ em x ≈ 10 - NSF, Re = 0.01, ε = 0.25, Wi = 1, β = 1 9 . 4. Resultados Numéricos 66 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 y 0 0.2 0.4 0.6 0.8 1 1.2 1.4 u 7k = 0 7k = 0:18 7k = 0:45 7k = 0:7 7k = 45 (a) 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 y -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 v 7k = 0 7k = 0:18 7k = 0:45 7k = 0:7 7k = 45 (b) 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 y 0 5 10 15 20 25 30 35 40 6̂ 7k = 0 7k = 0:18 7k = 0:45 7k = 0:7 7k = 45 (c) 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 y -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 7̂ 7k = 0 7k = 0:18 7k = 0:45 7k = 0:7 7k = 45 (d) 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 y 0 1 2 3 4 5 6 8̂ 7k = 0 7k = 0:18 7k = 0:45 7k = 0:7 7k = 45 (e) 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 y 0 500 1000 1500 2000 2500 3000 p 7k = 0 7k = 0:18 7k = 0:45 7k = 0:7 7k = 45 (f) Figura 4.19: Perfis das propriedades do escoamento com a variação de k̄ em x = 20 - NSF, Re = 0.01, ε = 0.25, Wi = 1, β = 1 9 . 4. Resultados Numéricos 67 0 1 2 3 4 5 6 7 8 y 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 u 7k = 0 7k = 0:18 7k = 0:45 7k = 0:7 7k = 45 (a) 0 1 2 3 4 5 6 7 8 y -8 -6 -4 -2 0 2 4 6 8 v #10-4 7k = 0 7k = 0:18 7k = 0:45 7k = 0:7 7k = 45 (b) 0 1 2 3 4 5 6 7 8 y 0.99 1 1.01 1.02 1.03 1.04 1.05 1.06 1.07 6̂ 7k = 0 7k = 0:18 7k = 0:45 7k = 0:7 7k = 45 (c) 0 1 2 3 4 5 6 7 8 y -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 7̂ 7k = 0 7k = 0:18 7k = 0:45 7k = 0:7 7k = 45 (d) 0 1 2 3 4 5 6 7 8 y 0.9985 0.999 0.9995 1 1.0005 1.001 1.0015 8̂ 7k = 0 7k = 0:18 7k = 0:45 7k = 0:7 7k = 45 (e) 0 1 2 3 4 5 6 7 8 y 0 5 10 15 20 25 30 35 40 45 50 p 7k = 0 7k = 0:18 7k = 0:45 7k = 0:7 7k = 45 (f) Figura 4.20: Perfis das propriedades do escoamento com a variação de k̄ em x ≈ 30 - NSF, Re = 0.01, ε = 0.25, Wi = 1, β = 1 9 . 4. Resultados Numéricos 68 4.2 Contração 4:1 As simulações nessa seção são realizadas numa contração de comprimento 40L ilus- trada pela figura abaixo e nas mesmas malhas detalhadas na Tabela 4.1. A Figura 4.22 ilustra a estrutura da malha não-uniforme. 20L 20L 1L 3L in flo w ou tfl owr XR θ Figura 4.21: Representação da contração 4:1, regiões de formação de vórtices e singulari- dade onde é feita a análise assintótica. Figura 4.22: Regiões mais refinadas da malha não-uniforme. 4. Resultados Numéricos 69 4.2.1 Fluido Newtoniano Nas simulações iniciais é considerado um fluido Newtoniano (β = 1) e os parâmetros adimensionais: • Comprimento característico: L = 1; • Número de Reynolds: Re = 0.01; • Passo temporal: ∆t = 10−4; • Inflow : perfil de velocidade desenvolvido, com u U = 0.375; • Coeficiente de escorregamento: k̄ = 0, 0.45, 1, 10, 100. Para a validação, são considerados resultados como tamanhos dos vórtices, XR, e comportamento assintótico das componentes da velocidade, do tensor do solvente e pressão próximo à singularidade para o caso no-slip. Os resultados obtidos são comparados com os presentes em [16], Figura 4.23, e mostram que, à medida que o valor da constante de fricção k̄ é aumentado, o tamanho do vórtice diminui até que não haja recirculação (k̄ ≈ 50); os valores obtidos neste trabalho são organizados na Tabela 4.4. Comparando esses valores com os produzidos por Ferrás [16], vê-se que os valores obtidos estão coerentes. 10-3 10-2 10-1 100 101 102 7k 0 0.5 1 1.5 X R Ferr4as [16] Nossos resultados - M2 Figura 4.23: Resultados do caso Newtoniano obtidos por [16] comparados aos deste tra- balho. Tabela 4.4: Tamanhos de vórtice, XR - Newtoniano, contração. k̄ M1 M2 0.0 1.465579 1.525190 0.45 1.164632 1.231791 1.0 0.976286 0.992115 10 0.369544 0.257082 100 0.0 0.0 4. Resultados Numéricos 70 O artigo apresenta ainda resultados do comportamento assintótico de algumas pro- priedades próximo à singularidade. Esses resultados teóricos são os mesmos de (4.1) e os resultados apresentados na Figura 4.25 confirmam as predições teóricas, com exceção da pressão, do caso no-slip para as propriedades do escoamento, isto é, ao se aproximarem da quina (r → 0), seus valores se aproximam da inclinação das respectivas retas. Isso não acontece para outros valores de k̄, conforme a Figura 4.26; os valores das propriedades se afastam da reta relativa a k̄ = 0. Ainda com simulações em que a constante de fricção k̄ assume diferentes valores, são analisados os perfis das propriedades do escoamento e os cortes são feitos em três regiões da geometria: x ≈ 10, x = 20 e x ≈ 30, como na Figura 4.24. 20 20 1 3 x ≈ 10 x = 20 x ≈ 30 y = 4 Figura 4.24: Cortes horizontal e vertical na geometria (L = 1). As Figuras 4.27, 4.28 e 4.29 exibem de forma clara a influência da condição slip nos perfis das propriedades do escoamento nas três regiões da geometria. Voltando a atenção à componente u da velocidade, é bastante evidente seu aumento na parede à medida que k̄ também aumenta, o que era esperado. No corte em x = 20, em particular, o aumento do valor de u também ocorre porém sem grandes variações com os diferentes valores do coeficiente de escorregamento. Em termos da pressão, uma diminuição em seu valor acontece nos três cortes estudados com a variação de k̄, atingindo os menores valores em x ≈ 30. Outro detalhe relevante é a variação da componente u da velocidade ao longo da geometria da contração na altura y = 4. A Figura 4.30a ilustra essa variação enquanto a Figura 4.30b apresenta a mesma análise para a pressão. Observa-se que aumentar k̄ faz com que, no centro do canal, ambas as propriedades tenham seus valores reduzidos. 4. Resultados Numéricos 71 -2.5 -2 -1.5 -1 -0.5 0 log r -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 lo g ju j 7k = 0!M1 7k = 0!M2 slope= 0:5445 (a) -2.5 -2 -1.5 -1 -0.5 0 log r -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1 -0.9 -0.8 -0.7 lo g jv j 7k = 0!M1 7k = 0!M2 slope= 0:5445 (b) -3 -2.5 -2 -1.5 -1 -0.5 0 log r -2 -1.5 -1 -0.5 0 0.5 1 lo g jT x x j 7k = 0!M1 7k = 0!M2 slope= !0:4555 (c) -3 -2.5 -2 -1.5 -1 -0.5 0 log r -2 -1.5 -1 -0.5 0 0.5 lo g jT x y j 7k = 0!M1 7k = 0!M2 slope= !0:4555 (d) -3 -2.5 -2 -1.5 -1 -0.5 0 log r -2 -1.5 -1 -0.5 0 0.5 1 lo g jT y y j 7k = 0!M1 7k = 0!M2 slope= !0:4555 (e) -3 -2.5 -2 -1.5 -1 -0.5 0 log r 3.778 3.779 3.78 3.781 3.782 3.783 3.784 3.785 3.786 lo g jp j 7k = 0!M1 7k = 0!M2 slope= !0:4555 (f) Figura 4.25: Comportamento assintótico da velocidade, das componentes do tensor Ts e da pressão próximo à quina - Re = 0.01. 4. Resultados Numéricos 72 -2.5 -2 -1.5 -1 -0.5 0 log r -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 lo g ju j 7k = 0 7k = 0:45 7k = 1 7k = 10 7k = 100 slope= 0:5445 (a) -2.5 -2 -1.5 -1 -0.5 0 log r -3.5 -3 -2.5 -2 -1.5 -1 -0.5 lo g jv j 7k = 0 7k = 0:45 7k = 1 7k = 10 7k = 100 slope= 0:5445 (b) -3 -2.5 -2 -1.5 -1 -0.5 0 log r -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 lo g jT x x j 7k = 0 7k = 0:45 7k = 1 7k = 10 7k = 100 slope= !0:4555 (c) -3 -2.5 -2 -1.5 -1 -0.5 0 log r -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 lo g jT x y j 7k = 0 7k = 0:45 7k = 1 7k = 10 7k = 100 slope= !0:4555 (d) -3 -2.5 -2 -1.5 -1 -0.5 0 log r -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 lo g jT y y j 7k = 0 7k = 0:45 7k = 1 7k = 10 7k = 100 slope= !0:4555 (e) -3 -2.5 -2 -1.5 -1 -0.5 0 log r 0.5 1 1.5 2 2.5 3 3.5 4 lo g jp j 7k = 0 7k = 0:45 7k = 1 7k = 10 7k = 100 slope= !0:4555 (f) Figura 4.26: Comportamento assintótico da velocidade, das componentes do tensor Ts e da pressão próximo à quina - Re = 0.01. 4. Resultados Numéricos 73 0 1 2 3 4 5 6 7 8 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 (a) 0 1 2 3 4 5 6 7 8 -1.5 -1 -0.5 0 0.5 1 1.5 10-3 (b) 0 1 2 3 4 5 6 7 8 -1.5 -1 -0.5 0 0.5 1 1.5 10-3 (c) 0 1 2 3 4 5 6 7 8 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 (d) 0 1 2 3 4 5 6 7 8 -1.5 -1 -0.5 0 0.5 1 1.5 10-3 (e) 0 1 2 3 4 5 6 7 8 0 1000 2000 3000 4000 5000 6000 7000 (f) Figura 4.27: Perfis das propriedades do escoamento com a variação de k̄ em x ≈ 10 - Re = 0.01. 4. Resultados Numéricos 74 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 (a) 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 (b) 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 -6 -5 -4 -3 -2 -1 0 1 (c) 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 -1.5 -1 -0.5 0 0.5 1 1.5 (d) 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 -1 0 1 2 3 4 5 6 (e) 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 -3000 -2000 -1000 0 1000 2000 3000 4000 5000 6000 7000 (f) Figura 4.28: Perfis das propriedades do escoamento com a variação de k̄ em x = 20 - Re = 0.01. 4. Resultados Numéricos 75 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 0 0.5 1 1.5 (a) 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 (b) 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 10-12 (c) 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 -3 -2 -1 0 1 2 3 (d) 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 10-12 (e) 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 0 500 1000 1500 2000 2500 3000 3500 (f) Figura 4.29: Perfis das propriedades do escoamento com a variação de k̄ em x ≈ 30 - Re = 0.01. 4. Resultados Numéricos 76 0 5 10 15 20 25 30 35 40 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 (a) 0 5 10 15 20 25 30 35 40 0 1000 2000 3000 4000 5000 6000 7000 (b) Figura 4.30: Perfil horizontal de (a) u e (b) p em y = 4 - Re = 0.01. 4.2.2 Oldroyd-B Os próximos resultados dizem respeito à simulação com o fluido Oldroyd-B, utilizando as formulações CSF e NSF, com os seguintes parâmetros adimensionais: • Comprimento característico: L = 1; • Razão das viscosidades: β = 1 2 ; • Tolerância da formulação NSF: tolNSF=10−6; • Número de Reynolds: Re = 0.01; • Número de Weissenberg: Wi = 1; • Passo temporal: ∆tCSF = 5× 10−5, ∆tNSF = 10−5; • Inflow : perfis de velocidade e de T desenvolvidos, com u U = 0.375; • Coeficiente de escorregamento: k̄ = 0, 0.1, 0.18, 0.5. O principal objetivo deste teste é comparar a precisão de ambas as formulações na captura dos valores dos tensores e da pressão. São comparados os resultados assintóticos de T xx, T xy, T yy e p produzidos pelas formulações bem como o comportamento das variáveis naturais λ, µ e ν na malha M1. Os artigos de referência [6] e [18] detalham os cálculos para a obtenção das proporcionalidades da velocidade u e componentes do tensor T. Segundo esses autores, o comportamento assintótico respeita u ∝ r 5 9 , T ∝ r− 2 3 , (4.4) e as variáveis naturais, conforme [11], são dadas por λ ∝ r−1.8310, µ ∝ r−0.28651, ν ∝ r1.2580, (4.5) válidos somente para o caso em que k̄ = 0. A Figura 4.31 mostra que, de fato, a NSF aproxima melhor os valores das componentes de T perto da singularidade e, inclusive, consideravelmente, a pressão conforme r → 0. Já para as componentes da velocidade, a diferença é menos significativa. 4. Resultados Numéricos 77 -2.5 -2 -1.5 -1 -0.5 0 log r -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 lo g ju j 7k = 0 - CSF 7k = 0 - NSF slope= 5=9 (a) -2.5 -2 -1.5 -1 -0.5 0 log r -1.4 -1.3 -1.2 -1.1 -1 -0.9 -0.8 -0.7 -0.6 lo g jv j 7k = 0 - CSF 7k = 0 - NSF slope= 5=9 (b) -2.5 -2 -1.5 -1 -0.5 0 log r -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 lo g jT x x j 7k = 0 - CSF 7k = 0 - NSF slope= !2=3 (c) -2.5 -2 -1.5 -1 -0.5 0 log r -2 -1.5 -1 -0.5 0 0.5 1 1.5 lo g jT x y j 7k = 0 - CSF 7k = 0 - NSF slope= !2=3 (d) -2.5 -2 -1.5 -1 -0.5 0 log r -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 lo g jT y y j 7k = 0 - CSF 7k = 0 - NSF slope= !2=3 (e) -2.5 -2 -1.5 -1 -0.5 0 log r 3.75 3.8 3.85 3.9 3.95 4 lo g jp j 7k = 0 - CSF 7k = 0 - NSF slope= !2=3 (f) Figura 4.31: Comparação do comportamento assintótico das propriedades do escoamento utilizando as formulações CSF e NSF na malha M1 - Re = 0.01, Wi = 1 e β = 1 2 . 4. Resultados Numéricos 78 Em seguida são apresentados resultados assintóticos com os valores k̄ = 0, 0.1, 0.18 e 0.5 produzidos também por ambas as formulações com a observação de que a inclinação presente nas Figuras 4.32 e 4.33 é referente a k̄ = 0. -2.5 -2 -1.5 -1 -0.5 0 log r -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 lo g ju j 7k = 0 7k = 0:1 7k = 0:18 7k = 0:5 slope= 5=9 (a) -2.5 -2 -1.5 -1 -0.5 0 log r -2 -1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 lo g jv j 7k = 0 7k = 0:1 7k = 0:18 7k = 0:5 slope= 5=9 (b) -2.5 -2 -1.5 -1 -0.5 0 log r -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 lo g jT x x j 7k = 0 7k = 0:1 7k = 0:18 7k = 0:5 slope= !2=3 (c) -2.5 -2 -1.5 -1 -0.5 0 log r -2.5 -2 -1.5 -1 -0.5 0 0.5 1 lo g jT x y j 7k = 0 7k = 0:1 7k = 0:18 7k = 0:5 slope= !2=3 (d) -2.5 -2 -1.5 -1 -0.5 0 log r -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 lo g jT y y j 7k = 0 7k = 0:1 7k = 0:18 7k = 0:5 slope= !2=3 (e) -2.5 -2 -1.5 -1 -0.5 0 log r 3.4 3.45 3.5 3.55 3.6 3.65 3.7 3.75 3.8 3.85 lo g jp j 7k = 0 7k = 0:1 7k = 0:18 7k = 0:5 slope= !2=3 (f) Figura 4.32: Propriedades do escoamento ao se aproximarem da quina da contração - CSF, Re = 0.01, Wi = 1 e β = 1 2 . 4. Resultados Numéricos 79 -2.5 -2 -1.5 -1 -0.5 0 log r -1 -0.8 -0.6 -0.4 -0.2 0 0.2 lo g ju j 7k = 0 7k = 0:1 7k = 0:18 7k = 0:5 slope= 5=9 (a) -2.5 -2 -1.5 -1 -0.5 0 log r -2 -1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 lo g jv j 7k = 0 7k = 0:1 7k = 0:18 7k = 0:5 slope= 5=9 (b) -2.5 -2 -1.5 -1 -0.5 0 log r 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 lo g j6 j 7k = 0 7k = 0:1 7k = 0:18 7k = 0:5 slope= !1:8310 (c) -2.5 -2 -1.5 -1 -0.5 0 log r -2 -1.5 -1 -0.5 0 0.5 lo g j7 j 7k = 0 7k = 0:1 7k = 0:18 7k = 0:5 slope= !0:28651 (d) -2.5 -2 -1.5 -1 -0.5 0 log r -3 -2.5 -2 -1.5 -1 -0.5 0 lo g j8 j 7k = 0 7k = 0:1 7k = 0:18 7k = 0:5 slope= 1:2580 (e) -2.5 -2 -1.5 -1 -0.5 0 log r 3.4 3.5 3.6 3.7 3.8 3.9 4 lo g jp j 7k = 0 7k = 0:1 7k = 0:18 7k = 0:5 slope= !2=3 (f) Figura 4.33: Propriedades do escoamento ao se aproximarem da quina da contração - NSF, Re = 0.01, Wi = 1 e β = 1 2 . 4. Resultados Numéricos 80 Abaixo, na Figura 4.34, são ilustrados os vórtices com tamanho em redução com o au- mento do coeficiente de escorregamento. Para estes resultados, considerou-se a formulação NSF e os valores k̄ = 0, 0.1, 0.18, e 0.5. (a) (b) (c) (d) Figura 4.34: Redução da dimensão dos vórtices com o aumento de k̄ - Oldroyd-B, NSF, M1. 4. Resultados Numéricos 81 4.2.3 PTT Considerando agora o modelo constitutivo PTT, as simulações são executadas com os parâmetros adimensionais: • Comprimento característico: L = 1; • Número de Reynolds: Re = 0.01; • Número de Weissenberg: Wi = 1; • ε = 0.25; • Razão das viscosidades: β = 1 9 (CSF) e β = 1 2 (CSF e NSF); • Passo temporal: ∆tCSF = 10−4, ∆tNSF = 10−5; • Inflow : perfil de velocidade e de T desenvolvidos com u U = 0.375; • Coeficiente de escorregamento: k̄ = 0, 0.18, 0.27, 0.36, 0.45, 1, 4.5, 45. Para validar os resultados, o trabalho de Ferrás [16] é utilizado como referência na comparação do tamanho do vórtice que se forma no canto superior da contração, para alguns valores de k̄. Considerando a formulação cartesiana, os resultados são apresentados na tabela abaixo e indicam boa concordância com a referência. Tabela 4.5: Tamanhos de vórtice, XR - PTT, contração. k̄ M1 M2 Ferrás [16] 0.0 1.551008 1.525190 1.5180 0.18 1.465579 1.446187 1.4794 0.27 1.465579 1.446187 1.4779 0.36 1.465579 1.446187 1.4806 0.45 1.465579 1.446187 1.4848 4.5 1.551008 - 1.5516 45 1.551008 - 1.5667 O estudo seguinte é sobre o comportamento assintótico das propriedades do escoa- mento próximo da quina da contração e, para isso, são utilizadas ambas as malhas M1 e M2 e alguns valores de k̄. Desta forma, busca-se comparar os resultados numéricos obtidos com os resultados teóricos de [27] dados por (4.2). Esses resultados são apresentados nas Figuras 4.35 e 4.36 em que as retas são referentes ao caso no-slip. Vê-se que, para este caso, os resultados numéricos estão em conformidade com os teóricos a menos da pressão, que não acompanha a inclinação esperada. Logo após, são exibidos os resultados das simulações com β = 1 2 com o intuito de avaliar as diferenças entre os resultados produzidos por ambas as formulações para o tensor polimérico. Inicialmente é medido o comprimento dos vórtices para diferentes coeficientes k̄ considerando-se as formulações CSF e NSF e as malhas M1 e M2, cujos resultados são apresentados na Tabela 4.6 e na Figura 4.37. Esses resultados indicam que a formulação NSF tem vórtices menos duradouros, ou seja, os vórtices tendem a diminuir mais rapidamente com o aumento de k̄ do que os produzidos pela CSF. Nas Figuras 4.38 e 4.39 é exibida a variação assintótica das propriedades do escoamento onde observa-se que o aumento de k̄ tende a afastar as propr