Estruturas coerentes e modelos reduzidos para o escoamento ao redor de um cilindro no regime.. por Iago de Carvalho Barbeiro - Versão HTML

ATENÇÃO: Esta é apenas uma visualização em HTML e alguns elementos como links e números de página podem estar incorretos.
Faça o download do livro em PDF, ePub, Kindle para obter uma versão completa.

UNIVERSIDADE DE SÃO PAULO

ESCOLA POLITÉCNICA

IAGO DE CARVALHO BARBEIRO

Estruturas coerentes e modelos reduzidos para o escoamento

ao redor de um cilindro no regime bidimensional periódico

São Paulo

2012

IAGO DE CARVALHO BARBEIRO

Estruturas coerentes e modelos reduzidos para o escoamento

ao redor de um cilindro no regime bidimensional periódico

Tese apresentada à Escola Politécnica da

Universidade de São Paulo para obtenção do

título de Doutor em Engenharia.

São Paulo

2012

IAGO DE CARVALHO BARBEIRO

Estruturas coerentes e modelos reduzidos para o escoamento

ao redor de um cilindro no regime bidimensional periódico

Tese apresentada à Escola Politécnica da

Universidade de São Paulo para obtenção do

título de Doutor em Engenharia.

Área de Concentração:

Engenharia Mecânica de Energia e Fluidos

Orientador:

Prof. Dr. Julio Romano Meneghini

São Paulo

2012

Este exemplar foi revisado e alterado em relação à versão original, sob

responsabilidade única do autor e com a anuência de seu orientador.

São Paulo, 30 de abril de 2012

Assinatura do autor

Assinatura do orientador

FICHA CATALOGRÁFICA

Barbeiro, Iago de Carvalho

Estruturas coerentes e modelos reduzidos para o escoamen-

to ao redor de um cilindro no regime bidimensional periódico /

I.C. Barbeiro. -- São Paulo, 2012.

59 p.

Tese (Doutorado) - Escola Politécnica da Universidade de

São Paulo. Departamento de Engenharia Mecânica.

1. Dinâmica dos fluidos 2. Modelos matemáticos 3. Estabili-

dade de sistemas I. Universidade de São Paulo. Escola Politéc-

nica. Departamento de Engenharia Mecânica II. t.

AGRADECIMENTOS

Ao professor Meneghini, pelo convite, pelo grande apoio e principalmente pela amizade

ao longo desses anos. Ao professor Aranha, que me trouxe o desafio desta tese, pelas

batalhas compartilhadas na compreensão dos problemas e na busca pelas questões.

À minha mãe e à minha irmã Fernanda queridas, pelo porto aconchegante sempre

posto ao viajante quase sempre ausente. Ao meu pai, pela confiança e pelo suporte ao

longo desses caminhos. Ao meu irmão Lincoln, companheiro nesses quase vinte e cinco

anos, por me tentar ao risco das alegrias simples. Aos pequenos, Bianca e Matheus, pelas

brincadeiras e pela alegria.

Aos grandes amigos do 54E: Freitas, Finamore, Furlaneto, Roque Jr., Fred, Luís,

Doug e Ricardo, pela amizade infalível e pelos anos de churrascos e risadas.

Aos companheiros do NDF: Roque Jr. (de novo), pelos tantos socorros (não)

técnicos; Pedrão, pela camaradagem; Ivan & Monica, pelos passeios e conversas; Gioria,

pela prontidão e precisão; Brunoc, pelo primeiro código de CFD e pela revisão cuida-

dosa; Pjabardo, pelo C++ e pelos brados; Pepe, pelas aulas de boxe; Baliño, também

pelos churrascos; Fsaltara, pela primeira orientação; Clóvis, pelas discussões e recomen-

dações; Lauro, pela atenção nas respostas; Fernanda, pela paciência; Ivone, pelo capricho;

Stergios, Adson, Franzini, Nemoto, Provasi, Alfredo e Chaves, pelas conversas e cafés.

À Fernanda Canavêz, por ter acreditado no nosso mergulho, pelo afeto e, é claro,

pelas vírgulas, gerúndios e outros realocados ou extraídos em prol da leitura deste texto.

Aos amigos de motocicletas: Marcio & Dó, Chico Landi e Tosi, pelos passeios,

experiências, socorros e histórias que completaram esses anos com muito mais que 2 rodas.

À FAPESP, pela bolsa de doutorado.

À FINEP e à Petrobras, pelos financiamentos que possibilitaram esta pesquisa.

ii

"En el fondo de un corredor, un no previsto muro me cerró

el paso, una remota luz cayó sobre mí. Alcé los ofuscados

ojos: en lo vertiginoso, en lo altíssimo, vi un círculo de cielo

tan azul que pudo parecerme de púrpura. Unos peldaños de

metal escalaban el muro. La fatiga me relajaba, pero subí, sólo

deteniéndome a veces para torpemente sollozar de felicidad.

Fui divisando capiteles y astrágalos, frontones triangulares y

bóvedas, confusas pompas del granito y del mármol. Así me

fue deparado ascender de la ciega región de negros laberintos

entretejidos a la resplandeciente Ciudad."

Jorge Luis Borges

iii

RESUMO

Esta tese trata o escoamento ao redor de um cilindro logo após a sua primeira insta-

bilidade, dentro do seu regime bidimensional periódico. A abordagem é principalmente

teórica, passa por experimentos e culmina em uma importante parte numérica que com-

plementa a teoria com evidências e ilustrações.

As principais contribuições são a análise sobre a composição modal da solução

dentro do regime periódico e o método desenvolvido para identificar autovetores de uma

linearização da equação de Navier-Stokes presentes em uma dada solução. As bases com-

postas pelos autovetores identificados servem para a projeção da equação de Navier-Stokes

e dão a essência dos modelos reduzidos deste estudo.

A aplicação numérica apresentada para Re = 60 traz duas iterações do processo,

com duas bases de autovetores de dimensões 12 e 24. Os modelos reduzidos são nume-

ricamente estáveis e a sua integração apresenta custo várias ordens mais baixo que o da

simulação numérica completa. As séries temporais das coordenadas e as bases de auto-

vetores possibilitam a recomposição do escoamento e a sua comparação com a simulação

numérica de referência. A análise de aderência foi baseada nas médias temporais, nos

valores de Strouhal e na estrutura dos harmônicos.

Ambos modelos reduzidos têm correspondência próxima com o comportamento

assintótico do escoamento e a tendência convergente das iterações é clara. As simetrias

espaciais e temporais dos harmônicos são facilmente identificadas na estrutura dos mode-

los, de forma que as bases construídas podem ser entendidas como conjuntos de estruturas

coerentes do fenômeno.

Palavras-chave: escoamento ao redor do cilindro, estruturas coerentes,

dinâmica linearizada, dinâmica não linear e modelos reduzidos.

iv

ABSTRACT

This thesis concerns the flow past a cylinder just after its first bifurcation, within its

two-dimensional periodic regime. The approach is mainly theoretical, goes through expe-

riments and is concluded by an important numerical part which complements the theory

with evidences and illustrations.

The main contributions are the analysis concerning the modal composition of

the solution within the periodic regime and a method to identify eigenvectors of some

linearizaton of the Navier-Stokes equation participating on a given solution. The bases

spanned by the identified eigenvectors are employed in the projection of the Navier-Stokes

equation and are central to the reduced models of this study.

The numerical results for Re = 60 present two iterations of the process, with

two bases of dimensions 12 e 24. The reduced models are numerically stable and their

