PGMEC · mostra, nos dois modelos, ... 4.11 Comparação entre erros obtidos para malhas de 501 pontos, ... resultado comparado com o obtido por Aktas [7]

PGMECPÓS-GRADUAÇÃO EM ENGENHARIA MECÂNICAESCOLA DE ENGENHARIAUNIVERSIDADE FEDERAL FLUMINENSE

Dissertação de Mestrado

MODELANDO PERDAS DE PRESSÃO

DEVIDO À IMPEDÂNCIA DE PAREDES

NA TRANSFERÊNCIA DE CALOR

TERMOACÚSTICA

HEITOR HERCULANO DE BARROS

DEZEMBRO DE 2014

HEITOR HERCULANO DE BARROS

MODELANDO PERDAS DE PRESSÃO DEVIDO ÀIMPEDÂNCIA DE PAREDES NATRANSFERÊNCIA DE CALOR

TERMOACÚSTICA

Dissertação apresentada ao Programa de Pós-graduação em Engenharia Mecânica da UFFcomo parte dos requisitos para a obtenção do tí-tulo de Mestre em Ciências em Engenharia Me-cânica

Orientador(es): Leonardo Santos de Brito Alves, Ph.D. (PGMEC/UFF)

UNIVERSIDADE FEDERAL FLUMINENSE

NITERÓI, DEZEMBRO DE 2014

MODELANDO PERDAS DE PRESSÃO DEVIDO ÀIMPEDÂNCIA DE PAREDES NATRANSFERÊNCIA DE CALOR

TERMOACÚSTICA

Esta dissertação foi julgada adequada para a obtenção do título de

MESTRE EM ENGENHARIA MECÂNICA

na área de concentração de Termociências, e aprovada em sua forma finalpela Banca Examinadora formada pelos membros abaixo:

Leonardo Santos de Brito Alves (Ph.D.)Universidade Federal Fluminense – PGMEC/UFF

(Orientador)

Leandro Alcoforado Sphaier (Ph.D.)Universidad Federal Fluminense – PGMEC/UFF

Daniel Rodriguez Alvarez (Ph.D.)Pontifícia Universidade Católica do Rio de Janeiro – PUC-Rio

Agradecimentos

Acima de tudo, agradeço à Deus, pela dádiva da vida, pela minha família e pelos

meus amigos.

Aos meus pais, José Maria e Maria Gildete, pelos ensinamentos e valores transmi-

tidos que me orientam a escolher entre o caminho certo e o errado e entre as coisas

fúteis e as essenciais.

À minha esposa Ellen Herculano, pela compreensão em relação ao tempo dispen-

sado para essa missão e pelo incentivo nos momentos mais difíceis.

Ao Professor Leonardo Alves, pela orientação precisa, competência e profissiona-

lismo.

Ao amigo e Pesquisador Renan Teixeira, pelo apoio fundamental na execução deste

trabalho.

Aos meus amigos do Laboratório de Mecânica Teórica e Aplicada (LMTA): Flávio,

Ricardo, Leandro, Rômulo, Jesus, entre outros.

Ao Programa de Pós-Graduação em Engenharia Mecânica da Universidade Federal

Fluminense, pela oportunidade de realizar o mestrado.

À Marinha do Brasil, pelo apoio e suporte recebido durante o período do mestrado.

iv

Resumo

Modelos hidrodinâmicos e termodinâmicos são utilizados para estudar o comporta-

mento de ondas termoacústicas dentro de uma cavidade unidimensional contendo um

gás perfeito ou um fluido supercrítico sob gravidade zero. Estas ondas termoacústicas

confinadas entre as paredes da cavidade se propagam e refletem diversas vezes indu-

zindo um rápido aquecimento do fluido fora da camada limite térmica, resultando num

aumento quase uniforme da temperatura do sistema. Este fenômeno de rápida rela-

xação da temperatura é conhecido como efeito pistão e tem sido extensivamente estu-

dado. O modelo hidrodinâmico utilizado consiste na solução das equações de Navier-

Stokes obtidas através do algoritmo de separação de fluxos com propriedade TVD para

cálculo dos termos invíscidos e o método Runge-Kutta SSP explícito para marchar no

tempo. Como o custo de simulação do modelo hidrodinâmico é extremamente ele-

vado, é necessário desenvolver um modelo termodinâmico. O modelo termodinâmico

original desenvolvido para a transferência de calor termoacústica é essencialmente a

equação do calor com um termo fonte proporcional a derivada temporal da temperatura

média. Este termo fonte modela o mecanismo de compressão adiabática responsável

pelo efeito pistão. Uma solução completamente analítica deste modelo é obtida através

da combinação da Técnica da Transformada Integral Generalizada e o Método da Ma-

triz Exponencial. Nenhum dos dois modelos apresenta resultados quantitativamente

semelhantes aos dados experimentais disponíveis, o que acreditamos ser causado pela

ausência de impedância nas paredes sólidas destes modelos. O objetivo principal deste

estudo é modelar as perdas de pressão nas paredes da cavidade durante a reflexão das

ondas termoacústicas devido a impedância acústica de parede. A inclusão deste efeito

mostra, nos dois modelos, que ele é realmente o responsável pela discordância na lite-

ratura entre resultados numéricos destes modelos e dados experimentais.

Palavras-chave: efeito pistão, ondas termoacústicas, fluidos compressíveis.

v

Abstract

Hydrodynamic and thermodynamic models have been utilized to study the behavior of

thermoacoustic waves inside a unidimensional cavity containing perfect gas or super-

critical fluid in zero gravity. These thermoacoustic waves entrapped between cavity

walls propagate and reflect several times inducing a rapid heating of the fluid beyond

the thermal boundary layer, resulting in a quasi-uniform increase of the system’s tem-

perature. The fast temperature relaxation phenomenon is known today as piston effect

and it has been extensively studied. The hydrodynamic model utilized consists of

the compressible unsteady Navier-Stokes equations solved through a flux-splitting al-

gorithm with a second order upwind method with TVD for convection term, second

order central difference conservative scheme for diffusion term and a second order

SSP Runge-Kutta for time marching. Once the computational cost of hydrodynamic

model is extremely high, it is necessary to develop a thermodynamic model. The ori-

ginal thermodynamic model developed for thermoacoustic heat transfer is essentially

the heat conduction equation with a source term proportional to the bulk temperature

time derivative. This source term models the adiabatic compression mechanism res-

ponsible for the piston effect. A fully analytical solution of this thermodynamic model

is obtained through a combination of the Generalized Integral Transform Technique

and the Matrix Exponential Method. None of the two models present quantity results

in according to available experimental results, which we believe occur due to the lack

of wall impedance in these models. The main purpose of this study is to model the

pressure losses that occurs when the thermoacoustic waves hit the walls by means of

a acoustic wall impedance. The inclusion of this effect, in both models, shows that it

is the responsible for the disagreement between the numerical results provided from

these models and the experimental data.

Keywords: piston effect, thermoacoustic waves, compressible fluids.

vi

Sumário

Agradecimentos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv

Resumo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v

Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii

Nomenclatura . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi

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

1.1 Motivação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.2 Revisão Bibliográfica . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

1.2.1 Estudos Teóricos . . . . . . . . . . . . . . . . . . . . . . . . . . 5

1.2.2 Estudos Experimentais . . . . . . . . . . . . . . . . . . . . . . . 11

1.2.3 Métodos Numéricos . . . . . . . . . . . . . . . . . . . . . . . . . 13

1.2.4 Técnica da Transformada Integral Generalizada . . . . . . . . . 19

1.3 Objetivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2. Modelos Matemáticos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

2.1 Modelo Hidrodinâmico . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

2.1.1 Equações Governantes . . . . . . . . . . . . . . . . . . . . . . . 23

2.1.2 Condições Iniciais e de Contorno . . . . . . . . . . . . . . . . . 25

2.1.3 Impedância Acústica . . . . . . . . . . . . . . . . . . . . . . . . 26

2.2 Modelo Termodinâmico . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

2.2.1 Equação do Calor com termo fonte de Compressão Adiabática 30

2.2.2 Pressão Termodinâmica . . . . . . . . . . . . . . . . . . . . . . 31


3. Métodos Numéricos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

3.1 Modelo Hidrodinâmico . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

3.1.1 Discretização Espacial . . . . . . . . . . . . . . . . . . . . . . . 35

vii

3.1.1.1 Fluxos Invíscidos . . . . . . . . . . . . . . . . . . . . 36

3.1.1.2 Fluxos Viscosos . . . . . . . . . . . . . . . . . . . . . 39

3.1.2 Discretização Temporal . . . . . . . . . . . . . . . . . . . . . . . 41


3.2 Modelo Termodinâmico . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

3.2.1 Sistema de Equações . . . . . . . . . . . . . . . . . . . . . . . . 43

3.2.2 Mudança para variáveis adimensionais e filtro . . . . . . . . . . 43

3.2.3 Transformação Integral . . . . . . . . . . . . . . . . . . . . . . . 45

4. Resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

4.1 Modelos Hidrodinâmicos . . . . . . . . . . . . . . . . . . . . . . . . . . 50

4.1.1 Equação de Burgers . . . . . . . . . . . . . . . . . . . . . . . . . 51

4.1.2 Tubo de Choque . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

4.1.3 Transferência de Calor Termoacústica . . . . . . . . . . . . . . 61

4.1.4 Efeito das condições de contorno para a pressão . . . . . . . . . 66

4.1.5 Efeito das condições de contorno para temperatura . . . . . . . 68

4.2 Modelos Termodinâmicos . . . . . . . . . . . . . . . . . . . . . . . . . . 71

4.2.1 Verificação Numérica . . . . . . . . . . . . . . . . . . . . . . . . 71

4.2.2 Análise do Efeito de Impedância das Paredes . . . . . . . . . . 72

4.3 Estimativa do valor da impedância acústica . . . . . . . . . . . . . . . . 78

5. Conclusões e Trabalhos Futuros . . . . . . . . . . . . . . . . . . . . . . . . 80

5.1 Conclusões . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

5.2 Trabalhos Futuros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

5.2.1 Evolução nas técnicas numéricas adotadas . . . . . . . . . . . . 81

5.2.2 Aspectos físicos e modelagem . . . . . . . . . . . . . . . . . . . 83

Referências . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

A. Computação Paralela com OpenMP Fortran API . . . . . . . . . . . . . . 90

A.1 Computação Paralela . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

viii

A.2 Breve Histórico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92

A.3 Filosofia Básica OpenMP . . . . . . . . . . . . . . . . . . . . . . . . . . 92

A.4 Paralelização do código Compressível . . . . . . . . . . . . . . . . . . . 93

A.5 Resultados de ganho com código paralelizado . . . . . . . . . . . . . . 95

B. Propriedades de Estabilidade Não-Linear . . . . . . . . . . . . . . . . . . 101

B.1 Monoticidade Preservada . . . . . . . . . . . . . . . . . . . . . . . . . . 101

B.2 TVD – Total Variation Diminishing . . . . . . . . . . . . . . . . . . . . 102

B.3 Amplitude Decrescente . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

B.4 Positividade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

B.5 Upwind Range Condition . . . . . . . . . . . . . . . . . . . . . . . . . . 107

B.6 Variação Total Limitada . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

B.7 Essencialmente Não-Oscilatória (ENO) . . . . . . . . . . . . . . . . . . 109

B.8 Resumo das condições de Estabilidade Não-Linear . . . . . . . . . . . 110

ix

Lista de Figuras

1.1 Representação esquemática da cavidade para estudo da transferência

de calor termoacústica. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.2 Refrigerador termoacústico [1]. . . . . . . . . . . . . . . . . . . . . . . . 4

1.3 Ondas termoacústicas para 0,25ta , 1,00ta , 1,50ta e 2,0ta [2]. . . . . . 7

1.4 Comparação entre as medições experimentais realizadas por Brown [3]

(círculos vazios) e os resultados obtidos numericamente por Huang e

Bau [4] (círculos cheios). . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.1 Impedância específica de contorno. . . . . . . . . . . . . . . . . . . . . . 27

2.2 Choque e reflexão da onda termoacústica. . . . . . . . . . . . . . . . . . 33

4.1 Deslocamento típico da descontinuidade na Equação de Burgers [5]. . 52

4.2 Comparação entre esquema de segunda ordem e solução exata da Equa-

ção de Burgers no tempo t = 2.0s. . . . . . . . . . . . . . . . . . . . . . 53

4.3 Comparação entre esquema SSPRK(2,2) e solução exata da Equação

de Burgers no tempo t = 2.0s. . . . . . . . . . . . . . . . . . . . . . . . . 53

4.4 Solução transiente produzida pelo RK(2,2) com passo no tempo ∆t =1.5∆x (SSP) e ∆t = 2.2∆x (Não-SSP) para o intervalo de 0 a 4s. . . . . 54

4.5 Representação esquemática do problema do tubo de choque e compor-

tamento das características ao longo do domínio [6] . . . . . . . . . . . 55

4.6 Comparação entre resultado numérico e solução analítica para veloci-

dade no tempo t = 10−2s. . . . . . . . . . . . . . . . . . . . . . . . . . . 58

4.7 Comparação entre resultado numérico e solução analítica para pressão

no tempo t = 10−2s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

4.8 Comparação entre resultado numérico e solução analítica para massa

específica no tempo t = 10−2s. . . . . . . . . . . . . . . . . . . . . . . . 58

4.9 Comparação entre resultado numérico e solução analítica para o nú-

mero de Mach no tempo t = 10−2s. . . . . . . . . . . . . . . . . . . . . . 59

x

4.10 Comparação entre resultado numérico e solução analítica para a velo-

cidade do som no tempo t = 10−2s. . . . . . . . . . . . . . . . . . . . . . 59

4.11 Comparação entre erros obtidos para malhas de 501 pontos, 1001 pon-

tos e 2001 pontos espaciais, problema do tubo de choque com tempo

físico final t f = 4×10−9s. . . . . . . . . . . . . . . . . . . . . . . . . . . 60

4.12 Ordem de erro espacial do código e ordens teóricas para comparação,

problema do tubo de choque com tempo físico final t f = 5×10−9s. . . 61

4.13 Comparação entre erros obtidos diferentes passos no tempo, problema

do tubo de choque com tempo físico final t f = 4×10−9s. . . . . . . . . 61

4.14 Comparação entre erros obtidos diferentes passos no tempo para o pro-

blema do tubo de choque com tempo físico final t f = 4×10−9s. . . . . 62

4.15 Comparação entre os resultados obtidos com o código e os mostrados

em Aktas [7], para quatro ondas termoacústicas nos tempos especifi-

cados na legenda. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

4.16 Comparação entre os resultados obtidos com o código e os mostra-

dos em Aktas [7] para temperatura ao longo da cavidade nos tempos

especificados na legenda. . . . . . . . . . . . . . . . . . . . . . . . . . . 63

4.17 Comparação entre os resultados obtidos com o código e os mostrados

em Aktas [7] para a massa específica ao longo da cavidade nos tempos


4.18 Resultados obtidos para a velocidade ao longo da cavidade nos tempos


4.19 Comportamento da pressão ao longo do tempo no centro da cavidade,

resultado comparado com o obtido por Aktas [7]. . . . . . . . . . . . . 66

4.20 Resultados da evolução da pressão ao longo do tempo com diferentes

valores de impedância acústica de parede Z . . . . . . . . . . . . . . . . 67

4.21 Tempo para atingir temperatura final na parede esquerda para os quatro

casos analisados. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

xi

4.22 Pressão na cavidade no tempo t = 2× 10−5s para diferentes tipos de

aquecimento. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

4.23 Ondas termoacústicas geradas pelo aquecimento com H = 0.020 e H =2.000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

4.24 Variação da pressão no ponto médio da cavidade para os quatro casos

analisados. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

4.25 Resultado do erro absoluto estimado P/P0 da solução com séries usando

N = 10 , 20, 30, 40 e 50. . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

4.26 Comportamento da pressão ao longo do tempo previsto pelo modelo

termodinâmico e o resultado experimental. . . . . . . . . . . . . . . . . 73

4.27 Resultados experimentais da pressão medidos por Brown [3] utilizando

nitrogênio e solução do código termodinâmico. . . . . . . . . . . . . . . 74

4.28 Resultados experimentais da pressão medidos por Brown [3] utilizando

ar padrão e solução do código termodinâmico. . . . . . . . . . . . . . . 75

4.29 Aquecimento da lâmina de níquel para tempo curto e longo, durante

experimento realizado por Lin e Farouk [8]. . . . . . . . . . . . . . . . . 76

4.30 Resultados experimentais da pressão medidos por Lin e Farouk [8]

utilizando ar padrão e solução do código termodinâmico. . . . . . . . . 77

4.31 Esquema e equipamentos que compõem um tubo de impedância [9]. . 78

A.1 Esquema de arquitetura com memória compartilhada e memória dis-

tribuída. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91

A.2 Sequência de execução de um código paralelizado com modelo Fork-

Join. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94

A.3 Speed-Up e Eficiência x Número de Processadores - Cluster 32 pro-

cessadores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

A.4 Speed-Up e Eficiência x Número de Processadores - Desktop . . . . . 97

B.1 Redução de um máximo global (a) e eliminação de um mínimo local

(b)[6]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

xii

B.2 Quadro esquemático das condições de Estabilidade Não-Linear [6]. . . 111

xiii

Lista de Tabelas

A.1 Tipo de execução por seção do código hidrodinâmico. . . . . . . . . . . 95

A.2 Resultado de Tempo de Execução e Speed-Up para o caso 1. . . . . . . 98




xiv

Nomenclatura

c velocidade do som

ρ densidade

T temperatura

Tc temperatura crítica

Tr temperatura reduzida do ponto crítico

u velocidade

cp calor específico a pressão constante

cv calor específico a volume constante

k condutividade térmica

e energia interna por unidade de massa

E energia total por unidade de massa

Ei vetor fluxo invíscido

Ev vetor fluxo viscoso

F filtro para temperatura adimensional

h entalpia por unidade de massa

H vetor termos fonte

k conditividade térmica

L comprimento da cavidade

M número de Mach

N número de termos do somatório da série que gera a solução

Ni norma

P pressão

PH pressão hidrodinâmica

PT pressão termodinâmica

xv

Q vetor de variáveis conservativas dependentes

Q vetor de variáveis primitivas dependentes

t tempo físico

ta tempo acústico

tD tempo de relaxamento térmico por difusão

tPE tempo de relaxamento do efeito pistão

x coordenada espacial

∆tE A passo do tempo do método de Euler avançado

Símbolos Gregos

α difusividade térmica

αP coeficiente de expansão térmica isobárica

βi autovalores

δi j delta de Kronecker

ηi coeficiente de transformação integral

γ razão entre calores específicos

κT compressibilidade isotérmica

µ viscosidade dinâmica

ψi autofunção

ψi autofunção normalizada

θ temperatura homogênea adimensional

θ vetor temperatura homogênea adimensional transformada

θi temperatura homogênea adimensional transformada

Θ temperatura adimensional

Θb temperatura adimensional global

ε coordenada espacial adimensional

υ condição de Courant-Friedrichs-Lewy (CFL)

xvi

Sobrescrito

0 estado inicial, parede direita

1 parede esquerda

b estado global

r reduzida do ponto crítico

c ponto crítico

D difusão térmica

H hidrodinâmico

i índice de somatório das séries

j índice de somatório das séries

P isobárico

PE efeito pistão

T isotérmico, termodinâmico

V isovolumétrico

E A relativo ao método de Euler avançado

xvii

Capítulo 1

Introdução

1.1 Motivação

O constante desenvolvimento dos computadores digitais de alta velocidade durante

o século XX tem causado uma grande mudança na forma de estudar e aplicar os princí-

pios da Mecânica dos Fluidos e de Transferência de Calor. Complexos problemas que

levariam anos de trabalho e análise em décadas anteriores, hoje podem ser avaliados

com um custo mínimo em poucos segundos de tempo computacional. Esta metodo-

logia, que consiste em resolver o sistema de equações diferenciais governantes com o

auxílio de um computador, quando aplicada à problemas de fenômenos de transporte,

é conhecida como Dinâmica dos Fluidos Computacional (denominação oriunda do idi-

oma inglês Computational Fluid Dynamics e normalmente abreviada por CFD) [10].

Tradicionalmente, métodos experimentais e teóricos eram usados para projetar equi-

pamentos. Hoje, embora experimentos continuem sendo extremamente importantes,

os mesmos são precedidos e direcionados a partir de resultados obtidos por meio de

simulações computacionais, principalmente quando envolvem geometrias complexas

e ambientes inóspitos [11].

Dentro deste contexto da Dinâmica dos Fluidos, os escoamentos compressíveis

exercem um fundamental papel, uma vez que praticamente todos os escoamentos de

interesse prático possuem algum grau de compressibilidade, e por vezes, esta compres-

1

sibilidade tem uma grande influência no seu comportamento. Simulações da circula-

ção de sangue no corpo humano, escoamento ao redor de uma aeronave e movimento

dos gases de combustão pelos dutos de saída de um carro de corrida são exemplos de

aplicações práticas [5].

Uma exemplo de escoamento compressível de interesse em muitos fenômenos na-

turais e aplicações industriais é a interação dinâmica fluido-acústica. Compressores

acústicos, levitadores acústicos e refrigerantes termoacústicos são exemplos de apli-

cações. Estas aplicações envolvem a formação de campos acústicos em domínios de

interesse, como uma cavidade, um sistema de armazenamento ou qualquer meio con-

finado. Embora a geração de ondas acústicas em fluidos possa ocorrer em decorrência

de efeitos térmicos, mecânicos ou como conseqüência do próprio escoamento, apenas

as ondas termoacústicas serão estudadas neste trabalho. Desta forma, este trabalho

foca especificamente nas ondas acústicas geradas por uma mudança repentina de tem-

peratura na fronteira do sistema.

Quando um gás, confinado em uma cavidade é submetido a um rápido aumento

Fig. 1.1: Representação esquemática da cavidade para estudo da transferência de calortermoacústica.

2

de temperatura na sua fronteira sólida, o volume fluido imediatamente vizinho é aque-

cido por condução e tende a expandir. Contudo, a expansão repentina do gás devido

a adição de energia é freada pela inércia do gás não perturbado e ocorre a indução de

uma onda de pressão denominada onda termoacústica [7]. O fenômeno físico descrito

e esquematizado na Figura 1.1 é denominado de transferência de calor termoacústica

devido ao fato de ser induzido termicamente, de consistir numa onda de pressão se pro-

pagando na velocidade acústica do fluido e de ocorrer transferência de calor durante a

propagação da região localmente comprimida. As ondas termoacústicas se manifestam

por um aumento de pressão que se propaga através do gás, acompanhada por perturba-

ções de temperatura, massa específica e velocidade. Ao atingir a extremidade oposta

ao aquecimento, esta onda termoacústica é refletida gerando um movimento cíclico

contínuo conhecido como efeito pistão [3]. Este fenômeno, inicialmente observado

experimentalmente, exerce influência no processo convectivo de gases industriais de

grande aplicação, como nitrogênio e hélio, mas se torna preponderante em fluidos em

estado supercrítico, que também são de grande interesse e possuem modernas aplica-

ções nos dias de hoje [12].

Interações termoacústicas são encontradas em muitos fenômenos naturais e aplica-

ções industriais, como motores e refrigeradores termoacústicos. Essas e outras apli-

cações envolvem campos acústicos gerados mecanicamente em cavidades preenchidas

com fluidos diversos. Uma onda sonora em um fluido compressível é normalmente

vista como um binário consistente entre oscilações de velocidade e pressão, onde os-

cilações de temperatura também estão presentes. Por exemplo, quando ondas sonoras

se deslocam ao longo de microcanais, transferência de calor oscilatória ocorre nas pa-

redes do canal. No caso de ondas sonoras intensas em gases pressurizadaos, a energia

acústica pode ser canalizada para produzir motores de potência, combustão pulsante,

bombas de calor, refrigeradores e separadores de mistura [8].

No caso de motores e refrigeradores termoacústicos, ocorre a interação entre uma

onda termoacústica contida em um tubo de ressonância e uma pilha. Um motor termo-

acústico absorve calor a alta temperatura e rejeita calor a baixa temperatura enquanto

3

produz energia acústica. Um refrigerador termoacústico, por sua vez, funciona de

maneira oposta, absorvendo calor a baixa temperatura e rejeitando calor a alta tempe-

ratura através do consumo de potência acústica [1]. O diagrama esquemático mostrado

na figura 1.2 exemplifica um modelo de um refrigerador termoacústico, que consiste

basicamente de um tubo ressonante reto, uma pilha de placas paralelas, trocadores de

calor e um transdutor de origem acústica (um auto-falante, por exemplo).

Fig. 1.2: Refrigerador termoacústico [1].

1.2 Revisão Bibliográfica

Nesta seção, os conceitos relevantes e últimos avanços existentes na literatura a

respeito dos fenômenos físicos e modelagem da transferência de calor termoacústica

em gases perfeitos são relacionados. O conhecimento atual a respeito de ondas termo-

acústicas em gases perfeitos e fluidos supercríticos são relacionados inicialmente. A

seção continua com o enfoque nos conhecimentos relativos ao conceito de impedância

acústica e com a revisão da bibliografia existente a respeito dos métodos numéricos

utilizados durante o desenvolvimento do trabalho de pesquisa. A seção é finalizada re-

