UNIVERSIDADE FEDERAL FLUMINENSETCE - Escola de EngenhariaTEM - Departamento de Engenharia Mecânica

PROJETO DE GRADUAÇÃO II'

&

$

%

Título do Projeto:

DECOMPOSIÇÕES PARA O TENSOR DE REYNOLDSE SUA APLICAÇÃO EM UM ESCOAMENTO DE

BOCAL CONVERGENTE-DIVERGENTE

'

&

$

%

Autor:

LUCAS FIRMINO MELO PEREIRA

'

&

$

%

Orientador:

RONEY LEON THOMPSON, D.Sc.

Data: 28 de Março de 2016

LUCAS FIRMINO MELO PEREIRA

DECOMPOSIÇÕES PARA O TENSOR DE REYNOLDS E SUAAPLICAÇÃO EM UM ESCOAMENTO DE BOCAL

CONVERGENTE-DIVERGENTE

Trabalho de Conclusão de Curso apresentado aoCurso de Engenharia Mecânica da Universidade Federal Flu-minense, como requisito parcial para obtenção do grau deEngenheiro Mecânico.

Orientador:

RONEY LEON THOMPSON, D.Sc.

Niterói

28 de Março de 2016

Ficha Catalográfica elaborada pela Biblioteca da Escola de Engenharia e Instituto de Computação da UFF

P436 Pereira, Lucas Firmino Melo

Decomposições para o tensor de Reynolds e sua aplicação em um

escoamento de bocal convergente-divergente / Lucas Firmino Melo

Pereira. – Niterói, RJ: [s.n.], 2016.

82 f.

Trabalho (Conclusão de Curso) – Departamento de Engenharia

Mecânica, Universidade Federal Fluminense, 2016.

Orientador: Roney Leon Thompson.

1. Turbulência. 2. Tensor de Reynolds. I. Título.

CDD 532.0257

···-UNIVERSIDADE FEDERAL FLUMINENSE

TCE - Escola de Engenharia

TEM - Departamento de Engenharia Mecânica

PROJETO DE GRADUAÇAO II

AVALIAÇAO FINAL DO TRABALHO

Título do Trabalho: DECOMPOSIÇÕES PARA O TENSOR DE REYNOLDS E SUA

APLICAÇÃO EM UM ESCOAMENTO DE BOCAL CONVERGENTE-DIVERGENTE

Parecer do Professor Orientador da Disciplina:

- Grau Final recebido pelos Relatórios de Acompanhamento:

- Grau atrjbuíc.lo ao grupo nos Seminários de Progresso:

Parecer do Professor Orientador: ._ . l Ü kvW u e& rlWMlVt1" RA ui, "'t rr ' o[, f t 0

Nome e Assinatura do Professor Orientador:

Prof.: Roncy Lcon Thompson. Assinatura:

Parecer Conclusivo da Banca Examinadora do Trabalho:

[2g Projeto Aprovado Sem Restrições

D Projeto Aprovado Com Restrições

Prazo concedido para cumprimento cfas exigências:

Discriminação das exigências e/ou observações adicionais:

----UNIVERSIDADE FEDERAL FLUMINENSE TCE - Escola de Engenharia TEM - Departamento de Engenharia Mecânica

PROJETO DE GRADUAÇAO II

AVALIAÇAO FINAL DO TRABALHO (continuação)

Aluno: Lucas Firmino Melo Pereira.

Composição da Banca Examinadora:

Prof.: Roney Leon Thompson, D.Se.

Prof.: Daniel Rodríguez Álvarez, D.Se.

Prof.: Raul Bernardo Yidal Pessolani, D.Se.

Data ele Defesa do Trabalho:

Departamento ele Engenharia Mecânica, 28/03/2016

Grau: 9. S

DEDICATÓRIA

Dedico este trabalho a meus pais Luiz e Valéria e minha irmã Luiza, que sempre acre-

ditaram na educação como sendo o caminho para o crescimento pessoal e cujo suporte foi

imprescindível para a conclusão deste curso.

vi

AGRADECIMENTOS

Agradeço ao professor Roney Thompson, que demonstrou grande paciência e dedicação

ao me orientar por um campo do conhecimento que está longe de ser considerado simples.

Agradeço aos colegas de iniciação científica: Matheus, Bernardo, Perroni, Bragança e

Madeira e ao amigo Fernando, com os quais horas de conversas proveitosas certamente cor-

roboraram para o que está aqui apresentado.

vii

RESUMO

No estudo da mecânica dos fluidos, uma das formas de se abordar o problema da turbu-

lência é através de modelagens RANS. No trabalho a ser desenvolvido, objetiva-se construir

bases não lineares para o tensor de Reynolds e verificar sua eficácia em um caso específico

de bocal convergente-divergente.

Palavras-Chave: Turbulência, RANS, Tensor de Reynolds

viii

ABSTRACT

In the study of fluid mechanics, one of the ways to approach the turbulence problem is

through RANS modeling. The main objective of this work is to build nonlinear basis to the

Reynolds Stress and verify their accuracy on a converging-diverging channel flow.

Key-Words: Turbulence, RANS, Reynolds Stress

ix

SUMÁRIO

1. INTRODUÇÃO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.1 Uma visão estatística . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

1.2 Cascata de Energia e Hipótese de Kolmogorov . . . . . . . . . . . . . . . . . . 4

2. EQUAÇÕES DE GOVERNO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.1 Hipótese do Contínuo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.2 Equações De Balanço . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.3 Abordagem Média . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.3.1 Hipótese de Viscosidade Turbulenta . . . . . . . . . . . . . . . . . . 13

2.3.2 Modelos Tradicionais . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

3. FORMULAÇÃO TEÓRICA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

3.1 O Tensor Não-Persistência a Deformação . . . . . . . . . . . . . . . . . . . . . 18

3.2 Objetividade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

3.3 Decomposições Tensoriais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

3.3.1 Decomposição Em Fase - Fora de Fase . . . . . . . . . . . . . . . . . 24

3.3.2 Decomposição Proporcional - Ortogonal . . . . . . . . . . . . . . . 26

3.4 Modelos propostos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

3.4.1 Modelo 1: B =αD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

3.4.2 Modelo 2: B =α01+α1D+α2D2 . . . . . . . . . . . . . . . . . . . . . 28

3.4.3 Modelo 3: B =α01+α1D+α2D2 +βP . . . . . . . . . . . . . . . . . . 30

3.4.4 Modelo 4: B =αD+βP . . . . . . . . . . . . . . . . . . . . . . . . . . 31

3.4.5 Modelo 5: B =αD+β01+β1P+β2P2 . . . . . . . . . . . . . . . . . . 31

3.4.6 Modelo 6: B =α01+α1D+α2D2 +β01+β1P+β2P2 . . . . . . . . . 32

x

xi

3.5 Adimensionalizações Propostas . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

4. O CASO EM ESTUDO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

5. RESULTADOS E DISCUSSÃO . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

6. CONCLUSÕES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

REFERÊNCIAS BIBLIOGRÁFICAS . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

7. APÊNDICE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

7.1 Equação de transporte para o tensor de Reynolds . . . . . . . . . . . . . . . . . 66

7.1.1 Equação de transporte da energia cinética turbulenta . . . . . . . 69

LISTA DE FIGURAS

1.1 Transferência de energia entre as diversas escalas da turbulência. Fonte:

Pope (2000) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

3.1 Plano π: representação geométrica de como a vorticidade relativa atua em

um plano de autodireções de D. Fonte: Thompson and Mendes (2005) . . . . 20

3.2 Representação do conceito de objetividade: Dois observadores analisando

um ponto P. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

3.3 Simplificação do conceito por trás do parâmetro R. . . . . . . . . . . . . . . . 27

4.1 Representação tridimensional da geometria estudada com ênfase para as es-

truturas turbulentas. Fonte: Laval and Marquillie (2009) . . . . . . . . . . . . 34

4.2 Imagem da malha utilizada com resolução de 1:16. Fonte: Marquillie et al.

(2011) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

4.3 Representação do perfil de velocidade u na entrada da geometria . . . . . . . 36

4.4 Razão entre o maior comprimento de cada volume da malha e a escala de

kolmogorov correspondente. Fonte: Laval and Marquillie (2009) . . . . . . . 36

5.1 Componentes da parcela deviatórica do tensor de Reynolds B em cada seção. 38

5.2 Perfil de velocidade u em cada seção. . . . . . . . . . . . . . . . . . . . . . . . 39

5.3 Componentes do tensor taxa de deformação D em cada seção. . . . . . . . . . 40

5.4 Parâmetro RI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

5.5 coeficiente Cµ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

5.6 Distribuição do coeficiente Cµ ao longo das respectivas seções. . . . . . . . . . 42

5.7 coeficiente CD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

5.8 Distribuição do coeficiente CD ao longo das respectivas seções. . . . . . . . . 44

xii

xiii

5.9 Componentes do tensor (αD) em cada seção. . . . . . . . . . . . . . . . . . . . 45

5.10 Parâmetro RI I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

5.11 Coeficientes Cµ0 e CD0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

5.12 Distribuição dos coeficientes Cµ0 e CD0 ao longo das respectivas seções. . . . 47

5.13 Coeficiente Cµ1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

5.14 Coeficiente CD1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

5.15 Coeficiente CD2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

5.16 Distribuição do coeficiente CD2 ao longo das respectivas seções. . . . . . . . . 50

5.17 Distribuição do coeficiente Cµ2 ao longo das respectivas seções. . . . . . . . . 51

5.18 Componentes do tensor (α01+α1D+α2D2) em cada seção. . . . . . . . . . . 52

5.19 Parâmetro RI I I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

5.20 Distribuição dos coeficientes CP ao longo das respectivas seções. . . . . . . . 54

5.21 Componentes do tensor (α01+α1D+α2D2 +βP) em cada seção. . . . . . . . 55

5.22 Parâmetro RIV . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

5.23 Componentes do tensor (αD+βP) em cada seção. . . . . . . . . . . . . . . . . 57

5.24 Parâmetro RV . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

5.25 Coeficiente CP2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

5.26 Distribuição dos coeficientes CP2 ao longo das respectivas seções. . . . . . . . 59

5.27 Parâmetro RV I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

5.28 Verificação da colaboração das anisotropias através do parâmetro R. . . . . . 60

NOMENCLATURA

ξ Escalar arbitrário.

η Escala de kolmogorov.

τη Escala temporal referente à η.

uη Escala de velocidade referente à η.

Re Número de Reynolds.

T Taxa de transferência de energia entre escalas.

ε Taxa de dissipação de energia.

K n Número de Knundsen.

Γ Função de mapeamento da trajetória de uma partícula de fluido.

L Gradiente de velocidade.

D Tensor taxa de deformação.

W Tensor vorticidade.

T Tensor representativo do campo de tensões do fluido.

T Parcela deviatórica do tensor das tensões.

ρ Massa específica do fluido

µ Coeficiente de viscosidade.

ν Viscosidade cinemática.

p Pressão.

g aceleração gravitacional.

u velocidade.

νt Viscosidade turbulenta.

u∗ Escala de velocidade utilizada na descrição da viscosidade turbulenta.

xiv

xv

l∗ Escala de comprimento utilizada na descrição da viscosidade turbulenta.

lm Comprimento de mistura

cD Constante arbitrária

Cµ Constante de valor 0,09 do modelo k-ε.

Cε1 Constante de valor 1,44 do modelo k-ε.

Cε2 Constante de valor 1,92 do modelo k-ε.

σε Constante de valor 1,3 do modelo k-ε.

σk Constante de valor 1 do modelo k-ε.

β Constante de valor 340 do modelo k-ω.

β∗ Constante de valor 9100 do modelo k-ω.

a Constante de valor 59 do modelo k-ω.

σ Constante de valor 0,5 do modelo k-ω.

σ∗ Constante de valor 0,5 do modelo k-ω.

k Energia cinética turbulenta.

B Tensor de Reynolds deviatórico.

P Tensor não-persistência a deformação.

ΩD Tensor taxa de rotação dos autovetores de D.

W Tensor vorticidade relativa

R Tensor ortogonal de mudança de base.

λ Autovalores arbitrários.

eDi Autovetor de D.

Imn Escalar representativo da magnitude de um filamento material.

x Coordenada espacial de um ponto arbitrário.

x∗ Coordenada espacial de um ponto arbitrário vista por um referencial O∗.

A Tensor arbitrário.

F Tensor arbitrário.

xvi

ϕF+A Parcela de A que comuta e não é ortogonal à F.

ϕF+−A Parcela de A que comuta e é ortogonal à F.

ϕF−A Parcela de A que não comuta e é ortogonal à F.

PFA Parcela proporcional da projeção de A em F .

PFA Parcela não proporcional da projeção de A em F .

ΦFA Parcela Em Fase da projeção de A em F .

ΦFA Parcela Em Fase da projeção de A em F .

R Parâmetro de performance do modelo sobre o tensor alvo.

1DD Tensor de quarta ordem de mudança de base.

M Tensor utilizado para verificar a colaboração das anisotropias do tensor de Reynolds.

1 INTRODUÇÃO

Em diversas situações de interesse na mecânica dos fluidos é possível se deparar com

escoamentos que apresentam as variáveis de interesse distribuídas de forma aparentemente

caótica e randômica, este padrão é característico de um escoamento turbulento. Tal fenô-

meno tem alto caráter difusivo, ou seja, há grande capacidade de se transportar propriedades

ao longo de direções que não são as principais do escoamento, notando-se o claro “efeito

de mistura”. Isto pode ser positivo, como na mistura de um combustível na câmara de com-

bustão de um motor. Trata-se de um fenômeno intrinsecamente transiente, ou seja, mesmo

que de forma mediana as propriedades convirjam para um patamar independente do tempo,

o valor integral das mesmas varia a cada instante, além disso o fenômeno é tridimensional,

rotacional e altamente dissipativo, isto é, precisa ser constantemente alimentado de modo a

evitar sua degeneração.

Tratar o escoamento como randômico pode soar de forma estranha a um primeiro mo-

mento, visto que acredita-se que o fenômeno seja determinístico e que possa ser descrito por

equações de governo. De fato estas equações descrevem o fenômeno, porém sua não lineari-

dade e sua extrema susceptibilidade a variações nas condições de contorno ou iniciais torna

a reprodução exata destas um ato extremamente improvável. Pequenos desvios nas condi-

ções de contorno ou de valor inicial ou mesmo erros numéricos no caso de uma abordagem

computacional são alimentados e podem provocar mudanças substanciais nos campos das

variáveis de interesse. Essa incapacidade de controlar as variáveis à tal nível de refino força

a se tomar, a primeiro momento, uma abordagem estatística.

O foco deste trabalho está nas modelagens RANS - Reynolds Averaged Navier-Stokes,

aonde assume-se uma média temporal sobre as equações de governo. Tal abordagem cria

um tensor chamado tensor de Reynolds, responsável por comportar a colaboração das flu-

tuações da velocidade para o campo médio de tensões do escoamento. Este tensor agrega

1

2

ao conjunto de equações já existentes um total de 6 novos termos e cria o famoso problema

de fechamento da turbulência. Como será visto mais adiante, uma forma muito usual de

modelá-lo, conhecida como hipótese de viscosidade turbulenta, já é matematicamente precá-

ria, uma vez que o tensor taxa de deformação é incapaz, na maioria das vezes, de representar

toda a complexidade do tensor de Reynolds.

Thompson et al. (2010) propõem algumas decomposições para este tensor além de uma

forma de avaliação da performance para cada modelo proposto. Neste trabalho serão vistas

bases que são função do tensor taxa de deformação e do tensor não-persistência à deforma-

ção, o motivo de sua escolha será visto mais adiante. Então, o presente trabalho permite,

além da verificação do desempenho, a avaliação de algumas das particularidades de cada

modelo.

Para a análise de tais bases, foi escolhido um escoamento estatísticamente bidimensional

sobre uma superfície levemente côncava. A simulação foi executada via DNS - Simulação

Numérica Direta - por Marquillie et al. (2011) em um número de Reynolds de Re=12600,

e os resultados usados para compor tais bases. O procedimento de determinação dos coe-

ficientes de cada base proposta necessita de resultados mais nobres do que o RANS pode

proporcionar, servindo como um “gabarito”. Haveria também a opção de fazer uso de uma

base LES - Large Eddy Simulation, porém como o LES faz uso de modelagens a partir de

uma certa escala, optou-se pela utilização de um DNS, que não faz uso de modelos e conse-

quentemente teria o conjunto de dados mais nobre frente a todas estas opções.

No capítulo seguinte, será feita uma revisão das equações de governo, da hipótese de

Boussinesq e de alguns dos modelos usualmente adotados para abordá-la. No capítulo 3,

irá se apresentar um importante tensor para os modelos que virão a ser propostos, além

da fundamentação matemática por trás das decomposições utilizadas e finalmente pode-se

mostrar como os coeficientes de cada modelo são obtidos. No quarto capítulo apresenta-se o

caso em estudo e no quinto capítulo serão apresentados os resultados.

3

1.1 Uma visão estatística

Como mencionado, mesmo as menores perturbações podem desviar os valores de dois

experimentos que, aparentemente sob mesmas condições de contorno e iniciais, deveriam

fornecer o mesmo resultado. Seria interessante a execução de um número suficientemente

grande destes experimentos e então a média da propriedade de interesse no mesmo ponto

e no mesmo tempo nos diversos experimentos entregaria um campo estatisticamente mais

confiável.

Performar a média de um grande número de experimentos independentes pode ser im-

praticável e uma solução possível é a de pensar em médias que ocorrem dentro de um único

experimento, como a média temporal e a média espacial dos campos.

Após observar-se o experimento durante um tempo suficientemente grande e consequen-

temente estando este em regime permanente1, a hipótese de ergodicidade afirma que a média

temporal tende à média amostral de diversos experimentos diferentes.

Seja uma variável aleatória contínua ξ= ξ(x, t ). A flutuação de seus valores no tempo em

torno de sua média pode ser escrita como:

ξ′(x, t ) = ξ(x, t )−ξ(x) (1.1)

Onde sua média temporal2 é definida como:

ξ(x) ≡ limt→∞

1

T

∫ T

0ξ(x, t )d t (1.2)

Combinando a equação (1.1) e a definição (1.2), é fácil ver que:

1 No regime permanente, as variáveis médias são invariantes quanto à mudança temporal.2 Neste texto, a média temporal será identificada através de uma barra sobrescrita à variável.

4

A média temporal de uma flutuação é nula.

ξ(x) = limt→∞

1

T

∫ T

0(ξ(x)+ξ′(x, t ))d t

ξ(x) = ξ(x) limt→∞

1

T

∫ T

0d t + lim

t→∞1

T

∫ T

0ξ′(x, t )d t︸ ︷︷ ︸

ξ′(x,t )

ξ′(x, t ) = 0 (1.3)

A média temporal da derivada espacial de uma propriedade é a derivada espacial da

média temporal desta propriedade.

limt→∞

1

T

∫ T

0

∂ξ(x, t )

∂xd t = ∂

∂xlim

t→∞1

T

∫ T

0ξ(x, t )d t

∂ξ(x, t )

∂x= ∂ξ(x, t )

∂x(1.4)

A média temporal entre o produto de uma flutuação da propriedade com seu valor médio

é zero.

ξ′(x, t ).ξ(x) = limt→∞

1

T

∫ T

0ξ′(x, t ).ξ(x)d t

= ξ(x) limt→∞

1

T

∫ T

0ξ′(x, t )d t = 0 (1.5)

1.2 Cascata de Energia e Hipótese de Kolmogorov

Fisicamente, o escoamento turbulento pode ser como visto uma série de estruturas3 com-

plexas que estão a cada instante consumindo energia e se degenerando. A ideia de que a

energia entra primeiro nas grandes escalas e então transfere-se das maiores para as menores

escalas foi originalmente proposta por Richardson na década de 1920, construindo assim a

ideia de uma cascata de energia.

Em 1941 Kolmogorov avançou neste conceito e postulou três hipóteses:

Em sua primeira hipótese, Kolmogorov afirma que para um número de Reynolds sufi-

3 Cada estrutura do escoamento tem sua própria escala espacial, temporal e consequentemente de velocidade,denotadas neste texto respectivamente por l ,τl , ul .

5

cientemente grande, as estatísticas das menores escalas do escoamento são isotrópicas, ou

seja, são invariantes quanto a mudança direcional. Isso quer dizer que mesmo que haja uma

tendência da turbulência por se apresentar mais expressivamente em uma direção, a transfe-

rência de energia entre as diversas escalas torna o caráter direcional indistiguível.

Em sua segunda hipótese, Kolmogorov afirma que, ainda em números de Reynolds sufi-

cientemente grandes, as menores escalas são uma função apenas da taxa de transferência de

energia ε para estas e da dissipação viscosa ν , então através de uma verificação dimensional

as famosas escalas de comprimento, temporal e de velocidade de Kolmogorov podem ser

definidas como:

η≡(ν3

ε

) 14 (1.6)

τη ≡(νε

) 12 (1.7)

uη ≡ η.τη =(ν.ε

) 14 (1.8)

Uma forma interessante de se enxergar as equações acima é verificando-se o número

de Reynolds, que pode ser entendido como uma razão de forças inerciais e forças viscosas,

assume o valor 1 para tal escala.

Reη =ρuηη

µ=

(ν.ε

) 14

.(ν3

ε

) 14

ν= 1 (1.9)

Seja l0, τl0 e U0 os parâmetros característicos da maior escala comportada pela geometria

do problema, e assumindo que no regime permanente a mesma quantidade de energia que

entra nas grandes escalas precisa ser dissipada à mesma taxa - T - para escalas menores. A

dissipação pode ser escrita como:

T = ε= U 3o

l0(1.10)

Assim, as microescalas podem ser escritas em função das maiores escalas, demonstrando

6

claramente que o aumento de Re provoca uma redução nas menores escalas.

η

l0=

(1

Re

) 34

(1.11)

τη

τ0=

(1

Re

) 12

(1.12)

uηu0

=( 1

Re

) 14 (1.13)

Por fim, em sua terceira hipótese, Kolmogorov afirma que em Re suficientemente altos,

há uma faixa de escalas suficientemente distantes das maiores escalas - por exemplo: l0 - e

suficientemente distantes das menores escalas -por exemplo: η - nas quais as estatísticas são

função apenas da taxa de dissipação ε e não mais da viscosidade ν. Estas hipóteses podem

ser sintetizadas em um diagrama simples:

Fig. 1.1: Transferência de energia entre as diversas escalas da turbulência. Fonte: Pope(2000)

Como apresentado por Pope (2000), pode-se ver que a escala de comprimento lE I -

Energy Range / Inertial Range - marca a divisão entre a produção de energia turbulenta,

que se dá em escalas superiores à esta e próximas à l0 e uma subregião inercial aonde a

energia transfere-se contínuamente até lD I - Dissipation Range / Inertial Range - que marca

o começo da ação da dissipação viscosa e a consequente perda de energia para a energia

interna.

2 EQUAÇÕES DE GOVERNO

2.1 Hipótese do Contínuo

Antes de prosseguir, é indispensável falar da hipótese de se abordar um meio como sendo

um contínuo de matéria. Esta abordagem traz algumas vantagens: pode-se deixar de refe-

renciar partículas individuais do escoamento e ganha-se a noção de volumes infinitesimais

por exemplo, trazendo também simplificações quanto às iterações destas partículas com sua

vizinhança. Para admitirmos tal hipótese é necessário que a menor escala do nosso esco-

amento seja suficientemente grande quando comparada com a escala molecular do fluido,

assim como a escala temporal. Viu-se que a menor escala é dada por η. e sendo λ o livre

caminho médio entre as moléculas que compõem o fluido em questão, a hipótese será tão

boa quanto menor for a razão:

K n = λ

η¿ 1 (2.1)

Pensando a matéria como um contínuo, não se fará referência à uma partícula específica

que é transportada ao longo do espaço pois este conceito não mais se aplica. Irá se pensar

em uma partícula de fluido que é referenciada por uma coordenada fixa X = (X1,X2,X3) em

um dado tempo t0 , tal coordenada é conhecida como coordenada material. Seja Γ(X, t ) a

função que realiza este mapeamento.

Γ(X, t0) = X

Γ(X, t ) = x

∂Γ(X, t )

∂t= v (2.2)

sendo x = (x1,x2,x3) uma coordenada de escolha arbitrária. Assim qualquer propriedade

7

8

do fluido pode ser descrita nos dois referenciais.

Λ(x, t ) =Λ(Γ(X, t ), t ) (2.3)

Logo, a derivada temporal no referencial material pode ser extendida para o referencial

espacial.

∂Λ(x, t )

∂t

∣∣∣X=X0

= ∂

∂t

(Λ(Γ(X, t ), t )

)∣∣∣X=X0

= ∂Γi (X, t ), t )

∂t

∂

∂xi

(Λ(Γ(X, t ), t )

)∣∣∣X=X0

+ ∂Λ(Γ(X, t ), t )

∂t

∣∣∣X=X0

= v(x, t ).∇Λ(x, t )+ ∂Λ(x, t )

∂t(2.4)

Este raciocínio pode ser extendido para tensores de ordem mais elevada, basta ser feito

componente à componente.

Ao ser advectado pelo escoamento, um filamento material dr = d x1,d x2,d x3 tem sua

evolução temporal dada da seguinte forma:

dr = ∂xk

∂xidxi ek (2.5)

dr = L.dr (2.6)

Aonde L é o grandiente de velocidade. Usualmente o gradiente de velocidade é decom-

posto de forma aditiva em uma parcela simétrica - Tensor taxa de deformação (D) e uma

antissimétrica - Tensor vorticidade (W).

L = D+W

D = 1

2

(L+LT

)(2.7)

W = 1

2

(L−LT

)

Um corpo material está sujeito a dois tipos de esforços: Superficiais e de campo. A soma

da resultante dos esforços de campo aplicados à um volume material Ω com a resultante

9

dos esforços superficiais aplicados à superfície ∂Ω deste volume material é igual a taxa de

variação do momento linear das partículas deste corpo material.

∂

∂t

ÑΩ

ρ(x, t )v(x, t )dV =ÑΩ

ρ(x, t )b(x,t)dV +Ï∂Ω

tdS (2.8)

Onde este vetor de esforços superficiais t pode ser escrito como o produto de um tensor de

segunda ordem que representa o campo de tensões do fluido e um versor normal à superfície

de interesse. t = T.n , assim aplica-se o teorema da divergência na integral de superfície e

pode-se reduzir a equação à uma forma diferencial. Conhecida como equação do movimento

de Cauchy - em notação diádica e indicial, respectivamente, como será utilizada adiante.

ρDv

Dt=∇.T+ρb (2.9)

ρai =∂Ti j

∂x j+ρbi (2.10)

2.2 Equações De Balanço

Ao se estudar o escoamento, trabalha-se com um conjunto de equações de balanço e de

continuidade que ao serem resolvidas permitem a obtenção das propriedades desejadas.

A hipótese de incompressibilidade que será adotada mais a frente afirma que a derivada

material da massa específica do fluido é nula. Iniciando-se com a conservação de massa,

como apresentado por Lai et al. (2009), é fácil mostrar que esta hipótese se reduz à:

Dm

Dt= D(ρdV )

Dt= dV

Dρ

Dt+ρDdV

DtDρ

Dt︸︷︷︸=0

+ρ∂vk

∂xk= 0

∂vk

∂xk= 0 (2.11)

Daqui em diante, a massa específica e será assumida como constante: ρ(x, t ) = ρ, assim

como faz Pope (2000). O tensor de tensões descrito na equação (2.9) é usualmente decom-

10

posto em uma parcela esférica que virá a ser conhecida como a pressão e outra deviatórica.

Ti j =−pδi j + Ti j (2.12)

Na concepção de um fluido newtoniano, a parcela deviatórica das tensões é linearmente

proporcional à taxa de deformação.

Ti j =µ(∂ui

∂x j+ ∂u j

∂xi− 1

3

∂uk

∂xkδi j

)(2.13)

vale ressaltar que a componente aqui referida como pressão não tem o mesmo sentido físico

que a pressão termodinâmica e deve ser entendida apenas como a componente isotrópica do

tensor T. Caso a velocidade do escoamento fosse à zero, esta seria a única componente que

restaria deste tensor.

Aplicando-se as considerações acima em (2.9), chega-se a:

ρ(∂ui

∂t+u j

∂ui

∂x j

)= ∂

∂x j

(−pδi j +µ

(∂ui

∂x j+ ∂u j

∂xi− 1

3

∂vk

∂xkδi j

))+ρgi (2.14)

admitindo a hipótese de incompressibilidade (2.11), e para efeitos de simplificação, rees-

crevendo o termo referente ao esforço de campo como o gradiente de seu potencial: ρgi =−ρ ∂

∂xig z, aonde g é a aceleração gravitacional e z a coordenada vertical na base cartesiana

pode-se assim incorporá-lo ao termo de pressão.

ρ(∂ui

∂t+u j

∂ui

∂x j

)= ∂

∂x j

(−pδi j +µ

(∂ui

∂x j+ ∂u j

∂xi

))(2.15)

Resolvendo-se o divergente, chega-se à forma mais conhecida deste balanço que são as

equações de Navier-Stokes:

ρ(∂ui

∂t+u j

∂ui

∂x j

)=− ∂

∂xip +µ

( ∂2ui

∂x j∂x j

)(2.16)

11

Tomando-se o divergente em (2.16):

ρ∂

∂xi

(∂ui

∂t+u j

∂ui

∂x j

)=− ∂

∂xi

∂

∂xip +µ ∂

∂xi

( ∂2ui

∂x j∂x j

)(2.17)

ρ( ∂∂t

∂ui

∂xi+ ∂

∂xi

(u j∂ui

∂x j

))=− ∂2

∂xi∂xip +µ

( ∂

∂x j

∂

∂x j

∂ui

∂xi

)(2.18)

ρ( ∂∂t

∂ui

∂xi+ ∂ui

∂x j

∂u j

∂xi+u j

∂

∂x j

∂ui

∂xi

)=− ∂2

∂xi∂xip +µ

( ∂

∂x j

∂

∂x j

∂ui

∂xi

)(2.19)

Por fim, adotando-se a hipótese de incompressibilidade, a equação 2.17 toma a forma de

uma equação de Poisson.

∂2

∂xi∂xip =−ρ∂ui

∂x j

∂u j

∂xi(2.20)

Esta equação, junto com as equações que surgirão de Navier-Stokes para cada direção, criam

um sistema de quatro equações e quatro incógnitas, que pode ser resolvido de forma aco-

plada. Uma das formas usuais de se enxergar a pressão é como o operador que corrige a

continuidade durante as iterações entre as equações acima.

2.3 Abordagem Média

O conceito de média temporal definido em (1.2) pode então ser aplicado na equação de

continuidade (2.11)

∂(ui +u′i )

∂xi= ∂(ui +u′

i )

∂xi= ∂ui

∂xi= 0 (2.21)

mostrando que a velocidade média continua a preservar esta característica, e consequente-

mente a flutuação de velocidade também:

∂(ui +u′i )

∂xi= ∂ui

∂xi+ ∂u′

i

∂xi= ∂u′

i

∂xi= 0 (2.22)

Aplicando o operador média temporal nas equações de balanço de quantidade de movi-

12

mento.

ρ(∂ui

∂t+u j

∂ui

∂x j

)= ∂

∂x j

(−pδi j +µ

(∂ui

∂x j+ ∂u j

∂xi

))(2.23)

ρ(∂(ui +u′

i )

∂t+u j

∂ui

∂x j

)= ∂

∂x j

(− (p +p ′)δi j +µ

(∂(ui +u′i )

∂x j+∂(u j +u′

j )

∂xi

))(2.24)

O termo de aceleração convectiva toma a seguinte forma:

u j∂ui

∂x j= ∂ui u j

∂x j= ∂ui u j

∂x j(2.25)

onde,

ui u j = (ui +u′i )(u j +u′

j ) = ui .u j +u′i u′

j (2.26)

logo, tem-se que,

u j∂ui

∂x j= u j

∂ui

∂x j+ui

∂u j

∂x j+∂u′

i u′j

∂x j= u j

∂ui

∂x j+∂u′

i u′j

∂x j(2.27)

Assim, isola-se o termo que contém divergente no lado esquerdo da equação e obtém-se

a conhecida equação de Reynolds:

(∂ui

∂t+u j

∂ui

∂x j

)= ∂

∂x j

(− p

ρδi j +ν

(∂ui

∂x j+ ∂u j

∂xi

)−u′

i u′j

)(2.28)

o último termo da equação (2.28) é conhecido como tensor de Reynolds, um tensor de se-

gunda ordem simétrico que pode aparecer na literatura de diversas formas, neste texto será

definido apenas como −u′i u′

j . Este é um tensor que acrescentará às equações 6 novas incóg-

nitas. Ao longo do tempo, diversos modelos surgiram com hipóteses de fechamento, aonde

propunham equações ou um conjunto delas que pudesse fechar o sistema.

13

2.3.1 Hipótese de Viscosidade Turbulenta

Na tentativa de se obter uma modelagem adequada para o tensor de Reynolds, surge a

hipótese de viscosidade turbulenta. Aonde modela-se o tensor de Reynolds em analogia à

parcela deviatórica do tensor das tensões para fluidos Newtonianos (2.13), porém diferente-

mente da hipótese das tensões viscosas que são fruto de iterações moleculares e que ocorrem

em escalas muito menores que a da menor estrutura do escoamento as flutuações de velo-

cidade são uma função da própria escala temporal e de comprimento da estrutura à qual

pertence o ponto analisado.

−u′i u′

j = νt

(1

2

(∂ui

∂x j+ ∂u j

∂xi

))− 2

3kδi j (2.29)

k é definido como a energia cinética turbulenta.

k = 1

2u′

k u′k (2.30)

O fato interessante a se observar é que esta proposta tenta tornar proporcional a parte

deviatórica do tensor de Reynolds ao tensor taxa de deformação. Tal abordagem é contestável

sob o ponto de vista matemático pois, como será visto nas propostas de novos modelos, o

tensor de Reynolds deveria ser perfeitamente proporcional com o tensor taxa de deformação

para que esta hipótese fosse perfeitamente eficiente.

Guardados os erros intrínsecos a esta modelagem, o problema passa então ser como

expressar o campo de νt da forma mais eficiente possível. Como dito anteriormente, o

tensor é uma função das próprias escalas em que se encontra e logo, como definido em Pope

(2000), tenta-se descrever este campo na forma:

νt = u∗.l∗ (2.31)

onde νt é proporcional a uma escala de velocidade e a uma escala de comprimento, que

por sua vez podem ser representados por outros campos escalares, nos chamados modelos

algébricos, ou modelados por equações de transporte, dando origem aos modelos de uma,

14

duas ou mais equações.

2.3.2 Modelos Tradicionais

A equação (2.29) pode ser escrita de forma simplificada como:

B = 2νt D (2.32)

aonde B será convencionado como o tensor de Reynolds deviatórico.

B =(−u′

i u′j +

1

3u′

k u′kδi j

)ei e j (2.33)

Dentro dos modelos algébricos, o modelo de comprimento de mistura estipula uma escala

de comprimento lm

l∗ = lm(x) (2.34)

a escala de velocidade pode então ser escrita como uma função deste comprimento e para

escoamentos entre placas paralelas e estatisticamente bidimensionais, como é o caso de um

vasto número de experimentos encontrados, esta velocidade pode ser descrita simplesmente

por,

u∗ = lm

∣∣∣∂u

∂y

∣∣∣ (2.35)

e a viscosidade turbulenta pode então ser escrita como:

νt = l 2m

∣∣∣∂u

∂y

∣∣∣ (2.36)

a escala de velocidade ainda pode ser escrita de forma genérica através do gradiente de

15

velocidade, representado, por exemplo, pela intensidade da taxa de deformação:

u∗ = lm

p2D : D (2.37)

νt = l 2m

p2D : D (2.38)

esta abordagem que aparentemente descreve a viscosidade turbulenta novamente transfere a

problemática de modelagem para o novo parâmetro lm , que de alguma forma ainda precisa

ser especificado. Como mencionado por Wilcox (1988), este parâmetro precisa ser estimado

empiricamente para cada caso e assim o modelo pode ser chamado de incompleto.

Uma outra forma de se propor a descrição desta escala de velocidade, como apresentada

por Pope (2000), foi dada por Kolmogorov e Prandtl, em 1942 e 1945 respectivamente. A

proposta consiste basicamente em mensurar esta escala como função da energia cinética

turbulenta - k.

u∗ = cp

k (2.39)

e neste caso a viscosidade turbulenta se apresenta da seguinte forma,

νt = lmcp

k (2.40)

este novo modelo baseado na energia cinética turbulenta não isenta a modelagem do com-

primento de mistura, além disso, exige também a especificação de k, que é usualmente feita

por sua equação de transporte:

∂k

∂t+u j

∂k

∂x j= ∂

∂x j

((ν+ νt

σk

) ∂k

∂x j

)+νt

(∂ui

∂x j+ ∂u j

∂xi

)∂ui

∂x j−ν∂u′

i

∂x j

∂u′i

∂x j(2.41)

o que faz este modelo pertencer a classe dos modelos de uma equação, uma vez que além do

conjunto de equações de governo já mencionadas, mais uma equação de transporte precisa

ser resolvida.

Um modelo proposto por Jones and Launder (1972) é capaz de descrever esta escala de

16

comprimento (2.34) com o auxílio de uma segunda propriedade, a taxa de dissipação - ε,

l∗ = cDk

32

ε(2.42)

e cD é uma constante. Novamente, uma equação de transporte pode ser resolvida para ε1.

∂ε

∂t+u j

∂ε

∂x j= ∂

∂x j

((ν+ νt

σε

) ∂ε∂x j

)+Cε1

ε

kνt

(∂ui

∂x j+ ∂u j

∂xi

)∂ui

∂x j−Cε2

ε2

k(2.43)

Esta modelagem constitui o famoso modelo k−ε. A necessidade de solução de duas equações

de transporte o classifica como um modelo a duas equações, fazendo com que este modelo

possa ser chamado de completo, uma vez que este par de equações de transporte é capaz de

eliminar a necessidade de se especificar previamente um comprimento.

Deste modo, pode-se reescrever a viscosidade turbulenta da seguinte maneira, aonde

tem-se o novo adimensional Cµ,

νt =Cµk2

ε(2.44)

As constantes do modelo, segundo Launder and Sharma (1974), são: Cµ = 0,09, Cε1 = 1,44,

Cε2 = 1,92, σε = 1,3, σk = 1,0.

Um segundo modelo a duas equações, apresentado por Wilcox (1988), faz uso de uma

propriedade conhecida como dissipação específica ω.

ω≡ ε

k(2.45)

aonde a viscosidade turbulenta pode ser reescrita como:

νt = k

ω(2.46)

com estas definições, pode-se tanto reescrever a equação (2.41) em função de ω, quanto

1 vale o comentário de que esta equação é obtida de forma empírica

17

utilizar a equação de transporte de ω no lugar de ε.

∂k

∂t+u j

∂k

∂x j= k

ω

(∂ui

∂x j+ ∂u j

∂xi

)∂ui

∂x j−β∗ωk + ∂

∂x j

((ν+σ∗ k

ω)

)∂k

∂x j(2.47)

∂ω

∂t+u j

∂ω

∂x j= a

(∂ui

∂x j+ ∂u j

∂xi

)∂ui

∂x j−βω2 + ∂

∂x j

((ν+σ k

ω)

)∂ω

∂x j(2.48)

este é o par de equações do modelo k −ω e suas constantes são: β = 340 , β∗ = 9

100 , a = 59 ,

σ= 12 e σ∗ = 1

2 .

3 FORMULAÇÃO TEÓRICA

3.1 O Tensor Não-Persistência a Deformação

Um importante tensor para os modelos que virão a ser propostos é o tensor não persistên-

cia à deformação. Essencialmente este tensor foi proposto na problemática da classificação

de escoamentos complexos.

Seja Imn um filamento material inicialmente alinhado com uma das direções principais

do tensor taxa de deformaçao D.

Imn = em .Den (3.1)

aonde em e en são direções associadas com as direções principais. E portanto sua varia-

ção temporal, vista por um observador que está na base dos autovetores de D pode ser escrita

como:

Imn = ˙em .Den + em .Den + em .D˙en

Imn = W.em .Den + em .Den + em .D.Wen (3.2)

Sendo:

DD =3∑

i=1λi eD

i eDi

D =3∑

i=1λi eD

i eDi +

3∑i=1

λi ˙eDi eD

i +3∑

i=1λi eD

i˙e

Di (3.3)

E definindo um tensor ΩD como o tensor capaz de descrever a taxa de variação de cada

18

19

autovetor de D.

˙eDi =ΩD eD

i

ΩD = ˙eDi eD

i (3.4)

E sendo ΩD antissimétrico, a equação (3.3), toma a seguinte forma:

D =3∑

i=1λi eD

i eDi +

3∑i=1

λiΩD eD

i eDi +

3∑i=1

λi eDi Ω

D eDi

D = D′+ΩD .D−D.ΩD (3.5)

sendo D′ a taxa de variação dos autovalores de D e (ΩD .D−D.ΩD) a taxa de rotação dos

autovetores.

Aplicando-se (3.5) em (3.2):

Imn = W.em .Den + em .(D′+ΩD .D−D.ΩD )en + em .D.Wen

Imn = em .(D′−W.D+ΩD .D−D.ΩD +D.W

)en (3.6)

Por fim, define-se o tensor vorticidade relativa W1, como sendo a diferença entre a taxa

de rotação promovida pelo escoamento livre e a taxa de rotação da base de D.

W = W−ΩD (3.7)

Chega-se à:

Imn = em .(D′+D.W−W.D

)en (3.8)

Aonde(D.W−W.D

)é chamado de tensor não-persistência à deformação P.

1 A barra superior em W não deve ser confundida com a notação adotada para a média temporal.

20

PD ≡(D.W−W.D

)D =

0 ω3(λ2 −λ1) ω2(λ1 −λ3)

ω3(λ2 −λ1) 0 ω1(λ3 −λ2)

ω2(λ1 −λ3) ω1(λ3 −λ2) 0

(3.9)

o símbolo referente ao tensor D sobrescrito ao tensor P é adotado como a notação para dizer

que o tensor P está na base de D.

Como exemplo, esteja Imn inicialmente alinhado com a direção principal de maior au-

tovalor de D, então o tensor P fornece uma medida de como este filamento tende a sair da

direção de maior esticamento. Uma representação geométrica que auxilia no entendimento

deste tensor, introduzida por Thompson and Mendes (2005), consiste no plano π. Um plano

que comporta duas das direções principais de D e é perpendicular à vorticidade relativa cor-

respondente.

Fig. 3.1: Plano π: representação geométrica de como a vorticidade relativa atua em um planode autodireções de D. Fonte: Thompson and Mendes (2005)

21

3.2 Objetividade

Um ponto importante para o desenvolvimento da teoria é o conceito de objetividade - ou

invariância euclidiana. De modo resumido, a ideia por trás da objetividade é a de que ao

se observar uma grandeza através de dois referenciais distintos, não sendo ambos necessa-

riamente inerciais, para que esta grandeza seja considerada objetiva ambos os observadores

precisam estar em consenso quanto ao evento físico observado.

Matematicamente a ideia acima se traduz da seguinte forma: Estando um observador em

um referencial ortogonal inercial O e seja O∗ um referencial que pode rotacionar e transladar

livremente no espaço. Seja x o vetor posição de uma partícula no referencial O e x∗ o

vetor posição no referencial O∗. Um mesmo ponto P no espaço pode ser descrito nos dois

referenciais através da seguinte transformação:

x∗ = R.x+c (3.10)

Na equação acima R é um tensor de mudança de base, capaz de efetuar, por exemplo,

uma rotação ou um espelhamento e c a distância entre as origens dos referenciais, como

ilustrado na figura (3.2):

Fig. 3.2: Representação do conceito de objetividade: Dois observadores analisando umponto P.

22

Para que uma grandeza seja considerada objetiva, ela deve satisfazer as seguintes rela-

ções:

α∗ =α (3.11)

f∗ = R.f (3.12)

A∗ = R.A.RT (3.13)

As equações (3.11), (3.12) e (3.13) retratam como um escalar, um vetor e um tensor

objetivo, respectivamente, se apresentam após a transformação.

Sendo assim, vale a análise sobre algumas grandezas:

A distância entre dois pontos é objetiva:

x∗A = R.xA +c

x∗B = R.xB +c

x∗A −x∗

B = R.(xA −xB ) (3.14)

A velocidade não é uma grandeza objetiva:

dx∗

d t= d

d t(R.x+c)

v∗ = x∗ = Rx+Rx+ c (3.15)

Assim como a aceleração:

a∗ = x∗ = Rx+2Rx+Rx (3.16)

O gradiente de velocidade é não objetivo:

L∗ = ∂v∗

∂x∗ = ∂v∗

∂x

∂x

∂x∗ = (R+RL)RT = RRT +RLRT (3.17)

23

fazendo uso do fato de que:

RRT = 1

d

dt

(RRT

)= 0 = RRT +RRT

RRT =−(RRT

)T(3.18)

sua decomposição em uma parte simétrica e outra antisimétrica toma a seguinte forma:

D∗ = 1

2

(RRT +RLRT +RRT +RLT RT

)= R

(1

2(L+LT )

)RT (3.19)

W∗ = 1

2

(RRT +RLRT −RRT −RLT RT

)= RRT +R

(1

2(L−LT )

)RT (3.20)

Mostra que D é objetivo e W não. Se D é objetivo, consequentemente seus autovetores

também são.

eD∗i = ReD

i (3.21)

˙eD∗i = ReD

i +R˙eDi (3.22)

˙eD∗i eD∗

i =ΩD∗ = RRT +RΩD RT (3.23)

E assim, a vorticidade relativa é objetiva.

W∗ = W∗−ΩD∗ = RRT +RWRT −

(RRT +RΩD RT

)= R(W∗−ΩD∗

)RT (3.24)

Conclui-se então que o tensor P é objetivo.

3.3 Decomposições Tensoriais

Sejam A e F dois tensores de segunda ordem e simétricos, é possível projetar A na base

dos autovetores de F . Tal projeção cria três partes:

A =ϕF+A +ϕF+−

A +ϕF−A (3.25)

24

ϕF+A é a parcela da projeção de A que comuta com F , isto quer dizer que partilham as

mesmas direções principais e podem ser chamados também de coaxiais, além disso, esta

parcela não é ortogonal à F .

ϕF+A .F =F .ϕF+

A (3.26)

ϕF−A é a parcela desta projeção que é ortogonal à F e não compartilha as mesmas dire-

ções, ou seja, seu produto interno2 é nulo:

<ϕF−A ,F >= 0 (3.27)

E por fim, há a possibilidade de haver um termo que é ao mesmo tempo coaxial e orto-

gonal, chamado de ϕF+−A . Assim:

<ϕF+−A ,F >= 0 (3.28)

ϕF+−A .F =F .ϕ

F+−A (3.29)

Estes três termos não criam uma única possibilidade de decomposição aditiva do tensor

A em F

3.3.1 Decomposição Em Fase - Fora de Fase

Seja a projeção de A na base F - em notação matricial :

AF =

AF

11 AF12 AF

13

AF21 AF

22 AF23

AF31 AF

32 AF33

(3.30)

2 Define-se como o produto interno de dois tensores A e F : < A,F > = tr (A.FT ) - Aonde tr(.) é o operadortraço, que pode ter sua notação abreviada para A :F

25

A decomposição dita Em Fase - Fora de Fase, assume a seguinte forma em notação matricial:

AF =

AF

11 0 0

0 AF22 0

0 0 AF33

︸ ︷︷ ︸

ΦFA

+

0 AF

12 AF13

AF21 0 AF

23

AF31 AF

32 0

︸ ︷︷ ︸

ΦFA

(3.31)

Aonde a primeira componente da equação acima - ΦFA - simboliza a parte em fase de A

com F e a segunda componente - ΦFA - a parte fora de fase.

Pode-se então buscar como que as parcelas de (3.25) se relacionam com as partes em

fase e fora de fase. Pode-se decompor a parte em fase na seguinte forma aditiva.

ΦFA =

AF

11 −α11 0 0

0 AF22 −α22 0

0 0 AF33 −α33

+

α11 0 0

0 α22 0

0 0 α33

(3.32)

Então, procura-se o termo ϕF+A na decomposição acima, por exemplo, na primeira par-

cela.

Tem-se que A :F = ϕF+A :F :

(AF11 −α11)λ1 + (AF

22 −α22)λ2 + (AF33 −α33)λ3 = AF

11λ1 + AF22λ2 + AF

33λ3 (3.33)

3∑i=1

AFi iλi −

3∑i=1

αi iλi︸ ︷︷ ︸0

=3∑

i=1AF

i iλi (3.34)

Aonde λ representa os autovalores de F .

Tem-se também que ϕF+A :ϕ

F+−A = 0 , e por isso:

(AF11 −α11)α11 + (AF

22 −α22)α22 + (AF33 −α33)α33 = 0 (3.35)

26

Da condição colocada em (3.34), obtém-se o seguinte sistema:

AF11 −α11 = ξλ1

AF22 −α22 = ξλ2 (3.36)

AF33 −α33 = ξλ3

Por fim, aplicando este conjunto de equações ao que foi dito em (3.33) tem-se:

ξλ21 +ξλ2

2 +ξλ23 = A :F

ξ(F :F ) = A :F

ξ= A :FF :F (3.37)

Para uma demonstração mais completa, incluindo a prova do porque a parcela ϕF+−A não

está na componente fora de fase favor recorrer à Thompson et al. (2010).

3.3.2 Decomposição Proporcional - Ortogonal

A segunda possibilidade que se cria é a de uma parcela PFA proporcional ao tensor F e

uma parcela PFA complementar.

PFA = ξF (3.38)

PFA = A−ξF (3.39)

Por fim, as duas possibilidades diferem essencialmente em como se aloca o termo que é

tanto coaxial à F quanto ortogonal ao mesmo.

PFA ,PF

A = ϕF+A ,ϕ

F+−A +ϕF−

A (3.40)

27

ΦFA,ΦF

A = ϕF+A +ϕF+−

A ,ϕF−A (3.41)

Fica claro portanto, que ambas as decomposições são capazes de descrever, em sua tota-

lidade o tensor de interesse. A que será utilizada nos modelos é a Em fase - Fora de Fase.

Em uma tentativa de se descrever a performance do modelo, pode-se pensar na simples

razão de magnitudes entre o que está sendo modelado e o tensor alvo.

‖modelo‖‖B‖ (3.42)

o que já define uma medida de representatividade entre modelo e tensor alvo, porém ao

pensar na parcela não capturada ortogonal ao que é modelado, e. Não é possível compor

uma relação complementar entre os dois.

‖modelo‖‖B‖ + ‖e‖

‖B‖ > 1 (3.43)

observando-se a simplificação geométrica abaixo, é possível repensar esta medida de repre-

sentatividade.

Fig. 3.3: Simplificação do conceito por trás do parâmetro R.

28

α+θ = π

2

arccos(‖modelo‖

‖B‖)+arccos

( ‖e‖‖B‖

)= π

2(3.44)

Pode-se definir então um parâmetro R capaz de quantificar o quanto um dado modelo é

capaz de representar o tensor alvo.

R=(1− 2

πarccos

(‖modelo‖‖B‖

))(3.45)

aonde R varia de 0 a 1. Desta forma, a magnitude da parcela não capturada Re pode ser

obtida da seguinte forma.

Re = 1−R (3.46)

3.4 Modelos propostos

Os modelos que virão a ser propostos trabalham sobre o tensor de Reynolds deviatórico

(2.33), serão propostos um total de seis modelos.

3.4.1 Modelo 1: B =αD

O primeiro modelo a ser analisado será um modelo linear com D. Esta construção é

similar à hipótese de boussinesq (2.29). E o coeficiente α pode ser obtido da mesma forma

em que (3.37).

3.4.2 Modelo 2: B =α01+α1D+α2D2

O segundo modelo à ser proposto é um modelo que captura as correlações não lineares

entre B e D. Pelo teorema de Cayley-Hamilton3, é possível verificar que qualquer polinômio

3 O teorema de Cayley-Hamilton afirma que a matriz satisfaz sua própria equação característica. De formasimplificada, para uma matriz A 3x3 - que é o foco deste texto.A3 − tr (A)A2 + 1

2

(tr (A)2 − tr (A2)

)A−det (A)1 = 0

29

em D pode ter sua ordem reduzida a 2. Com o auxílio de um tensor de quarta ordem 1DD ,

como definido em Thompson (2008).

1DD ≡3∑

i=1eD

i eDi eD

i eDi (3.47)

É possível encontrar a parte em fase de B com D.

1DD : B =3∑

i=1eD

i eDi eD

i eDi : Bpq ep eq

3∑i=1

eDi eD

i Bpq (eDi eD

i ) : (ep eq )

3∑i=1

eDi eD

i BpqδDi pδ

Di q =

3∑i=1

B Di i eD

i eDi

1DD : B =ΦDB (3.48)

e portanto,

ΦDB =α01+α1D+α2D2 (3.49)

onde 1 é o tensor identidade - 3x3, a obtenção dos coeficientes consiste em colocar a equação

(3.48) na base de D. Assim, constroi-se um sistema entre a diagonal principal de BD - que é

coaxial à D - e a matriz de Vondermonde que surgirá à direita.

B D

11

B D22

B D33

=

1 λ1 λ2

1

1 λ2 λ22

1 λ3 λ23

.

α0

α1

α2

(3.50)

Uma vez que D é um tensor simétrico, deve conter autovalores reais e o sistema acima

deve ser de fácil solução. Além disso, é natural esperar que o parâmetro definido em (3.45)

seja superior ao mesmo parâmetro no modelo 1.

RI I ≥RI (3.51)

30

Uma observação interessante sobre este modelo é a de que nem todos os termos de α01+α1D+α2D2 são necessariamente deviatóricos. Ao se tomar o traço deste modelo,

B =α01+α1D+α2D2

tr (B)︸ ︷︷ ︸=0

= tr (α01)+ tr (α1D)︸ ︷︷ ︸α1tr (D)=0

+tr (α2D2)

α2 =− 3α0

tr (D2)(3.52)

o que mostra que estes dois coeficientes estão correlacionados, além disso, pode-se decompor

D2 em uma parte esférica e outra deviatórica

D2 = D2es f +D2

dev (3.53)

Sendo

D2es f =

1

3tr (D2)1 (3.54)

as equações (3.52) e (3.53) quando aplicadas no modelo permitem sua redução para:

B =α1D+α2D2dev (3.55)

Conclui-se portanto que os termos não deviatóricos devem se aniquilar. Outra observação

é que tensores de base estatísticamente bidimensionais, como será o caso em estudo, não

contém informação na terceira direção, restando à identidade este papel para cumprir. O

tensor D2 é bidimensional, porém sua decomposição esférica-deviatórica, cria dois tensores

tridimensionais, aonde D2dev é capaz de suprir esta informação para o modelo.

3.4.3 Modelo 3: B =α01+α1D+α2D2 +βP

No terceiro modelo, ao se adicionar o tensor P ganha-se a capacidade de capturar, pelo

menos em parte, a ortogonalidade entre B e D como apresentado em (3.9). O procedimento

é muito similar à como funciona o modelo 2, com a diferença de que, após a parte em fase

31

ser capturada, pode-se então encontrar β da mesma forma que (3.37):

B−ΦDB = ΦD

B (3.56)

β= ΦDB : P

P : P(3.57)

Continuando o raciocínio referente à melhora progressiva das bases, é natural esperar que:

RI I I ≥RI I ≥RI (3.58)

3.4.4 Modelo 4: B =αD+βP

Este quarto modelo é função apenas das partes proporcionais à D e à P. Se torna portanto

uma referência do quanto que a captura da não linearidade acrescenta ao modelo.

O procedimento de obtenção dos coeficientes também é simples:

B =αD+βP

α= B : D

D : D(3.59)

β= (B−αD) : P

P : P(3.60)

Neste caso não fica clara a relação de RIV com os demais.

3.4.5 Modelo 5: B =αD+β01+β1P+β2P2

De forma análoga ao modelo 2, este polinômio só precisa ir até a ordem dois em P. Tal

modelo propõem a captura completa da parte de B em fase com P. O procedimento numérico

consiste na obtenção de α como se faz no modelo 1 e posterior obtenção da parte em fase

com P da mesma forma que no modelo 2. Mais uma vez, não se pode inferir como RV irá se

comportar frente aos demais.

32

3.4.6 Modelo 6: B =α01+α1D+α2D2 +β01+β1P+β2P2

Por fim, chega-se ao último e teoricamente o mais completo modelo. Com esta base deve

ser possível a captura completa de B. O procedimento numérico é similar ao modelo 3, e

após a parte em fase ser subtraída de B, a parte fora de fase é levada à base de P e então seus

coeficientes obtidos por mais um sistema simples.

Por fim, é esperado que:

RV I ≥RI I I ≥RI I ≥RI (3.61)

3.5 Adimensionalizações Propostas

Os coeficientes obtidos de cada modelo serão apresentados através de duas adimensio-

nalizações distintas, como apresentado em Alves (2014). Uma primeira com k e ε e uma

segunda com k, γ e p. O parâmetro γ pode ser visto como uma medida de intensidade do

tensor taxa de deformação para escoamentos entre placas paralelas e o parâmetro p é uma

medida de intensidade do tensor não-persistência P.

γ=p

2D : D (3.62)

p =p

P : P (3.63)

Assim, pode-se adimensionalizar os coeficientes descritos nos modelos propostos, to-

mando como base o primeiro e sexto modelo, respectivamente.

B =Cµk2

εD (3.64)

B =Cµ0k1+Cµ1k2

εD+Cµ2

k3

ε2D2 +Cβ0k1+Cβ1

k3

ε2P+Cβ2

k5

ε4P2 (3.65)

33

B =CDk

γD (3.66)

B =CD0k1+CD1k

γD+CD2

k

γ2D2 +CP0k1+CP1

k

pP+CP2

k

p2P2 (3.67)

4 O CASO EM ESTUDO

Para a obtenção destas bases será utilizada uma simulação DNS. Este tipo de simulação

não trabalha com modelos e é capaz de capturar todas as escalas do escoamento, porém

como visto em (1.11) é fácil perceber que o número de pontos em uma malha tridimensional

é da ordem de Re94 e consequentemente seu custo computacional a torna impraticável para

aplicações que exigem uma rápida resposta, tornando sua aplicabilidade restrita ao estudo

científico e a geometrias relativamente simples. Em contrapartida, é justamente esta alta

qualidade de resultados a torna interessante para o estudo de cada modelo proposto.

O caso em questão é dito um bocal convergente-divergente, executado com um algoritmo

que admite incompressibilidade. Que também pode ser pensado como um pequeno "morro",

servindo de obstáculo ao escoamento. As coordenadas são apresentadas adimensionalmente

e o canal tem altura 2 , comprimento de aproximadamente 12,56 (exatos 4π) e largura π (na

direção Z). Número de Reynolds de Reτ = 617 ou Re = 12600. Outros dados mais detalhados

podem ser obtidos em Marquillie et al. (2011). O procedimento numérico é melhor descrito

em Marquillie et al. (2008), que apesar de trabalhar com um Re mais baixo, utiliza uma

geometria similar.

Fig. 4.1: Representação tridimensional da geometria estudada com ênfase para as estruturasturbulentas. Fonte: Laval and Marquillie (2009)

34

35

A imagem (4.1) mostra a geometria em estudo com a representação dos vórtices baseada

no critério de isosuperfícies de Q para o valor 300, adimensionalizadas por uτ e metade da al-

tura do bocal. Este critério é não objetivo e caracteriza os vórtices como sendo isosuperfícies

do tensor,

Q = 1

2

(‖W‖2 −‖D‖2

)(4.1)

A simulação foi executada com aproximadamente 511 milhões de nós, 2304 x 385 x 576

divisões nas direções X,Y e Z respectivamente. Na apresentação dos resultados, é feito um

corte em um plano na direção Z, devido à simetria do escoamento, como o plano da figura

(4.2). Isso resulta em um total de 887040 (385 x 2304) nós resultantes.

Fig. 4.2: Imagem da malha utilizada com resolução de 1:16. Fonte: Marquillie et al. (2011)

Tanto a superfície inferior quanto a superior não permitem a passagem de fluido, atuando

como parede. O ponto de maior elevação da superfície inferior se localiza na coordenada de

abscissa x = 5.22 e está a y = 0.67 de altura. O sentido do escoamento se dá da esquerda para

a direita da figura (4.2) e a figura (4.3) representa o perfil de entrada.

São fornecidos dois conjuntos de dados, um arquivo de campos médios, que contém, para

cada par de coordenadas do plano, informações como: o gradiente de velocidade médio, as

velocidades médias, o tensor de reynolds. Um segundo arquivo traz todos os valores de in-

teresse para a equação de transporte do tensor de Reynolds, que não foram utilizados neste

problema, com exceção ao tensor taxa de dissipação ε cujo traço é utilizado nas adimensio-

nalizações propostas.

Uma razão entre entre as maiores distâncias entre os volumes vizinhos da malha e a

escala de kolmogorov (1.6) é mostrada na figura (4.4)

36

0 0.5 1 1.5

u

0

1

2

y

Fig. 4.3: Representação do perfil de velocidade u na entrada da geometria

Fig. 4.4: Razão entre o maior comprimento de cada volume da malha e a escala de kolmo-gorov correspondente. Fonte: Laval and Marquillie (2009)

Por fim, vale um comentário sobre a obtenção numérica do tensor P, a malha estudada

é estruturada, porém não cartesiana e não uniforme. A base dos autovetores do tensor taxa

de deformação D gira, isto é: ΩD existe e é diferente de zero. O procedimento então é

subtrair de (DW−WD) a parcela, (DΩD −ΩD D), como visto em (3.6). Como mostrado por

Thompson et al. (2010) a variação temporal desta parcela não precisa ser calculada, basta

ser feita uma decomposição Em Fase - Fora de Fase sobre a derivada material de D, que em

regime permanente se resume à:

D = ∂D

∂t︸︷︷︸=0

+u.∇D =ΦDD + ΦD

D (4.2)

onde segundo a equação (3.5), é justamente a parte fora de fase que é de interesse. Então,

(∇D) foi obtido através de uma interpolação da malha e sucessivamente P.

5 RESULTADOS E DISCUSSÃO

Foram escolhidas duas formas de apresentação dos campos de interesse, uma através de

superfícies de nível bidimensionais - quando aplicável - e para uma melhor interpretação de

alguns campos, também foram feitos cortes unidimensionais em seções transversais à direção

do escoamento. Um total de seis seções foram escolhidas e a partir de uma seção localizada

no ponto mais alto da superfície inferior, tomam-se duas seções igualmente espaçadas para

trás, em x = 1,74 e x = 3,48 respectivamente e três seções igualmente espaçadas para a frente

em x = 7,05, x = 8,89, x = 10,73, respectivamente. Além disso, também são avaliados os

valores de cada componente do tensor alvo em comparação com cada componente do tensor

modelado. O que permite verificar algumas vantagens e desvantagens de cada modelo.

Primeiramente vale observar a distribuição dos elementos do tensor alvo B, o perfil de

velocidade e o tensor taxa de deformação - mais especificamente onde suas componentes

chegam a zero - uma vez que este está presente em todas as bases em estudo. Sendo os

tensores B e D simétricos, apenas as componentes

B =

Bxx Bx y Bxz

By x By y By z

Bzx Bz y Bzz

(5.1)

Bxx , Bx y , Bxz , By y , By z e Bzz precisam ser representadas e para o tensor D, apenas as

componentes Dxx , Dx y , D y y , uma vez que além de ser simétrico, este é estatisticamente

bidimensional.

D =

Dxx Dx y 0

D y x D y y 0

0 0 0

(5.2)

37

38

Componentes de B

0 0.5 1 1.5 2.0

y

−0.01

−0.005

0

0.005

x = 1, 74

0.18 0.5 1 1.5 2.0

y

−0.005

0

0.005

x = 3, 48

0.67 1 1.5 2.0

y

−0.01

0

x = 5, 22

0.33 1 1.5 2.0

y

−0.05

−0.025

0

0.025

x = 7, 05

0 0.5 1 1.5 2.0

y

−0.01

0

0.01

x = 8, 89

0 0.5 1 1.5 2.0

y

−0.005

0

0.005

x = 10, 73

Fig. 5.1: Componentes da parcela deviatórica do tensor de Reynolds B em cada seção.

39

O perfil de velocidade u

0 1 2

u

0

0.5

1

1.5

2.0y

x = 1, 74

0.18 0.5 1 1.5 2.0

u

0

1

2

y

x = 3, 48

0 1 2

u

0.67

1

1.5

2.0

y

x = 5, 22

0 1 2

u

0.33

1

1.5

2.0

y

x = 7, 05

0 1 2

u

0

1

2

y

x = 8, 89

0 1 2

u

0

1

2

y

x = 10, 73

Fig. 5.2: Perfil de velocidade u em cada seção.

40

Componentes de D

0 0.5 1 1.5 2.0

y

−1

0

1

x = 1, 74

0.18 0.5 1 1.5 2.0

y

−1

0

1

x = 3, 48

0.67 1 1.5 2.0

y

−1

0

1

x = 5, 22

0.33 1 1.5 2.0

y

−1

0

1

x = 7, 05

0 0.5 1 1.5 2.0

y

−1

0

1

x = 8, 89

0 0.5 1 1.5 2.0

y

−1

0

1

x = 10, 73

Fig. 5.3: Componentes do tensor taxa de deformação D em cada seção.

41

Observando-se o perfil de velocidade, figuras (5.2) e as componentes de D, figuras (5.3).

pode-se perceber que Dx y - curva de cor laranja - vai à zero justamente na inversão do perfil

de velocidade correspondente e que as componentes Dxx e D y y são pequenas na região

pré-bump1, onde o perfil de velocidade tende a ser viscométrico. Obervando a região pós-

bump, na seção x = 10,73, pode-se perceber que o perfil ainda não se encontra perfeitamente

desenvolvido.

Para o primeiro modelo - B =αD, o parâmetro RI apresenta a seguinte característica:

Fig. 5.4: Parâmetro RI .

Esta é a melhor captura que a hipótese de Boussinesq (2.32) pode apresentar. É possível

observar uma a tendência de formação de uma região de baixa captura na região pré-bump

do campo de RI . Esta tendência de má captura na região se dá pela inversão do perfil de

velocidade e consequente queda expressiva da magnitude de D, uma vez que as componentes

Dxx e D y y são praticamente irrelevantes, fazendo com que o coeficiente divirja localmente -

figura (5.6) - e a captura seja fraca.

Quanto ao coeficiente Cµ,

Fig. 5.5: coeficiente Cµ.

Observa-se claramente que o coeficiente está longe de ser a constante 0.09 usualmente

adotada no modelo k-ε, além disso há uma região negativa. Os valores elevados do gráfico1 encontra-se na literatura o termo "bump"para referir-se à região convexa inferior.

42

acima se dão devido à região próxima à parede - fato evidenciado na figura (5.6), uma vez

que este coeficiente é adimensionalizado com k2

ε, que tende à crescer quando próximo à

parede.

Cµ

0 0.5 1 1.5 2.0

y

0

0.5

1

x = 1, 74

0.18 0.5 1 1.5 2.0

y

0

0.5

1

x = 3, 48

0.67 1 1.5 2.0

y

−2

−1

0

1

x = 5, 22

0.33 1 1.5 2.0

y

0

0.5

1

x = 7, 05

0 0.5 1 1.5 2.0

y

0

0.5

1

x = 8, 89

0 0.5 1 1.5 2.0

y

0

0.5

1

x = 10, 73

Fig. 5.6: Distribuição do coeficiente Cµ ao longo das respectivas seções.

43

Olhando-se para a figura (5.6) pode-se verificar que seu padrão é correlato com a com-

ponente Dx y , na região pré-bump, o coeficiente diverge em função de Dx y ser o único termo

expressivo. Na seção x = 5,22 percebe-se a tendência de desprendimento da camada limite

e consequentemente um pico no valor do coeficiente. Finalmente, na região pós-bump o

coeficiente se comporta bem em x = 7,05 e x = 8,89 já que, apesar de continuar havendo

inversão do perfil de velocidade, os termos Dxx e D y y não são mais insignificantes. Por fim,

em x = 10,73 a mesma situação da primeira seção volta a acontecer, com o coeficiente vindo

à divergir.

A distribuição bidimensional do coeficiente CD ,

Fig. 5.7: coeficiente CD .

demonstra que uma adimensionalização baseada em k - γ exibe um comportamento mais

bem definido, servindo como um parâmetro de análise melhor do que Cµ, e isto se dá devido

a dois motivos: por mais que a magnitude de D tenda à zero, a razão Dγ

não deixa o coefici-

ente apresentar saltos expressivos nessa região, além disso, mesmo onde a energia cinética

turbulenta k tende a zero, aqui a mesma é de ordem 1, o que diminui sua colaboração.

44

CD

0 0.5 1 1.5 2.0

y

−1

0

1

x = 1, 74

0.18 0.5 1 1.5 2.0

y

−1

0

1

x = 3, 48

0.67 1 1.5 2.0

y

−1

0

1

x = 5, 22

0.33 1 1.5 2.0

y

−1

0

1

x = 7, 05

0 0.5 1 1.5 2.0

y

−1

0

1

x = 8, 89

0 0.5 1 1.5 2.0

y

−1

0

1

x = 10, 73

Fig. 5.8: Distribuição do coeficiente CD ao longo das respectivas seções.

45

Componentes de (αD)

0 0.5 1 1.5 2.0

y

−10−3

0

10−3

0.002

x = 1, 74

0.18 0.5 1 1.5 2.0

y

0

0.0025

x = 3, 48

0.67 1 1.5 2.0

y

−10−3

0

10−3

x = 5, 22

0.33 1 1.5 2.0

y

0

0.01

x = 7, 05

0 0.5 1 1.5 2.0

y

0

0.005

x = 8, 89

0 0.5 1 1.5 2.0

y

−0.0025

0

0.0025

0.005

x = 10, 73

Fig. 5.9: Componentes do tensor (αD) em cada seção.

quanto à representação componente a componente deste tensor, verifica-se a falha na

captura da tridimensionalidade do fenômeno, a magnitude da componente x y tende a ser

exacerbada - como já era em D, tanto na região pré-bump quanto na região posterior. Além

disso, um segundo ponto que pode ser tido como negativo é que, ao se expressar o modelo

como um múltiplo de D, tem-se a incompressibilidade2 se apresentando na forma da simetria

2 Neste caso, a incompressibilidade fica sendo apenas ∂u∂x + ∂v

∂y = 0, sendo u e v as componentes da velocidadenas direções x e y respectivamente.

46

entre as componentes xx e y y , o que reduz ainda mais a liberdade do modelo em expressar

o tensor alvo.

No segundo modelo - B = α01 +α1D +α2D2, o parâmetro RI I se apresenta seguinte

maneira,

Fig. 5.10: Parâmetro RI I .

É possível verificar que a previsão de RI I > RI se confirma para todos os pontos do

domínio. Apesar do avanço, este campo retrata o limite do que se pode extrair com uma base

exclusivamente constituída por D.

Os coeficientes Cµ0 e CD0, que apresentam a mesma adimensionalização. Portanto tem a

mesma representação,

Fig. 5.11: Coeficientes Cµ0 e CD0.

47

Cµ0 e CD0

0 0.5 1 1.5 2.0

y

−0.5

0

0.5

x = 1, 74

0.18 0.5 1 1.5 2.0

y

−0.5

0

0.5

x = 3, 48

0.67 1 1.5 2.0

y

−0.5

0

0.5

x = 5, 22

0.33 1 1.5 2.0

y

−0.5

0

0.5

x = 7, 05

0 0.5 1 1.5 2.0

y

−0.5

0

0.5

x = 8, 89

0 0.5 1 1.5 2.0

y

−0.5

0

0.5

x = 10, 73

Fig. 5.12: Distribuição dos coeficientes Cµ0 e CD0 ao longo das respectivas seções.

é possível verificar que estas componentes são exatamente iguais ao elemento Bzz esca-

lanadas pela energia cinética turbulenta, o que é intuitivo uma vez que a identidade é o único

tensor capaz de igualar, em um primeiro momento, a terceira direção de B.

48

Os coeficientes Cµ1 e CD1,

Fig. 5.13: Coeficiente Cµ1.

Fig. 5.14: Coeficiente CD1.

Sobre o campo de α1 representado nas figuras (5.13) e (5.14), respectivamente por Cµ1 e

CD1 verifica-se que este é igual à α do modelo 1 -independentemente da adimensionalização,

o que condiz com a teoria:

B =α01+α1D+α2D2

B : D =α01 : D︸ ︷︷ ︸=0

+α1D : D+α2D2 : D

B : D

D : D=α1

D : D

D : D+α2

tr (D3)

D : D

α∣∣∣modelo 1

=α1 +α2tr (D3)

D : D(5.3)

Segundo a o teorema de Cayley-Hamilton para o tensor D:

D3 − tr (D)︸ ︷︷ ︸=0

D2 + 1

2

(tr (D)2︸ ︷︷ ︸

=0

−tr (D2))D−det (D)1 = 0 (5.4)

49

tomando-se o traço da equação acima:

tr (D3) = 3det (D) (5.5)

sendo o tensor D bidimensional, como é o caso analisado, seu determinante é nulo, logo

tr (D3) = 0. Portanto, retornando-se à equação (5.3), chega-se à conclusão de que:

α∣∣∣modelo 1

=α1

∣∣∣modelo 2

(5.6)

O coeficiente Cµ2, por ser adimensionalizado com ε2

k3 , já apresenta valores da ordem de

1011 quando próximo à parede, o que torna sua visualização em superfícies de nível pouco

aproveitável. Um valor mais interessante é o de CD2.

Fig. 5.15: Coeficiente CD2.

Mais uma vez, pode-se verificar como a adimensionalização baseada na intensidade da

taxa de deformação constrói um resultado mais bem comportado na parede.

Uma verificação interessante que pode ser feita neste ponto é a de comparar as figuras

(5.12) e (5.16) observando uma notável simetria entre os gráficos. Isto pode ser verificado

retornando-se à equação (3.52) aonde pode-se verificar a clara atuação destes componentes

em tentar tornar o modelo deviatórico em concordância com lado esquerdo da equação B.

50

CD2

0 0.5 1 1.5 2.0

y

−2

−1

0

x = 1, 74

0.18 0.5 1 1.5 2.0

y

−2

−1

0

1

x = 3, 48

0.67 1 1.5 2.0

y

−2.5

0

x = 5, 22

0.33 1 1.5 2.0

y

−2.5

0

2.5

x = 7, 05

0 0.5 1 1.5 2.0

y

−1

0

1

2

x = 8, 89

0 0.5 1 1.5 2.0

y

−1

0

1

x = 10, 73

Fig. 5.16: Distribuição do coeficiente CD2 ao longo das respectivas seções.

51

Cµ2

0 0.5 1 1.5 2.0

y

−2

−1

0

x = 1, 74

0.18 0.5 1 1.5 2.0

y

−2

−1

0

x = 3, 48

0.67 1 1.5 2.0

y

−2

−1

0

x = 5, 22

0.33 1 1.5 2.0

y

−2

−1

0

x = 7, 05

0 0.5 1 1.5 2.0

y

−2

−1

0

x = 8, 89

0 0.5 1 1.5 2.0

y

−2

−1

0

x = 10, 73

Fig. 5.17: Distribuição do coeficiente Cµ2 ao longo das respectivas seções.

Este gráfico replica um pouco do que foi dito para Cµ. Onde o perfil é viscométrico

- ou próximo disso - tende-se a ter um coeficiente mal comportado, além do agravante de

que neste caso a adimensionalização é feita com ε2

k3 , o que intensifica as regiões em que o

coeficiente tende a divergir.

Finalmente pode-se verificar componente a componente deste segundo modelo,

52

Componentes de (α01+α1D+α2D2)

0 0.5 1 1.5 2.0

y

−0.0025

0

0.0025

0.005

x = 1, 74

0.18 0.5 1 1.5 2.0

y

−0.002

0

0.002

0.004

x = 3, 48

0.67 1 1.5 2.0

y

−0.0025

0

0.0025

0.005

x = 5, 22

0.33 1 1.5 2.0

y

0

0.01

x = 7, 05

0 0.5 1 1.5 2.0

y

0

0.005

x = 8, 89

0 0.5 1 1.5 2.0

y

−0.0025

0

0.0025

0.005

x = 10, 73

Fig. 5.18: Componentes do tensor (α01+α1D+α2D2) em cada seção.

Contrariamente ao apresentado no primeiro modelo - figura (5.9), este segundo modelo

já começa a apresentar uma captura da terceira direção, fato evidenciado pela componente

zz, que diferentemente do que ocorre para o modelo de boussinesq, não é nula. A presença

do tensor identidade 1 no modelo 2 constrói a noção de tridimensionalidade para o modelo,

aonde, por mais que este venha a ser aniquilado conforme a equação (3.52) atesta, a tridimen-

sionalidade continua presente na parcela deviatórica do tensor D2 como afirmava a equação

53

(3.54). Além disso, não há mais a necessidade de simetria dada pela continuidade, dando

mais liberdade aos termos xx e y y .

O terceiro modelo - B = α01+α1D+α2D2 +βP - apresenta uma melhora significativa

em termos de captura - RI I I - frente aos modelos anteriores. Corroborando com o esperado

RI I I ≥RI I ≥RI . Há praticamente a captura completa do tensor em toda região do domínio,

havendo apenas uma região central onde o índice cai para 0,945.

Fig. 5.19: Parâmetro RI I I .

Este é um indicativo de nítido de que a inserção de um termo ortogonal à base expandida

em D produziu melhoras.

Pela própria construção do modelo, os campos dos coeficientes Cµ0, Cµ1, Cµ2 e CD0,

CD1, CD2 são iguais aos do modelo 2, não havendo necessidade de representá-los. O campo

bidimensional de Cβ será omitido pelos mesmos motivos de Cµ2. As superfícies de nível de

CP também não agregam muito à interpretação uma vez que uma série de saltos descontínuos

em todas a seções amostradas é observada.

54

CP

0 0.5 1 1.5 2.0

y

−1

0

1

x = 1, 74

0.18 0.5 1 1.5 2.0

y

−1

0

1

x = 3, 48

0.67 1 1.5 2.0

y

−1

0

1

x = 5, 22

0.33 1 1.5 2.0

y

−1

0

1

x = 7, 05

0 0.5 1 1.5 2.0

y

−1

0

1

x = 8, 89

0 0.5 1 1.5 2.0

y

−1

0

1

x = 10, 73

Fig. 5.20: Distribuição dos coeficientes CP ao longo das respectivas seções.

55

Componentes de (α01+α1D+α2D2 +βP)

0 0.5 1 1.5 2.0

y

−0.01

−0.005

0

0.005

x = 1, 74

0.18 0.5 1 1.5 2.0

y

−0.005

0

0.005

x = 3, 48

0.67 1 1.5 2.0

y

−0.01

0

x = 5, 22

0.33 1 1.5 2.0

y

−0.05

−0.025

0

0.025

x = 7, 05

0 0.5 1 1.5 2.0

y

−0.01

0

0.01

x = 8, 89

0 0.5 1 1.5 2.0

y

−0.005

0

0.005

x = 10, 73

Fig. 5.21: Componentes do tensor (α01+α1D+α2D2 +βP) em cada seção.

Assim como descrito pelo parâmetro RI I I a distribuição de cada componente deste mo-

delo é praticamente idêntica à do tensor alvo, figura (5.1).

56

O quarto modelo - B = αD+βP, que deve dar uma medida da relevância do termo não

linear apresenta a seguinte captura - RIV ,

Fig. 5.22: Parâmetro RIV .

O parâmetro RIV se comportou melhor do que RI I e RI , mostrando que o acréscimo de

uma parcela ortogonal via tensor P colabora muito mais do que uma expansão sobre D.

Devido à própria construção deste modelo os coeficientes replicam as características de

α do modelo 1 e β no modelo 3.

Finalmente, observa-se que uma base desprovida de 1 é incapaz de reconstruir a tridimen-

sionalidade no modelo, uma vez que P é também estatísticamente bidimensional e apesar de

o parâmetro RIV se portar melhor do que RI I , há esta defasagem intrínseca ao modelo 4.

As componentes xx,x y e y y são quase representadas em sua totalidade quando comparadas

com a figura (5.1), o que é interessante pois a simples soma de dois tensores ortogonais se

mostrou capaz de produzir um avanço grande.

57

Componentes de (αD+βP)

0 0.5 1 1.5 2.0

y

−0.005

0

0.005

0.01

x = 1, 74

0.18 0.5 1 1.5 2.0

y

−0.005

0

0.005

x = 3, 48

0.67 1 1.5 2.0

y

−0.01

0

0.01

x = 5, 22

0.33 1 1.5 2.0

y

−0.025

0

0.025

x = 7, 05

0 0.5 1 1.5 2.0

y

−0.005

0

0.005

0.01

x = 8, 89

0 0.5 1 1.5 2.0

y

−0.005

0

0.005

x = 10, 73

Fig. 5.23: Componentes do tensor (αD+βP) em cada seção.

58

No quinto Modelo - B =αD+β01+β1P+β2P2, o campo de RV apresentou uma captura

completa o que é um resultado superior às expectativas iniciais.

Fig. 5.24: Parâmetro RV .

Depois de observar como a base do quarto modelo foi capaz de descrever componente a

componente - com exceção da terceira direção - torna-se natural pensar que, ao acrescentar a

capacidade deste modelo em capturar esta terceira direção, ou seja, incrementar a identidade

1 à esta base - como foi visto nos segundo e terceiro modelos - o resultado observado em RV

torna-se compreensível.

No quinto modelo, α se apresenta igual ao modelo 1 pela própria construção, β0 se

apresenta igual à Cµ0/CD0 uma vez que é o coeficiente que multiplica a identidade e deve

capturar, em primeiro momento a terceira direção Bzz e finalmente β se comporta como no

modelo 3.

Restando o campo de β2 a ser representado, neste caso por sua adimensionalização CP2,

Fig. 5.25: Coeficiente CP2.

59

CP2

0 0.5 1 1.5 2.0

y

−1

−0.5

0

x = 1, 74

0.18 0.5 1 1.5 2.0

y

−1

−0.5

0

0.5

x = 3, 48

0.67 1 1.5 2.0

y

−1

0

x = 5, 22

0.33 1 1.5 2.0

y

0

0.01

0.02

x = 7, 05

0 0.5 1 1.5 2.0

y

−0.5

0

0.5

1

x = 8, 89

0 0.5 1 1.5 2.0

y

−0.5

0

0.5

x = 10, 73

Fig. 5.26: Distribuição dos coeficientes CP2 ao longo das respectivas seções.

Na distribuição de CP2 pode-se verificar um padrão muito similar com o registrado na

distribuição de CD2 - ou seja, simetria com a figura (5.12) - indicando que neste caso, esta é

a componente responsável tornar o modelo deviatórico.

60

No sexto e último modelo - B = α01+α1D +α2D2 +β01+β1P +β2P2, onde deveria

haver, em tese, a maior captura possível do tensor de Reynolds, o campo de RV I se expressa

da seguinte maneira,

Fig. 5.27: Parâmetro RV I .

indicando a impossibilidade da captura desta região central. A primeira hipótese foi a

de que isto se dá devido ao aumento dos termos anisotrópicos com a terceira direção, que

não são capturados por campos médios bidimensionais. Para se verificar esta ideia, pode-se

definir, por exemplo, um tensor M,

M=

0 0 Bxz

0 0 By z

Bzx Bz y 0

(5.7)

assim, pode-se verificar o impacto deste tensor dentro do tensor cheio B com o parâmetro

RM.

Fig. 5.28: Verificação da colaboração das anisotropias através do parâmetro R.

exibindo exatamente a relação de complementaridade. Mostrando que estes termos são

a razão do modelo não capturar esta região.

Além disso, como se pôde observar na figura (5.19) tudo o que poderia ser capturado

61

do tensor de Reynolds deviatórico já foi pela base 1,D,D2 e P, logo, o coeficiente β2, que

multiplica P2 não tem relevância - Se apresenta como algo da ordem de 10−7.

6 CONCLUSÕES

O modelo de Boussinesq verificou-se como o mais simples entre todos os apresentados,

em que a tridimensionalidade não é capturada, além da limitação dada pela relação de in-

compressibilidade que correlaciona alguns dos poucos termos que sobram para a construção

do modelo. O segundo modelo mostra um ganho evidente, não sofrendo de nenhum des-

ses dois problemas mencionados, porém representa um limitador para qualquer modelo não

linear que venha a ser construido sobre a base D. É no terceiro modelo que a vantagem

de se inserir um termo ortogonal fica evidente e a captura de B é praticamente absoluta. O

quarto modelo mostrou que o incremento de uma parcela ortogonal à D já é capaz de recons-

truir quase em sua totalidade as componentes do tensor de Reynolds deviatórico no plano.

Seguindo esta linha de raciocínio, o quinto modelo apresenta-se como uma solução ainda

melhor, uma vez que incorpora a terceira direção. E finalmente o sexto e mais completo

modelo apresenta o limitante do que pode ser escrito nas bases de D e P, sendo os termos

anisotrópicos com a terceira os resposáveis pela não captura absoluta do tensor B.

Neste ponto, fica evidenciado que a tentativa de se tratar os coeficientes como constantes

- especialmente se tratando de um escoamento relativamente complexo como o que foi ana-

lisado - é uma falha grave. Os coeficientes ainda precisam de uma descrição mais elaborada.

Para a realização do trabalho aqui apresentado, foram construídos uma série de algorit-

mos que permitiram a obtenção, em primeiro momento, do tensor P e dos coeficientes de

cada modelo, assim os algoritmos são capazes de extrair estes dados para uma malha bidi-

mensional, estruturada, não cartesiana e não uniforme. Servindo, por exemplo, para explorar

os dados necessários para a modelagem da equação de transporte para o tensor de Reynolds,

que não está no escopo deste trabalho.

Como sugestão para trabalhos futuros, considera-se portanto a tentativa de se descrever

os coeficientes, que pode vir, por exemplo, através da associação destes com propriedades

62

63

de tensores cinemáticos. O que pode levar à construção de um modelo RANS para ser

executado e verificado. Ou mesmo a possibilidade de se investigar a modelagem da equação

de transporte do tensor de Reynolds para este caso de bocal, que não foi abordada neste

texto.

REFERÊNCIAS BIBLIOGRÁFICAS

Felipe A. V. De Brangaça Alves. O uso de base de dados dns para construção de modelos

rans utilizando decomposições tensoriais e coeficientes de bases normalizadas. Master’s

thesis, UNIVERSIDADE FEDERAL FLUMINENSE, Julho 2014.

W.P. Jones and B.E. Launder. The prediction of laminarization with a two-equation model

of turbulence. International journal of heat and mass transfer, 15(2):301–314, 1972.

W Michael Lai, David H Rubin, and Erhard Krempl. Introduction to Continuum Mechanics.

Butterworth-Heinemann, 4a edition, 2009. ISBN 9780750685603.

B.E. Launder and B.I. Sharma. Application of the energy-dissipation model of turbulence

to the calculation of flow near a spinning disc. Letters in heat and mass transfer, 1:131–137,

1974.

Jean-Philippe Laval and Matthieu Marquillie. Direct numerical simulations of converging–

diverging channel flow. Progress in Wall Turbulence: Understanding and Modeling, pages

203–209, 2009.

Matthieu Marquillie, Jean-Philippe Laval, and Rostislav Dolganov. Direct numerical simu-

lation of a separated channel flow with a smooth profile. Journal of Turbulence, page 23,

Fevereiro 2008.

Matthieu Marquillie, Uwe Ehrenstein, and Jean-Philippe Laval. Instability of streaks in

wall turbulence with adverse pressure gradient. Journal of Fluid Mechanics, 681:205–240,

2011.

Stephen B. Pope. Turbulent Flows. Cambridge University Press, 1a edition, 2000. ISBN

9780521598866.

64

65

Roney L. Thompson. Some perspectives on the dynamic history of a material element.

International Journal of Engineering Science, 46:224–249, 2008.

Roney L. Thompson and Paulo R. Souza Mendes. Persistence of straining and flow classi-

fication. International Journal of Engineering Science, 43:79–105, 2005.

Roney L Thompson, Gilmar Mompean, and Laurent Thais. A methodology to quantify

the nonlinearity of the reynolds stress tensor. Journal of Turbulence, 11(33):1–27, Agosto

2010.

David C Wilcox. Reassessment of the scale-determining equation for advanced turbulence

models. AIAA journal, 26(11):1299–1310, November 1988.

7 APÊNDICE

7.1 Equação de transporte para o tensor de Reynolds

Uma equação de transporte para a flutuação de velocidade pode ser obtida subtraindo-se

da equação da velocidade instantânea (2.15) a equação de transporte da velocidade média

(2.28), respectivamente:

∂ui

∂t+u j

∂ui

∂x j=− 1

ρ

∂

∂xip +ν ∂

∂x j

((∂ui

∂x j+ ∂u j

∂xi

))(∂ui

∂t+u j

∂ui

∂x j

)=− 1

ρ

∂p

∂xi+ ∂

∂x j

(ν(∂ui

∂x j+ ∂u j

∂xi

)−u′

i u′j

)

onde a diferença do termo de aceleração convectiva pode ser reescrita como:

ui u j −ui u j =(ui +u′i )(u j +u′

j )−ui u j

= ui u j +ui u′j +u j u′

i +u′i u′

j −ui u j

= ui u′j +u j u′

i +u′i u′

j

exceto pela parcela u j u′i , o restante pode ser jogado para a direita da equação como um

termo de divergente, definindo-se uma equação de transporte para a flutuação de velocidade.

∂u′i

∂t+u j

∂u′i

∂x j=− 1

ρ

∂p ′

∂xi+ ∂

∂x j

(T ′

i j −u′i u′

j −ui u′j −u′

i u′j

)(7.1)

onde T ′i j = ν

(∂u′

i∂x j

+ ∂u′j

∂xi

).

Para que esta equação comece a assumir a forma de uma equação de transporte para o

tensor de Reynolds, precisa-se de um segundo índice livre, que pode ser obtido pela multi-

plicação da equação (7.1) por uma nova flutuação de velocidade u′k , além de se executar uma

66

67

nova média temporal.

u′k

∂u′i

∂t+u j u′

k

∂u′i

∂x j=− 1

ρu′

k

∂p ′

∂xi+u′

k

∂

∂x j

(T ′

i j −u′i u′

j −ui u′j −u′

i u′j

)

u′k

∂u′i

∂t+u j u′

k

∂u′i

∂x j=− 1

ρu′

k

∂p ′

∂xi+u′

k

∂

∂x jT ′

i j −u′k

∂

∂x ju′

i u′j︸ ︷︷ ︸

=0

+

+−u′k

∂

∂x jui u′

j −u′k

∂

∂x ju′

i u′j (7.2)

onde o segundo termo dentro do divergente deixa de existir segundo o que foi afirmado em

(1.5). observando-se como exemplo o primeiro termo - termo transiente - com o auxílio da

regra da cadeia,

u′k

∂u′i

∂t= ∂u′

i u′k

∂t−u′

i

∂u′k

∂t(7.3)

pode-se perceber que o tensor de Reynolds pode ser obtido caso produza-se uma segunda

equação análoga à (7.2) com os índices livres permutados

u′i

∂u′k

∂t+u j u′

i

∂u′k

∂x j=− 1

ρu′

i

∂p ′

∂xk+u′

i

∂

∂x jT ′

k j

−u′i

∂

∂x juk u′

j −u′i

∂

∂x ju′

k u′j (7.4)

e então estas duas equações podem ser somadas.

Usando-se a regra da cadeia, termo a termo a equação final se constrói,

Termo transiente:

u′k

∂u′i

∂t+u′

i

∂u′k

∂t= ∂u′

i u′k

∂t(7.5)

68

Aceleração convectiva:

u j u′k

∂u′i

∂x j+u j u′

i

∂u′k

∂x j= u j

∂u′i u′

k

∂x j(7.6)

O primeiro termo do lado direito, que contém as flutuações de pressão:

− 1

ρu′

k

∂p ′

∂xi− 1

ρu′

i

∂p ′

∂xk=− 1

ρ

(∂p ′u′

k

∂xi+ ∂p ′u′

i

∂xk

)+ 1

ρp ′

(∂u′

k

∂xi+ ∂u′

i

∂xk

)(7.7)

o primeiro termo do lado direito é usualmente transformado em um termo de divergente com

o auxílio do delta de kroenecker.

− 1

ρ

(∂p ′u′

k

∂xi+ ∂p ′u′

i

∂xk

)=− 1

ρ

∂

∂x j

(p ′u′

iδk j +p ′u′kδi j

)(7.8)

O segundo termo, parcela deviatórica das tensões:

u′k

∂

∂x jT ′

i j +u′i

∂

∂x jT ′

k j =∂

∂x j

(u′

k T ′i j +u′

i T ′k j

)−

(T ′

i j

∂u′k

∂x j+ T ′

k j

∂u′i

∂x j

)(7.9)

e reescrevendo o termo de tensões deviatóricas como: T ′i j = 2νD ′

i j , aonde D ′i j = 1

2

(∂u′

i∂x j

+ ∂u′j

∂xi

),

u′k

∂

∂x jT ′

i j +u′i

∂

∂x jT ′

k j = 2ν∂

∂x j

(u′

k D ′i j +u′

i D ′k j

)−2ν

(D ′

i j

∂u′k

∂x j+D ′

k j

∂u′i

∂x j

)(7.10)

O terceiro termo, usando a continuidade:

−u′k

∂

∂x jui u′

j −u′i

∂

∂x juk u′

j =−u′k u′

j

∂ui

∂x j−u′

i u′j

∂uk

∂x j(7.11)

e finalmente o último termo:

−u′i

∂u′k u′

j

∂x j−u′

k

∂u′i u′

j

∂x j=−

∂u′i u′

k u′j

∂x j(7.12)

69

Finalmente, a soma das equações (7.2) e (7.4) pode ser reescrita como:

∂u′i u′

k

∂t+u j

∂u′i u′

k

∂x j=− 1

ρ

∂

∂x j

(p ′u′

iδk j +p ′u′kδi j

)+ 1

ρp ′

(∂u′

k

∂xi+ ∂u′

i

∂xk

)+

+2ν∂

∂x j

(u′

k D ′i j +u′

i D ′k j

)−2ν

(D ′

i j

∂u′k

∂x j+D ′

k j

∂u′i

∂x j

)+

−u′k u′

j

∂ui

∂x j−u′

i u′j

∂uk

∂x j−∂u′

i u′k u′

j

∂x j(7.13)

a equação acima é a equação de transporte para o tensor de Reynolds, é usual ainda unir

todos os termos que contém divergente em um único,

∂u′i u′

k

∂t+u j

∂u′i u′

k

∂x j= ∂

∂x j

(− 1

ρ

(p ′u′

iδk j +p ′u′kδi j

)−u′

i u′k u′

j+

+2ν

(u′

k D ′i j +u′

i D ′k j

))−2ν

(D ′

i j

∂u′k

∂x j+D ′

k j

∂u′i

∂x j

)+

−u′k u′

j

∂ui

∂x j−u′

i u′j

∂uk

∂x j+ 1

ρp ′

(∂u′

k

∂xi+ ∂u′

i

∂xk

)(7.14)

7.1.1 Equação de transporte da energia cinética turbulenta

Para a obtenção da equação de transporte da energia cinética turbulenta, seguindo sua

definição (2.30) basta efetuar a contração indicial da equação (7.14) - ou aplicar o operador

traço na mesma,

∂u′i u′

i

∂t+u j

∂u′i u′

i

∂x j= ∂

∂x j

(− 2

ρ

(p ′u′

iδi j

)−u′

i u′i u′

j +4ν

(u′

i D ′i j

))+

−4ν

(D ′

i j

∂u′k

∂x j

)−2u′

i u′j

∂ui

∂x j(7.15)

70

após algumas manipulações com o terceiro e quarto termos, e usando a definição (2.30) para

k, chega-se à:

∂k

∂t+u j

∂k

∂x j= ∂

∂x j

(− 1

ρ

(p ′u′

iδi j

)− 1

2u′

i u′i u′

j +ν∂k

∂x j

)−ν

(∂u′

i

∂x j

∂u′i

∂x j

)+

−u′i u′

j

∂ui

∂x j(7.16)

assim, uma série de termos ainda permanece em aberto, devendo ser modelados para que esta

equação seja utilizável. Para se chegar na equação (2.41) assume-se a hipótese de que o fluxo

de energia seja contrário ao gradiente de energia cinética - encontrada na literatura como:

gradient-diffusion hypothesis, similar ao que é feito para o fluxo de calor, por exemplo.

− 1

ρ

(p ′u′

iδi j

)− 1

2u′

i u′i u′

j =νt

σk

∂k

∂x j(7.17)

além disso, reescrevendo o último termo com a própria hipótese de Boussinesq (2.32), chega-

-se a forma da equação (2.41)

∂k

∂t+u j

∂k

∂x j= ∂

∂x j

((ν+ νt

σk

) ∂k

∂x j

)+νt

(∂ui

∂x j+ ∂u j

∂xi

)∂ui

∂x j−ν∂u′

i

∂x j

∂u′i

∂x j(7.18)