integration is many orders less costly than that of the full simulation. The time series of

the modal coordinates and the eigenvectors bases allow the recomposition of the flow and

its comparison with the full simulation results. The convergence analysis was based on

the time averages, the Strouhal number values and the harmonic structure.

Both reduced models have close correspondence with the asymptotic behavior of

the flow and the convergent trend of the iterations is clear. The space and time symmetries

of the harmonics have a simple representation within the structure of the models, therefore

the identified bases can be understood as sets of coherent structures of the phenomenon.

Keywords: flow past a cylinder, coherent structures, linearized dyna-

mics, nonlinear dynamics and reduced models.

v

LISTA DE FIGURAS

Figura - 2.1 Extraído de Van Dyke (1982): desenvolvimento das bolhas estacioná-

rias. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

Figura - 2.2 Extraído de Prandtl (1934): esteira de von Kármán com Re = 250: (a)

câmera fixa ao referencial do cilindro, (b) câmera com a velocidade do

escoamento incidente. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

Figura - 2.3 Harmônico zero (média temporal): PIV & SIM.

. . . . . . . . . . . . . . . . . . . . 16

Figura - 2.4 Primeiro harmônico: PIV & SIM.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

Figura - 2.5 Segundo harmônico: PIV & SIM.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

Figura - 2.6 Terceiro harmônico: PIV & SIM.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

Figura - 4.1 Espectro do operador discreto linearizado na solução estacionária do

escoamento na bifurcação, em Re = 46.05, extraído de Lopez, Meneghini

e Saltara (2008). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

Figura - 6.1 Malha computacional com 1386 vértices, 4104 arestas e 2718 elementos.

45

Figura - 6.2 Solução estacionária para Re = 60: contornos de magnitude de veloci-

dade e linhas de corrente definindo as bolhas recirculantes. . . . . . . . . . . . 45

Figura - 6.3 U(t∗)=Us+U�(t∗): contornos de magnitude de velocidade. . . . . . . . . . . 46

Figura - 6.4 Espectro reduzido para Re = 60 com 12 autovalores, sendo 6 associa-

dos a autovetores simétricos (•) e mais 6 associados a autovetores anti-

simétricos (◦). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

vi

Figura - 6.5 Espectro reduzido para Re = 60 com 24 autovalores, sendo 12 asso-

ciados a autovetores simétricos (•) e mais 12 associados a autovetores

anti-simétricos (◦). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

Figura - 6.6 Linhas de corrente com as bolhas recirculantes da solução estacionária,

em laranja, e da média temporal da simulação de referência, em branco.

49

Figura - 6.7 Médias temporais das linhas de corrente com as bolhas recirculantes da

simulação de referência, em branco, do primeiro modelo reduzido (n =

12), em amarelo, e do segundo modelo reduzido, em vermelho (n = 24).

50

Figura - 6.8 Identificação dos harmônicos a que pertencem os autovetores. . . . . . . . 51

Figura - 6.9 Espectro reduzido para Re = 60, o mesmo da figura 6.4, agora identifi-

cando a qual dos harmônicos pertence cada modo. . . . . . . . . . . . . . . . . . . . 52

vii

SUMÁRIO

1 Introdução. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

2 O escoamento e a abordagem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.1 Observações experimentais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.2 Equações do movimento fluido . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

2.3 Modos e estruturas coerentes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.4 Decomposição em harmônicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

3 Discretização das equações . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.1 Discretização espacial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.2 Condições de contorno . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

3.3 Discretização temporal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

3.4 Solução estacionária . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

3.5 Malha computacional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

3.6 Código computacional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

4 Dinâmica do escoamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

4.1 Dinâmica linearizada . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

4.2 Dinâmica não linear. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

viii

5 Reconstrução modal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

5.1 Simetria dos autovetores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

5.2 Análise dos harmônicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

5.3 Identificação dos autovetores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

5.4 Projeção não linear . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

6 Aplicação numérica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

6.1 Modelo discreto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

6.2 Bases e espectros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

6.3 Modelos reduzidos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

7 Conclusões e trabalhos futuros. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

Referências . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

Apêndice A - Lista de publicações . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

A.1 Publicações em Revistas Internacionais (2) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

A.2 Publicações em Congressos Internacionais (9) . . . . . . . . . . . . . . . . . . . . . . . . . . 58

ix

1

1

INTRODUÇÃO

Esta tese desenvolve um estudo fundamental sobre o escoamento ao redor de um cilindro

logo após a sua primeira instabilidade, no início do regime bidimensional periódico. A

abordagem é principalmente teórica, passa por experimentos e culmina em uma impor-

tante parte numérica que complementa a teoria com evidências e ilustrações.

O escoamento ao redor do cilindro é um importante protótipo dos escoamentos ao

redor de corpos rombudos e a sua compreensão ainda desafia várias áreas da engenharia

e das ciências fundamentais. Na engenharia uma grande preocupação vem do seu caráter

oscilatório que deve ser levado em conta na concepção de estruturas que são projetadas

para resistir à ação de correntes de ar ou água. Na física e na matemática ainda existem

muitas questões cercando a sequência de regimes que levam este escoamento do seu estado

estacionário aos turbulentos.

Na engenharia nacional, os maiores interesses vêm da indústria petrolífera, cuja

produção é principalmente oriunda de poços submarinos. O problema está nas tubulações

que conectam os poços às plataformas flutuantes: elas podem passar dos 2.000 metros

de comprimento e devem resistir constantemente às fortes correntes da nossa costa sob a

ação do fenômeno de vibrações induzidas por vórtices. Esse fenômeno está relacionado

ao caráter oscilatório dos escoamentos ao redor de corpos rombudos que, nesse caso, é

capaz de forçar vibrações com amplitudes da ordem do diâmetro dos tubos, preocupando

os projetistas no que diz respeito à resistência das estruturas à fadiga. Nesse sentido, é

esperado que a melhor compreensão dos fenômenos envolvidos gere ferramentas de cálculo

mais sofisticadas que aumentem a segurança e a economia dos projetos.

2

Na vista da física e da matemática, os interesses focam na natureza dos compor-

tamentos dinâmicos observados com o aumento do seu parâmetro, o número de Reynolds1.

Tais comportamentos já foram, em grande parte, caracterizados experimentalmente e tam-

bém reproduzidos numericamente, mas ainda carecem de bases teóricas consistentes para

o estudo dos mecanismos envolvidos. Matematicamente, as mudanças de comportamento

surgem de bifurcações que são responsáveis pelo aumento da complexidade dinâmica dos

escoamentos, no que é chamado por Ruelle (1989) de rota para a turbulência.

Uma abordagem comum para o estudo dessas bifurcações parte de modelos redu-

zidos construídos através da projeção da equação de Navier-Stokes sobre pequenas bases

de funções ou modos. Um exemplo clássico de modelo reduzido é o de Lorenz (1963),

que trata numericamente a questão da sensibilidade de sistemas meteorológicos com rela-

ção às suas condições iniciais. O modelo de Lorenz é uma simplificação de um problema

convectivo que leva em conta apenas os três modos mais relevantes de uma discretização

em modos ortogonais e que, apesar de bastante simplificado, é capaz de reproduzir os

comportamentos caóticos observados em problemas meteorológicos reais.

Os modos também são chamados de estruturas coerentes e determinam a essência

dos modelos reduzidos. A sua definição matemática não é fixa e pode vir de várias fontes.

Os modos desta tese são autovetores de uma linearização da equação de Navier-Stokes e o