latando o desenvolvimento da Técnica da Transformada Integral Generalizada e aplica-

ções desta técnica em problemas de natureza próxima aos fenômenos aqui estudados.

4

1.2.1 Estudos Teóricos

O problema de ondas de pressão induzidas termicamente, denominadas como on-

das termoacústicas, em uma massa finita de gás perfeito quiescente sujeita a uma mu-

dança de temperatura em sua parede sólida foi analiticamente estudada por Trilling

[13]. Nesse trabalho foi determinado como a intensidade do som depende da variação

de temperatura na parede. As equações de escoamento compressível unidimensional

foram linearizadas e obtidas utilizando a técnica da transformada de Laplace. Além

disso, um modelo simplificado, que consiste basicamente de uma equação hiperbólica

para a condução modelando o movimento das ondas termoacústicas, foi comparado

com as equações de Navier-Stokes unidimensionais e as limitações estudadas.

Larkin [14] realizou o primeiro estudo numérico do movimento de fluido induzido

termicamente em um meio como resultado de um brusco aquecimento. Nessa aná-

lise unidimensional, a natureza acústica do fluido em movimento foi observado, mas

devido a uma alta difusão numérica, o comportamento das ondas de pressão não foi

completamente capturado. Ozoe et al. [15] também realizaram estudos numéricos

de ondas termoacústicas unidimensionais e bidimensionais em uma região confinada.

Esses estudos computacionais descreveram soluções de diferenças finitas da equações

compressíveis de Navier-Stokes para um gás com propriedades termo-físicas indepen-

dentes da temperatura. As soluções foram obtidas empregando esquemas upwind de

primeira ordem para resolver as equações governantes e, como consequência, todos os

resultados também demonstraram substanciais efeitos de difusão numérica.

Huang e Bau [4, 16, 17] obtiveram uma classe geral de soluções para o problema

de ondas termoacústicas usando o método da transformada de Laplace com inversão

numérica para equações modeladas como uma onda linear. Para mudanças bruscas

e graduais de temperatura na fronteira, estudando a geração e transmissão de ondas

termoacústicas planas para um gás com Pr = 3/4, campos assintóticos foram deriva-

dos para pressão, velocidade, temperatura e para o fluxo de calor na parede. Nesses

trabalhos, as equações não lineares foram resolvidas numericamente e os efeitos da

não-linearidade foram avaliados na natureza e características das ondas termoacústicas

5

geradas. O estudo também se concentrou em determinar as condições sobre as quais a

aproximação linear é adequada. Os autores concluíram que o método da transformada

de Laplace possui a vantagem de ser livre de problemas de natureza numérica como

dissipação artificial e resolução insuficiente da malha, resultando numa boa concor-

dância entre os resultados obtidos numericamente e os existentes experimentais para

tempos curtos. Por outro lado, uma comparação entre as soluções lineares e não-

lineares mostrou que para muitas aplicações práticas a aproximação linear é adequada,

e as soluções não lineares requerem uma enorme quantidade de tempo computacional

para ser produzida e, por isso, é restrita a análises de tempo muito curtos.

Brown e Churchill [18] investigaram ondas termoacústicas geradas por aqueci-

mento rápido. Os resultados foram computados para passos discretos de aquecimento

da parede e para uma taxa exponencial de aquecimento. As previsões numéricas tive-

ram uma boa concordância qualitativa com medidas experimentais em termos de forma

da onda, amplitude e taxa de decaimento da onda em termos de pressão. Comparações

quantitativas não foram possíveis, uma vez que as dimensões reais da cavidade em

que se realizou o experimento tornaram os tempos de simulação extremamente longos

inviabilizando uma simulação computacional correspondente.

Farouk et al. [2] estudaram numericamente o comportamento das ondas termo-

acústicas em uma cavidade bidimensional e investigaram como estas ondas podem

ser utilizadas como um efetivo mecanismo de remoção de calor. As características

das ondas termoacústicas bidimensionais foram investigadas tanto para aquecimento

quanto para resfriamento instantâneo e não-uniforme das paredes. Quando as paredes

verticais foram aquecidas instantaneamente ou gradualmente, as ondas induziram es-

coamentos dentro da cavidade. A figura 1.3 exemplifica quatro ondas termoacústicas

se deslocando ao longo da cavidade em quatro tempos específicos escritos em função

do tempo acústico (ta = L/c).

Aktas e Farouk [19] estudaram o problema termoacústico considerando proprie-

dades variáveis em função da temperatura numa cavidade bidimensional com aqueci-

mento uniforme e não-uniforme. Foram realizadas comparações entre os resultados

6

Fig. 1.3: Ondas termoacústicas para 0,25ta , 1,00ta , 1,50ta e 2,0ta [2].

obtidos em modelos unidimensionais e bidimensionais. Neste trabalho também foi ve-

rificada a interação entre ondas termoacústicas bidimensionais e escoamentos dentro

da cavidade uniformemente ou não-uniformemente aquecida em uma das suas frontei-

ras.

Zappoli et al. [20] investigaram numericamente o comportamento de uma cavi-

dade contendo gás carbônico em estado supercrítico, sujeito à gravidade zero e a um

aquecimento na sua fronteira. Neste trabalho, foi utilizado um modelo hidrodinâmico

baseado nas equações de Navier-Stokes na sua forma completa, não-linear e unidi-

mensional. Utilizando escalas de tempos curtos (acústicos) e longos (difusivos) foram

realizados experimentos numéricos, onde se verificou o fenômeno conhecido como

efeito pistão, ocorrendo a transformação de energia térmica em energia cinética em

uma camada limite em expansão. Ao invés dos processos difusivos se tornarem mais

lentos, comportamento esperado para fluidos supercríticos, foi observada uma rápida

homogeneização da temperatura em apenas 1% do tempo de relaxamento difusivo.

Nakano [21] estudou teoricamente e numericamente a transferência de calor e

massa de uma mistura binária de fluidos supercríticos ao longo da linha pseudo-crítica.

Neste estudo, utilizando uma mistura de ar supercrítico artificial com composição de

79% de nitrogênio e 21% de oxigênio, foi investigado o efeito pistão, o efeito de ter-

moforese e a interação entre eles. Um modelo hidrodinâmico derivado das equações

7

da dinâmica termo-fluida foi utilizado, no qual a compressibilidade do fluido, a tempe-

ratura, a pressão e a dependência da entropia em função das concentrações da mistura

são levadas em consideração. As equações governantes foram resolvidas através do

Método das Diferenças Finitas. O autor concluiu que o efeito pistão foi o responsá-

vel pela propagação da energia térmica ao longo do domínio e que as mudanças de

concentração ocorreram devido ao efeito de termoforese. Neste trabalho, também foi

possível estimar a razão de difusão térmica, que torna possível a correlação direta entre

o gradiente de temperatura e o gradiente de concentração.

Zhang e Shen [22] utilizaram o modelo hidrodinâmico para estudar a propagação e

reflexão de ondas termoacústicas próximas ao ponto crítico líquido-gás. Foram reali-

zadas investigações numéricas da resposta acústica à perturbações térmicas de acordo

com os experimentos realizados por Miura et al. [23] com gás carbônico supercrítico.

Os resultados obtidos mostram uma boa concordância com os dados experimentais

para tempos curtos. Diferentes características do efeito pistão foram analisadas neste

estudo. As análises revelaram as diferentes respostas do sistema ao aquecimento con-

tínuo ou por impulso em fluidos em estado próximo ao ponto crítico e o papel da

divergência da viscosidade média num processo identificado como efeito pistão secun-

dário.

Zhang e Shen [24], agora utilizando nitrogênio supercrítico, conduziram experi-

mentos numéricos para avaliar o comportamento termomecânico em resposta as dife-

rentes formas de perturbação térmica como aquecimento instantâneo ou gradual. Vá-

rios aspectos como geração, propagação, padrões de reflexão e evolução no tempo das

ondas termoacústicas foram investigados através do modelo hidrodinâmico baseado

nas equações de Navier-Stokes. Os resultados obtidos demonstraram que, essencial-

mente, dois tipos de ondas são geradas ao ajustar-se a taxa de aquecimento da fron-

teira. Impondo um aumento de temperatura extremamente rápido na fronteira, um

pulso acústico é emitido no fluido, no qual a energia é rapidamente dissipada devido

à atenuação viscosa e numerosas reflexões de fronteira. O aquecimento gradual, por

outro lado, emite ondas compressíveis que carregam energia e momento da fronteira

8

aquecida para o interior da cavidade sofrendo menor dissipação viscosa ao longo do

domínio. Foi demonstrado também que a transição do processo acústico para o pro-

cesso difusivo (mais lento) é também fortemente influenciado pela taxa de variação da

temperatura na fronteira.

Zhang e Shen [25] investigaram as emissões acústicas para diferentes temperaturas

próximas da linha crítica isocórica sob perturbações lineares e não-lineares de tempe-

ratura. Foram analisadas as respostas acústicas resolvendo as equações governantes

do modelo hidrodinâmico. Foi demonstrado que um padrão de onda termoacústica

homogênea domina o caso linear, largamente independente da temperatura reduzida.

Enquanto que, sob perturbação nao-linear, uma variação mínima em torno da tempe-

ratura reduzida pode levar a severas deformações da frente de onda.

Onuki e Ferrel [26] aplicaram um modelo termodinâmico baseado na equação de

calor com termo fonte de compressão adiabática para modelar o rápido equilíbrio adi-

abático observado experimentalmente em um fluido próximo ao ponto crítico líquido-

vapor. Eles observaram que a mudança de entropia ao longo do fluido leva a uma

mudança de pressão e, consequentemente, a temperatura adiabática muda ao longo

da massa fluida dentro da cavidade. Também foi observada que esta correlação en-

tre temperatura e pressão, após um curto intervalo de tempo, leva a uma mais rápida

homogeneização da temperatura que a resultante do processo de difusão pura. Dessa

forma, foi exposto o modelo termodinâmico que negligenciava o efeito da gravidade e

o escoamento dentro da cavidade, porém foi consistente com resultados experimentais

para temperatura ao longo do tempo.

Boukari et al. [27] mostraram como a compressibilidade extrema de fluidos puros

em estado próximo ao ponto crítico afetam significativamente a resposta dinâmica da

temperatura do sistema através de processos adiabáticos. Os autores utilizaram o mo-

delo termodinâmico para descrever o problema na ausência da gravidade e quantificar

os efeitos correlatos da temperatura e pressão através de soluções numéricas unidimen-

sionais. Os resultados obtidos neste trabalho utilizando xenon crítico aquecido entre

10 e 20mK acima do ponto crítico previram uma relaxação da temperatura de 99% em

9

menos de 5s. Outra conclusão do trabalho foi que o resfriamento adiabático é mais

rápido quando o fluido está cada vez mais próximo do ponto crítico.

Carlès e Dadzie [28] investigaram o efeito pistão em fluidos próximos ao ponto crí-

tico utilizando o modelo termodinâmico. O efeito pistão foi descrito como um binário

termomecânico acelerando a relaxação térmica próximo ao ponto crítico, sendo for-

temente influenciado pela divergência crítica da viscosidade média. Foram realizadas

análises teóricas próximas ao ponto crítico, sugerindo que o efeito pistão é governado

por dois tempos característicos que possuem comportamentos distintos à medida que o

estado termodinâmico se aproxima do ponto crítico. Enquanto o tempo característico

clássico do efeito pistão tende a zero no ponto crítico, o segundo tempo característico

tente a infinito. Os autores atribuíram o tempo característico clássico como responsável

pela relaxação da temperatura e o segundo tempo característico como responsável pelo

enfraquecimento do efeito pistão, consequência da divergência na viscosidade média

próximo ao ponto crítico. Por fim, é proposta uma estratégia experimental alternativa

para cálculo da viscosidade média no ponto crítico.

Alves [29] considerou a transferência de calor transiente em um fluido compres-

sível sob microgravidade próximo ao ponto crítico. A partir do modelo termodinâ-

mico, assumindo as propriedades dos fluidos como função da temperatura e pressão,

foi capturado o efeito pistão não-linear. O autor demonstrou que as condições iniciais

representam o fator mais importante no processo, uma vez que a variação das propri-

edades tem efeito menos relevante na transferência de calor termoacústica. Foi mos-

trado ainda, que a hipótese de propriedades constantes pode ser utilizada para modelar

o efeito pistão sob uma ampla faixa de temperaturas e pressões.

Nikolayev et al. [30] investigaram a transferência de calor em um fluido puro su-

jeito a uma temperatura ligeiramente acima da temperatura crítica. Os resultado do

modelo termodinâmico foram comparados com resultados de simulações diretas do

modelo hidrodinâmico para o gás carbônico e hexafluoreto de enxofre. A equação

de estado de Van der Waals foi utilizada para descrever os dois fluidos. O modelo

termodinâmico ofereceu resultados em concordância com o modelo hidrodinâmico,

10

provendo ainda um excelente ganho em termos de tempo computacional. O modelo

termodinâmico previu a evolução completa do processo de transferência de calor em

minutos, enquanto que as simulações do modelo hidrodinâmico levaram semanas de

tempo computacional. O modelo termodinâmico não requer a avaliação das variáveis

em cada célula do domínio computacional, garantindo uma melhor eficiência em re-

lação ao modelo hidrodinâmico. Esta vantagem de custo computacional possibilita a

utilização deste modelo para simulações bidimensionais ou tridimensionais.

Teixeira e Alves [31] realizaram um estudo comparativo utilizando o modelo hi-

drodinâmico e termodinâmico do problema de transferência de calor em uma cavidade

contendo fluido supercrítico à gravidade zero. Um método com pré-condicionamento

para baixo Mach foi empregado para gerar a solução numérica do modelo hidrodi-

nâmico, evitando a necessidade de resolver escalas de tempo acústicas. Um modelo

apropriado da evolução do efeito pistão no pseudo-tempo é incluindo para gerar uma

descrição fisicamente correta da evolução no tempo físico. A solução analítica do mo-

delo termodinâmico foi obtida através da Técnica da Transformada Integral Generali-

zada e o Método da Matriz Exponencial. Foi demonstrado que ambos modelos geram

resultados graficamente idênticos, mas apenas se as propriedades dos fluidos impostas

e a equação de estado satisfazem a relação generalizada de Mayer e a definição de

velocidade do som.

1.2.2 Estudos Experimentais

A geração de ondas termoacústicas em gases foi estudada experimentalmente por

poucos investigadores. Parang e Salah-Eddine [32] investigaram o fenômeno de con-

vecção termoacústica em um cilindro contendo ar em condições normais e em ambi-

entes com gravidade reduzida. Em suas medições, a temperatura do ar foi registrada

devido a pequenas oscilações de amplitude, mesmo registrando o experimento com um

sistema de baixa taxa de medições, a temperatura do ar aumentou mais rapidamente

que nos resultados numéricos para condução pura. Nesse estudo, nenhuma medição

de pressão foi reportada pelos autores.

11

Medidas experimentais da pressão geradas por aquecimento rápido foram registra-

das pela primeira vez por Brown [3]. Em seus experimentos, o processo de aqueci-

mento instantâneo foi obtido através de um circuito elétrico resistor-capacitor (RC). O

experimento foi realizado com o gás nitrogênio, hélio e ar padrão. As medidas de pres-

são obtidas mostram claramente a geração de ondas termoacústicas por aquecimento

instantâneo de parede previsto em teoria.

Lin e Farouk [8] também realizaram estudos experimentais utilizando um circuito

elétrico RC conforme descrito por Brown [3], porém utilizando um cilindro de me-

nor dimensão e somente para o ar padrão. Neste trabalho foi realizado uma análise

comparativa entre as medidas experimentais e os resultados numéricos para a evolução

da pressão ao longo do tempo. Para tempos curtos (até 0.002 segundos) os resultados

numéricos tiveram uma boa concordância, porém, as previsões numéricas não mais

concordam com as medidas experimentais para tempos mais longos. Os autores le-

vantaram a hipótese da divergência ter ocorrido devido à presença de vazamentos de

gás na câmara e de amortecimento na parede, que foi considerada rígida durante o

experimento.

Fig. 1.4: Comparação entre as medições experimentais realizadas por Brown [3] (cír-culos vazios) e os resultados obtidos numericamente por Huang e Bau [4](círculos cheios).

Huang e Bau [4] compararam seus resultados numéricos com as medições experi-

12

mentais de pressão ao longo do tempo obtidas por Brown [3]. Também foram obser-

vadas divergências para a evolução da pressão e amplitude da onda termoacústica para

tempos longos, conforme mostrado na figura 1.4, extraída do próprio artigo, onde os

círculos preenchidos representam os resultados numéricos e os círculos vazios repre-

sentam resultados experimentais. Os autores sugeriram que as diferenças para tempos

longos são consequência de falta de sensibilidade do microfone utilizado, vazamento

de gás devido ao aparato não está perfeitamente selado, fronteira das paredes sólidas

não serem rígidas, do fenômeno não ser estritamente unidimensional e dos efeitos ad-

versos de não-linearidade serem mais significantes após algumas reflexões da onda.

Nakano e Shiraishi [33] investigaram como a combinação entre alta compressibi-

lidade térmica e baixa difusão térmica afetam a propagação da energia térmica. Neste

experimento, o efeito pistão no nitrogênio supercrítico foi investigado. Uma vez que

a convecção natural devido à gravidade sob condições terrestres interfere no efeito

pistão, no intuito de suprimir este efeito, foi adicionado calor a partir do topo da cé-

lula experimental. Desta forma, conseguiu-se obter o perfil de temperatura esperado

resultante do efeito pistão próximo à linha pseudo-crítica.

Miura et al. [23] investigaram experimentalmente mudanças adiabáticas no gás

carbônico próximo ao ponto crítico em escalas de tempo acústico utilizando interfe-

rometria ultrasensível. Pequenos fluxos de calor aplicados na fronteira produziram

pulsos acústicos com amplitude da ordem de 10µs, que se propagam ao longo da cé-

lula interagindo com as paredes. Foi observado que o pulso se intensificava à medida

que o estado termodinâmico se aproximava do ponto crítico. Foi estudada também a

emissão e propagação das ondas acústicas no fluido a partir de diferentes aquecimentos

e como o calor aplicado é transformado em trabalho mecânico.

1.2.3 Métodos Numéricos

Como observado anteriormente nas figuras 1.3 e 1.4, a solução numérica do modelo

hidrodinâmico requer a utilização de esquemas capazes de resolver de maneira acurada

os gradientes severos encontrados tanto na camada limite térmica quanto na frente

13

de onda acústica. Estes métodos são discutidos em detalhes a seguir, primeiro para

resolução especial e, depois, para a temporal.

Lax [34] mostrou que para Leis de Conservação, a variação total, denotada por TV,

das soluções fisicamente possíveis não aumentam com o tempo. Onde TV é dada por:

T V =∫ L

0

∣∣∣∂u

∂x

∣∣∣d x (1.1)

e TV para um caso discreto, temos:

T V =∑j

∣∣∣u j+1 −u j

∣∣∣ (1.2)

Um método numérico é dito possuir Variação Total Decrescente, ou é dito possuir a

propriedade de estabilidade TVD, se:

T V (un+1) ≤ T V (un) (1.3)

Harten [35] provou que:

1. Um esquema monotônico é TVD, ou seja, sua Variação Total diminui ao longo

do tempo;

2. Um esquema TVD possui Monoticidade Preservada.

Então, se esquemas TVD de alta ordem podem ser construídos, estes esquemas te-

rão monoticidade preservada. Uma discussão mais elaborada se encontra no apêndice.

Laney [6] afirma que, no estudo de Leis de Conservação Escalar em domínios espa-

ciais infinitos, dizer que as soluções preservam a propriedade de monoticidade significa

dizer que, se as condições iniciais u(x,0) aumentam a monoticidade, a solução u(x, t )

aumenta a monoticidade para qualquer tempo t , e se as condições iniciais diminuem

a monoticidade, a solução terá a monoticidade decrescente para qualquer tempo t . A

propriedade chamada Monoticidade Preservada é uma consequência da preservação da

forma de onda, da destruição da forma de onda e da criação da forma da onda.

Segundo Pletcher et al. [11], a idéia central ao construir um esquema TVD é

14

tentar desenvolver um método de alta ordem que evite oscilações e exiba propriedades

similares aos esquemas monótonos. Para tais esquemas, a solução é de primeira ordem

próximo a descontinuidades e alta ordem em regiões suaves. A transição para alta

ordem é feita com o uso de limitadores de inclinação das variáveis dependentes ou

limitadores de fluxo.

Boris e Book [36], durante o desenvolvimento do algoritmo Flux-corrected Trans-

port, estudaram formas de manter a estabilidade em regiões de choque, evitar oscila-

ções numéricas espúrias e evitar a criação de novos máximos e mínimos ao longo da

solução. Como estes objetivos, os autores desenvolveram o primeiro limitador não-

linear adicionando uma diferença limitada entre fluxos de primeira e segunda ordem

para prevenir oscilações associadas a esquemas de segunda ordem.

O trabalho de Van Leer [37] trata do desenvolvimento de esquemas de diferenças

monotônicas de segunda ordem de acurácia através da inclusão de termos não-lineares

que avaliam o comportamento da solução na região em análise. A metodologia baseia-

se em limitadores da extrapolação das variáveis dependentes para células na fronteira

no intuito de prevenir saltos na solução.

Roe [38] propôs um esquema de diferenças utilizando limitadores de fluxo que ofe-

rece perfis agudos na solução da equação da advecção linear. O trabalho apresentou

a unificação de diversos esquemas de segunda ordem TVD que haviam sido propos-

tos independentemente por diversos autores. O trabalho permitiu visualizar como os

diversos esquemas estão relacionados entre si.

Sweby [39] obteve esquemas de diferenças explícitas com alta resolução e livre

de oscilações (TVD) adicionando limitadores de fluxos anti-difusivos ao esquema de

primeira ordem. Sweby também investigou as fronteiras de aplicação dos limitado-

res desenvolvidos e comparou estes limites teoricamente e numericamente com outros

limitadores existentes.

Tanto Sweby [39], quanto Roe [38] desenvolveram limitadores baseados em mode-

los TVD do esquema Lax-Wendroff, enquanto Harten [40] aplicou o conceito para um

fluxo modificado buscando o mesmo efeito dos seus antecessores. Harten apresentou

15

um novo método numérico para cálculo de soluções descontínuas de leis de conser-

vação hiperbólicas. O método de fluxo modificado de Harten consiste num primeiro

estágio, onde um esquema híbrido de diferenças finitas é utilizado juntamente com um

método de primeiro ordem não-oscilatório para prover uma variação monotônica da

solução próximo às descontinuidades, e um segundo estágio, onde uma compressão

artificial é aplicada às transições agudas nas descontinuidades.

É importante compreender que apenas métodos de diferenças centradas ou atrasa-

das podem ser feitos TVD com uma adequada modificação usando limitadores. Limi-

tadores podem ser escritos em termos de variações sucessivas de variáveis dependentes

ou de fluxos. Porém, por simplicidade, eles normalmente são escritos em função de

uma razão das variáveis calculada no ponto em análise [11].

Laney [6] resumiu as propriedades do TVD como se segue:

• A propriedade TVD está intimamente ligada a positividade, que é uma forte

condição de estabilidade não-linear;

• TVD implica em monoticidade preservada. Isto é desejável em circunstâncias

onde a preservação de monoticidade é muito fraca mas indesejável onde a pre-

servação de monoticidade é muito forte, uma vez que o TVD pode deixá-la ainda

mais forte.

• TVD tende a causar erros de clipping nos extremos, ou seja, alguns extremos são

aparados para que sua amplitude não cresça e cause saltos. Na prática, a maio-

ria dos métodos TVD cortam os extremos localizados entre primeira e segunda

ordem de acurácia.

• Em teoria, métodos TVD podem permitir oscilações, porém raramente acontece

na prática. TVD provê um limite superior nas oscilações, eliminando, no pior

caso, instabilidade com crescimento ilimitado.

• Apesar do incoviniente de aparar alguns extremos, TVD é um elemento pratica-

mente indispensável na teoria de estabilidade não-linear moderna e está presente

na maioria dos códigos atuais.

16

Estas propriedades não-lineares para discretização espacial são válidas para es-

quemas monotônicos e, desta forma, mantidas quando a marcha no tempo utiliza um

esquema de Euler, seja ele implícito ou explícito. Contudo, este método possui ape-

nas primeira ordem no tempo. Sua extensão para ordens maiores no tempo usando

esquemas tradicionais fazem com que estas propriedades sejam perdidas.

Discretizações temporais de alta ordem com Forte Preservação de Estabilidade

(termo vindo do inglês Strong Stability Preserving - SSP) foram desenvolvidas para

garantir a manutenção destas propriedades de estabilidade não-linear, necessárias na

solução numérica de equações diferenciais parciais com soluções possuindo fortes gra-

dientes ou descontinuidades. Métodos de multi-passo e multi-estágio SSP explícitos

tem sido empregados em conjunto com uma grande variedade de discretizações espa-

ciais, incluindo métodos Galerkin, métodos ENO (Essentially Non-Oscillatory), méto-

dos WENO (Weighted Essentially Non-Oscillatory) e métodos de diferenças espectrais

[41].

A ideia dos métodos SSP é assumir que a discretização temporal de Euler avançada