diferencial técnico aqui está no método de identificação dos autovetores mais relevantes,

que foi desenvolvido durante este projeto e mostrou bons resultados na representação da

primeira bifurcação e do regime periódico do escoamento.

O capítulo 2 deste texto começa com uma descrição do escoamento através de observações experimentais, onde são identificados os seus principais regimes dinâmicos

e simetrias, e segue com o início da abordagem matemática do problema. Ainda nesse

capítulo é introduzido o conceito estruturas coerentes que compõem o escoamento e a

primeira análise do problema vem com a decomposição de séries temporais do regime

periódico em harmônicos de Fourier.

1Definido no próximo capítulo.

3

No capítulo 3, é descrita a aplicação do método dos elementos finitos na discretização espacial do problema. São também apresentados o método de cálculo da solução

estacionária e o esquema de marcha temporal empregados. Essa é a estrutura que possi-

bilita a reprodução numérica do fenômeno e que dá subsídios para a sua análise dinâmica.

O tratamento do problema dinâmico é feito no capítulo 4 e começa pela análise da equação de evolução linearizada na solução estacionária, onde é discutida a estrutura

gerada pelo seu espectro e sua base de autovetores. A não linearidade é então agregada

a essa estrutura e possibilita o acoplamento dos autovetores, aumentando a complexi-

dade das soluções possíveis e viabilizando a estabilidade do ciclo-limite observado nos

experimentos e nas simulações numéricas.

O capítulo 5 aborda a composição modal da solução periódica e traz o método desenvolvido para identificaar os autovetores mais relevantes da solução periódica estudada.

O método é apresentado com o apoio da teoria dos subespaços de Krylov e é baseado

na equação de evolução linearizada associada a uma amostra tirada de uma simulação

numérica do escoamento. A simetria demonstrada para os autovetores da linearização é

incluída no método e, nesse ponto, são obtidos dois problemas característicos reduzidos:

um relacionado aos autovetores simétricos e outro aos antissimétricos.

Os resultados numéricos apresentados no capítulo 6 têm Re = 60 e incluem as bases de autovetores identificados juntamente com os espectros reduzidos e a análise dos

respectivos modelos obtidos através da projeção da equação de evolução nas bases. A

comparação entre os modelos reduzidos e a simulação de referência é feita com base no

comportamento assintótico das soluções e mostrou boa aderência no que confere à média

temporal e à frequência dominante do escoamento.

4

2

O ESCOAMENTO E A ABORDAGEM

2.1 Observações experimentais

O escoamento ao redor de um cilindro é um sistema dinâmico autônomo cujo parâmetro é

o seu número de Reynolds1, adimensional amplamente utilizado em mecânica de fluidos na caracterização de escoamentos viscosos internos e externos. A sua complexidade dinâmica

cresce com o número de Reynolds, seguindo uma sequência de bifurcações que o levam

do seu estado estacionário estável aos regimes turbulentos. Nesta seção são apresentadas

algumas tradicionais evidências experimentais que ilustram essa sequência.

Os experimentos de Coutanceau e Bouard (1977a) e Coutanceau e Bouard (1977b)

indicam que o escoamento estacionário é observado até Re ≈ 40. No início, o escoamento

é super viscoso e mantém aderidas ao cilindro as suas linhas de corrente. Este padrão

segue até Re ≈ 5, quando acontece a separação da camada limite. Dessa separação tem

origem o par de bolhas recirculantes que ficam coladas à jusante do cilindro e caracterizam

o restante da faixa de regime estacionário. Durante todo o regime estacionário, a simetria

do domínio e das condições de contorno também se mostra no escoamento, que é espelhado

pelo plano contendo a direção da velocidade ao longe e o eixo do cilindro, como se vê na

figura 2.1.

Esta também mostra o desenvolvimento das bolhas recirculantes, cujo compri-

mento cresce linearmente com o número de Reynolds até o surgimento da bifurcação para

o regime transiente, em Re ≈ 40. Nesse ponto surge uma perturbação oscilatória cres-

1Calculado a partir da velocidade incidente do escoamento U∞, do diâmetro do cilindro D e da

viscosidade cinemática do fluido ν: Re = U∞D .

ν

index-17_1.png

index-17_2.png

5

(a) Re=13.1

(b) Re=26

Figura 2.1: Extraído de Van Dyke (1982): desenvolvimento das bolhas estacionárias.

cente, que domina o escoamento e se desenvolve até uma amplitude de saturação, como

em um ciclo-limite. Diminuindo o nível de perturbações e aumentando a precisão do

experimento, é possível retardar a bifurcação, mas não indefinidamente. Esse processo se

torna inviável depois de apenas algumas unidades de número de Reynolds e o escoamento

entra definitivamente no regime transiente.

A bifurcação muda completamente a figura do escoamento e entra em cena o

fenômeno de formação e desprendimento alternado de vórtices. O resultado é a conhecida

esteira de von Kármán, que segue até os mais altos regimes. As duas fotografias da figura

2.2 trazem dois pontos de vista da esteira oscilatória para Re = 250. Na figura 2.2(a), a câmera está fixa ao referencial do cilindro e captura o escoamento incidente perturbado

pela sua presença. Já na figura 2.2(b), a câmera tem a mesma velocidade do escoamento incidente, eliminando-o e deixando em evidência a perturbação composta pelos turbilhões

de von Kármán.

Nas fotografias das figuras 2.1 e 2.2, os escoamentos carregam partículas refletoras capazes de reproduzir no negativo as trajetórias realizadas na janela de tempo da exposi-

ção fotográfica. Nos escoamentos estacionários a exposição pode ser longa e os traços de

trajetória correspondem a segmentos das linhas de corrente do escoamento. Já nos tran-

sitórios, quanto menor for o tempo de exposição melhor será a aproximação das linhas de

corrente instantâneas. A orientação e o comprimento dos traços dão uma aproximação do

index-18_1.jpg

index-18_2.jpg

6

(a)

(b)

Figura 2.2: Extraído de Prandtl (1934): esteira de von Kármán com Re = 250: (a) câmera fixa ao referencial do cilindro, (b) câmera com a velocidade do escoamento incidente.

campo de velocidades local, bastando levar em conta o tempo de exposição da fotografia,

lembrando que aqui também o erro cresce com esse tempo.

Uma característica marcante desse escoamento é a sua frequência dominante fs,

sentida nas oscilações da força de sustentação e observada nos campos de velocidade. O

seu tratamento é feito com base no número de Strouhal, um adimensional muito utilizado

na análise de escoamentos oscilatórios que é obtido com o auxílio de uma dimensão espacial

e uma velocidade características do escoamento:

f

St = sD ,

(2.1)

U∞

onde D é o diâmetro do cilindro e U∞ é o módulo da velocidade do escoamento incidente.

No estudo experimental de Roshko (1953) são identificados três diferentes regimes dentro da faixa transiente do escoamento, marcando o aumento da sua complexidade

dinâmica com o número de Reynolds. Junto com a primeira bifurcação instala-se um

regime oscilatório regular, identificado principalmente pela sua esteira bidimensional e

pelas frequências bem definidas que a compõem, sempre múltiplas naturais da frequência

de Strouhal. Este é o regime bidimensional periódico, também chamado simplesmente

de regime estável, o qual segue com mudanças contínuas de frequência e comprimento

de onda adimensionais até Re ≈ 150. Nessa vizinhança são observadas as primeiras

7

flutuações turbulentas. Elas trazem as primeiras tridimensionalidades do escoamento

e também agem sobre a sua dinâmica fundamental. A partir desse ponto, a frequência

dominante deixa de ser regular e passa a ser identificada em espectros mais largos, através

dos seus picos. Este é o regime de transição que segue até Re ≈ 300, quando as flutuações

turbulentas se instalam de vez na esteira, dando início ao regime irregular, observado nesse

estudo até Re ≈ 10.000.

O estudo de Roshko (1953) também foi pioneiro na identificação de duas estruturas predominantes na parcela oscilatória do campo de velocidade. Através da análise

das séries temporais de velocidade, medidas em vários pontos da esteira com o auxílio

de anemômetros, foram observadas duas estruturas: uma pulsando com a frequência de

Strouhal, e outra, de menor intensidade, pulsando duas vezes mais. A comparação entre

medições feitas em pontos espelhados pelo plano médio do escoamento permitiram uma

importante constatação sobre a distribuição espacial da intensidade das duas estruturas:

a primeira é anti-simétrica em relação ao plano médio do escoamento, enquanto a segunda

é simétrica.

Esse padrão de simetria relaciona as duas estruturas identificadas, que represen-

tam os dois primeiros harmônicos da série de Fourier do campo de velocidade, com as

duas componentes da força aplicada pelo fluido sobre o cilindro: a integral de um campo

de tensões anti-simétrico sobre a superfície do cilindro só pode ter resultante na direção

normal ao plano de simetria, ou seja, na direção da força de sustentação, enquanto a inte-

gral de um campo de tensões simétrico só pode ter resultante na direção contida no plano

de simetria e ortogonal ao cilindro, ou seja, na direção da força de arrasto. Por fim, essa

relação é confirmada pelo fato da força de sustentação do cilindro oscilar principalmente

com a frequência do primeiro harmônico e a força de arrasto com a frequência do segundo

harmônico.

Esta seção não exaure as observações experimentais desse escoamento, princi-

palmente no que diz respeito a suas tridimensionalidades, mas foca nas indicações mais

8

simples e importantes que caracterizam os diferentes regimes dinâmicos identificados. En-

tre estas se destaca a frequência dominante do escoamento, cuja regularidade pode ser

aferida dos espectros das medições temporais e é determinante na identificação dos regi-

mes. Uma referência bastante completa sobre esse escoamento é Zdravkovich (1997).

2.2 Equações do movimento fluido

O comportamento de um fluido viscoso é descrito matematicamente pelas equações de

Navier-Stokes, que são a aplicação ao movimento fluido das três equações de balanço da

mecânica: de massa, quantidade de movimento e energia. Considerando o escoamento in-

compressível e a massa específica do fluido uniforme, é possível simplificar as três equações

e ainda desacoplar a equação de energia das outras duas. Esta pode ser então eliminada

tomando o escoamento como isotérmico.

Os parâmetros principais do escoamento ao redor de um cilindro circular são o

seu diâmetro D, a velocidade do escoamento ao longe U∞ e a viscosidade cinemática do

fluido ν. Utilizando esses parâmetros, é possível adimensionalizar as equações deixando

em evidência o adimensional que caracteriza esse escoamento, o número de Reynolds

U

Re =

∞D .

(2.2)

ν

O presente estudo toma as equações de Navier-Stokes em sua forma cartesiana e

bidimensional, sendo o domínio fluido contido no plano da seção transversal do cilindro. O

eixo x está sempre alinhado com o escoamento ao longe, não perturbado, enquanto o eixo

y se encontra na direção transversal ao escoamento, ortogonal ao eixo do cilindro e ao eixo

x. Seguindo as simplificações, estas são as equações do movimento fluido consideradas:

1

∂tu = −(u.∇)u +

∇2u − ∇p

Re

(2.3)

∇.u = 0,

9

onde u = ui + vj é o campo de velocidade bidimensional e p é o campo escalar de pressão.

A primeira equação é o balanço de quantidade de movimento, proveniente da

aplicação da segunda lei de Newton, e a segunda é o balanço de massa ou equação de

continuidade. O lado esquerdo da primeira equação representa a variação temporal do

campo de velocidade, enquanto o lado direito tem o termo convectivo e os dois termos de

esforços internos, viscosos e de pressão, respectivamente.

O termo convectivo vem da adoção do referencial Euleriano e representa a acele-

ração que a partícula fluida deve sofrer para continuar seguindo a sua linha de corrente.

Quando a soma dos esforços viscosos e de pressão não é capaz de anular o termo convec-

tivo, existe um desbalanço e o escoamento é transiente (∂tu �= �0).

Uma abordagem teórica geral é tratar o cilindro imerso em um domínio fluido

infinito com velocidade ao longe uniforme e unitária, eliminando completamente a influên-

cia das fronteiras externas. É possível levar esse tratamento também para a abordagem

numérica, considerando uma ou mais dimensões do domínio infinitas, como em Lavinas,

Barbeiro e Aranha (2007). Este trabalho optou pela utilização de domínios finitos, portanto, com todas fronteiras exteriores posicionadas a distâncias finitas do cilindro. No

entanto, é mantida a intenção de minimizar a influência das fronteiras externas e os do-

mínios finalmente utilizados são efetivamente extensos.

A primeira preocupação na definição de um domínio fluido finito é que ele conte-

nha as regiões de interesse do estudo, as quais, nesse caso, são as vizinhanças do cilindro

e a sua esteira. No que concerne à viabilidade do domínio para o cálculo, é necessário

que se tenha informações suficientes sobre a solução do problema nos seus contornos.

Isso porque, além de uma condição inicial, a equação de Navier-Stokes pede condições de

contorno.

O domínio considerado é bidimensional, plano e delimitado por dois contornos

fechados. O primeiro é interior ao domínio e define o contato entre o fluido e o cilindro.

Já o segundo é exterior ao domínio e o separa do escoamento ao longe, também chamado

10

de não perturbado. No contorno que define o contato do fluido com o cilindro é aplicada

a condição de aderência dos escoamentos viscosos, ao passo que no contorno exterior é

imposto o valor de velocidade ao longe, não perturbada, uniforme e unitária.

A aplicação da condição de aderência é simples e no caso do cilindro fixo recai em

impor velocidade nula ao longo de todo o contorno. A complicação da segunda condição

está no posicionamento do contorno exterior, que influencia o resultado nas proximidaes do

cilindro. Estudos experimentais de Coutanceau e Bouard (1977a) e numéricos de Lavinas,

Barbeiro e Aranha (2007) mostram que domínios reduzidos limitam o crescimento das bolhas recirculantes das soluções estacionárias, além de contribuir para o aumento da força

de arrasto sobre o cilindro. Tais estudos também mostram que aumentando o domínio

é possível chegar em uma dimensão finita a partir da qual os resultados variam muito

pouco e podem ser considerados semelhantes aos resultados obtidos com domínios muito

maiores. Esse efeito é conhecido como blocagem e ambos os estudos também observam

que a dimensão desejável para o domínio, no caso estacionário, cresce com o número de

Reynolds.

2.3 Modos e estruturas coerentes

O conceito de estrutura coerente é bastante comum na abordagem de fenômenos envol-

vendo meios contínuos. A sua noção intuitiva vem da observação desses fenômenos nos

quais padrões espaciais são identificados e associados a determinados comportamentos

temporais. Em essência, a constatação de evidências dessa natureza leva à intuição de

que existe algum mecanismo simplificado para o fenômeno que facilite a sua abordagem

e compreensão.

O tratamento matemático desses problemas é usualmente feito através de equa-