de primeira ordem utilizando o método de linhas é fortemente estável apenas sob uma

certa norma, quando o passo de tempo ∆t é adequadamente restrito. Então podemos

encontrar uma discretização temporal de alta ordem que mantenha esta forte estabili-

dade para a mesma norma, através de uma restrição imposta ao passo no tempo. A

classe de métodos de discretização temporal SSP de alta ordem para métodos de linhas

aproximados de equações diferenciais parciais foram desenvolvidos por Shu e Osher

[42, 43] e inicialmente denominadas Discretizações Temporais TVD. Esses métodos

preservam as propriedades de estabilidade de um Euler avançado para qualquer norma

ou semi-norma. Na verdade, uma vez que os argumentos de estabilidade são basea-

dos em decomposições convexas de métodos de alta ordem em termos do método de

Euler de primeira ordem, qualquer função convexa será preservada por discretizações

temporais de alta ordem [44].

Shu e Osher [43] consideraram o método SSP Runge Kutta para solução da equação

17

hiperbólica

ut + f (u)x = 0 (1.4)

onde a derivada espacial, f (u)x , é discretizada por uma diferença finita TVD descrita

como −L(u). No processo de discretização, a malha espacial de pontos é representada

por x j e a malha de pontos temporais por t n . A discretização espacial L(u) possui a

propriedade de, quando combinada com um Euler avançado de primeira ordem,

un+1 = un +∆tL(un) (1.5)

garante que a Variação Total (TV) de uma solução discreta un não aumenta ao longo

do tempo, isto é, a propriedade Variação Total Decrescente (TVD) se verifica:

T V (un+1) ≤ T V (un), T V (un) =∑j|un

j+1 −unj | (1.6)

Para manter as propriedades SSP, enquanto buscamos alta ordem de acurácia no tempo,

alteramos a restição para a condição CFL:

∆t ≤ υ∆tE A (1.7)

onde υ é uma constante que depende do esquema de marcha escolhido.

Evidências numéricas apresentadas por Gottlieb e Shu [45] demonstram que osci-

lações podem ocorrer quando usamos estabilidade linear em métodos de alta ordem

sem SSP, ainda que a discretização espacial seja TVD combinada com uma discretiza-

ção temporal Euler avançada de primeira ordem. Isto ilustra que é mais seguro utilizar

uma discretização temporal SSP para resolver problemas hiperbólicos, além disso, mé-

todos SSP tem uma garantia extra de estabilidade e, em muitos casos, não aumentam

o custo computacional.

Shu e Oshar [43] apresentaram cada estágio do método Runge-Kutta explícito na

forma das equações

u0 = un (1.8)

18

ui =i−1∑j=0

(αi , j u( j ) +∆tβi , j F (u( j ))) (1.9)

un+1 = u(s) (1.10)

para métodos Runge-Kutta, o requisito de consistência é dado por∑i−1

k=0αi ,k = 1. Se

todos os coeficientes são negativos, cada estágio do método Runge-Kutta pode ser

rearranjado em combinações convexas de passos de Euler avançado

||u(i )|| =∣∣∣∣∣∣ i−1∑

j=0

(αi , j u( j ) +∆tβi , j F (u( j ))

)∣∣∣∣∣∣≤ i−1∑j=0

αi , j

∣∣∣∣∣∣u( j ) +∆tβi , j

αi , jF (u( j ))

∣∣∣∣∣∣ (1.11)

Esta observação motivou o seguinte teorema proposto por Shu e Osher [43]:

"Se o método de Euler avançado aplicado a (1.4) é SSP sob a restrição de passo

no tempo ∆t ≤∆tE A, isto é, a equação (1.5) é valida, se αi , j ,βi , j ≥ 0, então a solução

obtida para o método Runge-Kutta satifaz a barreira SSP sob a condição de passo no

tempo descrita na equação (1.12)."

∆t ≥ υ(α,β)∆tE A onde υ(α,β) = mi ni , jβi , j

αi , j(1.12)

A partir do teorema é possível notar que o passo capaz de preservar a condição SSP

é um produto de dois fatores: do passo no tempo do Euler avançado ∆tE A, o qual de-

pende apenas da discretização espacial, e do coeficiente υ(α,β), o qual depende exclu-

sivamente da discretização temporal. O coeficiente SSP ótimo para todas as possíveis

representações dos métodos Runge-Kutta de segunda ordem de dois estágios explícito

na formulação de Shu-Osher é υ= 1 [45].

1.2.4 Técnica da Transformada Integral Generalizada

A Técnica da Transformada Integral Generalizada (GITT) [46, 47] lida com expan-

sões da solução buscada em termos de infinitas bases ortogonais de autofunções, man-

tendo o processo de solução sempre dentro de um domínio contínuo. No entanto, ao

contrário da Transformada Integral Clássica (CITT) [48], o método pode ser aplicado

a sistemas não transformáveis incluindo problemas não lineares, tornando o método

19

aplicável a um número virtualmente infinito de problemas. O sistema resultante é ge-

ralmente composto de um conjunto de equações diferenciais acopladas, que pode ser

facilmente resolvido por rotinas numéricas bem estabelecidas que permitem o controle

da precisão. Porém, como as séries infinitas devem ser truncadas para que qualquer

aplicação possível seja feita, um erro de truncamento está envolvido. Este erro de-

cresce quando o número de termos aumenta, e a solução converge para um valor final.

Devido à natureza da representação em séries, a estimativa do erro pode ser facilmente

obtida, que resulta em um melhor controle do erro global da solução. A desvantagem

associada a essa abordagem é a necessidade de uma manipulação analítica mais elabo-

rada. No entanto, este esforço pode ser consideravelmente minimizado com o emprego

de computação simbólica [49].

Özisik e Murray [50], em 1974, aplicaram pela primeira vez a teoria da transfor-

mada integral em problemas de difusão com condições de contorno variáveis com a

posição e com o tempo. O problema proposto se caracterizava pela presença de ter-

mos não transformáveis pela CITT, que mesmo assim foram inseridos na fórmula de

inversão, o que resultou em um sistema infinito de equações diferenciais ordinárias

de primeira ordem para o potencial transformado. Para a recuperação do potencial

original, foram obtidas soluções aproximadas deste sistema de equações diferenciais,

nascendo aí a natureza híbrida analítico-numérica do procedimento adotado e dando

os primeiros passos na Técnica Transformada Integral Generalizada.

Outro estudo que também contribuiu consideravelmente para a evolução da teoria

da transformada integral foi publicada por Mikhailov [51]. Esse trabalho solucionou

problemas de difusão com coeficientes dependentes do tempo. O que gerou termos

que também não eram transformáveis pela CITT. Usando um problema de autovalor

auxiliar, dependendo do tempo e aplicando o mesmo procedimento usado por Özisik

e Murray [50], Mikhailov obteve um sistema infinito de equações diferenciais com

coeficientes variáveis para o potencial de transformação.

Ambos os trabalhos, Özisik e Murray [50] e Mikhailov [51], criaram as condições

necessárias para o desenvolvimento de uma nova metodologia capaz de resolver pro-

20

blemas de difusão, até então, insolúveis pelas técnicas clássicas, estabelecendo assim

os princípios da Técnica da Transformada Integral Generalizada - GITT.

Devido ao processo evolutivo encontrado no método de transformada integral,

Cotta [46], apresentou uma revisão dos formalismos clássicos, que são agora esten-

didos com ênfase na resolução de problemas não lineares e fortemente acoplados, e

sugeriu técnicas para melhorar a eficiência da solução. Desde então, a GITT começou

a ser disseminada rapidamente e se destacou como uma poderosa ferramenta híbrida

para resolução de problemas difusivos e difusivo-convectivos. Mais tarde, um livro

foi desenvolvido especificamente para aplicações da GITT em problemas difusivos e

difusivos-convectivos [47].

Existe uma vasta literatura sobre GITT. Desta forma, iremos apresentar aqui apenas

os trabalhos publicados que utilizaram esta técnica para resolver o modelo termodinâ-

mico. O primeiro trabalho desta natureza foi publicado por Alves [52], que utilizou

com sucesso a Técnica da Transformada Integral Generalizada num modelo termodi-

nâmico transiente para transferência de calor em fluidos supercríticos sob microgravi-

dade. Este trabalho se concentrou nos efeitos que diferentes condições de contorno de

aquecimento possui no rápido aumento da temperatura média, onde todas as condições

testadas levam a mesma solução de regime permanente. Ficou demonstrado que um

aquecimento via imposição de uma temperatura prescrita mais alta no contorno leva

ao menor tempo de relaxação.

Alves [29] avaliou a dependência da pressão e da temperatura nas propriedades

dos fluidos supercríticos utilizando um modelo termodinâmico do efeito pistão e de-

monstrando que a técnica GITT pode ser aplicada a versões não-lineares de modelos

clássicos de transferência de calor transiente em fluidos altamente compressíveis pró-

ximos ao seu ponto crítico. O autor mostrou que o tempo de relaxação da temperatura

média não é afetado pela variação das propriedades do fluido que torna o modelo ter-

modinâmica não-linear.

Texeira e Alves [31, 53] compararam a solução do modelo termodinâmico via

GITT com a solução de um modelo hidrodinâmico onde a escala acústica não é re-

21

solvida mas sim modelada, usando um método numérico com pré-condicionamento

para baixo Mach ao simular as equações de Navier-Stokes. As duas soluções apenas

fornecem os mesmos resultados quando certas relações termodinâmicas fundamentais

são satisfeitas no uso do modelo hidrodinâmico.

1.3 Objetivos

Independente do modelo utilizado para simular este problema, os resultados nu-

méricos obtidos não estão de acordo com os dados experimentais disponíveis. Como

as ondas acústicas refletem centenas e até mesmo milhares de vezes nas paredes da

cavidade antes do problema atingir o regime permanente, acredita-se que a ausência

deste efeito nos modelos é a principal causa desta discordância.

Desta forma, o objetivo principal deste estudo é modelar as perdas de pressão nas

paredes de uma cavidade contendo um fluido compressível por meio de uma impe-

dância acústica de contorno. Isto será feito tanto nos modelos hidrodinâmico quanto

termodinâmico deste problema, uma vez que ambos podem ser acoplados de modo a

eliminar a necessidade de se simular as escalas acústicas do problema nos casos onde

existe escoamento e o modelo termodinâmico não pode ser utilizado em sua íntegra,

reduzindo drasticamente o custo computacional destas simulações, especialmente para

problema bidimensionais e tridimensionais.

Os objetivos secundários deste trabalho são:

• Comparar os resultados obtidos a partir dos modelos com resultados experimen-

tais disponíveis na literatura para gases ideais.

• Realizar um estudo comparativo da aplicabilidade, vantagens e desvantagens

entre os modelos hidrodinâmico e termodinâmico incluindo impedância acústica

no contorno.

• Estudar o mecanismo de transferência de calor termoacústica, avaliando suas

propriedades e as variáveis que influenciam no processo de geração e transmis-

são de ondas termoacústicas.

22

Capítulo 2

Modelos Matemáticos

2.1 Modelo Hidrodinâmico

2.1.1 Equações Governantes

Aproximações precisas são essenciais para simular as interações dos campos acús-

ticos com as paredes sólidas do sistema. Um campo acústico é formado por uma série

de regiões de altas e baixas pressões chamadas condensações e rarefações, respecti-

vamente. Por esta razão, um modelo matemático preciso deve ser capaz de descrever

estas regiões de compressões e expansões levando em consideração o comportamento

compressível da substância. O campo acústico no domínio unidimensional é carac-

terizado pelo campo de pressão P (x, t ), pela massa específica ρ(x, t ), e pelo campo

de temperatura T (x, t ). Os princípios de conservação da massa, quantidade de mo-

vimento linear e energia devem ser expressos considerando todas as características

compressíveis do escoamento, utilizando uma equação de estado para fornecer a rela-

ção funcional entre a pressão, massa específica e temperatura. Para o nosso problema

de transferência de calor termoacústica, será empregada a formulação compressível

unidimensional completa das equações de Navier-Stokes como equações governantes,

aqui apresentada na forma vetorial.

∂Q

∂t= ∂Ev

∂x− ∂Ei

∂x(2.1)

23

O vetor de varáveis conservadas dependentes Q , é descrito por:

Q = ρ, ρu, ρE

(2.2)

Os fluxos invíscidos, expressos pelo vetor Ei , são definidos como:

Ei =ρu, ρu2 +P,

(ρE +P

)u

(2.3)

E, por fim, os fluxos viscosos, expressos pelo vetor Ev , definidos como:

Ev =

0,4µ

3

∂u

∂x,

4µ

3

∂u

∂x+k

∂T

∂x

(2.4)

onde ρ, u, P , µ, k e T representam a massa específica, a velocidade, a pressão, a

viscosidade dinâmica, a condutividade térmica e a temperatura, respectivamente.

Para o caso de gases termicamente perfeitos, utilizaremos a equação de estado para

gases perfeitos, a equação de energia termodinâmica e o fato de que, para tais gases, a

entalpia depende apenas da temperatura:

P = ρRT , ρH = ρE +P e h = h0f +

∫ T

0CP d s (2.5)

onde H = h + u2

2 e E = e + u2

2 . Obtendo as seguintes relações:

∂P

∂ρ= (γ−1)

(1

2u2 +CP T −h

)(2.6)

∂P

∂(ρu)=−(γ−1)u (2.7)

∂P

∂(ρE)= (γ−1) (2.8)

onde R é a constante dos gases, h é a entalpia por unidade de massa, h0f é a entalpia de

formação ou entalpia de referência, CP é o calor específico a pressão constante e γ é a

razão entre calores específicos. A velocidade do som no gás c pode ser obtida a partir

24

de:

c2 =(∂P

∂ρ

)s= γP

ρ= (γ−1)CP T (2.9)

Por fim, para o caso de gases caloricamente perfeitos, CP é constante e a entalpia pode

ser obtida através de:

h =CP T (2.10)

2.1.2 Condições Iniciais e de Contorno

Complementando nossa formulação matemática, após expormos as equações de

Navier-Stokes, faz-se necessário definir as condições iniciais e de contorno para todas

as variáveis do nosso problema (pressão, temperatura, velocidade, energia e massa

específica). As condições iniciais se referem ao estado do gás na cavidade e suas

respectivas propriedades no tempo inicial. As condições de contorno descrevem as

perdas devido à impedância acústica que se deseja modelar baseando-se na variação

de pressão na paredes da cavidade.

As condições iniciais de pressão, velocidade e temperatura utilizadas podem ser

resumidas em:

P (x,0) = P0 (2.11)

u(x,0) = 0 (2.12)

T (x,0) = T0 (2.13)

A partir destas condições são calculadas as variáveis conservadas no instante t = 0.

As equações abaixo mostram as condições de contorno para temperatura, velo-

cidade e para pressão no caso limite com impedância infinita, o caso particular da

condição de contorno da pressão que será detalhado na próxima seção.

T (0, t ) = T1(t ) e T (L, t ) = T0 (2.14)

u(0, t ) = 0 e u(L, t ) = 0 (2.15)

25

∂P

∂x

∣∣∣x=0

= 0 e∂P

∂x

∣∣∣x=L

= 0 (2.16)

onde T1(t ) representa a temperatura na parede esquerda da cavidade, onde ocorre a

excitação térmica que perturba o sistema. Uma vez que o aumento de temperatura na

parede, na prática, é obtido pode meio de uma fonte externa (circuito elétrico, chama,

explosão), a temperatura na parede nao muda instantaneamente, variando ao longo do

tempo. No caso de aquecimento super rápido, que possa ser considerado aproximada-

mente instântaneo, T1 é uma constante.

2.1.3 Impedância Acústica

Impedância acústica é um importante parâmetro utilizado para avaliar transmissão

e reflexão de uma onda sonora na fronteira de meios de diferente composição. Nos

dias atuais, esta propriedade se relaciona com diversas aplicações práticas, como no

projeto de transdutores ultrassônicos, projeto acústico de ambientes, propagação de

ondas submarinas, etc. Muitas vezes, a impedância acústica específica, tratada sim-

plesmente como impedância, é usada para descrever a propriedade acústica intrínseca

dos materiais.

Rossing et al. [54] definem impedância acústica específica ou impedância acústica

Zs(ω) para uma superfície como:

Zs(ω) = P

vi n(2.17)

onde vi n representa a componente da velocidade do fluido incidindo sobre a superfície

sob consideração, P representa a amplitude complexa da onda de pressão e ω repre-

senta a velocidade angular da onda termoacústica. A unidade no Sistema Internacional

(SI) de Impedância Acústica é [Pa.sm ], denominado de Rayl, em homenagem ao Lorde

Rayleigh. O inverso da impedância acústica denomina-se Admitância Acústica.

Quando se avalia o campo sonoro na interface entre dois meios distintos, a com-

ponente apropriada da velocidade da onda de pressão é aquela normal à superfície,

originando uma impedância acústica específica normal à superfície, também chamada

26

impedância específica do contorno. A mesma é dependente da onda sonora incidente,

uma vez que a componente normal do vetor velocidade de partícula em um ponto qual-

quer da interface é influenciado não somente pela pressão sonora local, como também

pelas ondas que chegam de todos os outros pontos do meio excitado. Consequente-

mente, não é possível especificar uma única impedância de contorno sem levar em

consideração a amplitude e as distribuições de fase da onda incidente sobre a interface.

Portanto, uma superfície plana com uma impedância específica de contorno apresenta

uma resposta diferente dependendo da forma da onda incidente [9].

Fig. 2.1: Impedância específica de contorno.

Em muitos casos, superfícies de materiais na vizinhança de fluidos podem ser ca-

racterizados como localmente reagentes, isto significa que ZS é independente da natu-

reza do campo de pressão acústica na vizinhança, bem como de outros fatores do meio.

Em algumas situações particulares, é razoável assumir que existe uma relação linear

entre a componente da velocidade local normal da partícula e a pressão naquele ponto

de contato. Desta forma a superfície do material pode ser caracterizada em termos de

uma única impedância específica de contorno [55].

A hipótese de localmente reagente significa que a velocidade do material na su-

perfície não é afetada por pressões próximas ao ponto de contato. Para este trabalho,

onde o problema consiste em ondas termoacústicas confinadas em uma cavidade uni-

dimensional, a hipótese de localmente reagente é perfeitamente aplicável, uma vez

que a pressão ocorre somente na direção x. No estudo unidimensional, mostrado na

27

figura 2.1, a onda termoacústica atinge a superfície perpendicularmente e é refletida

na mesma direção. Para uma superfície estacionária e sem nenhuma resposta ativa a

pertubação da onda, a condição de contorno apropriada para a amplitude complexa P

que satisfaz a equação de Helmholtz é dada por:

iωρP =−ZS∆P ·n (2.18)

onde i é o número complexop−1, ω representa a velocidade angular da onda termoa-

cústica, ∆P é a variação entre a pressão incidente sobre a superfície e pressão refletida

pela mesma, e n é vetor unitário área apontando na direção normal, saindo da super-

fície do material para o fluido. Quando as ondas sonoras atravessam uma interface

entre dois materiais diferentes, se eles possuírem a mesma impedância acústica, não

ocorre reflexão e a onda é transmitida ao segundo material. Se existir diferença de

impedância acústica, é esta diferença que origina uma maior ou menor reflexão das

ondas sonoras. Uma superfície que é perfeitamente rígida refletindo completamente a

onda de pressão tem |ZS | =∞. No outro extremo, uma superfície que absorve a pres-

são possui |ZS | = 0. Isto é, por exemplo, o que se assume para superfície superior do

oceano quando se avalia a transmissão de ondas submarinas.

As condições de contorno nas paredes da cavidade relacionadas à impedância acús-

tica são obtidas a partir da equação (2.18). Para nossa cavidade unidimensional, pode-

mos escrever:

∆P ·n =±∂P

∂x(2.19)

onde o sinal ± será positivo ou negativo de acordo com o sentido do vetor área unitário

n em relação a parede onde ocorre a reflexão (+ para parede direita, - para parede

esquerda). Utilizando uma variação de pressão a partir da pressão inicial P0, isto é:

P = P −P0 (2.20)

28

Retomando a equação (2.18) e substituindo os resultados acima, obtemos:

iωρ(P −P0) =±ZS∂P

∂x(2.21)

Definindo Z∗ = ZSiωρL como impedância adimensional, chegamos a expressão para con-

dição de contorno nas paredes:

Z∗∂P

∂x± P −P0

L= 0 (2.22)

onde Z∗ varia de 0 a ∞. No intuito de trabalhar com uma impedância que varia de 0 a

1, podemos utilizar a seguinte transformação:

Z∗ = log1

1−Z=⇒ Z = 1−e−Z∗

(2.23)

Fixamos assim, para todas as variáveis do escoamento compressível, todas as con-

dições iniciais e de contorno do modelo hidrodinâmico.

2.2 Modelo Termodinâmico

Neste trabalho, além do modelo hidrodinâmico baseado nas equações de Navier-

Stokes, um outro modelo do problema de ondas termoacústicas confinadas dentro de

uma cavidade é utilizado. Este modelo é chamado de termodinâmico porque avalia

os efeitos de mudanças na temperatura, pressão e volume levando em consideração

as variações de calor e trabalho do sistema em escala macroscópica [56]. O modelo

termodinâmico tem sido usado com sucesso na modelagem da transferência de calor

termoacústica utilizando a equação da calor com um termo fonte proporcional a deri-

vada temporal da temperatura média do sistema [27, 57]. Este termo fonte modela o

mecanismo de compressão adiabática responsável pelo efeito pistão.

29

2.2.1 Equação do Calor com termo fonte de Compressão Adiabática

O modelo termodinâmico não considera o escoamento do fluido, tratando a cavi-

dade como um sistema e considera um balanço de energia na cavidade na forma:

ρCP∂T

∂t−ρ (

CP − CV) κT

αP

dPn

d t= ∂

∂x

(k∂T

∂x

)(2.24)

onde Pn representa a pressão real dentro da cavidade, κT é definido como compressi-

bilidade isotérmica e αP como coeficiente de expansão térmica isobárica. Assim como

todas as propriedades existentes na equação, κT e αP são avaliados no tempo inicial:

κT = 1

ρ

∂ρ

∂P

∣∣∣T0

, αP =− 1

ρ

∂ρ

∂T

∣∣∣P0

(2.25)

Esta equação (2.24) diz que toda a entropia gerada pelo fluido é dissipada via con-

dução de calor. A entropia por sinal, possui uma contribuição não desprezível da vari-

ação de pressão. Esta contribuição é normalmente desprezada na derivação tradicional

da equação do calor.

As condições iniciais e de contorno para a temperatura T descrita pela equação

(2.24) são:

T (x, t = 0) = T0 , T (x = 0, t ) = T1 e T (x = L, t ) = T0 (2.26)

sendo T1 a temperatura da parede esquerda e T0 a temperatura inicial do gás dentro da

cavidade e também da parede direita.

30

2.2.2 Pressão Termodinâmica

A dedução da expressão para a pressão termodinâmica PT pode ser feita a partir da

forma diferencial completa da massa específica ρ(PT ,T ):

dρ = ∂ρ

∂P

∣∣∣T0

dPT + ∂ρ

∂T

∣∣∣P0

dT (2.27)

que pode ser dividida por d t para tornar-se:

dρ

d t= ∂ρ

∂P

∣∣∣T0

dPT

d t+ ∂ρ

∂T

∣∣∣P0

dT

d t(2.28)

Substituindo as definições de κT e αT , obtemos:

dρ

d t= ρ0κT

dPT

d t−ρ0αP

dT

d t(2.29)

Integrando a equação acima pelo comprimento da cavidade unidimensional, temos:

∫ L

0

dρ

d td x =

∫ L

0ρ0κT

dPT

d td x −

∫ L

0ρ0αP

dT

d td x (2.30)

Uma vez que nosso problema consiste numa cavidade unidimensional fechada por

paredes impermeáveis, ocorre a conservação da massa no sistema. Desta forma:

∫ L

0ρ(x, t )d x =

∫ L

0ρ(x, t = 0)d x = constante (2.31)

A partir do resultado acima, podemos deduzir para o lado esquerdo da equação

(2.30):

∫ L

0

dρ

d td x = d

d t

∫ L

0ρd x = d

d t[constante] = 0 (2.32)

Pode ser obtida uma temperatura média da massa fluida contida na cavidade, esta

temperatura tem papel fundamental na modelagem do mecanismo de compressão adi-

abática. A temperatura média Tb(t ) é definida como:

31

Tb(t ) = 1

L

∫ L

0T (x, t )d x (2.33)

Retomando a equação (2.30), de posse do resultado (2.32), e utilizando a definição

para temperatura global do sistema (2.33), podemos deduzir a expressão final para

pressão termodinâmica:

0 =∫ L

0ρ0κT

dPT

d td x −

∫ L

0ρ0αP

dT

d td x

0 = ρ0κTdPT

d tL−ρ0αP

∫ L

0

dT

d td x

0 = ρ0κTdPT

d tL−ρ0αP

d

d t

∫ L

0T d x

0 =ρ0κTdPT

d tL−ρ0αPL

dTb

d t

dPT

d t= αP

κT

dTb

d t(2.34)

chegamos assim a uma expressão para derivada temporal da pressão termodinâmica,

onde foi assumido que a pressão termodinâmica é espacialmente uniforme, ou seja,