ções a derivadas parciais, onde o conceito de estrutura coerente refere-se precisamente

à separação das variáveis espaciais e temporais. Nesse contexto, as estruturas coerentes

materializam-se na forma de funções espaciais chamadas de modos que servem para a

11

construção de soluções para as equações. Os modos não têm definição matemática pres-

crita nem forma única, de maneira que conjuntos distintos de modos podem apresentar

performances semelhantes na aproximação de uma dada solução.

A construção de uma solução é feita compondo-se modos de um dado conjunto

ponderados por coordenadas modais que, nos problemas transientes, são funções tempo-

rais. É conveniente que, dado um conjunto de modos, haja apenas uma composição de

coordenadas equivalente a uma dada solução e para que isso seja sempre possível basta

que o conjunto de modos tenha a mesma dimensão do espaço gerado por eles, ou seja,

que ele seja uma base desse espaço. No caso dessa condição não ser evidente, é conveni-

ente e simples ortogonalizar o conjunto de modos chegando-se à sua dimensão e a uma

base modal associada. Por fim, é a projeção da dada equação a derivadas parciais nessa

base que dá origem ao respectivo modelo reduzido, que tem as coordenadas modais como

variáveis. Assim, uma solução aproximada com o auxílio de uma dada base modal toma

a seguinte forma:

n

u�(x, t) =

qi(t)ψi(x),

(2.4)

i=1

onde u�(x, t) é a solução aproximada, qi(t) são as coordenadas, ψi(x) são os modos e n é

a dimensão da base modal.

A abordagem mais tradicional para escoamentos oscilatórios considera modos

vindos da base de autovetores do operador de Navier-Stokes linearizado em torno de

uma solução estacionária ou periódica. Este é um caminho conservador, onde os modos

têm sentido matemático e físico bem estabelecidos, mas vem com uma contrapartida

importante: muitas vezes a grande dimensão dos problemas característicos envolvidos

inviabiliza o cálculo da base completa de autovetores, obrigando a redução do foco dos

estudos para pequenas regiões do espectro.

No caso da primeira bifurcação do escoamento ao redor de um cilindro, é comum

os estudos se limitarem à região instável do espectro, calculando o par de valores e vetores

característicos que a protagonizam, como em Jackson (1987). Com os elementos da análise 12

linear é possível determinar o tipo de bifurcação e também prever o comportamento a

curto prazo das soluções transientes que partem do estado estacionário. Mas a influência

do termo quadrático logo tira a solução do plano definido por esse par de autovetores,

levando-a para um estado periódico assintótico que é, de fato, muito distante da previsão

linear, dado que estável e com uma frequência muito distante da componente imaginária

do par de autovalores instáveis, como mostra Barkley (2006).

Uma alternativa para a não linearidade é a expansão da solução do problema

em séries assintóticas, considerando as primeiras interações não lineares dos autoveto-

res instáveis. O modelo de Landau e Lifshitz (1959), bastante conhecido na literatura, quando expandido até terceira ordem, reproduz o ciclo-limite do fenômeno e apresenta

uma correção para a frequência prevista pela componente imaginária do par de autovalo-

res instáveis. No entanto, mesmo após a primeira correção assintótica, a sua frequência

ainda é distante das observações experimentais, como é visto em Sipp e Lebedev (2007).

Duas técnicas vêm sendo bastante aplicadas na busca por estruturas coerentes

em escoamentos: a decomposição de Karhunen-Loève, também conhecida como POD, de

Proper Orthogonal Decomposition, e a decomposição em modos de Koopman. Ambas se

baseiam em observações do fenômeno, que podem vir tanto de medições experimentais

como de simulações numéricas e, por este motivo, são chamadas de empíricas, já que não

têm relação implícita com o operador de Navier-Stokes.

A decomposição de Karhunen-Loève é baseada em uma matriz de covariância, a

qual, nos estudos envolvendo escoamentos, conforme Holmes, Lumley e Berkooz (1988),

é gerada a partir de um conjunto de amostras de campos de velocidade. Em princípio

não há necessidade de organização temporal das amostras, mas é comum que sejam séries

temporais. Os modos da decomposição são os autovetores da matriz de covariância, sendo

que a relevância energética de cada um é dada pelo respectivo autovalor. Essa relevância

é utilizada como critério na construção de bases reduzidas, onde os autovetores menos

energéticos são descartados, como em Ma e Karniadakis (2002).

13

A teoria de Koopman fornece um operador linear capaz de reproduzir uma dada

série temporal, mesmo que esta seja fruto de um fenômeno não linear, como mostram

Mezic (2005) e Rowley et al. (2009). A não linearidade é compensada pela dimensão arbitrária do operador de Koopman, cujos autovetores, que no caso geral são complexos,

são os modos procurados. Os autovalores complexos associados aos modos de Koopman

fornecem a taxa de crescimento e a frequência com que eles participam da série temporal

em questão. O produto final é um problema linear capaz de emular a evolução do pro-

blema não linear, pelo menos no trecho da sua história representado pela série temporal

analisada. No caso de fenômenos periódicos, os modos de Koopman são equivalentes aos

harmônicos de Fourier da série temporal, enquanto os autovalores são puramente imagi-

nários e trazem as respectivas frequências dos harmônicos.

2.4 Decomposição em harmônicos

Partindo dos resultados de Roshko (1953) e do conceito das estruturas coerentes foi deci-dido que a decomposição do escoamento em seus harmônicos de Fourier traria elementos

importantes para a análise do regime periódico. É desse passo que surgem os primeiros

indícios sobre as funções espaciais que compõem a solução observada nos experimentos

e nas simulações numéricas. Os resultados que seguem também aparecem em Barbeiro,

Aranha e Meneghini (2007) e Korkischko et al. (2010).

As séries temporais dos campos de velocidade utilizadas nas decomposições apre-

sentadas aqui têm duas origens: simulações numéricas e experimentos em água. As simu-

lações numéricas foram computadas pelo código de elementos finitos construído no âmbito

deste projeto e descrito no próximo capítulo, enquanto as medições experimentais foram

realizadas no canal de água recirculante do Núcleo de Dinâmica e Fluidos (NDF), em

parceria com o então doutorando Ivan Korkischko e utilizando a técnica de PIV (Particle

Image Velocimetry, ver Adrian (1991)). O número de Reynolds para a comparação dos resultados foi fixado em 100, que é próximo do limite inferior imposto pelas instalações

14

experimentais. Em ambos os casos foram calculados desde o harmônico zero até o terceiro.

Nas simulações numéricas os harmônicos foram calculados em paralelo com a inte-

gração temporal. Contando com aproximadamente mil passos de tempo para cada período

de oscilação da esteira, as integrais dos harmônicos foram calculadas com grande precisão

utilizando o método dos trapézios sobre um intervalo de tempo equivalente a dez ciclos

de oscilação da esteira. Já o cálculo dos harmônicos experimentais foi feito a posteriori,

a partir das medições acumuladas pelo aparato de PIV. A capacidade do equipamento

possibilitou a medição de 345 campos de velocidade equispaçadamente distribuídos em

nove ciclos de oscilação da esteira. Nesse caso, os harmônicos foram calculados com o

emprego do método dos mínimos quadrados e, apesar de bem menos detalhadas, as séries

experimentais geraram resultados comparáveis aos numéricos.

As constatações mais notáveis dessa série de resultados vêm das relações de si-

metria espacial observadas nos harmônicos. Essas relações estão ligadas à geometria do