somente varia com o tempo. Esta hipótese simplificadora é baseada na magnitude

do tempo característico acústico, que é muito menor que os tempos característicos do

efeito pistão e da difusão térmica.

A massa específica pode ser obtida da equação de estado. No nosso caso, onde

todas as propriedades são consideradas constantes e são avaliadas no tempo inicial,

integrando a equação (2.27), obtemos:

ρ

ρ0= 1+κT (P −P0)−αP (T −T0) (2.35)

onde P0 e ρ0 representam a pressão inicial e a massa específica inicial, respectiva-

mente.

32


A impedância acústica nas paredes causa uma diminuição da pressão dentro da

cavidade. Desta forma a pressão termodinâmica PT , deduzida na subseção anterior,

não ocorre na prática, uma vez que existem perdas de amplitude quando as ondas

termoacústicas são refletidas nas paredes.

Fig. 2.2: Choque e reflexão da onda termoacústica.

Sendo ∆P1, a amplitude da onda termoacústica no primeiro choque com a parede

esquerda, assumindo uma perda de amplitude devido à impedância na parede Z , temos

∆P1 = Z ×∆PT (2.36)

Após sofrer uma pequena perda de amplitude devido a impedância da parede es-

querda, a onda termoacústica refletida (∆P1) viaja ao longo da cavidade e interage

agora com a parede direita:

∆P2 = Z ×∆P1 = Z 2 ×∆PT (2.37)

O processo é o mesmo para a terceira colisão entre a onda termoacústica e a parede

esquerda:

∆P3 = Z ×∆P2 = Z 3 ×∆PT (2.38)

O processo se repete milhares de vezes, até que a onda termoacústica seja dissi-

33

pada pelas forças viscosas agindo na cavidade e pela impedância existente nas paredes.

Sendo n o número de colisões entre a onda termoacústica e as paredes da cavidade,

podemos escrever:

∆Pn = Z n ×∆PT (2.39)

O tempo característico acústico ta representa fisicamente o tempo que a onda ter-

moacústica leva para atravessar a cavidade inteira viajando a uma velocidade c. Desta

forma num tempo arbitrário t , podemos assumir que a onda termoacústica em estudo

irá se chocar aproximadamente tta

vezes contra as paredes da cavidade, ou seja

n = t

ta(2.40)

Substituindo o resultado acima, a pressão real ou modificada Pn(t ) pode ser obtida

através da expressão:

Pn −P0 = Z t/ta (PT −P0) (2.41)

onde P0 representa a pressão no tempo inicial (t = 0) na cavidade, Z representa a impe-

dância acústica que varia entre 0 (absorção total da onda termoacústica) e 1 (reflexão

sem perdas da onda termoacústica) e ta representa o tempo característico acústico,

definido através da expressão ta = L/c.

Uma vez que a equação (2.24) está escrita em função da derivada temporal da

pressão modificada dPnd t , derivando a equação acima em relação ao tempo, obtemos:

dPn

d t= 1

talog[Z ]Z t/t a(PT −P0)+Z t/t a dPT

d t(2.42)

tornando completo o sistema de equações deste problema, que inclui as equações

(2.24) e (2.34). Notamos ainda que a pressão modificada se torna igual a pressão

termodinâmica quando Z = 1, recuperando assim o modelo termodinâmico tradicio-

nal.

34

Capítulo 3

Métodos Numéricos

3.1 Modelo Hidrodinâmico

3.1.1 Discretização Espacial

As equações do modelo hidrodinâmico na forma vetorial, conforme discutido no

capítulo anterior, consistem num termo transiente a esquerda da igualdade e dois ter-

mos espaciais a direita. Para avaliação dos termos com derivação espacial, utilizare-

mos dois diferentes tipos de discretizações, um para cada termo a direita da igualdade.

Na equação (3.1), o termo ∂Ev∂x é chamado de termo viscoso, enquanto que o segundo

termo ∂Ei∂x é chamdo de termo invíscido, visto que não envolve dissipações de qualquer

espécie.

∂Q

∂t= ∂Ev

∂x− ∂Ei

∂x(3.1)

As equações de Navier-Stokes, que modelam a conservação de quantidade de mo-

vimento linear nas equações de governo, consistem nas equações de Euler somadas

a um termo viscoso que envolve tensão de cisalhamento. Dessa forma, os métodos

numéricos utilizados no cálculo dos fluxos invíscidos das equações de Euler são per-

feitamente aplicáveis aos termos invíscidos das equações de Navier-Stokes. Para o

termo viscoso, que é mais simples, aproximações de diferenças finitas centradas são

35

aplicadas.

3.1.1.1 Fluxos Invíscidos

Para cálculo dos fluxos invíscidos, são utilizados o método da Divisão do Vetor

de Fluxo (Flux Vector Splitting - FVS). A idéia central por trás do método FVS é

dividir as contribuições do fluxo em componentes positivas e negativas, onde a divisão

baseia-se nos autovalores do sistema e no comportamento destes autovalores ao longo

do escoamento. O problema fundamental a ser resolvido é determinar o fluxo correto

nas faces do volume de controle na direção positiva e negativa.

Escrevendo o fluxo invíscido da equação de Navier-Stokes em duas parcelas de

fluxo, uma positiva e outra negativa, temos:

Ei =E+i +E−

i (3.2)

onde cada parcela é baseada na divisão das características da matriz Jacobiana

A= ∂Ei

∂Q=

0 1 0

∂P∂ρ

−u2 ∂P∂(ρu) +2u ∂P

∂(ρE)(∂P∂ρ

−H)

u ∂P∂(ρu) u +H

(1+ ∂P

∂(ρE)

)u

=A++A− (3.3)

sendo, portanto, a divisão destas características obtidas através da separação das velo-

cidades das ondas, que são os autovalores positivos e negativos:

λi =λ+i +λ−

i , com λ+i ≥ 0 e λ−

i ≤ 0 (3.4)

Os autovalores da matriz Jacobiana (3.3) são:

λ1 = u , λ2 =C1 −C2 e λ3 =C1 +C2 (3.5)

36

onde C1 e C2 são dados por:

C1 = u + 1

2

∂P

∂(ρu)+ u

2

∂P

∂(ρE)

C2 = 1

2

√(∂P

∂(ρu)+u

∂P

∂(ρE)

)2

+4

(∂P

∂ρ+u

∂P

∂(ρu)+H

∂P

∂(ρE)

)(3.6)

A matriz Jacobiana A e seus autovalores, que representam fisicamente as três veloci-

dades do som no fluido, se relacionam pela equação:

M−1 A M =Λ (3.7)

onde Λ é a matriz diagonal na qual seus elementos representam os três autovalores λi .

Já as colunas da matriz M são os autovalores da matriz A. As matrizes M e M−1 são

definidas como:

M =

1 − ρ

4C2

ρ4C2

u − ρ4C2

(C1 −C2) ρ4C2

(C1 +C2)

−(∂P∂ρ

+u ∂P∂(ρu)

)/∂P

∂(ρE) − ρ4C2

(C3 −uC2) ρ4C2

(C3 +uC2)

, (3.8)

M−1 = ρ

4C4

ρ

C2

(C3−uC1

2

)ρu2C2

− ρ2C2

−C+5

C2

/∂ρ

∂(ρE)C+

6C2

/∂ρ

∂(ρE)u−C1

C2−1

−C−5

C2

/∂ρ

∂(ρE)C−

6C2

/∂ρ

∂(ρE)u−C1

C2+1

, (3.9)

onde C3, C4, C±5 e C±

6 são dados por:

C3 = H + u

2

(∂P

∂(ρu)+u

∂P

∂(ρE)

)(3.10)

C4 = ρ2

8C2

u2 +C3 −uC1 +

(∂P

∂ρ+u

∂P

∂(ρE)

)/ ∂ρ

∂(ρE)

, (3.11)

C±5 = (C1 ±C2)

(∂P

∂ρ+u

∂P

∂(ρu)

)+u(C3 ±uC2)

∂P

∂(ρE)e (3.12)

C±6 = ∂P

∂ρ+u

∂P

∂(ρu)+ (C3 ±uC2)

∂P

∂(ρE). (3.13)

37

Combinando as equações (3.3) e (3.7), chegamos a:

A+ = M Λ+ M−1 e A− = M Λ− M−1 (3.14)

O método FVS permite que as características positivas e negativas sejam separadas

através de uma apropriada discretização atrasada (upwind) e avançada (downwind)

utilizando uma formulação conservativa, uma vez que opera fluxos e não a velocidade

diretamente, como o método Wave Speed Splitting, por exemplo. O vetor fluxo, para

o termo invíscido das equações de governo, é uma função homogênea de ordem um,

logo podemos escrever:

Ei = ∂Ei

∂QQ = A Q (3.15)

e, consequentemente:

E±i = A± Q = (M Λ± M−1) Q (3.16)

A partir do resultado acima podemos obter a discretização espacial dos fluxos invís-

cidos através de uma aproximação de Diferenças Finitas atrasada/avançada de segunda

ordem com TVD utilizando um limitador de fluxo.

Desta forma, para o cálculo dos fluxos positivos, podemos escrever:

(∂E+

i

∂x

)i

= E+i (i )−E+

i (i −1)

d x+ φ+

1

2

E+i (i )−E+

i (i −1)

d x− φ+

2

2

E+i (i −1)−E+

i (i −2)

d x(3.17)

E para o cálculo dos fluxos negativos, temos:

(∂E−

i

∂x

)i= E−

i (i +1)−E−i (i )

d x+ φ−

1

2

E−i (i +1)−E−

i (i )

d x− φ−

2

2

E−i (i +2)−E−

i (i +1)

d x(3.18)

Nas equações (3.17) e (3.18) acima, φ1 e φ2 representam os limitadores de fluxo

conhecidos como superbee [11], dados pelas equações (3.19) e (3.20). No intuito

de garantir as condições necessárias para propriedade TVD, limitadores de fluxo são

normalmente escritos em termos de variações sucessivas das variáveis dependentes.

Estas variações, escritas sob a forma de razões r , devem ser calculadas para os dois

38

limitadores de forma a avaliar as variações no sentido positivo (r+1 e r+

2 ) e negativo (r−1

e r−2 ) conforme equações abaixo:

φ±1 = max[0,min(2r±

1 ,1),min(r±1 ,2)] (3.19)

φ±2 = max[0,min(2r±

2 ,1),min(r±2 ,2)] (3.20)

Desta forma, para os fluxos positivos, escrevemos:

r+1 = E+

i (i +1)−E+i (i )

E+i (i )−E+

i (i −1)(3.21)

r+2 = E+

i (i )−E+i (i −1)

E+i (i −1)−E+

i (i −2)(3.22)

e, similarmente, podemos escrever para os fluxos negativos:

r−1 = E−

i (i )−E−i (i −1)

E−i (i +1)−E−

i (i )(3.23)

r−2 = E−

i (i +1)−E−i (i )

E−i (i +2)−E−

i (i +1)(3.24)

3.1.1.2 Fluxos Viscosos

O termo viscoso da equação de governo consiste num termo relacionado com a

difusão térmica e em outro relacionado à viscosidade do fluido. A discretização dos

fluxos viscosos pode ser considerada simples, se comparada com a discretização dos

fluxos invíscidos mostrada na seção anterior. O termo viscoso consiste em:

Ev =

0

τxx

uτxx + qxx

(3.25)

onde a componente τxx , que refere-se à equação de conservação da quantidade de

movimento linear, pode ser escrita como:

τxx = 2µ

(∂u

∂x

)+λ

(∂u

∂x

)= 2µ

(∂u

∂x

)− 2

3µ

(∂u

∂x

)= 4µ

3

(∂u

∂x

)(3.26)

39

onde λ é a viscosidade global, cujo valor é dado por λ = −2/3µ de acordo com a

hipótese de Stokes. A componente uτxx + qxx , que refere-se à equação de conservação

da energia, pode ser escrita como:

uτxx + qxx = u4µ

3

(∂u

∂x

)+k

∂T

∂x(3.27)

Substituindo os resultados acima, obtemos:

Ev =

0

4µ3

(∂u∂x

)u 4µ

3

(∂u∂x

)+k ∂T

∂x

(3.28)

Aproximando as derivadas ∂u∂x e ∂T

∂x , por meio de uma discretização de Diferenças

Finitas Centradas de segunda ordem, escrevemos:

Ev(i ) =

0

4µi3

(ui+1−ui−12∆x

)ui

4µi3

(ui+1−ui−12∆x

)+kiTi+1−Ti−1

2∆x

(3.29)

onde o vetor Ev deve então ser calculado para todos os pontos de nossa malha unidi-

mensional.

Como estamos interessados no fluxo viscoso, que é o termo ∂Ev∂x da equação de

governo mostrado abaixo, faz-se necessária uma nova discretização centrada. Desta

forma, podemos escrever, para o fluxo do ponto i da malha unidimensional:

(∂Ev

∂x

)i= 1

2∆x

[Ev (i+1) −Ev (i−1)

](3.30)

Como calculamos os fluxos viscosos apenas entre os pontos 2 e N − 1 da nossa

malha, uma vez que as variáveis do nosso problema nos contornos são obtidas por meio

das condições de contorno, nao há problema de aplicar nossa discretização centrada ao

longo da malha.

40

3.1.2 Discretização Temporal

A discretização do termo transiente da equação de governo será realizada utilizando

a formulação de Shu-Osher do método Runge-Kutta explícito de segunda ordem com

dois estágios [45], conhecido na literatura como SSPRK(2,2), conforme deduzido no

capítulo introdutório. Este método, nesta formulação, tem cada estágio definido como:

1 estágio:

u(1) = un +∆tF (un) (3.31)

2 estágio:

u(2) =∆tF (u(1)) (3.32)

Atualização:

u(n+1) = 1

2un + 1

2u(1) + 1

2u(2) (3.33)

Partindo da equação de governo, podemos escrever:

∂Q

∂t= ∂Ev

∂x− ∂Ei

∂x= ∂Ev

∂x−

(∂E+

i

∂x− ∂E−

i

∂x

)(3.34)

Reescrevendo para a formulação Runge-Kutta, temos:

∂Q

∂t= F (Q) (3.35)

onde:

F (Q) = ∂Ev

∂x−

(∂E+

i

∂x− ∂E−

i

∂x

)(3.36)

de posse de F (Q), partimos para formulação Shu-Osher do SSPRK(2,2):

1 estágio:

Q(1) = Qn +∆t F (Qn) (3.37)

2 estágio:

Q(2) =∆t F (Q(1)) (3.38)

41

Atualização:

Q(n+1) = 1

2Qn + 1

2Q(1) + 1

2Q(2) (3.39)


Nesta subseção, é descrita a aproximação numérica utilizada para a condição de

contorno de impedância acústica utilizada nas paredes da cavidade e demostrada no

capítulo anterior.

Z∗∂P

∂x± P −P0

L= 0 (3.40)

Para a parede esquerda da cavidade (x = 0), utilizando a seguinte aproximação de

segunda ordem avançada para ∂P∂x :

∂P

∂x=−3

2P1 +2P2 − 1

2P3 (3.41)

E, substituindo a aproximação acima na equação 3.40, deduzimos:

Z

∆x

(− 3

2P1 +2P2 − 1

2P3

)− P1 −P0

L= 0 (3.42)

P1

(− ∆x

Z L− 3

2

)=−∆x

Z LP0 + 1

2P3 −2P2 (3.43)

P1 =

(∆xZ L P0 − 1

2 P3 +2P2

)(∆xZ L + 3

2

) (3.44)

Para a parede direita da cavidade (x = L), utilizando a seguinte aproximação de

segunda ordem atrasada para ∂P∂x :

∂P

∂x= 3

2PN x −2PN x−1 + 1

2PN x−2 (3.45)

E, utilizando o mesmo procedimento da parede oposta, substituímos a aproximação

acima na equação 3.40, podemos deduzir:

Z

∆x

(3

2PN x −2PN x−1 + 1

2PN x−2

)− P1 −P0

L= 0 (3.46)

42

PN x

(− ∆x

Z L− 3

2

)=−∆x

Z LP0 + 1

2PN x−2 −2PN x−1 (3.47)

PN x =

(∆xZ L P0 − 1

2 PN x−2 +2PN x−1

)(∆xZ L + 3

2

) (3.48)

Para simulação computacional do problema através das equações de governo dis-

cretizadas conforme demonstrado nesta seção, foi implementado um código na lingua-

gem Fortran 90 com precisão dupla e paralelizado. A especificação dos computadores

utilizados durante a execução do trabalho, estudos de custo computacional e resultados

obtidos com a paralelização estão especificados no apêndice A.

3.2 Modelo Termodinâmico

3.2.1 Sistema de Equações

Conforme mostrado na seção referente aos modelos matemáticos, o modelo termo-

dinâmico resulta no sistema acoplado de equações diferenciais:

ρCP∂T

∂t−ρ (

CP − CV) κT

αP

dPn

d t= ∂

∂x

(k∂T

∂x

)(3.49)

dPT

d t= αP

κT

dTb

d t(3.50)

dPn

d t= 1

talog[Z ]Z t/t a(PT −P0)+Z t/t a dPT

d t(3.51)

onde a temperatura média é definida como Tb(t ) = 1L

∫ L0 T (x, t )d x.

3.2.2 Mudança para variáveis adimensionais e filtro

Dispomos de um sistema de equações diferenciais acopladas formado pelas equa-

ções (3.49), (3.50) e (3.51). Antes de procedermos à solução do problema através da

Técnica da Transformada Integral Generalizada, é interessante reescrever o problema

de maneira adimensional, para isto as seguintes transformações são propostas:

ξ= x

L, τ= t

tD, θ = T −T0

T1 −T0, Π= P

P0(3.52)

43

α∗P =αP ·∆T , κ∗T = κT ·P0, δ= tD

ta(3.53)

onde tD = L2/α representa o tempo característico de difusão térmica, a difusividade

térmica é dada pela expressão α = k/(ρCP ) e a variação de temperatura é calculada

como ∆T = T1 −T0.

Após algumas operações algébricas, podemos escrever o sistema em termos das

variáveis adimensionais:

∂θ

∂τ= ∂2θ

∂ξ2+

(1− 1

γ

)κ∗Tα∗

P

dΠn

dτ(3.54)

dΠT

dτ= α∗

P

κ∗T

dθb

dτ(3.55)

dΠn

dτ= δ log[Z ]Z δτ(ΠT −1)+Z δτdΠT

dτ(3.56)

onde a temperatura global adimensional pode ser obtida através da equação:

θb(τ) =∫ 1

0θ(ξ,τ)dξ (3.57)

com condições iniciais e de contorno adimensionais dadas por:

θ(ξ,τ= 0) = 0 , θ(ξ= 0,τ) = 1 e θ(ξ= 1,τ) = 0 (3.58)

Um importante passo no procedimento de solução é fazer as condições de contorno

homogêneas antes de aplicar a transformada integral. Uma significante melhoria na

taxa de convergência das séries é atingida devido a este passo. Isto pode ser feito

introduzindo um filtro para a temperatura adimensional dado por:

θ(ξ,τ) =F(ξ)+Θ(ξ,τ) (3.59)

44

onde o filtro é definido como a solução de regime permanente:

F= 1−ξ (3.60)

A partir da equação (3.59) temos o valor do filtro médio e da temperatura adimen-

sional homogênea global respectivamente:

Fb = 1

2e θb(τ) =Fb +Θb = 1

2+Θb (3.61)

Podemos obter também as condições iniciais e de contorno para a nova variável Θ:

Θ(ξ,τ= 0) = ξ−1 , Θ(ξ= 0,τ) = 0 e Θ(ξ= 1,τ) = 0 (3.62)

O novo sistema, com variáveis adimensionais e incluindo o filtro na temperatura

pode ser escrito como:

∂Θ

∂τ= ∂2Θ

∂ξ2+

(1− 1

γ

)κ∗Tα∗

P

dΠn

dτ(3.63)

dΠT

dτ= α∗

P

κ∗T

dΘb

dτ(3.64)

dΠn

dτ= δ log[Z ]Z δτ(ΠT −1)+Z δτdΠT

dτ(3.65)

Θb(τ) =∫ 1

0Θ(ξ,τ)dξ (3.66)

3.2.3 Transformação Integral

Utilizando a Técnica da Transformada Integral Generalizada, as dependências es-

paciais e temporais do problema devem ser separadas, transformando o problema num

sistema de equações diferenciais ordinárias. Propomos assim uma solução para a equa-

ção (3.63) da seguinte forma:

Θ(ξ,τ) =N∑

i=1ψi (ξ)Θi (τ) (3.67)

45

onde a autofunção é obtida do problema de Sturm-Liouville

∂2ψi

∂ξ2+β2

i ψi (ξ) = 0 (3.68)

com as seguintes condições de contorno:

ψi (0) = 0 e ψi (1) = 0 (3.69)

chegamos às autofunções e aos autovalores, respectivamente:

ψi (ξ) = sin[βiξ] , βi = iπ (3.70)

onde i = 1,2, ...,∞. Como estas autofunções são ortogonais, a temperatura homogênea

transformada pode ser definida, a partir das equações dos autovalores, a relação:

Θi (τ) =∫ 1

0ψi (ξ)Θ(ξ,τ)dξ (3.71)

utilizando a norma:

Ni =∫ 1

0ψi (ξ)2dξ= 1

2(3.72)

temos a autofunção normalizada ψi =ψi /p

Ni .

Multiplicando a equação (3.63) pela autofunção ψi , integrando o resultado pelo

comprimento adimensional da cavidade e aplicando a transformação (3.71) temos:

dΘi

dτ=

∫ 1

0ψi (ξ)

∂2Θ

∂ξ2dξ+

(1− 1

γ

)κ∗Tα∗

P

dΠn

dτηi (3.73)

onde ηi representa o coeficiente de transformação integral, definido por simplicidade

como:

ηi =∫ 1

0ψi (ξ)dξ=p

2(1−cos[βi ]

βi

)(3.74)

46

Integrando por partes o primeiro termo do lado direito da equação (3.73), aplicando

as condições de contorno (3.62) e (3.69) e aplicando a transformação integral (3.71)

podemos obter:

∫ 1

0ψi (ξ)

∂2Θ

∂ξ2dξ=

∫ 1

0

d 2ψi

dξ2=−β2

i

∫ 1

0ψi (ξ)Θ(ξ,τ)dξ=−β2

iΘi (τ) (3.75)

Substituindo o resultado obtida acima na equação (3.73), chegamos:

dΘi

dτ=−β2

iΘi (τ)+(1− 1

γ

)κ∗Tα∗

P

dΠn

dτηi (3.76)

A equação transformada integral acima está sujeita a uma condição inicial transfor-

mada, obtida aplicando o mesmo procedimento de transformação na condição inicial

em (3.62), logo:

Θ(τ= 0) =∫ 1

0(ξ−1)ψi (ξ)dξ=−

p2

βi(3.77)

A temperatura média escrita em termos de Θ(τ):

Θb =∫ 1

0

N∑j=1

ψ j (ξ)θ j (τ)dξ=N∑

j=1

∫ 1

0ψ j (ξ)θ j (τ)dξ=

N∑j=1

η jΘ j (τ) (3.78)

Substituindo o resultado acima na equação da pressão termodinâmica (3.64) ,te-

mos:

ΠT = α∗P

κ∗T

N∑j=1

η jΘ j (3.79)

Chegamos finalmente a um sistema de equações formado pelas equações (3.76),

(3.65) e (3.79), descrevendo respectivamente a temperatura adimensional, a pressão

modificada pela impedância e a pressão termodinâmica. A pressão modificada e a

pressão termodinâmica estão sujeitas as seguintes condições iniciais:

Πn(τ= 0) = 1 , ΠT (τ= 0) = 1 (3.80)

47

uma vez que se assume que ambas são iguais a pressão inicial no tempo t = 0.

O sistema formado pelas equações (3.76), (3.65) e (3.79) possui um elevado nú-

mero de derivadas temporais, o que torna o processo de solução numérica lento. A

solução encontrada foi combinar as equações de forma a dispor os termos que pos-

suem derivadas temporais no lado esquerda da equação, e os outros termos no lado

direito, onde obtemos:

dΘi

dτ−

(1− 1

γ

)ηi Z δτ

N∑j=1

dΘ j

dτ=−β2

iΘi (τ)+(1− 1

γ0

)κ∗Tα∗

P

ηiδ log[Z ]Z δτ(ΠT −1) (3.81)

que pode ser escrito de uma forma mais compacta como:

N∑j=1

Ai , jdΘi

dτ= Bi , jΘi +Ci , j (ΠT −1) (3.82)

onde

Ai , j = δi , j −(1− 1

γ

)ηiη j Z δτ (3.83)

Bi , j =−β2i δi , j (3.84)

Ci , j =(1− 1

γ0

)κ∗Tα∗

P

ηiδ log[Z ]Z δτδi , j (3.85)

sendo Ai , j a matriz de coeficientes de transformação integral, com δi , j representando o

delta de Kronecker, Bi , j uma matriz contendo os termos que multiplicam a temperatura

e Bi , j a matriz de termos independentes da temperatura.

A matriz A é inversível, e sua matriz inversa é dada por:

Ai , j =(γ−1)

(∑Nk=1η[k]2 −η[i ]2

)Z δτ−γ

(γ−1)∑N

k=1η[k]2Z δτ−γ , para i = j (3.86)

Ai , j =− (γ−1)η[i ]η[ j ]Z δτ

(γ−1)∑N

k=1η[k]2Z δτ−γ , para i 6= j (3.87)

De posse da matriz inversa de A, o sistema modficado que será resolvido consiste

em

48

N∑j=1

dΘi

dτ= A−1

i , j [Bi , jΘi +Ci , j (ΠT −1)] (3.88)

No intuito de obter a solução do sistema formado pela equação acima e pelas equa-

ções (3.65) e (3.79), o somatório das séries foi truncado. Uma vez que as equações

(3.68) a (3.70) representam um sistema de autovalores da classe de Sturm-Liouville, a

série proposta como solução (3.67) é convergente. Então, a precisão da solução pro-

posta para o modelo termodinâmico depende inteiramente do número de termos da

série N . Logo, quanto maior for N , menor será o erro de truncamento introduzido

na solução do sistema. Por outro lado, aumentando o número de termos da série, au-

menta o custo computacional para se resolver numericamente o sistema de N equações

acopladas resultante do modelo.

Após obter-se a solução do sistema de equações acopladas para o vetorΘ, a solução

analítica para temperatura adimensional pode ser obtida:

θ(ξ,τ) = 1−ξ+N∑

i=1ψi (ξ)Θi (τ) (3.89)

Da mesma forma, obtemos também a temperatura global adimensional do sistema:

θb(τ) =Fb +Θb = 1

2+

N∑i=1

ηiΘi (τ) (3.90)

O modelo termodinâmico foi implementado no pacote de computação simbólica

Mathematica e o sistema foi resolvido algebricamente através do comando NDSolve.

Utilizando o processo descrito acima, o custo computacional é extremamente baixo e

a solução pode ser obtida em poucos minutos.

49

Capítulo 4

Resultados

Neste capítulo expomos e discutimos os resultados obtidos a partir dos códigos im-

plementados. Os resultados foram divididos entre aqueles obtidos através do modelo

hidrodinâmico e os obtidos através do modelo termodinâmico, esta divisão se faz ne-

cessária, uma vez que as considerações e hipóteses, principalmente para a variação da

pressão dentro da cavidade, são diferentes entre os dois métodos.

É importante ressaltar que o código hidrodinâmico foi inteiramente desenvolvido e

paralelizado pelo autor desta dissertação (ver apêndice), sem e com os efeitos de im-

pedância na parede. Por outro lado, partindo de um código termodinâmico tradicional

desenvolvido pelo orientador em outro trabalho, o autor incorporou a modelagem dos

efeitos de impedância.

4.1 Modelos Hidrodinâmicos

Através do modelo hidrodinâmico é possível realizar simulações de diversos tipos,

alterando os parâmetros de condições de contorno de forma rápida e prática. Devido a

esta facilidade este modelo foi utilizado tanto para modelar a impedância acústica de

parede, objetivo principal do trabalho, quanto para estudar a características das ondas

termoacústicas geradas a partir de diferentes taxas de mudança da temperatura de ex-

citação da cavidade contendo fluido compressível. Antes de discutir em detalhes estes

50

resultados, são reportadas as etapas de verificação e validação do código desenvolvido.

4.1.1 Equação de Burgers

No intuito de avaliar os esquemas TVD e SSP propostos, a Equação de Burgers in-

víscida foi utilizada como um dos casos testes. Esta equação possui um termo convec-

tivo não-linear e um termo transiente, similar às equações de interesse neste trabalho,

porém mais simples por ser uma equação escalar e não vetorial.

∂u

∂t+u

∂u

∂x= 0 (4.1)

A equação, escrita acima, pode ser vista como uma simples analogia das Equações

de Euler para o escoamento de um fluido invíscido. A equação de Burgers é conside-

rada uma equação da onda não-linear, onde cada ponto da frente de onda se propaga

com velocidade diferente. Uma consequência desta mudança na velocidade da onda

é a coalescência (ou dispersão) das características que leva a formação de ondas de

choque (ou rarefação) da Mecânica dos Fluidos. Outra vantagem de utilizar a equa-

ção de Burgers é que a mesma possui soluções analíticas exatas para diversos casos

onde aparecem ondas de choque e rarefação, podendo ser comparadas com resultados

obtidos numericamente.

A Equação de Burgers pode também ser escrita na forma conservativa como:

ut + f (u)x = 0 (4.2)

onde f (u) = u2

2 representa o fluxo invíscido.

A solução analítica exata da Equação de Burgers pode ser obtida usando o método

das características. Se a condição inicial for u(x, t = 0) = u0(x), a solução da equação

(4.1) ou (4.2) para t > 0 será obtida da equação transcendental:

u(x, t ) = u0(x −u(x, t )t ) (4.3)

onde x foi substituído por x −u(x, t )t na condição inicial u0(x). Este procedimento

51

pode ser utilizado para calcular qualquer solução desta equação.

Fig. 4.1: Deslocamento típico da descontinuidade na Equação de Burgers [5].

Para nossos testes com a equação de Burgers, a condições iniciais assumidas são:

u(x,0) = 0.5 se x = 0

u(x,0) = 0.0 se x > 0(4.4)

com as seguintes condições de contorno:

u(x = 0, t ) = 0.5

u(x = 1, t ) = 0.0(4.5)

Para comparar os resultados, foi implementado um esquema numérico utilizando

uma discretização espacial de segunda ordem TVD com limitador de fluxo, mesmo

procedimento utilizado na discretização dos fluxos invíscidos da equações governantes

do modelo hidrodinâmico, detalhados na seção 3.1.1.1. Para discretização temporal,

foi utilizada um esquema Runge-Kutta explícito de segunda ordem SSP (SSPRK(2,2)).

Porém na primeira análise foi utilizado um passo no tempo ∆t = 2.2∆x, onde não

temos a propriedade SSP, enquanto que na segunda análise, foi utilizado o passo no

tempo ∆t = 1.5∆x onde temos uma forte preservação de estabilidade (SSP), uma vez

que a barreira SSP é ∆t = 2.0∆x .

O esquema utilizado na primeira análise produz o resultado numérico mostrado na

figura 4.2, onde oscilações espúrias dominam na região do choque, mostrando uma

grande instabilidade numérica.

52

0.0 0.2 0.4 0.6 0.8 1.0

0.0

0.1

0.2

0.3

0.4

0.5

0.6

xHmL

Velo

cid

ade

HmsL

Solução Exata

Código RKH2,2L

Fig. 4.2: Comparação entre esquema de segunda ordem e solução exata da Equação deBurgers no tempo t = 2.0s.

Já a segunda análise, que respeita a barreira SSP, ofereceu os resultados mostrados

na figura 4.3. A solução numérica possui uma excelente concordância em relação a

solução analítica exata, o que serve de evidência da importância de utilizar esquemas

que possuem esta estabilidade numérica não-linear para controlar oscilações induzidas

por gradientes severos ou descontinuidades.

0.0 0.2 0.4 0.6 0.8 1.0

0.0

0.1

0.2

0.3

0.4

0.5

xHmL

Velo

cid

ade

HmsL

Solução Exata

Código SSPRKH2,2L

Fig. 4.3: Comparação entre esquema SSPRK(2,2) e solução exata da Equação de Bur-gers no tempo t = 2.0s.

Para finalizar o teste numérico utilizando a Equação de Burgers, foi feita uma aná-

lise da propagação da onda de choque no intervalo de 0 a 4s, mostrada na figura 4.5.

Verifica-se que as oscilações se acumulam ao longo do tempo quando se usa um passo

53

no tempo acima da barreira do SSP, que é de ∆t = 2.0∆x. Porém quando utilizamos

um passo no tempo abaixo desta barreira, verifica-se que resultado numérico descreve

precisamente o choque com excelente concordância com a solução analítica.

æææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææ

ààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààà

à

à

à

à

à

à

à

à

à

à

à

à

à

à

à

à

à

à

à

à

à

àààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààà

ìììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

òòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòò

ò

ò

òòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòò

ôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôô

ô

ô

ôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôô

ççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççç

ç

0.0 0.2 0.4 0.6 0.8 1.0

0.0

0.1

0.2

0.3

0.4

0.5

0.6

xHmL

Velo

cid

ade

HmsL

æ Não-SSP 1s

à Não-SSP 2s

ì Não-SSP 4s

ò SSP 1s

ô SSP 2s

ç SSP 4s

Fig. 4.4: Solução transiente produzida pelo RK(2,2) com passo no tempo ∆t = 1.5∆x(SSP) e ∆t = 2.2∆x (Não-SSP) para o intervalo de 0 a 4s.

4.1.2 Tubo de Choque

Suponha um tubo unidimensional contendo duas regiões preenchidas com um fluido

em repouso sujeito à diferentes temperaturas. Suponha ainda que as duas regiões es-

tão separadas por um diafragma rígido. Se o diafragma é instantaneamente removido,

o desbalanceamento de pressão causa um escoamento transiente unidimensional con-

tendo uma onda de choque que se move na direção de menor pressão, uma onda de

rarefação que se move na direção de maior pressão e uma linha de contato que separa

os fluidos de alta e baixa pressão. Este arranjo é conhecido como Tubo de Choque. O

escoamento no tubo de choque é mostrado na figura abaixo.

O escoamento no tubo de choque tem sempre velocidade inicial igual a zero. Re-

movendo esta restrição, o problema do tubo de choque torna-se o problema de Rie-

mann. Desta forma o tubo de choque é, na verdade, um caso especial do problema de

Riemann. Como o tubo de choque, o problema de Riemann também gera uma onda

de choque, uma onda de rarefação e uma linha de contato, contudo, diferentemente do

tubo de choque, uma ou duas dessas ondas podem estar ausentes.

54

Fig. 4.5: Representação esquemática do problema do tubo de choque e comportamentodas características ao longo do domínio [6]

O problema do tubo de choque é um caso interessante para testarmos o código

desenvolvido para este trabalho. Além do fato de possuir solução analítica exata, a

solução do tubo de choque é dita auto-preservativa, ou auto-similar. Isto é, a solução

se estende similarmente ao longo do espaço e do tempo mantendo sua forma, então

u(x, t1) e u(x, t2) são similares para quaisquer dois tempos t1 t2. Ou, dito de uma outra

forma, a solução depende apenas da simples variável xt e não de x e t separadamente.

Reduzindo o número de variáveis independentes, a auto-similaridade simplifica as téc-

nicas de solução e, frequentemente, leva a soluções analíticas, como no nosso caso do

tubo de choque. Outros exemplos de soluções auto-similares incluem problemas de

camada-limite laminar permanente, jatos, esteiras, e ondas de choque planares, cilín-

dricas e esféricas de Sedov [6].

A solução exata do tubo de choque é obtida separadamente para cada região do

domínio. Considerando inicialmente a região da onda de choque e definindo para a ve-

locidade à esquerda do choque como u2, e à direita como u1 = uR , podemos escrever:

c22

c21

= P2

P1

γ+1γ−1 + P2

P1

1+ γ+1γ−1

P2P1

(4.6)

55

u2 = u1 + c1

γ

P2P1

−1√γ+12γ

(P2P1

−1)+1

(4.7)

S = u1 + c1

√γ+1

2γ

(P2

P1−1

)+1 (4.8)

onde c1 e c2 representam a velocidade do som nas regiões 1 e 2, respectivamente, e

S representa a entropia. Consideremos agora a linha de contato. Definindo para a

velocidade à esquerda e à direita do contato como u3 e u2, podemos escrever:

u2 = u3 (4.9)

P2 = P3 (4.10)

Finalmente, consideramos uma onda de rarefação. Definindo para a velocidade à

esquerda e à direita da expansão como u4 = uL e u3, respectivamente, e c4, a velocidade

do som na região 4, temos:

u(x, t ) = 2

γ+1

(x

t+ γ−1

2u4 + c4

)(4.11)

c(x, t ) = u(x, t )− x

t= 2

γ+1

(x

t+ γ−1

2u4 + c4

)− x

t(4.12)

P = P4

(c

c4

)2γ/(γ−1)

(4.13)

Combinando os resultados mostrados para a região de choque, contato e rarefação

e após uma série de manipulações, chegamos a uma equação para a razão de pressões

P2/P1 através do choque em função da razão de pressões conhecidas P4/P1 = PL/PR .

u3 + 2c3

γ−1= u4 + 2c4

γ−1(4.14)

56

u3 = u4 + 2c4

γ−1

[1−

(P3

P4

)(γ−1)/2γ]

(4.15)

u2 = u4 + 2c4

γ−1

[1−

(P2

P4

)(γ−1)/2γ]= u4 + 2c4

γ−1

[1−

(P1

P4

P2

P1

)(γ−1)/2γ]

(4.16)

Resolvendo a equação 4.16 para P4/P1 temos:

P4

P1= P2

P1

[1+ γ−1

2c4(u4 −u2)

]−2γ/(γ−1)

(4.17)

onde finalmente chegamos no resultado desejado:

P4

P1= P2

P1

1+ γ−1

2c4

c4 −u1 − c1

γ

P2P1

−1√γ+12γ ( P2

P1−1)+1

(4.18)

A equação obtida (4.18) é uma equação implícita não-linear, que pode ser resolvida

utilizando o método de Newton ou da Secante. Desta forma dispomos da solução

analítica exata para comparar com os resultados numéricos obtidos com o código.

A malha utilizada no código de segunda ordem TVD no espaço e de segunda ordem

SSP no tempo possui 401 pontos e o passo no tempo utilizado foi ∆t = 1×10−6s. Fo-

ram feitas comparações entre os resultados oferecidos pelo código e a solução analítica

exata para a velocidade (figura 4.6), pressão (figura 4.7), massa específica (figura 4.8),

número de Mach (figura 4.9) e velocidade do som (figura 4.10). Os resultados produzi-

dos possuem uma excelente concordância para todas as variáveis, tanto nas regiões de

choque, quanto na rarefação e na linha de contato. Não foi verificada nenhuma espécie

de oscilação espúria ou instabilidade numérica.

A análise de acurácia dos resultados obtidos para as variáveis dependentes feita

pela superposição entre soluções numéricas e analíticas é limitada e não permite uma

avaliação criteriosa. Desta forma, se faz necessário um estudo mais criterioso, conhe-

cido como análise de ordem. No intuito de garantir que os esquemas implementados

57

ææææææææææææææææææææææææææææææææ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æææææææææææææææææææææææææææææ

æææææææææææææææææææææææààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààà

ààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààà

à

à

à

à

à

àààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààà

-10 -5 0 5 10

0

50

100

150

200

250

xHmL

Velo

cid

ad

eHm

sL

æ Solução Analítica à Resultado do Código

Fig. 4.6: Comparação entre resultado numérico e solução analítica para velocidade notempo t = 10−2s.


æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ


æææææææææææææææææææææææ

ààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààà

à

à

à

ààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààà

-10 -5 0 5 10

0

20 000

40 000

60 000

80 000

100 000

xHmL

Press

ão

HPa

L


Fig. 4.7: Comparação entre resultado numérico e solução analítica para pressão notempo t = 10−2s.


æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

ææææææææææææææææ

æææææææææææææ


ààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààà

à

à

àààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààà

-10 -5 0 5 10

0.0

0.2

0.4

0.6

0.8

1.0

xHmL

Ma

ssa

Esp

ecíf

ica

Hkg

m3

L


Fig. 4.8: Comparação entre resultado numérico e solução analítica para massa especí-fica no tempo t = 10−2s.

58


æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ



æææææææææææææææææææææææààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààà

àààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààà

à

àààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààà

à

à

à

à

à


-10 -5 0 5 10

0.0

0.2

0.4

0.6

0.8

xHmL

Nú

mero

de

Ma

ch


Fig. 4.9: Comparação entre resultado numérico e solução analítica para o número deMach no tempo t = 10−2s.


æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ




àààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààà

à

à

à

à

à

à

à

à

à

à

à

à

à

à

à

ààààààààààà

àààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààà

à

à

à

à

à


-10 -5 0 5 10

320

340

360

380

400

xHmL

Velo

cid

ad

ed

oso

mHm

sL


Fig. 4.10: Comparação entre resultado numérico e solução analítica para a velocidadedo som no tempo t = 10−2s.

atingem suas respectivas ordens de erro teóricas, faz-se necessário uma análise do

comportamento do erro com a variação de malha. Com este objetivo, foram realizadas

análises do erro variando a malha espacial e a malha temporal.

A análise de ordem foi feita utilizando o problema do Tubo do Choque, com um

tempo final de simulação de 4×10−8 segundos. O motivo de ser escolher um tempo fí-

sico final curto para a análise é minimizar os erros acumulados a cada passo no tempo,

de modo que o erro dominante seja espacial e não temporal. O erro foi calculado a par-

tir da diferença dos resultados numéricos para massa específica obtidos entre malhas

de diferentes espaçamentos.

Para a malha espacial, foram utilizadas quatro malhas: 501 pontos, 1001 pontos,

59

2001 pontos e 4001 pontos. A partir das quatro malhas foram calculados os erros

absolutos para três malhas, sempre subtraindo o resultado da malha mais grosseira pelo

resultado da mais refinada subsequente. Para a análise do erro espacial, foi utilizado

um passo no tempo ∆t = 4×10−9s. A figura 4.11 mostra o erro absoluto para as três

malhas utilizadas.

æ æ æ æ æ æ æ æ æ æ æ æ

æ

æ

æ

æ

æ

æ

æ

æ æ

æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æà à à à à à à à à à à à à à à à à à à à à à à à à à à

à à

à

à

à

à

à

à

à

à à

à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à àìììììììììììììììììììììììììììììììììììììììììììììììì

ì

ìììììì

ì

ì

ì

ìì

ìì

ì

ì

ì

ì

ì

ìì

ìì

ì

ì

ì

ìììììì

ì

ììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììì

-0.4 -0.2 0.0 0.2 0.4

0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

xHmL

Erro

na

ma

lha

esp

acia

lH1

0K

gm

3L

æ Erro 501 pts

à Erro 1001 pts

ì Erro 2001 pts

Fig. 4.11: Comparação entre erros obtidos para malhas de 501 pontos, 1001 pontos e2001 pontos espaciais, problema do tubo de choque com tempo físico finalt f = 4×10−9s.

Utilizando os resultados obtidos numericamente para o erro na discretização espa-

cial do esquema para diferentes ∆x, é possível calcular a ordem de erro espacial. A

partir da figura 4.12 verificamos que a ordem do erro espacial O(∆x2) se confirma para

o código implementado, uma vez que a inclinação da reta que representa o ordem do

código implementado (linha azul) é próxima à inclinação da reta de ordem dois (linha

amarela).

Para análise temporal foram utilizadas cinco diferentes malhas, com o seguintes

passos no tempo: ∆t = 8×10−7s, ∆t = 4×10−7s, ∆t = 2×10−7s, ∆t = 1×10−7s, ∆t = 5×10−8s. O erro absoluto foi calculado subtraindo a solução da malha mais grosseira pela

solução da malha mais refinada subsequente, desta forma, a partir das cinco malhas

utilizadas foi possível calcular o erro absoluto para quatro delas. O comportamento do

erro se verificou conforme o previsto, mostrado na figura 4.13.

Usando o mesmo procedimento adotado na avaliação de ordem espacial, também

confirmamos a ordem de erro temporal O(∆t 2) para o código implementado. A figura

60

0.010 0.020 0.0300.015

0.10

1.00

0.50

0.20

2.00

0.30

3.00

0.15

1.50

0.70

DxHmL

Erro

na

ma

lha

esp

acia

lH1

0K

gm

3L

Ordem numérica

Ordem teórica 1

Ordem teórica 2

Fig. 4.12: Ordem de erro espacial do código e ordens teóricas para comparação, pro-blema do tubo de choque com tempo físico final t f = 5×10−9s.


æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

ææ

ææææææææææææææàààààààààààààààààààààààààà

àà

à

à

à

àà

à

à

à

àà

à

à

à

àà

ààààààààààààààààààììììììììììììììììììììììììììììì

ìììì

ì

ì

ìììì

ììììììììììììììììììììììòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòò

òòòòòòòòòòòòòòòòòòòòòòòòòòò

-0.2 -0.1 0.0 0.1 0.2

0.0000

0.0002

0.0004

0.0006

0.0008

0.0010

0.0012

xHmL

Erro

na

Ma

lha

Tem

po

ral

HKg

m3

L

æ Erro Dt=8.´10-7

s à Erro Dt=4.´10-7

s

ì Erro Dt=2.´10-7

s ò Erro Dt=1.´10-7

s

Fig. 4.13: Comparação entre erros obtidos diferentes passos no tempo, problema dotubo de choque com tempo físico final t f = 4×10−9s.

4.14 compara a ordem de erro obtida com as curvas teóricas de ordem 1 e 2, onde cla-

ramente verificamos que a curva numérica possui a mesma inclinação da curva teórica

O(∆t 2).

4.1.3 Transferência de Calor Termoacústica

Com a versão invíscida do código verificada e validada, os termos viscosos são

incluídos para simulação do problema de transferência de calor termoacústica. Antes

61

0.500.20 0.30 0.70

0.00010

0.00050

0.00020

0.00030

0.00015

0.00070

DtH10 sL

Erro

na

ma

lha

tem

po

ra

lHK

gm

3L

Ordem numérica

Ordem teórica 1

Ordem teórica 2

Fig. 4.14: Comparação entre erros obtidos diferentes passos no tempo para o problemado tubo de choque com tempo físico final t f = 4×10−9s.

de incluir o efeito de impedância das paredes, a soluçao gerada pelo novo código é

comparada com resultados expostos em artigos/publicações científicas referenciadas.

As ondas termoacústicas geradas por um aumento instantâneo de temperatura fo-

ram analisadas inicialmente. A variação de temperatura na parede esquerda foi:

TL(t ) =

T0 para t ≤ 0

T0 +∆T para t > 0(4.19)

onde ∆T é o aumento de temperatura imposto na parede.

A figura 4.15 mostra a comparação para quatro ondas termoacústicas obtidas nos

tempos físicos t/ta = 0,25, t/ta = 1,0, t/ta = 1,5 e t/ta = 2,0. Onde ta = L/c é o

tempo acústico, que representa fisicamente o tempo necessário para a onda de pressão,

que viaja na velocidade do som c, atravessar o comprimento da cavidade L. Para este

caso, o valor do tempo acústico calculado foi ta = 2,83µs. Os resultados (linhas cheias)

são mostrados em conjunto com os resultados obtidos por Aktas [7] (linhas tracejadas).

O gráfico mostra uma boa concordância entre os resultados obtidos e os expostos

na referência [7]. Algumas pequenas diferenças já eram esperadas, uma vez que as

condições de contorno nas paredes não são as mesmas. É importante ressaltar que as

características peculiares da forma da onda, frente de onda aguda e a cauda alongada,

é similarmente capturada, mostrando que o método utilizado descreve a onda termoa-

62

æ ææ

ææ

æ

æ

æ

æ

æ

æ

æ

æ

æ æ æ æ æ

àà

àà

à

à

à

à

à

ì ì ì ì

ì

ì

ì

ìì

ì

ìì

ìì ì

ì ì ì ì

ò

ò

ò

ò

ò

òò

òò ò ò ò

0.0 0.2 0.4 0.6 0.8 1.0

102 000

103 000

104 000

105 000

xL

Press

ão

HPa

Lt=0.25ta

t=1.00ta

t=1.50ta

t=2.00ta

æ Aktas 0.25ta

à Aktas 1.00ta

ì Aktas 1.50ta

ò Aktas 2.00ta

Fig. 4.15: Comparação entre os resultados obtidos com o código e os mostrados emAktas [7], para quatro ondas termoacústicas nos tempos especificados nalegenda.

cústica corretamente. Para a análise mostrada na figura 4.15 foi utilizada TL = 400K,

pressão inicial na cavidade P0 = 101.325Pa, temperatura inicial T0 = 300K, compri-

mento da cavidade L = 1,0mm, uma malha de 701 pontos espaciais e um passo no

tempo ∆t = 1×10−9s.

æ

ææ æ

ææ

æ

æ

æ

æ

æ

æ

æ æ

à

à

àà

à

àà

àà

à

à

à

à

ì

ì

ìì ì ì ì

ì

ìì

ì

ì

ì

ì

ì

ìì

ì

ò

ò

ò

ò

ò

òò

òò

ò

ò

ò

0.0 0.2 0.4 0.6 0.8 1.0

300.0

300.5

301.0

301.5

302.0

302.5

303.0

xL

Tem

pera

tu

ra

HKL

Código 0.25ta

Código 1.00ta

Código 1.50ta

Código 2.00ta

æ Aktas 0.25ta

à Aktas 1.00ta

ì Aktas 1.50ta

ò Aktas 2.00ta

Fig. 4.16: Comparação entre os resultados obtidos com o código e os mostrados emAktas [7] para temperatura ao longo da cavidade nos tempos especificadosna legenda.

Variações na temperatura do gás ao longo da cavidade são mostradas na figura

4.16. O perfil de temperatura mostra uma tendência similar ao resultado mostrado para

63

a pressão. A figura demonstra claramente o efeito da onda termoacústica no transporte

de calor da extremidade para o interior da cavidade. Nas regiões cada vez mais pró-

ximas da parede esquerda, onde TL = 400K, a temperatura é significativamente mais