domínio, que é simétrico em relação à linha média do escoamento definida pelo eixo x e,

por isso, tem todos os seus pontos (x, y) associados a respectivos (x, −y). Essa propriedade

do domínio permite que qualquer função escalar f(x, y) nele definida seja decomposta em

uma parcela simétrica fs(x, y) e outra antissimétrica fa(x, y)

1

f s(x, y) =

(f (x, y) + f (x, −y))

2

(2.5)

1

f a(x, y) =

(f (x, y) − f(x, −y)) ,

2

que são complementares, dado que a sua soma leva de volta à função inicial. No caso de

funções vetoriais, tais como campos de velocidade, a relação de simetria faz a mímica de

um espelho posicionado na linha média do escoamento, associando os vetores f(x, y)2 aos seus reflexos f(x, −y). No caso em questão, de um plano bipartido por uma linha média,

a simetria entre vetores posicionados em pontos opostos é constatada pela igualdade das

componentes paralelas à linha média e pela simples oposição no sinal das componentes

2 f(x, y) = fx(x, y)i + fy(x, y)j.

15

ortogonais à linha. Isso quer dizer que, em um campo vetorial simétrico3, a componente paralela à linha é um campo escalar simétrico enquanto a componente ortogonal é um

campo escalar antissimétrico

f s(x, y) = f sx(x, y)i + fay(x, y)j

(2.6)

f a(x, y) = f ax(x, y)i + fsy(x, y)j,

onde as parcelas simétricas e antissimétricas de fx(x, y) e fy(x, y) vêm de 2.5.

Dados os campos escalares f(x, y) e g(x, y) com o produto escalar e a norma

usuais

�f(x, y), g(x, y)� =

f (x, y)g(x, y)dΩ

(2.7)

�f(x, y)� =

�f(x, y), f(x, y)�

onde Ω é o domínio em questão, é possível quantificar as parcelas simétrica e antissimétrica

de f(x, y) através da seguintes razões

�fs(x, y)�

�fa(x, y)�

Rs =

; R

(2.8)

�f(x, y)�

a = �f(x, y)�

lembrando que

Rs + Ra = 1,

(2.9)

pois

�fs(x, y), fa(x, y)� = 0

(2.10)

para qualquer f(x, y), por simetria. No que segue, utilizando essas razões, os campos são

chamados de simétricos quando Rs > 0.95 e antissimétricos quando Ra > 0.95.

3Em um campo vetorial antissimétrico, o quadro se inverte, de modo que a componente paralela à linha de simetria é um campo escalar antissimétrico e a ortogonal é um campo escalar simétrico, garantindo que as parcelas simétrica e antissimétrica sejam complementares na reconstrução de f(x, y).

index-28_1.jpg

index-28_2.jpg

index-28_3.png

index-28_4.png

16

A seguir são apresentados campos escalares das componentes de velocidade x e

y dos harmônicos de Fourier, nas figuras 2.3, 2.4, 2.5 e 2.6. Com exceção do termo zero, que é real puro, todos os outros são complexos4 e, por isso, são apresentados em suas partes reais e imaginárias r e i. Os contornos são apresentados em escalas de cinza e

os extremos da escala são dados logo abaixo de cada figura entre colchetes (o sentido da

escala é usual, do branco para o preto). Os diferentes harmônicos são apresentados em

figuras separadas, contendo duas fileiras horizontais de imagens cada um, com a primeira

referente aos experimentos (PIV) e a segunda referente às simulações numéricas (SIM). O

enquadramento é o mesmo em todos os gráficos, com o cilindro posicionado à montante

para enfatizar a região da esteira próxima, onde se localizam os máximos dos harmônicos.

(a) PIV - u0,r,x[-0.09;1.27] (b) PIV - u0,r,y[-0.19;0.33]

(c) SIM - u0,r,x[-0.17;1.26] (d) SIM - u0,r,y[-0.64;0.64]

Figura 2.3: Harmônico zero (média temporal): PIV & SIM.

A figura 2.3 traz o harmônico zero, que é simplesmente a média temporal do escoamento. A sua componente x deixa em evidência a esteira média paralela, que começa

logo atrás da região de recirculação e segue até grandes distâncias à jusante. Esse para-

lelismo é confirmado pelos contornos da componente y, que decaem logo após a região

de recirculação. Assim como a solução estacionária, a média temporal é um campo de

velocidade simétrico, com a sua componente escalar x simétrica e a y antissimétrica.

4Os harmônicos complexos são entidades oscilatórias definidas a menos de uma fase arbitrária θ, de modo que um harmônico u0n(x, y) = u0n,r(x, y) + iu0n,i(x, y) é sempre equivalente às suas rotações complexas uθn(x, y) = eiθu0n(x, y), para qualquer ângulo θ real. Os harmônicos das séries numéricas e experimentais foram rotacionados até chegarem em fases próximas com o objetivo de facilitar a sua comparação visual.

index-29_1.jpg

index-29_2.jpg

index-29_3.jpg

index-29_4.jpg

index-29_5.png

index-29_6.png

index-29_7.png

index-29_8.png

17

O primeiro harmônico, visto na figura 2.4, é o protagonista da dinâmica do problema, oscilando com a frequência de Strouhal e contendo a instabilidade geradora do

regime periódico. Como já havia sido notado por Roshko (1953), a principal componente da esteira de von Kármán é antissimétrica e lidera a hierarquia dos harmônicos com a

maior amplitude. O seu comprimento de onda local5 mostra o espaçamento entre os turbilhões, que não é uniforme, mas crescente à medida que os turbilhões são convectados.

Esse comportamento é consequência da relação de dispersão fixada pela frequência de

oscilação e a velocidade local da média temporal, que dita a velocidade de convecção dos

turbilhões. A antissimetria implica em resultantes de força sobre o cilindro que têm dire-

ção ortogonal à linha de simetria, esclarecendo o comportamento da força de sustentação

medida, que oscila predominantemente com frequência de Strouhal.

(a) PIV - u1,r,x[-0.28;0.27] (b) PIV - u1,i,x[-0.34;0.33] (c) PIV - u1,r,y[-0.49;0.55] (d) PIV - u1,i,y[-0.56;0.51]

(e) SIM - u1,r,x[-0.41;0.41] (f) SIM - u1,i,x[-0.41;0.41] (g) SIM - u1,r,y[-0.60;0.63] (h) SIM - u1,i,y[-0.56;0.64]

Figura 2.4: Primeiro harmônico: PIV & SIM.

A primeira parcela oscilatória simétrica vem com o segundo harmônico, visto na

figura 2.5. Os seus valores máximos situam-se mais à jusante e a sua distribuição espacial respeita a mesma relação de dispersão descrita para o caso do primeiro harmônico, sendo

que agora o comprimento de onda local deve ser a metade do primeiro, já que a sua

frequência é o dobro e a velocidade média local é a mesma. É do segundo harmônico que

vem a principal parcela oscilatória da força de arrasto medida no cilindro, pelo oposto da

relação de simetria da força resultante explicada no parágrafo anterior.

5Medido em um corte feito na linha média (y = 0) da componente y do primeiro harmônico.

index-30_1.jpg

index-30_2.jpg

index-30_3.jpg

index-30_4.jpg

index-30_5.png

index-30_6.png

index-30_7.png

index-30_8.png

index-30_9.jpg

index-30_10.jpg

index-30_11.jpg

index-30_12.jpg

index-30_13.png

index-30_14.png

index-30_15.png

index-30_16.png

18