alta comparada a região mais afastada devido a presença da camada limite térmica. Na

teoria clássica de transferência de calor, a região fora da camada limite térmica tem

um temperatura igual a temperatura inicial T0. Este resultado ilustra, conforme já ci-

tamos, que o efeito termoacústico transfere energia da camada limite para fora dela,

promovendo um rápido aquecimento do gás nesta região e pode ser considerado um

mecanismo adicional de transporte de calor.

æ

æ ææ

æ

æ

æ

æ

æ

æ æ

à

àà

àà

àà

à

à

à

à

ì

ìììì ì ì

ì

ì

ìì

ì

ì

ìì

ì ì ì

ì

ò

ò

ò

òòò

òò

ò òò

0.0 0.2 0.4 0.6 0.8 1.01.13

1.14

1.15

1.16

1.17

xL

Mass

aesp

ecíf

ica

Hkg

m3

L

Código 0.25ta

Código 1.00ta

Código 1.50ta

Código 2.00ta

æ Aktas 0.25ta

à Aktas 1.00ta

ì Aktas 1.50ta

ò Aktas 2.00ta

Fig. 4.17: Comparação entre os resultados obtidos com o código e os mostrados emAktas [7] para a massa específica ao longo da cavidade nos tempos especi-ficados na legenda.

O comportamento da massa específica é mostrado na figura 4.17. Próximo a parede

esquerda, a massa específica do gás é extremamente baixa, devido a expansão resul-

tante do aumento instantâneo de temperatura. Um efeito semelhante ocorre na parede

da direita, onde ocorre um aumento da massa específica devido ao resfriamento do gás

dentro da camada limite térmica induzida pela menor temperatura desta parede. As

magnitudes são menores nesta figura pois esta camada limite começa a ser formada

apenas após a primeira onda termoacústica incidir sobre esta parede.

O resultado para velocidade é mostrado na figura 4.18. É possível verificar que o

64



ææææææææææææææææææææææ

æææææææææææææææææ


ææææææææææææææææææææææææææææææææææææææææææææææææææææ

æ

æ

æ

æ

æ

æ

æææææææææ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

ææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææàààà

ààààààààààààààààààà

àààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààà

àààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààà

àààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààà

ààààààààààààààààààààààààààààààààààààààààààààààààààààààààààà

àààààààààààààààààààààààààààààààààààààààà

àààààààààààààààààààààààààààà

ààààààààààààààààààààà

ààààààààààààààà

ààààààààààààà

àààààààààààààààààààààààààààààààà

àààà

à

à

à

à

à

à

à

ààààààààààààààààìììììì

ììììììììììììììììììììì

ììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììì

ìììììì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ììììììììììììììììììì

ìììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììì

ììììììììììììì

ììììììììììììììììì

ìììììììììììììììììììììì

ìììììììììììììììììììììììììììì

ìììììììììììììììììììììììììììììììììììì

ììììììììììììììììììììììììììììììììììììììììììììììì

ììììììììììììììììììììììììììììììììììììììììììììììììììììììììì

ììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììòòòòòòò

òòòòòòòòòò

òòòòòò

òòòòòòòòòòòòòòòòòòò

òòòòòòòòòòòòòòòòòòòòòòòòòòòòòòò

òòòòòòòòòòòòò

òòòòòòòòòòòòòòòòòò

òòòòòòòòòòòòòòòòòòòòòòòò

òòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòò

òòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòò

òòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòò

òòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòò

òòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòò

òòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòò

òòòòòòòòòòòòòòòòò

0.0 0.2 0.4 0.6 0.8 1.0

-2

0

2

4

6

xL

Velo

cid

ad

eHm

sL æ t=0.25ta

à t=1.00ta

ì t=1.50ta

ò t=2.00ta

Fig. 4.18: Resultados obtidos para a velocidade ao longo da cavidade nos tempos es-pecificados na legenda.

efeito pistão gera um campo de velocidades não desprezível dentro da cavidade.

A variação temporal da pressão no centro da cavidade (x/L = 0.5) é também com-

parada na figura 4.19. A pressão no centro permanece constante e igual ao seu valor

inicial até a chegada da primeira onda termoacústica. Então, a pressão cresce repenti-

namente com o passar da primeira frente de onda aguda, mas decresce gradualmente

com o passar da primeira cauda. Um comportamento similar é observado todas as

vezes que a onda termoacústica passa por esta região, porém as oscilações de pres-

são eventualmente decaem devido a dissipação viscosa e difusão de quantidade de

movimento e calor dentro da cavidade. Nosso trabalho tem como objetivo modelar

condições de contorno de pressão, porém, para efeitos comparativos, foi modelada a

mesma condição de contorno utilizada por [7]. A solução utilizando a condição de con-

torno de pressão (CC(ρ)) mostrada com linhas vermelhas nesta figura possui uma boa

concordância com o obtido por Aktas [7], mostrada com linhas verdes. Com relação

à condição de contorno de pressão (CC(P)), as diferenças aumentam a cada reflexão,

devido às diferentes aproximações utilizadas para o contorno.

65

0 5 10 15

1.000

1.005

1.010

1.015

1.020

1.025

tta

PP

0 Solução CCHPLSolução CCHΡLAktas

Fig. 4.19: Comportamento da pressão ao longo do tempo no centro da cavidade, resul-tado comparado com o obtido por Aktas [7].

4.1.4 Efeito das condições de contorno para a pressão

Para os resultados anteriores, a impedância acústica foi desprezada (Z = 1), assu-

mindo as paredes da cavidade como perfeitamente rígidas e sem perda alguma durante

a reflexão da onda termoacústica. Esse valor foi assumido este valor devido ao fato

dos casos analisados anteriormente terem por objetivo validar o código implementado

para simulação de ondas termoacústicas. Uma vez que os artigos usados para com-

paração não utilizam impedância acústica na modelagem, a parede foi considerava

perfeitamente rígida. Contudo, apesar deles não incluírem explicitamente perdas por

impedância da parede, condições de contorno mais sofisticadas incluem perdas por

outros efeitos [19, 22].

Infelizmente, os dados experimentais existentes na literatura foram gerados em ca-

vidades de tamanho relativamente grandes, fazendo com que as malhas necessárias

para simular a mesmas condições numericamente utilizem um número muito elevado

de pontos para os recursos computacionais disponíveis, apesar do código estar parale-

lizado (ver apêndice A). Desta forma, para ilustrar o efeito da impedância na simulação

do modelo hidrodinâmico, utilizamos uma cavidade fictícia com apenas 1mm de com-

primento. É por esta razão que as oscilações observadas na figura 4.20 possuem uma

frequência tão alta. Esta também é a principal razão pela qual a maioria absoluta dos

66

artigos publicados na literatura estudando ondas termoacústicas são restritos a proble-

mas unidimensionais. Por isso o modelo termodinâmico é tão importante, uma vez que

ele reduz o tempo de simulação drasticamente quando incorporado aos modelos hidro-

dinâmicos, tornando desnecessário simular as escalas acústicas. Consequentemente,

a descoberta da importância que a impedância da parede tem nas simulações somente

terá aplicação prática se este efeito por incorporado no modelo termodinâmico. Isto

será feito na próxima seção.

0 100 200 300 400 500 600

1.00

1.02

1.04

1.06

1.08

1.10

1.12

1.14

tta

PP

0

Z=1.000000

Z=0.999999

Z=0.950000

Z=0.900000

Fig. 4.20: Resultados da evolução da pressão ao longo do tempo com diferentes valoresde impedância acústica de parede Z .

A figura 4.20 mostra o resultado evolução da pressão ao longo do tempo utilizando

quatro valores de impedância acústica, mostrando que a modelagem da impedância

acústica desenvolvida neste trabalho tem o efeito desejado de simular uma absorção

da onda termoacústica no contato com a parede da cavidade. A medida que o valor de

Z decresce, ocorre consequentemente uma maior perda de pressão a cada contato da

onda termoacústica com a parede da cavidade. A partir do gráfico podemos verificar

que a medida que aumentamos esta perda, aumenta a tendência da pressão no sistema

retornar para pressão inicial P0. Este fenômeno é coerente com a observação expe-

rimental, conforme gráfico mostrado na figura 1.4 discutida no capítulo introdutório

deste trabalho.

67

4.1.5 Efeito das condições de contorno para temperatura

Esta seção tem por objetivo analisar o efeito da taxa de variação da temperatura

da parede na geração e transmissão das ondas termoacústicas. Na prática é impossível

variar a temperatura da parede da cavidade instantaneamente. No máximo, através de

um circuito com grande eficiência, é possível atingir uma diferença de 100 graus em

alguns milissegundos. Para isto, serão utilizados quatro casos testes, onde a parede da

cavidade é aquecida de acordo com a expressão proposta por Shen e Zhang [24].

TL = f (t ) = Ti +∆T [1−exp( −t

H ta

)] (4.20)

H = 0.002

H = 0.020

H = 0.200

H = 2.000

A cavidade é inicialmente preenchida com nitrogênio a temperatura T0 = 300K e

pressão P0 = 101325Pa. A figura 4.21 mostra as diferentes taxas de mudança da tempe-

ratura na fronteira esquerda da cavidade. Os resultados desta seção foram produzidos

utilizando uma malha de 501 pontos numa cavidade de 1mm utilizando um passo no

tempo ∆t = 1×10−9s.

æ

æææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææ

à

à

à

à

à

ààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààà

ì

ì

ì

ì

ì

ì

ì

ì

ììììì

ììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììì

ò

ò

ò

ò

ò

ò

ò

ò

ò

ò

ò

ò

òò

òò

òò

òò

òò

òò

òò

òòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòò

0 2 4 6 8 10

0.0

0.2

0.4

0.6

0.8

1.0

tta

T-

T0

T1

-T

0

æ H=0.002

à H=0.02

ì H=0.2

ò H=2.0

Fig. 4.21: Tempo para atingir temperatura final na parede esquerda para os quatro ca-sos analisados.

68

A partir da análise das figuras 4.21 e 4.22, podemos a verificar diferentes formas

de onda termoacústica. Para H = 0.002 e H = 0.020 observamos que a onda mantém

o formato observado no aquecimento instantâneo. Já com H = 0.200, não observamos

mais a cauda da onda. A principal diferença deste caso em relação aos dois primeiros

é que a temperatura final da parede só é atingida quando t ∼ ta , enquanto nos casos

anteriores isso ocorre com t ¿ ta . No último caso (H = 2.000) a parede ainda está

sendo aquecida e, por isso, a amplitude da onda é mínima. De qualquer forma, ela

ainda tem o mesmo formato de onda mostrada com H = 0.200.

0.0 0.2 0.4 0.6 0.8 1.0

1.000

1.005

1.010

1.015

1.020

1.025

1.030

xL

PP

0

H=0.002

H=0.02

H=0.2

H=2.0

Fig. 4.22: Pressão na cavidade no tempo t = 2×10−5s para diferentes tipos de aqueci-mento.

A figura 4.23 mostra a evolução da pressão na cavidade em tempos sucessivos para

H = 0.020 e H = 2.000. É possível visualizar também as diferentes características

das ondas geradas. A onda gerada pelo aquecimento com H = 0.020 possui maior

amplitude para os quatro tempos avaliados. Com H = 2.000 não ocorre a formação

de uma cauda e a onda é praticamente quadrada. Para estes dois casos ocorre uma

aumento da pressão média na cavidade a medida que a onda percorre a cavidade ao

longo do tempo sendo refletida pelas paredes, porém, como mostrado na figura, o

aquecimento rápido gera uma maior pressão na cavidade para um mesmo tempo.

A pressão no ponto médio da cavidade, para cada um dos perfis de temperatura

analisados através dos quatro casos teste, são mostrados na figura 4.24. A partir dessa

figura é possível verificar que o aquecimento instantâneo (linha azul - H = 0.002) gera

uma onda termoacústica de maior amplitude, porém a mesma é amortecida até a che-

69


ææææææææ

ææææææ

ææææææææææææææææææææææææææææææ

æ

æ

æ

æ

ææææææ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

æ

ææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææ

ààààààààààààààààààààààààà

àààààààààààààààààà

ààààààààààààààà

àààààààààààà

àààààààààà

ààààààààà

ààààààà

àààààà

àààààààààààààààààààààààààààà

à

à

àààààààà

àà

à

à

à

à

à

à

à

à

à

à

à

à

àààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààà

ììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììì

ìììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììì

ìììììììììììììììììììììììììììììììììììììììì

ìììììììììììììì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ì

ìììììììììììììì

ìììììììììììììììììììììììììììììììììììììììììììììììììììììììììììì

òòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòò

òòòòòòòòòòòòòòòòòòòòòò

òòòòòòòòò

ò

ò

ò

ò

ò

ò

ò

ò

ò

ò

ò

ò

ò

òòòòòòòòòòòòòòòò

òòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòò

ôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôô

çççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççççç

ááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááá

ááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááááá

ííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííí

íííííííííííííííííííííííííííííííííííííííííííííííí

ííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííííí

0.0 0.2 0.4 0.6 0.8 1.0

1.000

1.005

1.010

1.015

1.020

1.025

xL

PP

0

æ t=0.7ΜsHH=0.02Là t=1.4ΜsHH=0.02Lì t=3.6ΜsHH=0.02Lò t=4.2ΜsHH=0.02Lô t=0.7ΜsHH=2.0Lç t=1.4ΜsHH=2.0Lá t=3.6ΜsHH=2.0Lí t=4.2ΜsHH=2.0L

Fig. 4.23: Ondas termoacústicas geradas pelo aquecimento com H = 0.020 e H =2.000.

gada de um próximo pico de pressão. Os casos com H = 0.020 e H = 0.200 temos

sucessivas diminuições da amplitude máxima da onda, porém, a pressão média da

cavidade se estabiliza no mesmo valor. O único caso apresenta um comportamento

diferente em relação a pressão é o aquecimento com H = 2.000, onde verificamos que

não ocorre pico e a pressão aumenta a uma taxa quase constante.

0 1 2 3 4 5 6 7

1.000

1.005

1.010

1.015

1.020

1.025

1.030

tta

PP

0

H=0.002

H=0.02

H=0.2

H=2.0

Fig. 4.24: Variação da pressão no ponto médio da cavidade para os quatro casos anali-sados.

70

4.2 Modelos Termodinâmicos

O modelo termodinâmico clássico, por ser um modelo simplificado que possui

solução analítica, oferece resultados com um baixo custo computacional. A versão

simplificada deste modelo desenvolvida durante este trabalho, que inclui a impedância

acústica de parede no termo fonte de compressão adiabática da equação do calor, foi

implementado no ambiente de computação simbólica Mathematica. Mesmo utilizando

um número de termos alto, de forma a garantir a convergência da solução em série, o

Mathematica efetua os cálculos do modelo e resolve o sistema de equações acopladas

em algumas dezenas de segundos. Isto, sem dúvida, consiste numa grande vantagem,

uma vez que o modelo hidrodinâmico estudado necessita de um tempo da ordem de

dias para disponibilizar um resultado para nosso problema de ondas termoacústicas

numa cavidade. Desta forma, comparações com dados experimentais são agora possí-

veis com os recursos computacionais disponíveis.

4.2.1 Verificação Numérica

Antes de avaliar os resultados oferecidos pelo modelo termodinâmico, o número

ótimo de termos N das soluções em série do modelo foi avaliado. Uma vez que,

conforme citado no capítulo de modelos matemáticos, um filtro foi utilizado na solução

transformada pelo método GITT no intuito de se acelerar a convergência, é esperado

um número relativamente baixo de termos necessários para se garantir a convergência.

Foram feitas análises utilizando um número de termos N = 10, 20, 30, 40 e 50.

Conforme mostrado na figura 4.25, um número de termos N = 30, já garante um erro

absoluto na pressão da ordem de 10−3, o que já é consideravelmente baixo e sufici-

ente para nosso estudo. Uma vez que os fluidos aqui estudados são considerados gases

perfeitos e um filtro foi utilizado para se obter a solução transformada, um número

alto de termo realmente não se faz necessário. Para outros estudos, utilizando fluidos

supercríticos por exemplo, faz necessário um número maior de termos, conforme ex-

posto por Alves [29]. Como o filtro é uma solução de regime permanente, ele acelera

mais a convergência para tempos elevados, quando a solução está próxima do regime

71

permanente. Consequentemente, ele não é tão eficiente para tempos curtos, por isso

esta análise de convergência precisa ser feita para tempos curtos. Este aspecto é confir-

mado pela figura 4.25, onde podemos verificar uma maior amplitude do erro absoluto

para tempos baixos, enquanto que para tempos longos, o erro absoluto é baixo, mesmo

para um número pequeno de termos das soluções em série. A análise foi realizada

utilizando como gás o ar padrão (γ= 1.402) e o valor da impedância acústica utilizado

foi Z = 0.9999.

æ

æ

æ

æ

æ

æ

æ

æ

ææææææææ

ææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææ

à

à

à

à

àààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààààà

ì

ì

ììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììììò

òòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòòô

ôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôôô

0.00000 0.00005 0.00010 0.00015 0.00020

0.0000

0.0002

0.0004

0.0006

0.0008

0.0010

TempoHsL

DP

P0

æ N=10

à N=20

ì N=30

ò N=40

ô N=50

Fig. 4.25: Resultado do erro absoluto estimado P/P0 da solução com séries usandoN = 10 , 20, 30, 40 e 50.

4.2.2 Análise do Efeito de Impedância das Paredes

Tendo exposto o estudo a respeito do número de termos a serem utilizados nas so-

luções em série do modelo termodinâmico, já é possível analisar o comportamento das

soluções oferecidas pelo modelo. Foram utilizados três experimentos com geração e

transmissão de ondas termoacústicas com fluidos compressíveis existentes na litera-

tura, no intuito de comparar com as soluções obtidas em nossa simulação do modelo

termodinâmico incluindo impedância acústica.

Em seu experimento, Brown [3] utilizou um tubo plástico montado verticalmente,

de comprimento L = 0,62m, com diâmetro interno de 0,02m, fechado nas extremida-

des, preenchido inicialmente com nitrogênio e aquecido a partir de baixo com uma

72

pequena folha fina de metal sustentada por uma seção de isolamento cerâmico. A

folha metálica foi aquecida rapidamente através um circuito RC. Após um segundo o

circuito foi desligado, deixando o filme resfriar perdendo calor para o ambiente interno

do tubo. A onda de pressão induzida foi medida através de um microfone localizado

a uma distância de 0,305m da folha metálica. A temperatura máxima obtida na folha

metálica foi de T1 = 480K, o que significa um ∆T = 180K, uma vez que a tempera-

tura inicial foi de T0 = 300K. A folha metálica atingiu esta temperatura num tempo

∆t = 0.0001s.

Foi observado que o tempo do crescimento e decaimento da pressão dentro da ca-

vidade previsto pelo modelo termodinâmico e o obtido experimentalmente estão defa-

sados no tempo por um fator Ω= 1/100. Conforme pode ser visto no gráfico da direita

na figura 4.26, o crescimento da pressão e posterior decaimento ocorrem em 0.025s no

experimento. Já os resultados obtidos pelo modelo termodinâmico prevêem a ocorrên-

cia do fenômeno em 2.5s, como mostrado no gráfico da esquerda na figura 4.26. No

intuito de possibilitar uma comparação qualitativa, a escala de tempo do modelo ter-

modinâmico foi ajustada multiplicando o tempo pelo fatorΩ. Utilizando esta correção,

os resultados são coerentes com os dados experimentais existentes, mostrando que esta

diferença na escala de tempo têm origem no modelo termodinâmico tradicional.

Fig. 4.26: Comportamento da pressão ao longo do tempo previsto pelo modelo termo-dinâmico e o resultado experimental.

O código termodinâmico desenvolvido aqui, assume condições de contorno de tem-

peratura de primeiro tipo, isto é, temperaturas prescritas nas paredes. Desta forma, no

73

intuito de simular o experimento realizado por Brown [3] relatado acima, nós assu-

miremos que o aquecimento foi instantâneo e a parede de baixo do cilindro estava a

temperatura T1.

Modelo Termodinâmico Z=0.9926Resultado ExperimentalModelo Termodinâmico Z=1.0

0.000 0.005 0.010 0.015 0.020 0.025Tempo H sL0

50

100

150

200

250Pressão H Pa L

0.015 0.025Tempo H s L

500

1000Pressão H Pa L

Pressão Termodinâmica

Fig. 4.27: Resultados experimentais da pressão medidos por Brown [3] utilizando ni-trogênio e solução do código termodinâmico.

No gráfico 4.27, podemos verificar que, uma vez que a pressão é assumida uni-

forme ao longo da cavidade para o modelo termodinâmico, não é possível descrever

as oscilações de alta frequência observadas experimentalmente. A pressão descrita no

modelo consiste em uma pressão média dentro da cavidade. Como já foi dito anteri-

ormente, as oscilações mostradas na curva experimental ocorrem devido ao aumento

de pressão na região do microfone durante a passagem da onda termoacústica. Após

a passagem é verificada uma queda, e um novo aumento só ocorre após a reflexão da

onda termoacústica na parede.

O resultados obtidos para pressão no experimento foram reproduzidos apenas até

o tempo t = 0.029s, pois não haviam dados para tempos maiores. O aspecto mais

importante a ser observado na figura 4.27 é que o modelo termodinâmico com impe-

dância acústica previu corretamente o decaimento da pressão e a amplitude da pressão

média dentro da cavidade. O valor da impedância utilizado para gerar o gráfico foi

Z = 0.99261101, o que significa uma reflexão de 99,261101% da amplitude de onda

termoacústica incidente contra as paredes da cavidade. Este valor foi obtido manual-

74

mente, sem o uso de uma rotina de otimização.

O segundo gráfico mostrado na figura 4.27 representa a pressão termodinâmica

(Z = 1) dentro da cavidade. Isto é, a pressão transiente ao longo da cavidade caso

o efeito de impedância acústica da parede não fosse computado. Pode-se verificar

claramente que a pressão na cavidade, sem a impedância, sobe para aproximadamente

1000Pa no mesmo intervalo de tempo do experimento.


0.000 0.002 0.004 0.006 0.008 0.010 0.012 0.014Tempo H sL0

50

100

150

200

250Pressão H Pa L

0.006 0.014Tempo H s L

250

500Pressão H Pa L


Fig. 4.28: Resultados experimentais da pressão medidos por Brown [3] utilizando arpadrão e solução do código termodinâmico.

Brown [3], em seu trabalho, relatou também um segundo experimento, similar ao

anterior, porém utilizando o ar padrão como fluido preenchendo o tubo. O resultado

do experimento, bem como nossos resultados teóricos estão mostrados na figura 4.28.

O resultados obtidos para pressão no experimento foram gerados somente até o tempo

t = 0.014s.

Conforme mostrado na figura 4.28, utilizando um valor para impedância acústica

no modelo termodinâmico de Z = 0.99281101, a previsão para a pressão média com a

cavidade contendo ar também possui uma boa concordância com o obtido experimen-

talmente por Brown [3]. A exemplo do que ocorreu no primeiro experimento mostrado,

tanto o decaimento da pressão, quanto a amplitude de pressão média dentro da cavi-

dade foram corretamente previstos. Na figura também é mostrado o gráfico da pressão

termodinâmica dentro da cavidade caso a impedância acústica não fosse incluída no

75

modelo, mostrando a importância da inclusão da impedância das paredes no modelo.

Lin e Farouk [8] também realizaram um experimento com ondas termoacústicas

em uma cavidade. Foi utilizado um tubo plástico de diâmetro interno de 38mm e

comprimento L = 201mm, isolado termicamente nas extremidades com discos espessos

de mica. Uma lâmina de níquel foi utilizada para efetuar o aquecimento do fluido

confinado dentro do tubo através de um circuito RC. O aquecimento neste experimento

atingiu uma temperatura máxima de 3430C, porém, a taxa de aquecimento não foi alta

e a temperatura máxima foi atingida após um intervalo de tempo ∆t = 0.02s e não

foi mantida, caindo exponencialmente logo após o pico, conforme mostrado na figura

4.29.

Fig. 4.29: Aquecimento da lâmina de níquel para tempo curto e longo, durante experi-mento realizado por Lin e Farouk [8].

Conforme mostrado na 4.29, o aquecimento da lâmina de níquel é relativamente

lento. O modelo termodinâmico assume temperatura T1 prescrita na parede esquerda

já em t = 0, o que significa um aquecimento instantâneo. Assumir um aquecimento de

0,02s como instantâneo pode parecer uma aproximação grosseira, porém, no intuito de

efetuar mais um estudo comparativo com o modelo termodinâmico, realizamos uma

simulação assumindo T1 = 343K, que foi o valor máximo atingindo pela lâmina de

níquel utilizada no experimento.

A figura 4.30 mostra o resultado experimental de Lin e Farouk [8] e a solução

obtida com o código termodinâmico utilizando um valor para impedância acústica

Z = 0.9987. A previsão para a pressão média obtida com o modelo possui uma boa con-

76


0.000 0.005 0.010 0.015 0.020Tempo H sL0

20

40

60

80

100

120

140

Pressão H Pa L

0.01 0.02Tempo H s L

150

300Pressão H Pa L