(a) PIV - u2,r,x[-0.07;0.07] (b) PIV - u2,i,x[-0.08;0.07] (c) PIV - u2,r,y[-0.09;0.10] (d) PIV - u2,i,y[-0.10;0.10]

(e) SIM - u2,r,x[-0.11;0.11] (f) SIM - u2,i,x[-0.12;0.12] (g) SIM - u2,r,y[-0.13;0.13] (h) SIM - u2,i,y[-0.13;0.13]

Figura 2.5: Segundo harmônico: PIV & SIM.

A antissimetria dos termos ímpares e a relação de dispersão também são consta-

tadas no terceiro harmônico, apresentado na figura 2.6 com o seu comprimento de onda de um terço do fundamental. Os contornos ruidosos são decorrentes de uma região de

sombra, isso porque o feixe de laser do equipamento de PIV, que incide ortogonalmente

à fronteira inferior (y = ymin), é obstruído pela presença do cilindro, o que prejudica a

qualidade da medição na região 0 < y < ymax e −0.5D < x < 0.5D.

(a) PIV - u3,r,x[-0.03;0.05] (b) PIV - u3,i,x[-0.03;0.03] (c) PIV - u3,r,y[-0.05;0.05] (d) PIV - u3,i,y[-0.05;0.05]

(e) SIM - u3,r,x[-0.05;0.05] (f) SIM - u3,i,x[-0.05;0.05] (g) SIM - u3,r,y[-0.10;0.10] (h) SIM - u3,i,y[-0.10;0.10]

Figura 2.6: Terceiro harmônico: PIV & SIM.

As constatações das simetrias temporais e espaciais são bastante fortes e ilustram

a estrutura da composição modal da solução, que será retomada no capítulo 5. Por fim, a notável aderência dos resultados numéricos aos experimentais é um importante argumento

para a validação do modelo discreto que servirá de base para o resto da análise.

19

3

DISCRETIZAÇÃO DAS EQUAÇÕES

3.1 Discretização espacial

A discretização espacial das equações de Navier-Stokes é feita através do Método de

Elementos Finitos, possibilitando o tratamento numérico do problema. O domínio é

dividido em elementos triangulares que são agrupados em malhas não estruturadas e a

formulação adotada é padrão. Utilizam-se polinômios de segundo grau para as velocidades

e de primeiro grau para as das pressões, sempre garantindo a continuidade das funções

entre os elementos. Essa configuração também é conhecida como elemento de Taylor-Hood

e satisfaz a condição div-stability, como mostram Taylor e Hood (1973), Gunzburger (1989)

e Brenner (2008).

Seguindo as referências do parágrafo anterior, as equações de quantidade de mo-

vimento e de continuidade passam pelos passos usuais da formulação fraca e tomam a

forma integral, também semelhante à apresentada por Foias O. Manley e Temam (2004):

� �

∂u.δu dΩ = − (((u.∇)u).δu)dΩ

∂t

1

+

(∇u.∇δu + ∇v.∇δv) dΩ −

(p (∇.δu)) dΩ

Re Ω

(3.1)

1

((∇u.n) δu) d∂Ω +

(p (n.δu)) d∂Ω

Re ∂Ω

∂Ω

((∇.u) δp) dΩ = 0,

onde a velocidade δu = δui + δvj e a pressão δp são as respectivas variáveis virtuais do

método das potências virtuais, Ω é o domínio fluido bidimensional, ∂Ω é a fronteira do

domínio e n é a normal exterior à fronteira. A integração por partes diminui a ordem do

20

termo de dissipação viscosa e também faz com que apareçam as integrais das tensões na

fronteira, tanto viscosas quanto de pressão.

O passo seguinte é a substituição das variáveis contínuas de velocidade e pres-

são, físicas e virtuais, pelas análogas discretas, dadas pelas funções de forma ponderadas

pelos valores nodais. Sendo um esquema de Galerkin, as variáveis físicas e virtuais são

representadas pelo mesmo tipo de função de forma.

O valor de cada componente de velocidade é aproximado no interior de cada

elemento de malha triangular pela soma de seis funções de forma quadráticas, ponderadas

por seis valores nodais. Nesse tipo de elemento, os valores nodais coincidem com os valores

da função aproximada nos nós, localizadas nos vértices dos triângulos e nos pontos médios

de suas arestas. Para garantir a continuidade do campo de velocidade, elementos vizinhos

compartilham tanto nós de vértices como de arestas, de forma que o número total de nós

de cada componente de velocidade para uma dada malha é dado pela soma do seu número

de vértices com o seu número de arestas. Sendo assim, o número de nós utilizados na

representação de um campo de velocidade bidimensional é simplesmente o dobro desse

valor:

nU = 2(nvértices + narestas),

(3.2)

onde nvértices é o número de vértices da malha e narestas é o seu número de arestas.

A pressão é aproximada no interior de cada elemento de malha triangular pela

soma de três funções de forma lineares, ponderadas por três valores nodais. Aqui também

os valores nodais coincidem com os valores da função aproximada nos nós, localizadas

apenas nos vértices dos triângulos. Os campos de pressão também devem ser contínuos e

assim o número total de nós para uma dada malha coincide com o seu número de vértices,

ou seja:

nP = nvértices.

(3.3)

21

Começando o tratamento discreto do problema, seguem os vetores de valores

nodais de velocidade

U = [U1, . . . , Un ]

U

(3.4)

δU = [δU1, . . . , δUn ],

U

onde os índices ímpares carregam os valores nodais da componente x e os pares da com-

ponente y, na mesma ordem de nós. Logo abaixo são apresentados os vetores dos valores

nodais de pressão

P = [P1, . . . , Pn ]

P

(3.5)

δP = [δP1, . . . , δPn ].

P

Após a introdução das variáveis discretas, a equação algébrica final da forma

fraca fica composta predominantemente de termos bilineares nos valores nodais, físicos e

virtuais, com coeficientes na forma de produtos das funções de forma e de suas derivadas

integrados nos elementos da malha. As exceções são o termo convectivo, quadrático

na velocidade física e linear na virtual, e as forças de vínculo que aqui incorporam a

dependência na velocidade física

� �

∂u

δUt.M. ˙

U =

.δu dΩ

∂t

δUt.N(U).U = −

(((u.∇) u) .δu) dΩ

δUt.D.U =

(∇u.∇δu + ∇v.∇δv) dΩ

Ω�

δUt.R.P = −

(p (∇.δu)) dΩ

(3.6)

Ω �

1

δUt.Fv = −

((∇u.n) δu) d∂Ω

Re ∂Ω

δUt.Fp =

(p (n.δu)) d∂Ω

∂Ω

δPt.Rt.U =

((∇.u) δp) dΩ,

onde os vetores ˙U, Fv e Fp têm o mesmo formato do vetor U e carregam, respectivamente,

22

os valores nodais de ∂tu, as forças de vínculo viscosas e as de pressão

˙

U = [ ˙

U1, . . . , ˙

Un ]

U

F

(3.7)

v = [Fv,1, . . . , Fv,n ]

U

Fp = [Fp,1, . . . , Fp,n ].

U

Concatenando o par de vetores das variáveis físicas U e P e o par de vetores das

variáveis virtuais δU e δP, é possível remontar a equação 3.1 na seguinte forma matricial:

t 

  

t 

  

t 

˙

1

δU M 0 U

δU  D + N(U) R

U

δU

F

Re

  

  v + Fp

 

   = 

 

   + 

 

 , (3.8)

δP

0

0

P