Fig. 4.30: Resultados experimentais da pressão medidos por Lin e Farouk [8] utili-zando ar padrão e solução do código termodinâmico.

cordância até o tempo t = 0.01s. O crescimento da pressão também foi corretamente

previsto durante a primeira metade do experimento. Porém, após o tempo t = 0.01s,

apesar da amplitude entre o experimento e a solução modelada estarem concordantes,

o modelo começa a prever o retorno da pressão à pressão inicial P0, enquanto que a

amplitude da pressão obtida no experimento continua com um leve crescimento. Este

resultado era esperado, uma vez que, no intervalo de comparação entre os dados expe-

rimentais e a solução do modelo, a temperatura gerada pelo circuito RC continua su-

bindo. Desta forma, ainda ocorre uma geração de ondas termoacústicas que, conforme

descrito nos capítulos anteriores, aparecem devido a qualquer variação de temperatura.

Como o modelo termodinâmico desenvolvido assume temperatura T1, desde t = 0, as

ondas termoacústicas são geradas apenas nos instantes iniciais e depois se dissipam

devido às forças viscosas do fluido e, principalmente, a impedância das paredes. A fi-

gura 4.30 também mostra o gráfico da pressão termodinâmica dentro da cavidade caso

a impedância acústica não fosse incluída no modelo, mostrando que a impedância é

fundamental, mesmo quando o aquecimento da parede não é devidamente modelado.

77

4.3 Estimativa do valor da impedância acústica

Para avaliação do valor real da impedância acústica de uma superfície faz-se ne-

cessário a execução de um experimento utilizando um tubo de impedância contendo

uma amostra do material que compõe a superfície. O experimento consiste em um

onda incidente e uma onda refletida propagando-se ao longo do eixo de um tubo cilín-

drico com a superfície teste na extremidade do tubo [54]. Um transdutor (auto-falante)

gera estas ondas na outra extremidade do tubo. O esquema pode ser visto na figura

4.31, onde os itens 1 e 2 indicam as duas posições ocupadas por microfones/sensores

que realizam as medições do campo de pressão, o item 3 indica a amostra de material

acústico, o item 4 uma tampa rígida e o item 5 um alto-falante ou transdutor.

Fig. 4.31: Esquema e equipamentos que compõem um tubo de impedância [9].

Segundo Rossing et al.[54], a expressão para o coeficiente de reflexão ℜ utilizando

o arranjo do tubo de impedância mostrado na figura 4.31 é dada por:

ℜ= Z2 −Z1

Z2 +Z1(4.21)

A impedância específica Zn do meio pode ser definida como:

Zn = ρcn (4.22)

onde o valor da impedância específica é função direta da massa específica e da veloci-

dade do som no meio.

78

Uma vez que os resultados experimentais produzidos por Lin e Farouk [8] e Brown

[3] foram realizados em tubos cilíndricos de PVC, podemos assumir, no intuito de ob-

ter uma estimativa do valor real da nossa impedância acústica, nossa cavidade como

um tubo de impedância. No nosso modelo termodinâmico, a impedância acústica Z

consiste na mesma grandeza física ℜ, uma vez que significa a parcela de onda termoa-

cústica refletida após o contato com a parede da cavidade.

Utilizando Z1 = 4,286×104 para o ar e Z2 = 3,27×106 para o PVC [9], obtemos

através da equação 4.21 o valor ℜ = 0,97413. O valor utilizado no modelo termodi-

nâmico para ajustar-se a solução aos resultados dos experimentos realizados por Lin e

Farouk [8] com ar foi Z = 0,9987. Para os experimentos realizados por Brown [3] uti-

lizando ar, o valor da impedância usada no nosso modelo foi Z = 0,99281101. Os dois

valores utilizados estão acima da estimativa do valor da impedância realizada aqui, o

que significa que assumimos uma perda na parede menor que a estimada utilizando

valores tabelados para o material PVC, (tubo) e para o ar, gás que preenche a cavi-

dade. Caso o valor utilizado da impedância acústica para prever o comportamento da

pressão ao longo do tempo nos experimentos fosse menor que o estimado, poderíamos

estar atribuindo perdas de outros fenômenos à impedância acústica. Porém, como o

valor foi mais alto, podemos afirmar que o efeito da impedância acústica é provavel-

mente o único responsável pelas divergências entre os resultados experimentais e as

previsões numéricas. Esta é a provável razão do novo modelo incluindo impedância

ter gerado uma boa concordância com os dados experimentais, o que não ocorria com

a versão do modelo sem impedância das paredes.

79

Capítulo 5

Conclusões e Trabalhos Futuros

5.1 Conclusões

Os resultados dos estudos mostrados no capítulo anterior nos permitem tirar algu-

mas conclusões a respeito da modelagem da impedância acústica de parede na trans-

ferência de calor termoacústica e das características dos modelos utilizados durante o

estudo.

• A impedância acústica de parede, para o problema de ondas termoacústicas con-

finadas numa cavidade com gases ideias, é relevante na previsão do comporta-

mento da pressão do sistema ao longo do tempo. A partir da modelagem da

impedância acústica utilizando os modelos hidrodinâmicos e termodinâmicos

do problema, os resultados numéricos, diferente do que tem sido mostrado na

literatura recente, oferecem excelentes previsões quando comparados com os

resultados experimentais. Os resultados obtidos mostram que os modelos utili-

zados em outros trabalhos para modelagem do efeito pistão se limitam a prever o

comportamento da pressão dentro da cavidade para temos curtos, onde os efeitos

da impedância ainda não se acumularam devido ao pequeno intervalo de tempo

percorrido.

• O modelo hidrodinâmico do problema, utilizando as Equações de Navier-Stokes,

é flexível e, através da modelagem da condição de contorno incluindo impedân-

80

cia acústica de parede, pode oferecer excelentes previsões. Contudo, possui em

alto custo computacional e mesmo utilizando a tecnologia de computação para-

lela aqui implementada é necessária uma boa estrutura de processamento numé-

rico para se obter as soluções em tempos longos de simulação, que incluem o

tempo característico de difusão térmica.

• O modelo termodinâmico do problema, apesar de assumir algumas simplifica-

ções, mostrou-se aplicável e ofereceu resultados com boa concordância qualita-

tiva com os experimentos, além de possuir um custo computacional baixo, por

possuir uma solução analítica. Este trabalho representou o primeiro passo no de-

senvolvimento deste modelo que, após ter seu desenvolvimento completo, será

essencial para simular os efeitos termoacústicos em modelos hidrodinâmicos.

• A geração e transmissão de ondas termoacústicas relaciona-se diretamente com

a velocidade de aquecimento da parede. Quando a parede atinge a tempera-

tura final de aquecimento antes do intervalo de um tempo característico acústico

(t/ta ¿ 1) a onda termoacústica possui uma frente aguda e uma cauda longa.

Caso a parede sofra um aquecimento lento e só atinga a temperatura final de

aquecimento após o tempo t/ta = 1, não há a formação de uma cauda e a onda

possui uma forma mais quadrada.

5.2 Trabalhos Futuros

Refletindo sobre os estudos e conclusões da investigação numérica da transferência

de calor termoacústica, vislumbramos pontos com necessidade de aprofundamento de

pesquisa, metodologias a serem adotadas, bem como outros aspectos aqui não avalia-

dos, mas que são continuações lógicas induzidas pelos resultados obtidos.

5.2.1 Evolução nas técnicas numéricas adotadas

O problema de Transferência de Calor Termoacústica envolve escalas de tempo

característico bem diferentes, uma vez que os efeitos acústicos possuem um tempo ca-

81

racterístico da ordem de 10−4 segundos, os fenômenos de difusão térmica ocorrem num

tempo da ordem de 103 segundos. Desta forma, para capturar o efeito termoacústico

faz-se necessário simular o problema até o regime permanente, controlado pela difu-

são térmica, porém capturando a propagação de ondas acústicas. A consequência deste

extenso tempo físico de simulação marchando num passo de tempo curto são enormes

tempos computacionais, da ordem de meses e anos. Faz-se necessário então a exis-

tência e disponibilidade de uma potente estrutura de computadores aliada a métodos

eficientes de computação paralela. Uma opção, já vislumbrada pelo grupo de pesquisa

no qual este trabalho é inserido, é a utilização de computação paralela com memória

distribuída GPU (Graphics Processing Unit). Essa seria uma opção futura ao OpenMP

Fortran utilizado neste trabalho, que mostrou bons resultados, porém exige computa-

dores com um elevado número de processadores compartilhando a mesma memória.

De posse de uma estrutura computacional mais poderosa será possível simular tempos

longos com o código hidrodinâmico, tarefa não realizada aqui por possuir um tempo

completamente inviável e incompatível com a duração de um curso de mestrado.

O esquema utilizado na discretização espacial dos fluxos invíscidos neste traba-

lho foi o Flux Vector Splitting - FVS com propriedade TVD e limitador de fluxo. O

esquema é eficiente e capta as ondas de choque sem oscilações expúrias ou instabi-

lidades conforme proposto. Porém possui o inconveniente de, na região do choque,

diminuir a ordem de erro para a unidade. Um esquema Weno (Weighted Essentially

Non-Oscillatory), em fase de implementação no grupo de pesquisa, aproxima um po-

linômio de maiores ordens sem pagar o preço dessa queda de ordem na região da

descontinuidade. A utilização de um esquema desse tipo é uma evolução natural, que

pode ser aplicada ao problema discutido neste trabalho e será implementada em breve.

Isto deve possibilitar, em tese, uma diminuição da densidade da malha espacial adotada

neste trabalho.

O modelo termodinâmico para o problema de transferência de calor termoacústica

ainda precisa ser melhorado para prever, com precisão, o comportamento da pressão

ao longo do tempo. A pressão média da cavidade e o decaimento da pressão estão qua-

82

litativamente de acordo com os experimentos, porém a avaliação da escala de tempo

precisa ser melhorada. O trabalho atual representa apenas o primeiro passo no desen-

volvimento deste modelo.

5.2.2 Aspectos físicos e modelagem

Um extensão natural do trabalho aqui realizado, que já é objeto de estudo de um

outro membro do grupo, é a avaliação do efeito da variação não-linear das propriedades

do fluido o aquecimento na extremidade da cavidade. No trabalho atual, as proprieda-

des do fluido são avaliadas no tempo T0 e seu valor não muda ao longo do tempo, o

que é uma aproximação, principalmente quando lidamos com fluidos supercríticos.

Outra extensão natural deste trabalho, é a utilização de um modelo bidimensional,

onde efeitos da gravidade passam a ser considerados. Surge então um fenômeno con-

vectivo adicional, conhecido como convecção natural, causado pela gravidade devido

a diferença entre massas específicas do fluido à diferentes alturas dentro da cavidade

bidimensional. Este problema envolve maior sofisticação na modelagem física e ma-

temática e pode ser utilizado como tema de uma tese de doutorado de um membro do

grupo.

O valor da impedância acústica aqui utilizado foi variado no intuito de se propor

um valor que ajustasse os resultados numéricos aos experimentais, avaliando se este

efeito de impedância ajusta as previsões numéricas com os resultados experimentais.

Faz-se necessário então, uma vez que o valor da impedância específica ZS é função

da natureza da interface entre os meios e da velocidade do som na parede e no fluido,

um estudo focado no valor real desta grandeza física, de forma a verificar a consistên-

cia entre os valores utilizados nas simulações termoacústicas e o valor obtido por um

estudo específico experimental para impedância.

83

Referências

[1] Yiqiang Lin. Acoustic Wave Induced Convection and Transport in Gases under

Normal and Micro-gravity Conditions. PhD thesis, Drexel University, 2007.

[2] Bakhtier Farouk, Elaine S. Oran, e Toru Fusegi. Numerical study of thermoa-

coustic waves in an enclosure. Physics of Fluids, 12(5):1052–1061, May 2000.

[3] M. A. Brown. Thermally Induced Pressure Wave in a Gas: Experimental Ob-

servation and Theoretical Prediction of Thermoacoustic Convection. PhD thesis,

University of Pennsylvania, Philadelphia, PA, 1992.

[4] Y. Huang e H. H. Bau. Thermoaocustic waves in a semi-infinite medium. Inter-

national Journal of Heat and Mass Transfer, 38(8):1329–1345, 1995.

[5] John D. Anderson. Computational Fluid Dynamics: The Basics with Applicati-

ons. Mc Graw Hill Series in Mechanical Engineering. Mc Graw Hill Education,

indian edition, 2013.

[6] Culbert B. Laney. Computational Gasdynamics. Cambridge University Press,

1998.

[7] Murat K. Aktas. Thermoacoustically Induced and Acoustically Driven Flows and

Heat Transfer in Enclosures. PhD thesis, Drextel University, 2004.

[8] Yiqiang Lin e Bakhtier Farouk. Experimental and numerical studies of thermally

induced acoustic waves in an enclosure. Journal of Thermophysics and Heat

Transfer, 22(1):105–114, January-March 2008.

[9] Leonardo Ferreira Lopes. Uso de materiais porosos em filtros acústicos. Master’s

thesis, Universidade Federal de Santa Catarina, 2006.

[10] Abdulnaser Sayma. Computational Fluid Dynamics. Ventus Publishing ApS,

2009.

84

[11] Richard H. Pletcher, John C. Tannehill, e Dale A. Anderson. Computational Fluid

Mechanics and Heat Transfer. Series in Computational and Physical Process in

Mechanical and Thermal Sciences. CRC Press, third edition, 2013.

[12] P. Carlès. A brief review of the thermophysical properties of supercritical fluids.

The Journal of Supercritical Fluids, 53(1):2–11, 2010.

[13] L. Trilling. On Thermally Induced Sound Fields. Group report. M.I.T. Fluid

Dynamics Research Group, 1954.

[14] B. K. Larkin. Heat flow to a confined fluid in zero gravity. Progress in Astronau-

tics and Aeronautics, 20:819–832, 1967.

[15] H Ozoe, N Sato, e S. W. Churchill. The effect of various parameters on thermo-

acoustic convection. Chemical Engineering Communications, 5(1-4):203–221,

1980.

[16] Y. Huang e H. H. Bau. Thermoacoustic convection. in heat transfer in microgra-

vity. Heat Transfer Division Newsletter, ASME HTD-269:13–22, 1993.

[17] Y. Huang e H. H. Bau. Thermoaocustic waves in a confined medium. Internati-

onal Journal of Heat and Mass Transfer, 40(2):407–419, 1997.

[18] M. A. Brown e S.W. Churchill. Finite-difference computation of the wave motion

generated in a gas by a rapid increase in the bounding temperature. Computers

and Chemical Engineering, 23(3):357–376, 1999.

[19] Murat K. Aktas e B. Farouk. Numerical simulation of developing natural convec-

tion in an enclosure due a rapid heating. International Journal of Heat and Mass

Transfer, 46(12):2253–2261, 2003.

[20] B. Zappoli, D. Bailly, Y. Garrabos, B.L. Neindre, P. Guenoun, e D. Beysens.

Anomalous heat transport by the piston effect in supercritical fluids under zero

gravity. Physical Review, A(41):2264–2267, 1990.

85

[21] A. Nakano. Studies on piston and soret effect in a binary mixture supercritical

fluid. International Journal of Heat and Mass Transfer, 50:4678–4687, 2007.

[22] B. Shen e P. Zhang. Thermoacoustic wave propagation and reflection near the

liquid-gas critical point. Physical Review E, 79, 2009.

[23] Y. Miura, S. Yoshihara, M. Ohnishi, K. Honda, M. Matsumoto, J. Kawai, M. Ishi-

kawa, H. Kobayashi, e A. Onuki. Hidh-speed observation of the piston effect near

the gas-liquid critical point. Physical Review E, (74), 2006.

[24] B. Shen e P. Zhang. On the transition from thermoacoustic convection to diffusion

in a near-critical fluid. International Journal of Heat and Mass Transfer, 53:

4832–4843, 2010.

[25] B. Shen e P. Zhang. Thermoacoustic waves along the critical isochore. Phys. Rev.

E, 83(1):011–115, 2011.

[26] A. Onuki, H. Hao, e R. A. Ferrell. Fast adiabatic equilibration in a single-

component fluid near the liquid-vapor critical point. Physical Review A, 41(4):

2256–2259, 1990.

[27] H. Boukari, Shaumeyer J.N., Briggs M.E., e Gammon R.W. Critical speeding up

in pure fluids. Physical Review A, A(41):2260–2263, 1990.

[28] P. Carlès e K. Dadzie. Two typical time scales for piston effect. Physical Review

E, (71), 2005.

[29] L.S. de B. Alves. An integral transform solution for unsteady compressible heat

transfer in fluids near their thermodynamic critical point. Thermal Science, 17

(3):673–686, 2013.

[30] V. S. Nikolayev, A. Dejoan, Y. Garrabos, e D. Beysens. Fast heat transfer calcu-

lations in supercritical fluids versus hydrodynamic approach. Physical Review E,

(67), 2003.

86

[31] Patricia C. Teixeira e L.S. de B. Alves. Modeling supercritical heat transfer

in compressible fluids. International Journal of Thermal Science, 88:267–278,

2015.

[32] M. Parang e A. Salah-Eddine. Thermoacoustic convection heat-transfer pheno-

menon. AIAA Journal, 22(7):1020–1022, 1984.

[33] A. Nakano e M. Shiraishi. Piston effect in supercritical nitrogen around the

pseudo-critical line. International Communications in Heat and Mass Transfer,

32:1152–1164, 2005.

[34] Peter D. Lax. Hyperbolic systems of convervation laws and mathematical theory

of shock waves. Society for Industrial and Applied Mathematics, 1973.

[35] A. Harten. High order schemes for hyperbolic conservation laws. Journal of

Computational Physics, 49:357–385, 1983.

[36] Boris J. e Book D. Flux-corrected transport: I shasta, a fluid transport algorithm

that works. Journal of Computational Physics, 11:38–69, 1973.

[37] B. van Leer. Towards the ultimate conservation difference scheme, ii: monotoni-

city and convervation combined in a second-order scheme. Journal of Computa-

tional Physics, 14:361–370, 1974.

[38] P. L. Roe. Some contributions to the modeling of discontinuous flows. In Proce-

edings of 1983 AMS-SIAM Summer Seminar on Large Scale Computing in Fluid

Mechanics, volume 22, pages 163–193, 1985.

[39] P. K. Sweby. High resolution schemes using flux limiter for hyperbolic conser-

vation laws. SIAM J. Num. Anal., 21:995–1011, 1984.

[40] A. Harten. The artificial compression method for computation of shocks and

contact discontinuities: Iii. self-adjusting hybrid schemes. Math. Comput., 32:

363–389, 1978.

87

[41] Sigal Gottlieb, David Ketcheson, e Chi-Wang Shu. Strong Stability Preserving

Runge-Kutta and Multistep Time Discretizations. World Scientific, 2011.

[42] C. W. ans Osher S. Shu. Efficient implementation of essentially non-oscillatory

shock-capturing schemes. Journal of Computational Physics, 77:439–471, 1998.

[43] C. W. Shu. Total-variation-diminishing time discretizations. SIAM J. Sci. Stat.

Comput., 9:1073–1084, 1988.

[44] Sigal Gottlieb. On high order strong stability preserving runge-kutta and multi

step time discretizations. Journal of Scientific Computing, 25(1/2):105–127,

2005.

[45] Gottlieb S. e Shu C. W. Total variation diminishing runge-kutta schemes. Math.

Comput., 67:73–85, 1998.

[46] R. M. Cotta. Integral Transforms in Computational Heat and Fluid Flow. CRC

Press, 1993.

[47] R. M. Cotta. The Integral Transform Method in Thermal and Fluids Science and

Engineering. Begell House, Inc., 1998.

[48] M. Necati Ozisik. Heat Conduction. John Wiley and Sons Inc, second edition,

1993.

[49] S. Wolfram. The Mathematica Book. Wolfram Media/Cambridge University

Press, 5th edition edition, 2003.

[50] M. N. Özisik e R. L. Murray. On the solution of linear diffusion problems with

variable boundary condition parameters. Journal of Heat Transfer, 98c:48–51,

1974.

[51] M. D. Mikhailov. General solutions of the heat equation in finite regions. Inter-

national Journal of Engineering Science, 7:577–591, 1972.

88

[52] Leonardo S. de B. Alves. Unsteady compressible heat transfer in supercritical

fluids. In Proceedings of 14th Brazilian Congress of Thermal Sciences and En-

gineering, volume ENCIT2012-347, pages 1–9, 2012.

[53] P. C. Texeira e Leonardo S. de B. Alves. Piston effect characteristic time depen-

dence on equation of state model choice. In Proceedings of Int. Symp. on Adv. in

Comp. Heat Transfer, volume CHT-12-EN02, pages 1–11, 2012.

[54] Rossing T.D. et al. Spring handbook of acoustics. Springer Science and Business

Media, 2007.

[55] Carlos Spa, Adan Garriga, e Jose Escolano. Impedance boundary conditions for

pseudo-spectral time-domain methods in room acoustics. Applied Acoustics, 71

(5):402 – 410, 2010.

[56] R.E. Sonntag, C. Borgnakke, e G.J. Van Wylen. Fundamentals of Thermodyna-

mics. Number v. 7. Wiley, 2010.

[57] A. Onuki e R. A. Ferrell. Adiabatic heating effect near the gas-liquid critical

point. Physica, A(164):245–264, 1990.

[58] Miguel Hermanns. Parallel Programming in Fortran 95 using OpenMP. SAE -

DMT - Universidad Politécnica de Madrid, 2002.

89

Apêndice A

Computação Paralela com OpenMP Fortran API

A.1 Computação Paralela

Na necessidade de maior poder de processamento, os desenvolvedores de compu-

tadores começaram a pensar numa forma de utilizar os computadores existentes de

maneira conjunta. A partir desta ideia, se originou a computação paralela e iniciou-se

um novo campo para programadores e pesquisadores [5].

Nos dias atuais, cluster de computadores em paralelo são muito comuns em labo-

ratórios de pesquisa bem como em empresas ao redor do mundo sendo amplamente

utilizados em cálculos complexos como simulação de explosões atômicas, de quebra

de proteínas e de fluxos turbulentos [5].

O desafio então no campo de Computação Paralela é desenvolver códigos capazes

de utilizar a capacidade das máquinas disponíveis no intuito de resolver complexos

problemas em menor tempo. Mas programação paralela não é uma tarefa fácil, uma

vez que uma grande variedade de arquiteturas existem e o ganho com a paralelização

entre um código e outro pode variar consideravelmente [58]. Principalmente duas

famílias de arquiteturas podem ser identificadas:

• Memória Compartilhada: as máquinas em paralelo são construídas de forma

que seus grupos de processadores tem acesso a uma memória comum a todos.

A tecnologia conhecida como SMP é baseada nesta arquitetura, onde SMP é a

abreviação para Symmetric Multi Processing[58].

• Memória Distribuída: nestas máquinas paralelas cada processador tem sua pró-

pria memória privada e a informação e trocada entre os processadores por meio

de mensagens. O nome de cluster é comumente usado para este tipo de monta-

90

gem de computadores[58].

Fig. A.1: Esquema de arquitetura com memória compartilhada e memória distribuída.

A figura A.1 esquematiza as duas arquiteturas e a utilização de memória por cada

uma delas. Cada uma das duas famílias tem suas vantagens e desvantagens e a atual

programação paralela tenta explorar suas vantagens focando em apenas uma dessas

arquiteturas. Nos últimos anos, uma nova indústria foi criada com o objetivo de servir

como uma boa base para o desenvolvimento de programas paralelizados em máquinas

de memória compartilhada: OpenMP.

OpenMP é uma interface de programação (Aplication Program Interface - API),

portátil, baseada no modelo de programação paralela de memória compartilhada para

arquiteturas de múltiplos processadores. É composto por três componentes básicos:

• Diretivas de Compilação;

• Biblioteca de Execução;

• Variáveis de Ambiente.

OpenMP está disponível para uso com os compiladores C/C++ e Fortran, podendo ser

executado em ambientes Unix e Windows (Sistemas Multithreads). Definido e mantido

por um grupo composto na maioria por empresas de hardware e software, denominado

como OpenMP ARB (Architecture Review Board).

91

A.2 Breve Histórico

No início da década de 90, fabricantes de máquinas de memória compartilhada e

similar, desenvolveram extensões do compilador Fortran com um conjunto especial de

instruções denominadas de diretivas de execução que permitiam:

• Usuário, programador Fortran (programação serial), adicionar instruções para

identificar loops que poderiam ser paralelizados;

• O compilador passa a ser responsável pela paralelização automática desses loops

por entre os processadores do ambiente SMP.

O primeiro padrão, ANSI X3H5, para testes foi lançado em 1994, mas nunca foi

adotado devido ao grande uso e interesse, na época, por máquinas de memória distri-

buída (clusters).

Um novo padrão foi desenvolvido em 1997 a partir do padrão anterior, quando as

arquiteturas de memória compartilhada se tornaram mais comum.

Em 28 de Outubro de 1997 foi divulgado e disponibilizado o OpenMP Fortran API

e no final de 1998 foi disponibilizado o OpenMP C/C++ API.

A.3 Filosofia Básica OpenMP

A programação paralela Fortran OpenMP segue algumas filosofias básicas que se

enquadram em conceitos da Ciência da Computação, e que merecem ser citados aqui,

no intuito de proporcionar um maior entendimento da sistemática de paralelização de