δP

Rt

0

P

δP

0

e, por fim, eliminando as variáveis virtuais δU e δP e substituindo

1

K(Re, U) =

D + N(U),

(3.9)

Re

chega-se a:

  

  

˙

M 0 U

K(Re, U) R U

Fv + Fp

   = 

   + 

 .

(3.10)

0

0

P

Rt

0

P

0

Nessa forma, o problema numérico ainda apresenta três complexidades importan-

tes: a aplicação das condições de contorno, a derivada temporal de ˙U e a não linearidade

do termo convectivo N(U). Todas são tratadas ainda neste capítulo, nas próximas seções.

3.2 Condições de contorno

A abordagem das condições de contorno é padrão, sendo que cada grau de liberdade

posicionado na fronteira deve receber alguma informação externa prévia, na forma de uma

condição do tipo natural ou do tipo essencial. O esquema de imposição das condições é

23

o mesmo, tanto no cálculo das soluções transientes quanto no das estacionárias, diferindo

apenas em relação ao tipo da condição.

As integrais das tensões no contorno definem as forças de vínculo e são armaze-

nadas nos vetores Fv e Fp, sendo que os valores Fv,i e Fp,i (referentes a um dado grau de

liberdade i) só participam efetivamente das equações dinâmicas dos graus de liberdade de

contorno, pois valem zero nas equações dos graus de liberdade internos.

As condições do tipo natural são as que trazem informações para o cálculo das

integrais de tensões referentes a certos graus de liberdade de contorno, possibilitando o

cálculo dos valores Fv,i e Fp,i, que completam as equações desses graus de liberdade. A sua

implementação é direta, bastando calcular as integrais referentes aos graus de liberdade

envolvidos e preencher as respectivas posições do vetor forçante do sistema.

Já as condições do tipo essencial são as que impõem valores a certos graus de

liberdade de contorno. Nesse caso, não existe preocupação com o cálculo das integrais

de contorno, uma vez que as equações dinâmicas dos graus de liberdade que recebem

condições essenciais são descartadas, dando espaço às equações de vínculo. A implemen-

tação desse tipo de condição se resume a zerar as linhas da matriz referentes aos graus

de liberdade em questão, colocando 1 apenas nas posições que coincidem com a diagonal

principal e substituindo os valores da forçante pelos valores a serem impostos.

Uma vez resolvido o sistema, as forças de vínculo dos graus de liberdade que

receberam condições essenciais podem ser calculadas através das respectivas equações

dinâmicas que haviam sido descartadas.

3.3 Discretização temporal

Nessa fase do processo, logo após a discretização espacial das equações, é feita a discreti-

zação em relação ao tempo. Um esquema para o cálculo da derivada temporal do campo

de velocidades é proposto e o termo ˙U sai da equação, dando lugar para uma expressão

24

contendo campos de velocidades de tempos vizinhos. A variável contínua de tempo t é

substituída pela discreta tn ou simplesmente n, seguindo:

tn = n∆t

Un = U(tn) = U(n∆t)

(3.11)

P n = P (tn) = P (n∆t),

onde ∆t é o passo de tempo.

Aqui foi empregada uma variação do método de Euler, ver Deuflhard e Bornemann

(2002), que aproxima a derivada temporal da velocidade em primeira ordem

˙

Un+1 − Un

Un+1 ≈

.

(3.12)

∆t

Introduzindo o esquema na equação 3.10 e substituindo o termo convectivo por uma aproximação também de primeira ordem, chega-se ao seguinte esquema de marcha,

que é semi-implícito

 

1

 D + ˆ

N(Un) − 1 M R

Un+1

N(Un)Un − 1 MUn − Fn

Re

∆t

 

∆t

v − Fn

p

 

 = 

 , (3.13)

Rt

0

P n+1

0

onde:

ˆ

N(Un)Un+1 = N(Un)Un+1 + N(Un+1)Un.

3.4 Solução estacionária

A estratégia escolhida para o cálculo da solução estacionária foi a imposição de ˙U = 0

seguida da resolução do problema não linear decorrente via método de ponto fixo. A

informação de simetria observada nos experimentos é utilizada e o problema é resolvido

em apenas metade do domínio, com a condição de simetria aplicada na fronteira que

corresponde à linha média do escoamento, como em Fornberg (1984).

Por fim, o esquema iterativo decorrente do método de ponto fixo exige a recons-

25

trução e a resolução do seguinte sistema linear a cada passo k:

 

1

 D + ˆ

N(Uk) R

Uk+1

N(Uk)Uk − Fk

Re

 

v − Fk

p

 

 = 

 ,

(3.14)

Rt

0

P k+1

0

onde:

ˆ

N(Uk)Uk+1 = N(Uk)Uk+1 + N(Uk+1)Uk.

3.5 Malha computacional

O código de elementos finitos desenvolvido utiliza elementos triangulares organizados em

malhas não estruturadas bidimensionais. As vantagens desse tipo de discretização são

a flexibilidade em geometrias complexas e a facilidade na construção de gradientes de

refinamento elevados, bastante valorizada em problemas envolvendo camadas limite.

As malhas deste estudo são construídas utilizando o código computacional

BAMG, de Bidimensional Anisotropic Mesh Generator, que é um gerador de malhas bi-

dimensional, anisotrópico e adaptativo. O código é parte de um projeto desenvolvido por

um grupo de pesquisa do INRIA, na França, e sua disponibilidade para fins de pesquisa é

livre. Informações mais detalhadas sobre o método empregado no código são encontradas

em Alauzet e Frey (2003a), Alauzet e Frey (2003b), Borouchaki e Fray (1998), Castro-Díaz

F. Hecht e Pironneau (1997).

Para que a simetria do problema contínuo seja mantida, a malha também deve

ser simétrica, não apenas em seu contorno mas em toda a sua extensão, elemento a

elemento. Para garantir essa simetria, apenas uma metade da malha é construída pelo

gerador e passada para o código de elementos finitos. A segunda metade é adicionada na

fase de pré-processamento e a relação de nós e elementos simétricos é guardada para as

análises posteriores. Essa relação é de grande importância para o método a ser apresentado

e possibilita a decomposição de qualquer campo, escalar ou vetorial, em componentes

26

simétricas e antissimétricas (sempre em relação ao plano médio do escoamento).

3.6 Código computacional

O código computacional que calcula as matrizes elementares de 3.6 e monta o problema

3.10 é baseado nas bibliotecas Getfem++ e Gmm++, de Pommier e Renard (2011), de-senvolvidas em linguagem C++ e distribuídas livremente para fins de pesquisa. Enquanto

Getfem++ é um biblioteca de elementos finitos flexível e eficiente que possibilita o cálculo

das integrais elementares da formulação fraca 3.6 e o gerenciamento de dados geométricos e de graus de liberdade do problema; Gmm++ é uma biblioteca de matrizes esparsas que

auxilia na criação das matrizes e no seu gerenciamento, incluindo funções eficientes para

as operações básicas.

A cada passo os sistemas 3.13 e 3.14 são montados com a ajuda das bibliotecas Getfem++ e Gmm++ e resolvidos com a utilização das rotinas do UMFPACK, que são

de grande eficiência e também distribuídas livremente para fins de pesquisa, ver Davis

(2002) e Davis e Duff (1999).

Por fim, os problemas característicos reduzidos descritos no capítulo 5 são resolvidos pelas funções da biblioteca computacional LAPACK, de Anderson et al. (1999).

27

4

DINÂMICA DO ESCOAMENTO