um código inicialmente serial.

• Paralelismo baseado em trilhas: Processo de memória compartilhada que con-

siste em threads, que representa uma unidade mínima capaz de processar uma

sequência de comandos.

• Paralelismo explícito: O OpenMP é um recurso de programação paralela que

permite ao programador total controle sobre a paralelização do código.

92

• Modelo Fork-Join: todos os programas OpenMP começam com uma trilha mes-

tre até a primeira bifurcação (Fork), onde diversas trilhas são criadas e os outros

processadores começam a ser executados. Os comandos do programa inseri-

dos na região paralela serão executados pelas diversas trilhas criadas. Quando

o conjunto de trilhas completam a execução dos comandos na região paralela,

as trilhas são sincronizadas e finalizadas (Join), permanecendo apenas a trilha

“mestre”. A figura A.2 esquematiza este processo.

• Paralelismo baseado em diretivas de compilação: todo paralelismo do OpenMP

é especificado através de diretivas de compilação que criam trilhas, sincronizam,

inciam e terminam seções paralelas ou serias do código .

• Paralelismo recursivo: permite que uma região paralela exista em outras regiões

que podem ser executadas em paralelo, e assim sucessivamente. Este aspecto

depende inteiramente da implementação OpenMP.

• Trilhas dinâmicas: é possível, durante a execução do programa alterar o número

de trilhas que serão criados.

• Variáveis classificáveis: como OpenMP é baseado no modelo de programação de

memória compartilhada, existem dois tipos básicos de variáveis: variáveis glo-

bais e variáveis locais ou privadas. A maioria das variáveis são compartilhadas

entre trilhas através de variáveis globais, porém é necessário ser avaliar criteri-

osamente a definição de uma variável para que as informações sejam trocadas

entre

A.4 Paralelização do código Compressível

Antes de inciar o processo de paralelização do código compressível, foi feito um

estudo do custo computacional envolvido com a execução de cada grupo de tarefas.

Foram utilizadas variáveis em determinadas linhas do código que recebiam o horá-

rio do sistema ao longo da sequência de execução. A partir da subtração entre duas

93

Fig. A.2: Sequência de execução de um código paralelizado com modelo Fork-Join.

varáveis localizadas em diferentes pontos do código, é possível quantificar o tempo

execução da seção entre as duas linhas. Uma análise preliminar utilizando uma malha

de 10.000 pontos com apenas dois passos no tempo e apenas um estágio do método

Runge-Kutta foi feita para se confirmar as estimativas de custo por seção do código.

Cada grupo de tarefas foi definido da seguinte forma:

• Inicialização - Consiste na seção do código que inicializa as variáveis, aloca ma-

trizes e vetores, obtém informações como número de processadores disponíveis

para execução, abre arquivos que serão preenchidos, calcula a malha espacial

e calcula as condições iniciais do problema. Este grupo de tarefas consumiu

0,456% do tempo de execução do código.

• Cálculo dos Fluxos Invíscidos – Consiste na discretização espacial dos fluxos

positivos e negativos ao longo de toda malha a cada passo no tempo. Este grupo

de tarefas consumiu 98,615% do tempo de execução do código.

• Cálculo dos Fluxos Viscosos – Consiste na discretização espacial dos fluxos

viscosos ao longo de toda malha a cada passo no tempo. Este grupo de tarefas

consumiu 0,473% do tempo de execução do código.

94

• Discretização Temporal – Consiste basicamente no cálculo da solução no passo

posterior a partir dos dados de fluxo calculados na seção anterior. Este grupo de

tarefas consumiu 0,189% do tempo de execução do código.

• Gravação de resultados nos arquivos e finalização. Este grupo de tarefas consu-

miu 0,267% do tempo de execução do código.

Como já era esperado, praticamente todo o custo computacional para execução do

nosso código compressível deve-se ao cálculo da discretização espacial dos fluxos in-

víscidos. Desta forma, esta tarefa foi completamente paralelizada, enquanto que outras

tarefas foram parcialmente paralelizadas. Desta forma tivemos a seguinte situação a

respeito da paralelização das seções do código:

Tab. A.1: Tipo de execução por seção do código hidrodinâmico.

Seção do código Tipo de Execução

Inicialização SerialCálculo dos fluxos viscosos e invíscidos Paralelo

Discretização temporal Runge-Kutta (2 estágios) ParaleloCálculo das condições de contorno Serial

Atualização da solução e das variáveis primitivas ParaleloGravação dos arquivos Serial

A.5 Resultados de ganho com código paralelizado

Para executar os códigos dispomos de um computador pessoal AMD FX-6300 3,5

GHz com 6 núcleos e de um cluster com 32 processadores Intel Core 2-Duo cada. Foi

feito um estudo avaliando o tempo de execução variando o número de processadores e

comparando com o tempo de execução serial. O ganho máximo (Speed-Up) obtido foi

de 12,79 vezes o tempo de execução serial para o cluster de 32 processadores. Para o

computador pessoal FX-6300, o ganho máximo obtido foi de 3,55 vezes o tempo de

execução serial. As tabelas A.2, A.3, A.4 e A.5 mostram os resultados obtidos para

quatro casos distintos:

95

• Caso 1: Execução com uma malha de 401 pontos, com tempo físico final de

1E-5 segundos no cluster.


1E-5 segundos no cluster.


1E-6 segundos no computador pessoal.


1E-6 segundos no computador pessoal.

Através das figuras A.3 e A.4 é possível analisar o comportamento do Speed-Up e

da eficiência, que é calculada como a razão entre o ganho obtido com a paralelização

com n processadores e o próprio número n.

Fig. A.3: Speed-Up e Eficiência x Número de Processadores - Cluster 32 processado-res

96

Fig. A.4: Speed-Up e Eficiência x Número de Processadores - Desktop

97

Tab. A.2: Resultado de Tempo de Execução e Speed-Up para o caso 1.

Num.Processadores Tempo Execução(s) Speed-Up(×) Eficiência

1 276,99 1 100%2 141,80 1,95 97,56%3 96,65 2,86 95,43%4 72,97 3,79 94,80%5 61,55 4,50 89,91%6 54,13 5,11 85,19%7 48,62 5,69 81,30%8 43,21 6,40 80,04%9 41,26 6,71 74,51%

10 35,99 7,69 76,88%11 35,33 7,83 71,20%12 33,32 8,30 69,20%13 35,89 7,71 59,30%14 32,97 8,39 59,94%15 32,94 8,40 56,00%16 25,79 10,73 67,06%17 36,70 7,54 44,35%18 33,45 8,27 45,95%19 33,60 8,23 43,34%20 30,87 8,96 44,82%21 31,42 8,81 41,93%22 29,77 9,29 42,25%23 32,73 8,45 36,76%24 34,48 8,02 33,44%25 25,53 10,43 41,71%26 29,99 9,23 35,48%27 36,08 7,67 38,41%28 28,80 9,61 34,41%29 36,31 7,62 26,28%30 32,57 8,50 28,32%31 34,43 8,04 25,93%32 34,09 8,12 25,36%

98



1 704,78 1 100,00%2 358,15 1,97 98,39%3 244,15 2,89 96,22%4 184,98 3,81 95,25%5 154,42 4,56 91,28%6 133,75 5,27 87,82%7 122,81 5,74 81,98%8 107,14 6,58 82,23%9 97,79 7,21 80,08%

10 88,35 7,98 79,77%11 83,60 8,43 76,64%12 79,46 8,87 73,91%13 75,09 9,39 72,20%14 78,21 9,01 64,37%15 73,74 9,56 63,72%16 62,18 11,33 70,84%17 87,00 8,10 47,65%18 82,46 8,55 47,48%19 82,13 8,58 45,16%20 75,20 9,37 46,86%21 76,87 9,17 43,66%22 72,38 9,74 44,26%23 81,23 8,68 37,72%24 77,08 9,14 38,10%25 85,71 8,22 32,89%26 70,89 9,94 38,24%27 79,16 8,90 32,98%28 90,90 7,75 27,69%29 59,38 11,87 40,93%30 65,10 10,83 36,09%31 84,00 8,39 27,07%32 55,11 12,79 39,96%

99



1 64,18 1,00 100,00%2 34,43 1,86 93,20%3 24,49 2,62 87,35%4 21,21 3,03 75,64%5 18,09 3,55 70,95%6 23,65 2,71 45,23%



1 133,87 1,00 100,00%2 82,74 1,62 80,90%3 60,4 2,22 73,88%4 48,37 2,77 69,19%5 42,304 3,16 63,29%6 47,536 2,82 46,94%

100

Apêndice B

Propriedades de Estabilidade Não-Linear

Conforme visto no capítulo 1, a presença de Estabilidade Não-Linear é indispensá-

vel nos esquemas numéricos para problemas que envolvem choques e contatos, fenô-

menos físicos que induzem grandes oscilações espúrias nos resultados em vez de so-

luções monotônicas e estáveis. Desta forma, este apêndice tem por objetivo expor os

conceitos ligados às propriedades de Estabilidade Não-Linear. A teoria moderna de

Estabilidade Não-Linear inclui propriedades que podem ser dividades em dois grupos:

aquelas que são consideradas difíceis de provar ou induzir num determinado esquema

numérico, e aquelas relativamente mais simples de serem introduzidas e, consequente-

mente, mais populares e mais práticas. Como representantes do primeiro grupo temos

a Monoticidade Preservada, TVD e Amplitude Decrescente. Do segundo grupo, pode-

mos citar a Positividade e Upwind Range Condition [6].

B.1 Monoticidade Preservada

No estudo de Leis de Conservação Escalar em domínios espaciais infinitos, di-

zer que as soluções preservam a propriedade de monoticidade significa dizer que, se

as condições iniciais u(x,0) aumentam a monoticidade, a solução u(x, t ) também au-

menta a monoticidade para qualquer tempo t , e se as condições iniciais decaem a

monoticidade, a solução terá a monoticidade decrescente para qualquer tempo t . A

propriedade chamada Monoticidade Preservada é uma consequência da preservação

da forma de onda, da destruição da forma de onda e da criação da forma da onda [6].

Suponha que um método numérico possui a propriedade Monoticidade Preservada.

Então, se as condições iniciais são monotônicas, estes métodos não permitem oscila-

ções e, deste modo, são estáveis. Contudo, de uma forma geral, Monoticidade Preser-

101

vada é uma condição de estabilidade pobre pelas seguintes razões:

• Monoticidade Preservada não produz estabilidade em soluções não-monotônicas;

• A maioria das tentativas de forçar Monoticidade Preservada termina por for-

çar outras fortes condições de estabilidade não-linear, como a positividade, por

exemplo;

• Monoticidade Preservada não permite pequenas oscilações benignas nas solu-

ções monotônicas. Pequenas oscilações não necessariamente indicam séria ins-

tabilidade, e na tentativa de purgar todos os erros oscilatórios pode-se causar

ainda maiores erros não-oscilatórios.

Por fim, Monoticidade Preservada é muito fraca em algumas condições e muito

forte em outras. A condição de Monoticidade Preservada, por ser uma condição de es-

tabilidade não-linear, se aplica também à métodos lineares. Métodos de monoticidade

preservada linear possuem primeira ordem de acurácia no máximo (Teorema de Godu-

nov). Devido a isto, este teorema é frequentemente usado para justificar métodos ine-

rentemente não-lineares, os quais são não-lineares até mesmo quando as equações de

governo são lineares. Contudo, o teorema de Godunov é possivelmente melhor usado

como argumento contra a a condição Monoticidade Preservada, um critério mais rela-

xado em relação às pequenas oscilações pode reduzir ou evitar o sacrifício na ordem

de acurácia do método. Mas, apesar de suas desvantagens, Monoticidade Preservada é

um elemento estabelecido na teoria moderna de estabilidade não-linear [6].

B.2 TVD – Total Variation Diminishing

Como sua maior desvantagem, a Monoticidade Preservada falha ao ser aplicada

em soluções não-monotônicas. Em essência, TVD é o menor passo possível além da

monoticidade preservada, isto é, a mais fraca condição possível que implica em preser-

vação de monoticidade e ainda se aplica a estabilidade de ambas soluções monotônicas

e não-monotônicas [6].

102

De uma maneira mais restrita, para a maioria das funções vistas nos escoamentos

compressíveis, a variação total da solução exata pode ser definida como se segue:

T V (u(•, t ) = max∞∑

i=−∞|u(xi+1, t )− (xi , t )| (B.1)

A variação total de uma função em um domínio infinito é a soma dos extremos,

máximo contado positivamente e mínimo contado negativamente. As duas fronteiras

infinitas são sempre extremos, e eles são contados uma só vez e todos os outros extre-

mos são contados duas vezes.

Baseando-se nos conceitos de preservação, criação e destruição da forma de onda,

sabe-se que máximos não crescem, mínimos não decrescem, e nenhum novo máximo

ou mínimo são criados em soluções de leis de conservação escalares. Resumindo, ne-

nhum dos fenômenos que causam aumento na variação total podem acontecer na solu-

ção de leis de conservação escalares. Então, soluções de leis de conservação escalares

são Total Variation Diminishing (TVD) e pode-se escrever para todo t2 ≥ t1:

T V (u(•, t2)) ≤ T V (u(•, t1)) (B.2)

Seguem alguns aspectos relevantes sobre a propriedade TVD:

• A maioria das tentativas de forçar TVD, assim como ocorre na propriedade Mo-

noticidade Preservada, também acabam por forçar outras condições de estabili-

dade não-linear, como a positividade.

• TVD implica em monoticidade preservada. Isto é desejável em circunstâncias

onde a preservação de monoticidade é muito fraca mas indesejável onde a pre-

servação de monoticidade é muito forte, uma vez que o TVD pode deixá-la ainda

mais forte.

• TVD tende a causar erros de clipping nos extremos, este fenômeno consiste na

perda de ordem nos extremos, esta perda é de segunda ordem ou até maior. Teori-

camente, clipping não ocorre em todos os extremos e pode ser apenas moderado

103

onde ele ocorre. Contudo, na prática, a maioria dos métodos TVD cortam todos

os extremos entre primeira e segunda ordem de acurácia.

• Em teoria, métodos TVD podem permitir oscilações, porém raramente acontece

na prática. TVD provê um limite superior nas oscilações, eliminando no pior

caso instabilidade com crescimento ilimitado.

Apesar do inconveniente de gerar perda de ordem nas regiôes de máximo e mínimo,

TVD é um elemento praticamente indispensável na teoria da estabilidade não-linear

corrente.

B.3 Amplitude Decrescente

TVD permite grande oscilações de impulso, pelo menos em algumas circunstân-

cias, a condição chamada Amplitude Decrescente (Range Diminishing) é comparati-

vamente forte e praticamente elimina as oscilações de impulso [6]. A solução exata de

leis de conservação escalar em domínios espaciais infinitos possui as seguintes propri-

edades:

1. Máximos não crescem;

2. Mínimos não decrescem;e

3. Nenhum extremo é criado ao longo do tempo.

Esta propriedade é consequência direta da preservação e destruição da forma de

onda e da criação da forma de onda monotônica. Isto implica que ambas amplitudes,

locais e globais, da solução exata contraem ao longo do tempo. Mais especificamente,

a amplitude local é mantida durante a preservação e criação da forma de onda e redu-

zida durante a destruição da forma de onda em choques. Esta propriedade é ilustrada

na figura 2. Suponha que uma aproximação numérica possua a propriedade Amplitude

Decrescente, o que significa que o máximo não aumenta, o mínimo não decresce e

104

nenhum novo máximo nem mínimo é criado. Então a aproximação numérica é obvia-

mente livre de qualquer oscilação de impulso e saltos. Infelizmente, a condição causa

uma perda de acurácia [6].

Fig. B.1: Redução de um máximo global (a) e eliminação de um mínimo local (b)[6].

Métodos que possuem a propriedade Amplitude Decrescente são formalmente, no

máximo, de segunda ordem, não sendo possível esta propriedade ser mantida para mai-

ores ordens, devido aos erros numéricos causados pelos clippings. Na grande maioria

dos casos, a ordem de acurácia destes métodos está entre um e dois, pouquíssimos mé-

todos Amplitude Decrescente atigem segunda ordem nos extremos. Seguem algumas

características da propriedade:

• Existem métodos na literatura que são capazes de formar Amplitude Decrescente

precisamente e diretamente.

• Amplitude Decrescente implica em TVD, o que é desejável quando num método

onde a propriedade TVD é fraca, porém indesejável, quando o método já possui

uma forte propriedade TVD.

• Amplitude Decrescente sempre causa clipping nos extremos. A ordem de acurá-

cia, mesmo de métodos de alta ordem, sempre caem para dois (no máximo) nos

extremos.

• Amplitude Decrescente elimina completamente qualquer oscilação expúria, sal-

tos e problemas nos extremos.

105

B.4 Positividade

A propriedade Positividade é baseada no método de divisão da velocidade de onda

em suas velocidades que se propagam na direção positiva e negativa:

∂u

∂t+ A+∂u

∂x− A+∂u

∂x= 0 (B.3)

A equação diferencial acima foi rescrita por Harten [35] utilizando uma notação

em termos de coeficientes

un+1i = un

i +C+i+1/2(un

i+1 −uni )−C−

i−1/2(uni −un

i−1) (B.4)

Sabemos que esse método é derivado naturalmente utilizando a ideia em que C±i+1/2 ≥

0, contudo sabemos que muitas diferentes formas que atendem a condição acima. Su-

ponhamos que um método numérico possa ser escrito na forma dividida de onda com

coeficientes não-negativos, tal que:

C+i+1/2 ≥ 0, C−

i+1/2 ≥ 0 (B.5)

C+i+1/2 + C−

i+1/2 ≤ 1 (B.6)

para todo i. Esta é a chamada condição de positividade, e foi primeiramente sugerida

por Harten [35]. A positividade se localiza entre duas condições de estabilidade não-

linear, TVD e Upwind Range Condition. Em particular, positividade implica em TVD

e Upwind Range Condition implica em positividade. Existem um mito que a positivi-

dade elimina oscilações espúrias e saltos, porém, ela apenas as limita, o que já é uma

consequência de interesse.

Seguem algumas caraterísticas da propriedade positividade:

• Positividade é relativamente fácil de induzir em um esquema, isto faz desta pro-

priedade, uma das mais amplamente citadas entre as condições de estabilidade

não-linear. Contudo, enquanto muitos artigos discutem positividade, a vasta mai-

oria deles, na verdade discutem a propriedade Upwind Range Condition, que

106

consiste num caso especial de positividade.

• Posivitidade impica em TVD, como ocorre em outros métodos já citados.

• Positividade causa erros simulares ao clipping nos extremos.

• Positividade também limita a ordem formal de acurácia entre um e dois.

• Positividade permite largas oscilações espúrias nos extremos.

B.5 Upwind Range Condition

Aos estudarmos as propriedades da velocidade de uma frente de onda u(x, t n+1)

verificaremos que elas são herdadas da frente de onda no tempo anterior u(x, t n). Isto

significa que se u(x, t n) possui uma certa amplitude de onda, definida pelo seu máximo

e seu mínimo, então u(x, t n+1) deve possuir a mesma amplitude. Escrevendo matema-

ticamente, a solução exata da lei de conservação escalar deve satisfazer a seguinte

condição:

minxi−1≤x≤xi

u(x, t n) ≤ u(xi , t n+1) ≤ maxxi−1≤x≤xi

u(x, t n) (B.7)

para 0 ≤λc(xi , t n+1) ≤ 1. E ainda, esta outra condição:

minxi≤x≤xi+1

u(x, t n) ≤ u(xi , t n+1) ≤ maxxi≤x≤xi+1

u(x, t n) (B.8)

para −1 ≤λc(xi , t n+1) ≤ 0.

Esta é a chamada Upwind Range Condition. As características da propriedade estão

listadas abaixo:

• Upwind Range Condition é relativamente fácil de induzir em um esquema, isto

faz desta propriedade, uma das mais amplamente citadas entre as condições de

estabilidade não-linear.

• Upwind Range Condition implica em positividade. Alguns autores a tratam

como um caso especial da propriedade positividade.

107

• Upwind Range Condition implica na condição de Amplitude Decrescente, ex-

ceto possívelmente nos pontos sônicos.

• Upwind Range Condition sempre causa quedas de ordem nos extremos, exceto

possívelmente nos pontos sônicos.

• Upwind Range Condition assume λ|c| ≤ 1, onde λ = ∆t∆x e c é a velocidade da

onda.

B.6 Variação Total Limitada

Todas as condições de estabilidade não-linear descritas até agora são consideradas

fortes, ao menos, sob certas circunstâncias. Digamos que a solução exata de uma lei

de conservação possua a seguinte propriedade:

T V (u(•, t ) ≤ M ≤∞ (B.9)

para todo t e para alguma constante M . Por exemplo, M pode ser a variação total da

condição inicial. Então a solução de leis de conservação escalar são ditas possuírem

Variação Total Limitada (Total Variation Bounded – TVB) se elas possuem um limite

M que não é excedido para todo t . TVB é, facilmente, a mais fraca condição de

estabilidade não-linear. Na verdade, TVB permite largas oscilações, prevendo apenas

que as oscilações expúrias não crescem ilimitadamente ao longo do tempo. Logo, TVB

simplesmente assegura que o método não irá divergir, ao menos de maneira oscilatória.

Características da propriedade TVB:

• As tentativas de induzir TVB, podem induzir uma outra propriedade de estabili-

dade não-linear.

• TVB é a mais fraca condição de estabilidade não-linear. Todas as outras condi-

ções de estabilidade implicam em TVB, exceto a Monoticidade Preservada, mas

TVB não implica em nenhuma das condições de estabilidade aqui vistas.

108

• TVB não força erros de clipping, como a maioria das condições de estabilidade

não-linear.

• TVB tolera largas oscilações espúrias, impondo apenas que as oscilações perma-

neçam limitadas. Embora TVB imponha um limite superior nas oscilações, este

limite pode ser mais alto que o desejado.

• Como a mais fraca condição constante na teoria de estabilidade não-linear mo-

derna, TVB é a condição mais comum aplicada a outras equações, que não sejam

Leis de Conservação Escalares.

B.7 Essencialmente Não-Oscilatória (ENO)

A partir da análise da equação B.2, podemos facilmente deduzir que:

T V (u(•, t2)) ≤ T V (u(•, t1))+O(∆xr ) (B.10)

para todo t2 ≥ t1, onde O(∆xr ) refere-se ao termo positivo de ordem r . Então, soluções

exatas de leis de conservação escalares em domínios infinitos são chamadas Essenci-

almente Não-Oscilatórias(ENO). Supondo que uma aproximação numérica herda a

propriedade ENO. Em outras palavras, estamos assumindo:

T V (um) ≤ T V (un)+O(∆xr ) (B.11)

para todo m ≥ n e qualqer valor arbitrário r > 0. Logo, a equação acima implica em:

T V (un+1) ≤ T V (un)+O(∆xr ) (B.12)

Uma vez que TVD implica em ENO, o qual, por sua vez implica em TVB. ENO,

no limite ∆x → 0 , torna-se TVD. A constante r pode ser a ordem de acurácia normal

do método. Então ENO permite aumentos de segunda ordem na variação total de uma

aproximação numérica e, em particular, aumentos de segunda ordem nos máximos,

109

e decréscimo de segunda ordem nos mínimos, eliminando o clipping tão comum no

TVD. Infelizmente, ENO teoricamente permite largas oscilações e saltos, mais que o

TVD e menos que o TVB, uma vez que que esses fenômenos podem crescer a uma

ordem r a cada passo no tempo. Além disso, como TVD e TVB, é extremamente

difícil provar diretamente que um método é ENO. Naturalmente, TVD, Amplitude

Decrescente, Positividade e Upwind range Condition, todos implicam em ENO.

B.8 Resumo das condições de Estabilidade Não-Linear

A figura B.2 esquematiza as condições de Estabilidade Não-Linear de maneira de-

crescente, isto é, da mais forte para a mais fraca a medida que descemos no fluxograma.

Uma linha conectando duas condições significa que a condição acima implica na con-

dição abaixo. Indo de baixo para cima, a força da condição de estabilidade aumenta,

tanto em termos de efeitos positivos de redução de oscilações espúrias e saltos, quanto

em termos de efeitos negativos reduzindo acurácia. Mais especificamente, quando a

força da condição de estabilidade aumenta, oscilações e saltos espúrios decrescem,

enquanto a acurácia formal diminui nos extremos ou ainda globalmente [6].

Por fim, vale a pena mencionar a expressão “método TVD” empregada atualmente

a um determinado método numérico, significa que

1. Envolve uma solução sensível obtida por meio de uma média de fluxos, ou rara-

mente, uma média da própria solução;

2. Usa uma solução sensível para forçar qualquer condição de estabilidade não-

linear que implica em TVD, como Upwind Range Condition, por exemplo;

3. Esta condição utilizada limita a ordem de acurácia nos extremos, usualmente

para ordem entre um e dois; e

4. Foi desenvolvido depois da invenção do termo TVD (por volta do ano de 1983).

110

Fig. B.2: Quadro esquemático das condições de Estabilidade Não-Linear [6].

111

Documents

PGMEC · mostra, nos dois modelos, ... 4.11 Comparação entre erros obtidos para malhas de 501 pontos, ... resultado comparado com o obtido por Aktas [7]