BR99H0135

"TRANSFORMAÇÃO INTEGRAL DAS EQUAÇÕES DE NAVIER-STOKES PARA ESCOAMENTO LAMINAR EM CANAIS DE

GEOMETRIA BIDIMENSIONAL ARBITRÁRIA"

Jesus Salvador Pérez Guerrero

TESE SUBMETIDA AO CORPO DOCENTE DA COORDENAÇÃO DOS

PROGRAMAS DE PÓS-GRADUAÇÃO EM ENGENHARIA DA UNIVERSIDADE

FEDERAL DO RIO DE JANEIRO COMO PARTE DOS REQUISITOS

NECESSÁRIOS PARA A OBTENÇÃO DO GRAU DE DOUTOR EM CIÊNCIAS EM

ENGENHARIA MECÂNICA.

Aprovada por:

/ ' . L •• • : --L _

Prof. Renato Machado Cotta, Ph.D.

Prof. Átila Pantaleão Silva Freire, Ph.D.

Prof. Mikhail Dimiter Mikhailov, Ph. D..

.JProf. Sergio H. Sphaier, Dr.-Ing.

Prof. Francesco Scofano NgloJD.Sc.

Prof. Carlos Antonio Cabral dos Santos, D.Sc.

RIO DE JANEIRO, RJ - BRASIL

OUTUBRO DE 1995

, 3 0 - 2 0

PAGES ARE MISSING IN THE

ORIGINAL DOCUMENT

Page 57

PEREZ GUERRERO, JESUS SALVADOR

Transformação Integral das Equações de Navier-Stokes paraEscoamento Laminar em Canais de Geometria BidimensionalArbitrária [Rio de Janeiro] 1995.xvi ,121 p. 29,7 cm (COPPE/UFRJ, D.Sc., EngenhariaMecânica, 1995)Tese - Universidade Federal do Rio de Janeiro, COPPE.1. Equações de Navier-Stokes.2.Escoamento em Canais de Geometria Arbitrária.3. Transformação Integral.I. COPPE/UFRJ II. Título (série)

Ill

Resumo da tese apresentada à COPPE/UFRJ como parte dos requisitos necessários para

obtenção do grau de Doutor em Ciências (D. Sc).

TRANSFORMAÇÃO INTEGRAL DAS EQUAÇÕES DE NAVIER-

STOKES PARA ESCOAMENTO LAMINAR EM CANAIS DE

GEOMETRIA BIDIMENSIONAL ARBITRÁRIA

Jesus Salvador Pérez Guerrero

OUTUBRO, 1995

Orientador: Renato Machado Cotta

Programa: Engenharia Mecânica

O desenvolvimento do escoamento laminar em canais de geometria arbitrária é

estudado, resolvendo-se as equações de Navier-Stokes na formulação de função corrente,

pela Técnica da Transformada Integral Generalizada (GITT).

A função corrente é expandida em autofünções obtidas ao considerar-se unicamente

os termos difusivos da formulação original. As equações de Navier-Stokes são

transformadas num sistema diferencial ordinário infinito, ao fazer-se uso das fórmulas de

transformação e inversão. As séries infinitas são truncadas para fins computacionais, de

acordo com um procedimento de controle automático de erro, resolvendo-se o sitema

diferencial ordinário com subrotinas bem estabelecidas, em bibliotecas matemáticas

amplamente disponíveis.

O clássico caso do desenvolvimento do escoamento entre placas paralelas é

inicialmente analisado, tanto para a condição de entrada uniforme e paralela, como

IV

irrotacional. O efeito do truncamento do comprimento do duto na precisão da solução

obtida é abordado. A convergência das soluções via GITT é apresentada e comparada com

resultados obtidos por diferenças finitas e elementos finitos para diferentes valores do

número de Reynolds.

O caso do escoamento sobre um degrau ("backward-facing step") é estudado a

seguir. Comparações com dados experimentais existentes na literatura indicam uma

excelente concordância. A covalidação numérica é estabelecida para um caso teste,

conseguindo-se perfeita concordância com os resultados considerados benchmark na

literatura recente. Os resultados mostram-se mais coerentes fisicamente que outros

conseguidos por métodos puramente numéricos, em particular para situações onde já se

identificam efeitos tridimensionais.

Por último, um caso teste de duto irregular é estudado e os resultados comparados

com os existentes na literatura, mostrando-se uma boa concordância, com excelentes taxas

de convergência do campo da função corrente ao longo de todo o canal, para diferentes

números de Reynolds.

Abstract of Thesis presented to COPPE/UFRJ as partial fulfilment of the

requirements for the degree of Doctor of Sciences (D. Sc).

INTEGRAL TRANSFORMATION OF THE NAVIER-STOKES

EQUATIONS FOR LAMINAR FLOW IN CHANNELS OF

ARBITRARY TWO-DIMENSIONAL GEOMETRY

Jesus Salvador Perez Guerrero

OCTOBER, 1995

Thesis Supervisor: Renato Machado Cotta

Department: Mechanical Engineering

Laminar developing flow in channels of arbitrary geometry was studied by solving

the Navier-Stokes equations in the streamfunction-only formulation through the

Generalized Integral Transform Technique (GITT).

The stream function is expanded in an infinite system based on eigenfunctions

obtained by considering solely the diffusive terms of the original formulation. The Navier-

Stokes equations are transformed into an infinite system of ordinary differential equations,

by using the transformation and inversion formulae. For computational purposes, the

infinite series is truncated, according to an automatic error control procedure. The ordinary

differential is solved through well-established scientific subroutines from widely available

mathematical libraries.

The classical problem of developing flow between parallel-plates is analysed first,

as for both uniform and irrotational inlet conditions. The effect of truncating the duct

length in the accuracy of the obtained solution is studied. A convergence analysis of the

VI

results obtained by the GITT is performed and compared with results obtained by finite

difference and finite element methods, for diferent values of Reynolds number.

The problem of flow over a backward-facing step then follows. Comparisons with

experimental results in the literature indicate an excellent agreement. The numerical

co-validation was established for a test case, and perfect agreement is reached against

results considered as benchmarks in the recent literature. The results were shown to be

physically more reasonable than others obtained by purely numerical methods, in particular

for situations where three-dimensional effects are identified. Finally, a test problem for an

irregular by shoped duct was studied and compared against results found in the literature,

with good agreement and excellent convergence rates for the streamfunction field along the

whole channel, for different values of Reynolds number.

Vil

A mi adorada hijita Estefanita.

A mi amada esposa y companera de

toda Ia vida Mariluz.

A mi adorada madre Alicia y mi

inolvidable padre Manuel.

VIII

AGRADECIMENTOS

Ao Prof. Renato M. Cotta pela orientação.

À CAPES, pelo suporte financeiro durante todo o doutorado.

Ao Núcleo de Atendimento de Computação em Alto Desempenho (NACAD/UFRJ)

pela utilização de seus recursos computacionais.

À minha amada esposa Mariluz, pela paciência e apoio durante todo o trabalho.

Aos meus irmãos Manuel e Miguel, pelo incentivo.

Aos meus grandes amigos Roseane Barcelos Santos e Ricardo Gondim dos Santos,

pela inestimável ajuda prestada neste trabalho.

A Alípio Castelo Branco Jr. e Luis Cláudio Gomes Pimentel pela colaboração.

A todo pessoal do Laboratório de Transmissão e Tecnologia do Calor (LTTC), pela

colaboração sempre atenta em todos os momentos.

Aos professores da COPPE/UFRJ pelos conhecimentos fornecidos em suas aulas.

Ao Dr. José Ananias e Sra. Ana Lúcia Figueira da Silva pela ajuda num momento

crítico.

IX

GLOSSÁRIO

- LETRAS LATINAS

a Posição para especificar a altura do degrau.

A jjk Coeficiente de transformação.

A o,jjç Coeficiente de transformação.

b Metade da separação entre placas paralelas

B yk Coeficiente de transformação.

b j , b2 Especifica separação das paredes do duto.

B<»jk Coeficiente de transformação.

c Parâmetro de compressão da escala.

C jjk Coeficiente de transformação.

Coojk Coeficiente de transformação.

D j : Coeficiente de transformação.

D h Diâmetro hidráulico (Dh = 4b).

F Filtro para caso de domínio irregular.

FOOjjk Coeficientes de transformação do duto irregular

FOljjk Coeficientes de transformação do duto irregular

F02jjk Coeficientes de transformação do duto irregular

F03jjk Coeficientes de transformação do duto irregular

F lop Coeficientes de transformação do duto irregular

Fl ljjk Coeficientes de transformação do duto irregular

F12;jk Coeficientes de transformação do duto irregular

FOy Coeficientes de transformação do duto irregular

Fly Coeficientes de transformação do duto irregular

F2jj Coeficientes de transformação do duto irregular

F3jj ' Coeficientes de transformação do duto irregular

fj Distribuição de velocidade axial na entrada do duto (caso dutos

regulares).

Distribuição de velocidade axial na saída do duto (caso dutos

irregulares).

Coeficientes de transformação do duto irregular

FCDj Coeficientes de transformação do duto irregular

Coeficientes de transformação do duto irregular

Coeficientes de transformação do duto irregular

Coeficientes de transformação do duto irregular

Coeficientes de transformação do duto irregular

Coeficientes de transformação do duto irregular

Subíndice de termo na série.

Subíndice de termo na série.

Subíndice de termo na série.

Valor da função-corrente nas paredes do duto.

Posição de truncamento do duto.

Número de termos de truncamento da série.

Norma.

Campo de pressões adimensional e dimensional respectivamente.

Coeficiente de balanço de massa.

Vazão por unidade de comprimento.

( A\Jb\Re Número de Reynolds Re = I (Caso do Duto Regular);

Re = — (Caso Geral do Duto Irregular),v J

Velocidade axial adimensional e dimensional respectivamente.

Velocidade média axial.

Velocidade transversal adimensional e dimensional respectivamente.

Posição de realocamento da condição de contorno transformada.

Função do sistema diferencial ordinário.

Coordenada vertical adimensional e dimensional respectivamente.

Autofunção.

Autofunção normalizada

Hl

H2

H3

i

j

k

L

N ,

N;

P ,

q

Q

Ü

Ü

Ü

,k2

Nt

,N*

P

u ,U

v ,

*o

X

y-

u

•V

*

y

XI

- LETRAS GREGAS

s Controle automático do erro global.

4 Variável Auxiliar.

r\ Variável de transformação do domínio.

|ij Autovalores.

fj.; Autovalores.

v Viscosidade.

p Densidade.

\|/ Função-corrente.

\\im Função-corrente do campo de velocidade completamente desenvolvido.

\j7j Campo Transformado, caso do duto domínio Írrregular.

co Vort ic idade.

-SÍMBOLOS MATEMÁTICOS:

V Gradiente.

V2 Laplaciano.

V 4 Operador biharmônico

xn

ÍNDICE DE FIGURAS

FIGURA 2 .1 -

FIGURA 2.2 -

FIGURA 3.1-

FIGURA 3.2-

FIGURA 3.3 -

FIGURA 3.4 -

FIGURA 4 .1 -

FIGURA 4.2 -

FIGURA 4.3 -

FIGURA 4.4 -

FIGURA 4.5 -

FIGURA 4.6 -

FIGURA 4.7 -

FIGURA 4.8 -

FIGURA 4.9 -

FIGURA 5.1-

FIGURA 5.2 -

FIGURA 6.1-

FIGURA 6.2.a,b

FIGURA 6.3 -

Representação geral do problema em duto regular 5

Especificação das condições de contorno 11

Tipos de condições de contorno na entrada do duto 25

Velocidade u (x, 0) ao longo do duto para condições de contorno

uniforme e irrotacional na entrada, comparada com resultados de

camada limite 34

Desenvolvimento do perfil de velocidades u (x, y), ao longo do duto

para condições de contorno uniforme e irrotacional na entrada

(Re = 300) 34

Desenvolvimento do perfil de velocidades v (x, y), com condição de

contorno uniforme e paralela na entrada do duto 35

Esquema do problema de escoamento sobre um degrau ("Backward

Facing Step") 37

Desenvolvimento da componente transversal da velocidade ao longo

do canal (Re = 438) 46

Desenvolvimento da componente transversal da velocidade ao longo

do canal (Re = 1374) 46

Isolinhas da Função Corrente (Re = 438) 47

Isolinhas da Função Corrente (Re = 1374) 47

Comparação dos Perfis da componente transversal da velocidade ao

longo do duto (x = 14 e x = 30 ; Re = 1600) 54

Desenvolvimento da velocidade transversal ao longo do duto

(Re = 1600) 54

Isolinhas da Função Corrente (Re = 1600) 55

Isolinhas da Função Corrente (Re = 1942) 55

Esquema do duto irregular 57

Geometria que identifica o filtro F (x, y) 61

Geometria do problema-teste: Canal com expansão gradual 83

Condições de contorno no duto com espansão gradual 84

Isolinhas da Função Corrente (Re = 10) 93

Xlll

FIGURA 6.4 - Isolinhas da Função Corrente (Re = 100) 93

FIGURA 6.5 - Vorticidade ao longo da parede do duto (Re =10) 98

FIGURA 6.6 - Vorticidade ao longo da parede do duto (Re =100) 98

FIGURA A. 1 - Geometria do Escoamento Completamente Desenvolvido entre

placas paralelas 109

XIV

ÍNDICE DE TABELAS

TABELA 3.1-

TABELA 3.2-

TABELA 3.3 -

TABELA 3.4-

TABELA 3.5.a,b-

TABELA 4.1-

TABELA 4.2 -

TABELA 4.3 -

TABELA 4.4 -

TABELA 4.5 -

TABELA 6.1-

TABELA 6.2 -

TABELA 6.3. a,b -

TABELA 6.4.a,b -

Convergência da velocidade no centro do duto, u (x, 0), para

Re=0. (u= 1, v = 0) 29

Convergência e comparação da velocidade no centro do duto,

u (x, 0), para Re = 40. ( condições de contorno na entrada u = 1,

v = 0) 30

Convergência e comparação da velocidade no centro do duto,

u (x, 0), para Re = 300. (condições de contorno na entrada u = 1,

v = 0) 30

Convergência e comparação da velocidade no centro do duto,

u (x, 0), para Re = 300. (condições de contorno na entrada u = 1,

o = 0) 31

Velocidade na linha de centro do duto, u (x, 0) , para diversas

posições de truncamento do duto 32

Convergência da componente longitudinal do campo de

velocidades (Re = 438 ) 44

Convergência da componente longitudinal do campo de

velocidades (Re= 1374) 45

Convergência da componente longitudinal do campo de velocidade

emx=14 (Re=1600) 49

Convergência da componente longitudinal do campo de velocidade

em x=30 (Re=1600) 50

Convergência da componente longitudinal do campo de velocidade

(Re = 1942) 53

Convergência da Função Correntei)/ (x/ xout,0.5) (Re = 1 0 ) 91

Convergência da Função Corrente \\f (x / xout,0.5) (Re =100 ) 92

Convergência da vorticidade na parede do duto (Re =10) 95

Convergência da vorticidade na parede do duto (Re = 100) 96

XV

ÍNDICE

Página

CAPÍTULO I - INTRODUÇÃO 1

LI - MOTIVAÇÃO E OBJETIVOS 1

1.2 - A TÉCNICA DA TRANSFORMADA INTEGRAL GENERALIZADA 4

CAPÍTULO II - ESCOAMENTO EM DUTOS REGULARES 5

II. 1 - DEFINIÇÃO DO PROBLEMA 5

II.2 - SOLUÇÃO VIA GITT 9

11.2.1 - Problema Auxiliar 10

11.2.2 - Par Transformada-Inversa 13

11.2.3 - Transformação Integral do Problema 15

11.2.4 - Algoritmo Computacional 18

CAPÍTULO III - ESCOAMENTO ENTRE PLACAS PARALELAS 24

111.1 - DESCRIÇÃO DO PROBLEMA 24

111.2 - REVISÃO BIBLIOGRÁFICA 27

111.3 - RESULTADOS 28

CAPÍTULO IV - ESCOAMENTO SOBRE UM DEGRAU 36

IV.l - DESCRIÇÃO DO PROBLEMA 36

IV.2 - REVISÃO BIBLIOGRÁFICA 40

IV.3 - RESULTADOS 43

CAPÍTULO V - ESCOAMENTO EM CANAIS DE GEOMETRIA ARBITRÁRIA. 56

V.l - DEFINIÇÃO DO PROBLEMA 56

V.2 - SOLUÇÃO VIA GITT 59

V.3 - O FILTRO F (x, y) 60

V.4 - PROBLEMA DE AUTOVALOR 63

V.4.1 - Propriedade de Ortogonalidade 65

V.5 - PAR TRANSFORMADA-INVERSA 65

XVI

V.6 - TRANSFORMAÇÃO INTEGRAL DO PROBLEMA 67

V.7 - ALGORITMO COMPUTACIONAL 77

CAPÍTULO VI - ESCOAMENTO NO INTERIOR DE UM DUTO COM EXPANSÃO

GRADUAL 81

VI. 1 - DEFINIÇÃO DO PROBLEMA 81

VI.2 - REVISÃO BIBLIOGRÁFICA 85

VI.3 - SOLUÇÃO DO PROBLEMA 86

VI.3.1 - Duto Truncado 86

VI.3.2 - Duto Infinito 88

VI.4 - RESULTADOS 89

CAPÍTULO VII - CONCLUSÕES 99

CAPÍTULO VIII - REFERÊNCIAS BIBLIOGRÁFICAS 101

APÊNDICE A - Escoamento completamente Desenvolvido entre Placas Paralelas 109

APÊNDICE B - Derivadas da Fórmula de Inversão 112

APÊNDICE C - Coeficientes do Sistema Diferencial Ordinário 114

I - INTRODUÇÃO

I.I - MOTIVAÇÃO E OBJETIVOS

O estudo do desenvolvimento hidrodinâmico de escoamentos laminares em dutos é

um tema de interesse permanente na engenharia pelas inúmeras aplicações em

equipamentos hidráulicos e térmicos, na busca pela determinação acurada de fatores de

atrito e coeficientes de transferência de calor ao longo do canal. A modelagem teórica deve

ser realizada a partir das equações de Navier-Stokes, ou através das equações simplificadas

de camada limite. A validade da aplicação neste último caso é circunscrita a situações de

número de Reynolds altos e em posições não muito próximas da entrada do duto. Este

aspecto é analisado no Capítulo III do presente trabalho.

A dificuldade na análise numérica do desenvolvimento de escoamentos em dutos,

está na adequada especificação de condições de contorno na entrada e saída do canal e na

necessidade de se resolver, com precisão controlada, um sistema diferencial parcial elíptico

e não linear, representado pelas equações de Navier-Stokes.

Soluções para problemas de escoamento dentro de dutos foram obtidas ainda em

meados dos anos 60 para o caso de placas paralelas, ampliando-se os estudos com a

evolução dos métodos numéricos para resolver equações diferenciais parciais e dos

computadores, tanto em velocidade quanto em memória. Ainda assim, a gama de análises

para escoamentos em canais não-regulares permanece restrita, ao introduzir regiões de

recirculação e adicionar a dificuldade do adequado mapeamento do domínio irregular.

Por outro lado, problemas convectivo-difusivos têm sido resolvidos com grande

sucesso pela Técnica da Transformada Integral Generalizada ( GITT ) [ 1 ], em inúmeras

aplicações.

Por sua concepção híbrida numérico-analítica, esta técnica tem as vantagens de

manter um controle sobre o erro relativo dos resultados, o qual é estabelecido a priori e

controlado automaticamente, e de ser uma ferramenta de simples implementação

computacional.

A ampla gama de aplicações do método, pode ser agrupada nas seguintes classes:

- Problemas com coeficientes variáveis nas condições de contorno [ 2-5 ]:

• Aplicações: problemas condutivos com número de Biot dependente do tempo,

convecção forçada em dutos aletados externamente, etc.

- Problemas com coeficientes variáveis na equação [ 6- 8 ]:

• Aplicações: análise transiente de aletas com dissipação dependente do tempo,

desenvolvimento simultâneo de velocidade e temperatura dentro de canais, etc,

- Problemas com Contornos Variáveis [ 9-15 ]:

• Aplicações: problemas de contorno móvel, tais como fusão/solidificação,

ablação e oxidação; domínio de forma irregular com respeito ao sistema de

coordenadas considerado, como no caso da transferência de calor em

trocadores com dutos triangulares, trapezoidais, hexagonais, etc.

- Problemas que Envolvem Problema Auxiliar Difícil [ 16-22 ]:

• Aplicações: convecção forçada interna em regimes transiente e periódico,

transferência de calor em escoamentos com condução axial, convecção em

dutos retangulares, problemas conjugados de transferência de calor, trocadores

de calor bitubulares, etc.

- Problemas Não- Lineares [23-41, 76-78 ]:

• Aplicações: condução de calor com condutividade térmica variável, condições

de contorno com troca radiante, solução das equações de camada limite para

regimes laminar e turbulento, e as equações de Navier-Stokes.

Com relação às equações de Navier-Stokes, a técnica está estendendo-se

gradativamente para problemas que envolvem escoamento laminar e turbulento e para as

aplicações em coordenadas cilíndricas.

Neste contexto, o objetivo do presente trabalho é estender a aplicabilidade da GITT

na solução das equações de Navier-Stokes para casos mais gerais, explorando a formulação

em função corrente. Desta forma, pretende-se abordar a solução do problema de

escoamento em desenvolvimento no interior de dutos bidimensionais com domínio regular

e irregular para o caso do escoamento incompressível e laminar em regime permanente,

abrindo-se assim perspectivas para aplicação da presente técnica em casos bem mais gerais

com dependência no tempo e em regime turbulento. Os problemas de escoamento

resolvidos no presente trabalho têm sido escolhidos, tanto pelo interesse físico que eles

apresentam, como também porque servem como excelentes testes de validação de métodos

computacionais.

O presente trabalho encontra-se inserido no escopo da cooperação técnica entre a

University College of Swansea, Institute for Numerical Methods in Engineering, e o

PEM/COPPE/UFRJ, sob fmancimento do British Council. Neste esforço mais amplo,

busca-se a co-validação de resultados de simulações via transformação integral e elementos

finitos, fazendo-se uso da interessante característica de controle automático de erro da

presente técnica, permitindo o estabelecimento de resultados benchmark para problemas -

teste de uso corrente na literatura, como aqueles aqui analisados.

1.2 - A TÉCNICA DA TRANSFORMADA INTEGRAL GENERALIZADA

A aplicação da GITT pode ser resumida nos seguintes passos:

a) Definição do problema auxiliar, com base nos termos difusivos da formulação

original.

b) Solução do problema auxiliar e obtenção das autofunções, autovalores, normas

e propriedade de ortogonalidade.

c) Desenvolvimento do par transformada-inversa.

d) Transformação integral do problema diferencial parcial num sistema diferencial

ordinário acoplado.

e) Truncamento do sistema diferencial ordinário infinito e solução numérica do

sistema diferencial resultante, para obtenção dos campos transformados.

f) Obtenção do potencial original, fazendo-se uso da fórmula de inversão.

A idéia básica na técnica generalizada, é relaxar a necessidade de encontrar-se uma

transformação integral exata, ou seja, que resulte num sistema diferencial transformado

não-acoplado. Assim, pode-se escolher um problema auxiliar (de autovalor) que seja

característico do problema original ou não, desenvolver o par transformada-inversa e

efetuar a transformação integral, chegando-se a um sistema diferencial ordinário infinito e

acoplado. Após truncamento em ordem suficientemente grande para a precisão requerida,

automaticamente selecionada durante o próprio processo de solução, o sistema diferencial

ordinário é resolvido numericamente por algoritmos bem estabelecidos, com controle

automático de erro, disponíveis em bibliotecas de subrotinas científicas. A fórmula

explícita de inversão fornece então uma representação analítica nas demais variáveis

independentes eliminadas pela transformação integral.

II - ESCOAMENTO EM DUTOS REGULARES

II. 1 - DEFINIÇÃO DO PROBLEMA

O problema de escoamento em duto regular é definido geometricamente por duas

placas planas infinitas e paralelas, entre as quais escoa um fluido Newtoniano com

condições de contorno especificadas na entrada e saída do canal (Figura 2.1).

Deseja-se encontrar as características do desenvolvimento do escoamento ao longo

do duto.

f.[y] f2[y*]

-b

FIGURA 2.1 - Representação geral do problema em duto regular.

O problema físico é formulado ao considerar-se:

- Escoamento bidimensional em regime laminar e estacionário.

- Escoamento incompressível.

- Propriedades do fluido constantes.

- Impermeabilidade e não deslizamento nas paredes do duto.

- Fluido newtoniano.

O escoamento dentro do duto é governado pelas equações de continuidade e de

Navier-Stokes, que em termos das variáveis primitivas adimensionais; são dadas por:

õx õy

õu. õu Õp 4 f õ V õ2n)

u

õx õy õx Re^õx õy

õv õv õp 4 j õ2v õ2v

õx õy õy Rel^õx õy J

onde os parâmetros adimensionais foram definidas a partir de:

(2.1.C)

(2.2.a,b)

* * *u = — , v = — , P = - ^ T (2.2.C, d, e)

U U ' V pU2

notando-se que "*" identifica as variáveis dimensionais. O número de Reynolds foi

definido a partir do diâmetro hidráulico e da velocidade média na entrada do duto:

UDh 4UbRe = = (2.2.f)

Como é sabido, uma das dificuldades na solução das equações (2.1) está no

desconhecimento do campo de pressões nos contornos.

Uma alternativa que permite evitar esta dificuldade, é fazer uso da função-corrente

\)/ (x,y), representada a partir do campo de velocidades, na forma:

| ( , y ) ôv|/(x,y)u(x ,y) = — , v(x,y) = - — (2.3.a,b)

Portanto, derivando a equação (2.1.b) em relação a "y" e a equação (2.1.c) em

relação a "x", com posterior subtração de ambas e tendo-se em conta (2.1.a), obtém-se a

seguinte equação não linear do tipo biharmônico para a função-corrente:

ôfy ô\\i Õ3\\i dy õ\ dy di|/ 4 4+ = V VI/ (Z.T-)

Re Yôy õx3 õy Õxdy2 õx õx2õy õx õy3 Re

válida para todo x > 0 , -1 < y < 1 e onde V4 representa o operador biharmônico:

õx ôx õy ( 2 - 5 )

As oito condições de contorno requeridas na solução da equação (2.4) são

especificadas ao considerar-se:

- Impermeabilidade e não-deslizamento nas paredes do duto:

u ( x , - l ) = 0 ; v ( x , - l ) = 0 ; x > 0 (2.6.a, b)

u ( x , l ) = 0 ; v (x , l ) = 0 ; x > 0 (2.6.c, d)

que ao ser representadas em termos da função-corrente tornam-se:

— = 0 (2.7.a, b)

^ = 0 (2.7.c,d)

onde k[ e k2 são constantes, que especificam o valor da função-corrente nas paredes.

- Perfil de velocidades conhecido na entrada do duto:

u(0,y) = f,(y) ; v(0,y) = 0 (2.8.a,b)

Usando a definição da função-corrente (2.3) e integrando-se diretamente:

j (Ç)d^ ; ^ ' y ) = 0 (2.9.a,b)

- Perfil de velocidades paralelo completamente desenvolvido para x —> oo :

; v(oo,y) = o (2.io.a,b)

ou em termos da função-corrente:

onde "q" é uma constante que garante o balanço de conservação da massa.

A constante k2 pode ser encontrada a partir de (2.9.a) ou (2.11.a). Fazendo-se uso

desta última:

k 2 = k t + 2 q (2.12)

Desta forma, o problema de escoamento em duto regular fica formulado apenas em

termos da função-corrente.

II.2 - SOLUÇÃO VIA GITT

Para facilitar a aplicação do método (GITT) e melhorar sua performance

computacional, é conveniente fazer uma homogeneização parcial das condições de

contorno. Para este fim, se considera o seguinte filtro:

onde ^ ( y ) representa a função-corrente na região completamente desenvolvida

M̂7 (°°>y)) e •t' (x>y) é o potencial a ser determinado.

Substituindo (2.13) em (2.4), e nas condições de contorno (2.7,9-11) obtém-se:

10

ô(j) cP_§_ dty õ3 (j) 3(j) a3 <|> S(() ô3 <|)

õy õx3 õy dxôy2 õx õx2õy ~ õx ôy3

| K, 5 ( | d i t / M a>j, = _4_ v

dy õxdy 2 dy 3 õx Re|

dy 5x 3 dy õxdy 2 dy 3 õx Re(2.14.a)

(2.14.b,c)5y

4õy

o, y)

(2.14.d,e)

(2.14.f,g)

Assim, o problema encontra-se reformulado para a solução da função auxiliar

<Kx,y).

II.2.1 - PROBLEMA AUXILIAR

Dado às características homogêneas do sistema (2.14) na direção "y", escolheremos

o problema de autovalor nessa coordenada.

= fi(y)

11

•• y

u=0 v=0

u = 0 v = 0

u = -q(i-y2)

x v = 0

jf,- l

©d

ax

¥

yY|/ = k, + 2q —— = 0

—y x

3 y3

y

ox

j— 1

J

y = 0 — = 0

X

— = 0

FIGURA 2.2 - Especificação das condições de contorno.

12

Um problema de autovalor associado à versão linear homogênea de (2.14),

detalhadamente discutido em [34] e já empregado na solução das equações de Navier-

Stokes via GITT [34-36,41] é definido como:

(2.15.a)

com condições de contorno:

(2.15.b,c)

(2.15.d,e)

onde Y[ (y) e ji; são as autofunções e autovalores, respectivamente, que satisfazem a

seguinte propriedade de ortogonalidade:

, i * j

e Nj é a norma ou integral de normalização.

A equação (2.15) pode ser analiticamente resolvida, encontrando-se uma

combinação linear de funções trigonométricas e hiperbólicas:

Y,(y) = i sinu^y sinhfj-jy, l — 2,4,0,...

(2.18.a-b)

Os autovalores são encontrados ao resolver-se as seguintes equações

transcendentais:

13

, , , ,

. . . . (2.19.a-b)i , 1 = 2,4,6,... v '

A norma, definida como:

i

Y2dy ; i = l , 2 , 3 , . . . (2.20)

apresenta o seguinte valor numérico:

Ni = 2 ; i = 1,2,3,... (2.21)

Para os próximos passos no uso da GITT, é conveniente normalizar a autofunção a

partir de:

Nf2

onde Yj (y) representa a autofunção normalizada.

II.2.2 - PAR TRANSFORMADA INVERSA

O uso da GITT baseia-se na idéia de que uma função pode ser representada como

uma expansão em autofunções, provenientes de um problema auxiliar, o qual tem

informação sobre os termos difusivos do problema original. Assim, consideramos que a

função (J) (x, y) pode ser representada como:

14

onde <(» i (x) é uma função desconhecida, que depende apenas de "x", representando os

coeficientes a serem determinados na expansão proposta.

Uma relação que permita definir fo (x), pode ser encontrada ao se fazer uso das

propriedades de ortogonalidade da autofunção Yj.

Operando ambos os membros de (2.23) com í Vy :

dy (2.24)J Yj(y)*(x,y)dy= J Yj (y)

e reordenando convenientemente:

J Yj(y)4>(x,y)dy= X <f>i(*) J % (y) Yj(y) dy (2.25)

Pelas propriedades de ortogonalidade que gozam as autofunções Yj , pode-se

facilmente deduzir que os termos do somatório em (2.25) tem um único valor diferente de

zero quando i = j , logo:

fc(x)= JYi(y)(j)(x,y)dy (2.26)

Esta última expressão define a transformação integral, indicando que o campo

original pode ser transformado pelo seu produto interno com a autofunção.

Denomina-se potencial transformado a <t>j(x), sendo a eq. (2.26) conhecida como

fórmula da transformada. A eq. (2.23), que por sua vez recupera o potencial original a

partir do conhecimento dos campos transformados, denomina-se fórmula da inversa.

15

II.2.3 - TRANSFORMAÇÃO INTEGRAL DO PROBLEMA

O processo de transformação integral das equações diferenciais parciais (2.14) num

sistema diferencial ordinário, é iniciado ao se fazer uso do operador J Yj (y) dy na

equação (2.14.a):

dv + fi ôy ax3 d y + J Yi dy õxôy

Re L -I

(2.27)

Neste ponto, é conveniente fazer uso da definição da inversa (2.23), que facilitará

reagrupar convenientemente os termos e aplicar as propriedades de ortogonalidade onde for

possível.

A utilização da inversa assegura a transformação de uma equação diferencial

parcial, inclusive onde os termos não forem transformáveis exatamente.

16

dY,dy

dy +hi d y Lfcf d x

dxL k = l dx2 dy

d y -

d y -

JY«L k = l

dy + 1 dy dxJ -Y: dy +

dx dy2 dy - h dx

y d2^ft dx2yft dx2 dy2 dy +

dy

(2.28)

Reordenando convenientemente e aplicando as propriedades de ortogonalidade

(2.16) e (2.17), no primeiro e terceiro termos do lado direito de (2.28):

dx4

00 00

ZZj=l k=l

C i j^ dx

ijk dx dx2

d ^ijk d

A

dx3C

dx U"'J dx

(2.29.a)

onde os coeficientes podem ser calculados analiticamente a partir das integrais:

17

l\ : :u —

d y-dy (2.29.b)

B-., = Y Y - d v (2.29.c)

'ijk - J iY i dy-dy (2.29.d)

1 A2

(2.29.e)

A»ij = i ^i^j-T^dydy

(2.29.Í)

ij J Y i dy2 dy-dy (2.29.g)

»u = J Y i Y J~x^ d y-i J dy-(2.29.h)

Note-se que a equação (2.29.a), apresenta também as mesmas características não-

lineares, as quais não são removidas na transformação integral. O mesmo procedimento de

transformação é aplicado nas condições de contorno (2.14.f, g, h, i), obtendo-se

respectivamente:

dy (2.30.a)

^j (0)

dx= 0 (2.30.b)

18

O») = O (2.30.c)

dx= O (2.30.d)

Desta forma, o problema diferencial parcial foi transformado num sistema de

equações diferenciais ordinárias não lineares, infinito e acoplado, com condições de

contorno em dois pontos.

II.2.4 - ALGORITMO COMPUTATIONAL

O sistema diferencial ordinário não-linear (2.29) deve ser resolvido a partir de

procedimentos numéricos devido à impossibilidade de se obter soluções analíticas.

Por outro lado, obter soluções exatas do sistema (2.29) implicaria em se resolver o

sistema infinito de equações, o que seria impossível de ser realizado sob o ponto de vista

computacional.

Uma característica importante do método é a garantia de convergência das soluções

para ordens crescentes de truncamento nas séries. Esta característica indica que é possível

obter soluções com um número de algarismos significativos "exatos" (convergidos) para

um determinado número de termos nas expansões. Por isso, diz-se que o método da

transformada integral é um método de precisão controlada, devido a esse controle

estabelecido na ordem de truncamento das expansões, que pode ser automaticamente

determinado durante o processo de solução, assemelhando bastante o presente

procedimento ao de uma solução puramente analítica.

Após truncamento do sistema infinito numa ordem N, a eq.(2.29) torna-se:

19

dx4 ar d<t>k

" A i ji * dx dx2 " A i j k d x * k

N

I °°ij dx3 ^ dx °°ij dx

(2.31)

Por outro lado, as condições de contorno (2.30.C, d) encontram-se especificadas no

infinito. Nos métodos numéricos clássicos, essa dificuldade é tradicionalmente resolvida

ao considerar-se as condições de contorno no infinito especificadas numa posição

suficientemente distante da origem, resolvendo-se assim o problema para essas condições

artificiais. Mas, como não é possível saber-se a priori, se tal aproximação pelo

truncamento do domínio infinito afeta os resultados finais dentro da precisão desejada,

torna-se indispensável resolver-se repetidamente o problema, considerando-se as condições

de contorno em outras posições e verificando se a convergência das soluções foi atingida

nas regiões de interesse.

Esta dificuldade é facilmente contornada pelo presente método, por uma

transformação da variável independente x, redefinindo o domínio de [ 0, oo ] a [ 0, 1 ].

A transformação do domínio aqui considerada foi:

-ex (2.32)

onde c é um parâmetro de compressão de escala.

Um aspecto interessante desta transformação é que ela permite comprimir a escala

onde os gradientes são menos significativos.

O uso da transformação do domínio é completado ao definir-se:

20

(2.33)

Assim, o problema a ser resolvido seria (2.31), usando-se a equação (2.32) em

conjunto com a (2.33), ficando as condições de contorno (2.30.C, d) especificadas em r) =

1.

A solução numérica deste sistema pode ser obtida através de algoritmos bem

estabelecidos e testados, disponíveis em bibliotecas de subrotinas matemáticas, como o

sistema IMSL [42]. Esta biblioteca possui a subrotina DBVPFD para a solução de

problemas diferenciais ordinários não lineares, com condições de contorno em dois pontos,

que apresentam comportamento rígido (stiff) .

Esta rotina resolve problemas do tipo:

x ' = f ( x , X ) , x e [ a ,b ] (2.34.a)

g[X(a),X(b)]=0 (2.34.b)

Seu algoritmo é baseado na rotina PASVA3 [43], e faz uma discretização sobre

uma malha não uniforme, que é escolhida adaptativamente para conseguir um erro local

aproximadamente o mesmo em qualquer posição, mantendo-se um controle de erro global

automático. O sistema algébrico não-linear resultante é resolvido mediante o método de

Newton generalizado.

Para fazer uso da DBVPFD é necessário reescrever o sistema de quarta ordem

como um sistema de primeira ordem:

Xj = fc (2.35.a)

21

dx

dx

(2.35.b)

(2.35.c)

X 3N+i " dx3 (2.35.d)

dX3N+i

dx dx1 (2.35.e)

e, auxiliando-se da regra da cadeia, é possível reescrever em termos do domínio

transformado r\:

(2.35.Í)

.fin'sendo —- uma função que depende apenas de r[ e do parâmetro c.Vdx;

O sistema diferencial ordinário pode ser reescrito como:

dX 3N + iN

dTJ

Re4

N Nv —/ , (k=l

- B; j k X N + j X 2 N + k - A j j k X N + j X k

X3N+j +B°oij

parai= 1, 2,..., N (2.36.a)

22

As outras 3N equações que completam o sistema diferencial de primeira ordem são

dadas por:

dX;(2.36.b)

com as condições de contorno:

Xi(0)=jYi JL

f,©dH Ty-V dy (2.36-c)

• = o (2.36.d)

(2.36.e)

dX t(l) = 0 (2.36.Í)

para i = 1, 2 , . . . , N.

Obtidos assim os potenciais transformados, o campo de função-corrente pode ser

calculado a partir da fórmula de inversão (2.23), considerando também (2.13). O campo de

velocidades, conseqüentemente, será obtido a partir de (2.3.a,b).

Como todas as tarefas computacionais intermediárias são realizadas dentro dos

limites de precisão prescritos pelo usuário, resta controlar a ordem de truncamento das

expansões e, em decorrência, a ordem do sistema diferencial ordinário, para se chegar ao

controle automático do erro global. Para tanto, faz-se uso da natureza analítica da presente

técnica, realizando-se um teste de convergência interno no algoritmo em cada posição onde

se deseja obter a solução, a partir da fórmula:

s = max( x , y )

23

N+AN

i=N+l

N+AN

(2.37)

incrementando-se em intervalos, o número de termos na expansão, até que o valor de s

satisfaça a precisão requerida em todo o domínio. Os valores obtidos com ordens N

menores, ainda não completamente convergidas, oferecem excelentes aproximações

iniciais no procedimento iterativo da subrotina DBVPFD para as ordens superiores, N +

AN.

24

III - ESCOAMENTO ENTRE PLACAS PARALELAS

III.l - DESCRIÇÃO DO PROBLEMA

Este clássico problema consiste no estudo do desenvolvimento do escoamento

laminar entre duas placas paralelas semi-infinitas, sendo que na entrada o escoamento

preenche completamente as placas paralelas e na saída (no infinito) o escoamento é

completamente desenvolvido.

Segundo o desenvolvido no Capítulo II, o problema estará completamente

formulado ao considerar-se as condições de contorno na entrada. Dois tipos de condições

de contorno são utilizadas com maior frequência na literatura, sendo que a mais usual

(Figura 3.1.a) especifica que a velocidade longitudinal é uniforme e constante, com

escoamento paralelo:

u = l , v = 0 (3.1.a,b)

Um outro tipo de condição de contorno especifica que o perfil de velocidades é

uniforme e irrotacional na entrada:

5vu = l , co = 0 | — = 0 | (3.2.a,b)

u = l

v=0

u = l

25

u = 0 v = 0

u = 0 v = 0

u = 0 v = 0

u = 0 v = 0

v = 0

v = 0

u=0 v=0

u = l

co=0

u = 0 v = 0

u = 0 v = 0 =0 v = 0

v = 0

FI GURA 3.1 - Tipos de condições de contorno na entrada do duto.

26

Ambas as condições de contorno (3.1) e (3.2) apresentam uma descontinuidade em

y = ± 1. Para evitar esta descontinuidade, um outro tipo de condição de contorno é também

usado, no qual o perfil de velocidades é idealizado como uniforme e paralelo em x = -oo

[44].

No presente trabalho se utilizará a condição (3.1) , principalmente por ter sido a

mais comumente usada em trabalhos anteriores, o que permitirá fazer comparações de

soluções via GITT com outros métodos.

Também se considerará a condição de contorno (3.2) , para compará-la com as

soluções a partir de (3.1) e estabelecer sua influência no desenvolvimento do perfil de

velocidades.

No caso das condições (3.1) , o perfil de velocidades na entrada e o coeficiente de

balanço q são dados por:

u = f 1 ( y ) = l (3.3)

q = l (3.4)

enquanto as condições de contorno na entrada, segundo (2.36.c,d), tornam-se:

X, (0 )= J ^ — - - y | d y (3.5.a)

dX;(0)- ^ = 0 ; i = l , 2 , 3 , . . . ,N (3.5.b)

No outro caso, com condições de contorno segundo as equações (3.2), o perfil de

velocidade f, (y) e o coeficiente de balanço q são os mesmos, visto que a condição de

contorno (3.5.a) se mantém inalterada. O efeito de irrotacionalidade é imposto através das

condições:

27

d Xj (0)= 0 ; i = l , 2 , 3 , . . . , N (3.6)

III.2 - REVISÃO BIBLIOGRÁFICA

O problema do desenvolvimento do escoamento laminar entre placas paralelas tem

sido abordado na literatura com uma frequência razoável.

A primeira solução numérica das equações de Navier-Stokes para o escoamento

entre placas paralelas foi feita por Wang e Longwell [45], publicada no clássico artigo de

1964. Estes autores formularam o problema em termos da vorticidade e função corrente,

especificando escoamento irrotacional na entrada (casos b, c da Fig. 3.1). Consideraram

uma transformação da coordenada axial para um tratamento mais adequado da condição de

contorno no infinito. A implementação computacional foi feita através do método de

diferenças finitas, reportando-se resultados convergidos até Re = 300. Mostraram ainda

que o desenvolvimento da velocidade longitudinal era caracterizado por uma concavidade

com um mínimo local no centro do duto.

Brandt e Gillis [46] formularam o problema com base na função corrente,

resultando em um problema do tipo biharmônico-não linear. Truncaram o domínio infinito

em distâncias muito afastadas da entrada para especificar condições de escoamento

completamente desenvolvido. Na entrada consideraram velocidade uniforme e paralela

(u = 1, v = 0) e, usando o método de diferenças finitas, obtiveram resultados até Re = 1000.

McDonald, Denny e Mills [47] empregaram uma formulação de vorticidade e

função corrente e consideraram dois tipos de condição de contorno na entrada: velocidade

uniforme e paralela, equações (3.1) e irrotacional, equações (3.2). O método empregado foi

diferenças finitas, usando uma malha de 21 x 201 para resolver o caso de Re = 300000. O

método de solução foi estendido até Re = 300 para o caso de escoamento irrotacional na

entrada, mas aparentemente com pouca precisão devido à malha grosseira utilizada (21 x

41). Os resultados para Re = 300 , apresentam uma boa concordância com Brandt e Gillis

28

[46] no caso de entrada uniforme e paralela, sendo que no caso irrotacional os resultados de

Wang e Longwell [45] encontram-se ligeiramente diferentes daqueles reportados em [47].

Narang e Krishnamoorthy [48] linearizaram os termos de inércia na equação de

Navier-Stokes e assumiram velocidade não uniforme na entrada do duto, truncando o

domínio infinito. Obtiveram uma representação analítica da função corrente, e resolveram

um sistema de 150 equações para determinar os coeficientes desta solução semianalítica,

obtendo resultados aproximados desde Re = 2 até 4000.

Comini e Del Guidice [49] empregaram o método de elementos finitos para uma

formulação em variáveis primitivas, considerando velocidade uniforme na entrada do

canal, e reportando resultados para Re = 40.

III.3 - RESULTADOS

Os cálculos foram realizados no VAX8810, e um erro relativo de 10"4 foi imposto

na solução do sistema diferencial ordinário para todos os casos, bem como no controle do

erro global do campo original.

Resultados benchmark são mostrados nas Tabelas (3.1 - 3.4) para Re = 0, 40 e 300,

que são os valores mais típicos encontrados na literatura.

Em posições próximas da entrada , a convergência, para a precisão desejada, é algo

mais lenta que em regiões mais afastadas da entrada. Este é um comportamento esperado e

típico da metodologia empregada. A função \y„ (y ) , usada como filtro, é o perfil

completamente desenvolvido de \\i (x, y); então, ao avançar o escoamento ao longo do

duto, \j/ (x, y) tende para \\i^ (y) , requerendo-se um número cada vez menor de termos

na expansão.

As quatro tabelas mostram claramente que o controle de erro é estabelecido numa

determinada posição, ao se variar o número de termos da série que representa o campo.

29

Outro comportamento relevante é o fato de que ao se aumentar o número de

Reynolds, maior é o número de termos necessários na expansão para convergência numa

precisão desejada. Esta tendência é devida da representação da função-corrente numa

expansão em autofunções, que provém de um problema de autovalor puramente difusivo,

extraído do caso Re —» 0.

Assim, os operadores difusivos vão cedendo sua importância para os termos

convectivos, que têm a característica de termo fonte, ao se incrementar o número de

Reynolds.

Nas Tabelas (3.2) - (3.4) verifica-se os resultados anteriores, reportados na

literatura, encontrados a partir da aplicação dos métodos de diferenças finitas e elementos

finitos. Note-se que ordens de truncamento relativamente baixas, já fornecem resultados

com desvio pouco significativo, em relação aos resultados obtidos por técnicas puramente

numéricas.

TABELA 3.1 - Convergência da velocidade no centro do duto, u (x, 0 ), para Re —» 0.

(condições de entrada: u = 1, v = 0)

N

3

5

7

9

11

13

x = 0.2

1.031

1.058

1.065

1.066

1.067

1.066

x = 0.4

1.193

1.198

1.198

1.198

1.198

1.198

x-0.6

1.319

1.320

1.321

1.321

1.321

1.321

x = 0.8

1.405

1.405

1.406

1.406

1.406

1.406

30

TABELA 3.2 - Convergência e comparação da velocidade no centro do duto, u (x, 0),

para Re = 40 (condições de contorno na entrada: u = 1, v = 0).

N

3

5

7

9

11

13

Ref. [46]

Ref. [49]

x = 0.2

0.977

1.012

1.020

1.022

1.022

1.022

1.0223

1.0243

x = 0.4

1.075

1.083

1.083

1.083

1.083

1.083

1.0849

1.0884

x=0.6

1.162

1.165

1.166

1.166

1.166

1.166

1.1693

1.1737

x=0.8

1.246

1.250

1.251

1.251

1.251

1.251

1.2535

1.2580

TABELA 3.3 - Convergência e comparação da velocidade no centro do duto,u(x,0), para

Re = 300 (condições de contorno na entrada: u = 1, v = 0).

3

5

7

9

11

13

15

17

Ref. [47]

;i..:::;X:^0j2pp:::::;

0.918

0.974

0.998

1.006

1.008

1.007

1.007

1.007

1.008

1.054

1.072

1.071

1.071

1.071

1.071

1.071

1.071

1.075

1.273

1.279

1.280

1.280

1.280

1.280

1.280

1.280

1.283

ww-mmÊMmm1.423

1.425

1.425

1.425

1.425

1.425

1.425

1.425

1.425

31

TABELA 3.4 - Convergência e comparação da velocidade no centro do duto, u(x,0),

para Re = 300 (condições de contorno na entrada: u = 1, co= 0).

N

5

7

9

11

13

15

17

19

Ref. [45]

Ref. [47]

x = 0.20833

1.022

1.041

1.047

1.048

1.049

1.050

1.051

1.052

1.0581

1.050

x = 0.8333

1.134

1.145

1.153

1.159

1.163

1.166

1.168

1.170

1.1880

1.170

x = 3.3333

1.316

1.324

1.328

1.331

1.333

1.335

1.336

1.337

1.3572

1.34

x=7.5

1.437

1.440

1.441

1.442

1.443

1.443

1.444

1.444

1.4509

1.44

Para o estudo do efeito do truncamento do domínio infinito, construiu-se um código

auxiliar, onde impõe-se a condição de contorno da saída em diversas posições truncadas do

duto. Os resultados dos testes realizados são mostrados nas Tabelas (3.5.a,b). Nos dois

casos, é observado que a posição do truncamento afeta as zonas próximas a saída, tendo

menor influência em posições mais próximas da entrada do canal.

Uma inadequada especificação dessa posição não produz problemas de divergência

na solução das equações diferenciais, mas deteriora a solução em zonas próximas à saída.

Por isso é necessário verificar a posição de truncamento em várias posições intermediárias,

a fim de se controlar a convergência dos resultados.

A incerteza da posição adequada do contorno é eliminada ao se considerar a

transformação proposta em (2.32), onde o problema sempre é resolvido num domínio finito

de [0,1].

32

Uma outra vantagem em se usar a transformação do domínio é a redução do número

de pontos na malha computacional, necessários para resolução do sistema diferencial

ordinário.

TABELA 3.5.a,b - Velocidade na linha de centro do duto, u(x,0), para diversas

posições de truncamento do duto.

mWSÊÊÊÊÊÊmí

0.2

0.6

1.0

1.4

1.8

2.0

1.022

1.166

1.323

1.421

1.480

1.5

llÉIIlillil1.022

1.166

1.322

1.417

1.463

1.475

í:i::i:::SSí;;iiiI|;:Sí6i:y;íí:t::::ííÍ

1.022

1.166

1.322

1.417

1.463

1.475

WÊ^&ÈÊSmmM1.022

1.166

1.322

1.417

1.463

1.475

(a) Re = 40

X

0.20833

0.8333

3.3333

7.5

L=8

1.007

1.071

1.280

1.437

L-Í0

1.007

1.071

1.280

1.425

L = 15

1.007

1.071

1.280

1.425

L ==00

1.007

1.071

1.280

1.425

(b) Re = 300

As Tabelas (3.3) e (3.4) mostram o comportamento da convergência para dois

diferentes casos de condições de contorno na entrada. O caso de entrada uniforme e

paralela apresenta uma tendência de convergência mais rápida que o caso de entrada

irrotacional.

33

Nota-se, por exemplo, na posição x = 0.8333, que apenas 7 termos foram

necessários para a convergência com 3 dígitos significativos ao se considerar entrada

uniforme e paralela, enquanto 19 termos foram requeridos ao se assumir entrada

irrotacional.

Comparações entre resultados obtidos a partir das equações de camada limite

[29,50] e Navier-Stokes (para as duas opções de condições de entrada) são mostradas na

Figura (3.2). Demonstra-se que soluções das equações de Navier-Stokes, com números de

Reynolds crescentes, se aproximam das soluções das equações da camada limite. Ressalta-

se, entretanto, as diferenças neste comportamento entre as duas situações de entrada no

canal.

A solução das equações de Navier-Stokes, para o caso de entrada com escoamento

irrotacional tem notadamente melhor concordância com a solução de camada limite,

inclusive para números de Reynolds relativamente baixos. Fica portanto evidente a

necessidade de ser bem caracterizada a condição de entrada real do escoamento a ser

simulado, em particular para baixos números de Reynolds, e sua grande influência na

decisão de adotar-se o modelo simplificado de camada limite, tendo em vista os resultados

acima.

Nas Figuras (3.3) e (3.4) são mostrados resultados típicos do desenvolvimento dos

perfis das componentes de velocidade, u e v.

Nota-se da Figura (3.3), que o perfil de velocidades, para o caso de entrada

irrotacional, se desenvolve mais rapidamente do que no caso de entrada uniforme e

paralela. Da Figura (3.4), observa-se a progressiva atenuação axial e transversal da

componente transversal da velocidade, no caso de entrada uniforme e paralela, a partir da

condição v = 0 na entrada, que volta a assumir essa distribuição no escoamento

completamente desenvolvido.

34

1.50-

1.40-

1.30-

1.20 —

1.10 —

1.00-

U = 1 , W =" 0 ;t> III U * 1 , V * 0

= 300

Re =600

Ro = 1200

O Camada Limite [29],[50]

0.0000\ i rr \ \ r r i 1111 i r i i i 1111 i i i TTTTT]

0.0001 0.0010 0.0100 0.10CX=x/(Dh - Re)

FIGURA 3.2 - Velocidade u (x,0) ao longo do duto para condições de contorno

uniforme e irrotacional na entrada, comparada com resultados de camada limite

u = 1 , v = 0

U = 1 , ( 0 = 0 Re = 3001.00

0.50 —

0 . 0 0 -

-0.50 —

-1.00

0 1 2 3 4 5 6 7 8 9 10X

FIGURA 3.3 - Desenvolvimento do perfil de velocidades u (x,y) ao longo do duto para

condições de contorno uniforme e irrotacional na entrada (Re = 300).

35

1.00-

0.80—I

0 .60 -

0.40 —

0 .20 -

0.00-

Re = 300

a

j .

= 7.5

= 3.33333

= 0.83333

= 208333

i n i i i i i I i i i i i i i i i i i i n i n i n n n i r t

-0.20 -0.15 -0.10v(x,y)

-0.05 0.00

FIGURA 3.4 - Desenvolvimento do perfil de velocidades v (x,y), com condição de

contorno uniforme e paralela na entrada do duto.

36

IV - ESCOAMENTO SOBRE UM DEGRAU

IV.l - DESCRIÇÃO DO PROBLEMA

Neste capítulo estuda-se o desenvolvimento do escoamento laminar entre placas

paralelas que abruptamente ampliam sua separação numa determinada posição (Figura

4.1.a). O objetivo deste clássico problema-teste em mecânica de fluidos computacional

("Backward Facing Step"), é avaliar se as soluções, provenientes do uso de um

determinado método numérico, representam as zonas de recirculação adequadamente e

estudar a influência que as condições de contorno específicas na saída do canal têm sobre

tal solução.

O problema do escoamento sobre um degrau como caso-teste tem sido simplificado

ao excluir-se do domínio de solução o duto de placas paralelas anterior ao degrau, e

assumindo-se que na entrada e na saída do segundo canal o escoamento é completamente

desenvolvido (Figura 4.1 .b).

Assim, o problema pode ser incluído dentro da classe de dutos regulares até aqui

analisados. Esta definição do problema do degrau, permite que se façam comparações de

soluções obtidas pela aplicação de diferentes métodos numéricos, bem como com

resultados experimentais. As condições de contorno na entrada podem então ser

especificadas como:

O ; - l < y < 0u(O,y) = (4.1.a, b)

[6(y-yz) ; 0 < y < l

v(0 ,y) = 0 ; - 1 < y < 1 (4.1.c)

u = 6 y - 6 y 2

u = 0v = 0

u = f,(y)v = 0

u = 0v = 0

37

(a)

u=0 v=0

u=0 v=0

(b)

u=0 v=0

u = 0 v=0

(c)

4'

u = -d-y2)

v = 0

FIGURA 4.1 - Esquema do problema de escoamento sobre um degrau ("Backward Facing

Step").

38

Considerar que na entrada do segundo canal o escoamento é completamente

desenvolvido é uma hipótese válida para comparações numéricas entre diferentes métodos.

Mas, nas situações reais, o perfil de velocidades no início do degrau não necessariamente

encontra-se completamente desenvolvido. Assim, uma alternativa viável para comparações

de resultados numéricos com medições experimentais é utilizar os valores medidos na

entrada do duto, interpolá-los adequadamente e adotar-se este perfil como condição de

contorno na simulação.

O estudo experimental que serve de referência às comparações que se seguem, é

devido a Denham e Patrick [51]. Neste caso os pontos medidos na entrada são interpolados

via cubic-splines. Considera-se também que o escoamento é paralelo, por não dispor-se de

nenhuma informação sobre a componente transversal da velocidade, v (Figura 4.1.c).

Portanto, as condições de contorno na entrada são dadas por:

f,(y)

v(0,y) = 0 ; - 1 < y < 1 (4.2.c)

onde f( (y) representa a curva que interpola os valores da componente u medidos na

entrada.

Sendo assim estabelecidas as condições na entrada do canal e sendo o escoamento

completamente desenvolvido na saída, é simples reescrever as condições de contorno

transformadas na entrada, segundo (2.14.f,g) e (2.36.c,d).

Para o perfil teórico dado pela eq. (4.1):

39

r f° rX i ( 0 ) = J Y,<KO,y)dy = J Y^CO,y)dy + J ?

- 1 - I O(4.3.a)

onde

4>(0,y) = i 3 1- - y - ~ ; 0< y < 1

4 ; 2

(4.3.b-c)

dX,(0)

dT)= 0 ; i = l,2,3,...,N (4.4)

e o coeficiente de balanço torna-se

q = (4.5)

No caso da comparação com os dados experimentais, o processo de transformação é

similar, expressando-se a integral acima nos dois subdomínios.

Assim :

= J Yi(j>(O,1

X,(0)= J Yi(j>(O,y)dy +-1 a

(4.6.a)

r (4.6.b)

- ^ - - q ; a < y < l

40

j (0)= 0 ; i = l,2,3,...,N (4.7)d 0 ; i l,2,3,...,N

onde o coeficiente q édeterminado através do balanço global de massa

Ah

IV.2 - REVISÃO BIBLIOGRÁFICA

O problema de escoamento sobre um degrau tem sido tratado com razoável

frequência na literatura tanto para regime laminar como para turbulento, sendo na mecânica

dos fluidos computacional um problema-teste que também motiva o interesse da pesquisa

em termos físicos.

Denham e Patrick [51], estudaram este problema experimentalmente. Mediram a

componente longitudinal, u, em diversas posições ao longo do duto, utilizando

anemometria laser com sensibilidade direcional. Conseguiram medições para Re = 438 até

1374 (de acordo à definição de Re no presente trabalho). Não reportaram vórtices

localizados na parede superior do duto, mas encontraram recirculações transversais na

região do vórtice primário. Compararam seus resultados com avaliações numéricas,

obtendo alguma concordância, mas concluíram que o comprimento da recirculação

estimado numericamente era maior que na situação experimental.

Campion-Renson e Crochet [52] utilizaram o método de elementos finitos com uma

formulação vorticidade-função corrente. Usaram elementos retangulares quadráticos e

obtiveram resultados até Re = 50.

41

Leschziner [53], utilizou três esquemas de diferenças finitas para o cálculo de

recirculações dentro de escoamentos em estado estacionário, considerando malhas grossas

e finas. Seus resultados, comparados com os de Denham e Patrick [51], mostraram

coerência, encontrando-se discrepância similar entre resultados teóricos e experimentais.

Armaly, Durst, Pereira e Schõnung [54], fizeram um estudo experimental amplo,

para número de Reynolds entre 140 e 16000, efetuaram comparações com resultados

teóricos. Usando anemometria laser, encontraram que o escoamento sobre degrau com

razão de aspecto 0.5 possuía uma estrutura laminar bidimensional para números de

Reynolds inferiores a 800. Acima deste valor o escoamento permanecia laminar, mas com

características tridimensionais, até Re = 2400. A partir deste valor de Re, o escoamento se

encontrava numa região de transição, tornando-se totalmente turbulento-bidimensional

acima de Re = 13200. Na faixa de números de Reynolds estudada, encontraram até três

tipos de vórtice. O vórtice primário permanecia para todos os valores de Reynolds, mas

variava o comprimento da recirculação, atingindo um máximo no limite do regime laminar

para transição, e mantendo-se quase constante no regime turbulento. O vórtice secundário

aparecia quando o escoamento era tipicamente tridimensional, desaparecendo quando o

escoamento atingia o regime turbulento. Já o vórtice terciário aparecia no início da

transição, desaparecendo em Re = 4600. Seus resultados concordaram com cálculos

teóricos, obtidos a partir do código TEACH [55], até Re = 800, concluindo-se que efeitos

tridimensionais, que foram detectados experimentalmente, aumentaram as discrepâncias

entre resultados teóricos e experimentais.

Thamgam e Knight[56], usaram uma formulação das equações de Navier-Stokes em

coordenadas generalizadas e em variáveis primitivas e utilizaram o método de volumes

finitos para sua solução. Fizeram um estudo do comprimento da recirculação ao se variar

a altura do degrau, para diferentes valores do número de Reynolds, desenvolvendo uma

correlação para predizer o comprimento do vórtice primário.

42

Gartling [57], motivado pelos Minisymposium sobre "Outflow Boundary

Conditions" organizado pela University College of Swansea em 1989 e 1991, quando se

definiram problemas-teste para se avaliar diferentes condições de contorno de saída,

resolveu o problema do degrau para Re = 1600, fornecendo resultados benchmark para

futuras comparações numéricas. Utilizou o método de elementos finitos, discretizando o

domínio com diferentes números e tamanhos de elementos do tipo biquadrático. A solução

numérica foi feita usando o código NACHOS II [58], o qual é baseado numa formulação de

Galerkin. A verificação dos resultados foi feita através do programa FIDAP [59], o qual é

também baseado no método de elementos finitos. Em todos os casos foi empregada a

formulação em variáveis primitivas das equações de Navier-Stokes.

Rogers e Kwak [60], discretizaram as equações de Navier-Stokes num esquema

upwind e usaram a formulação de compressibilidade artificial. Seus resultados tiveram boa

concordância com os de Armaly et ai. [54], mas distanciaram-se para altos números de

Reynolds.

Srinivasan e Rubin [61], recentemente, usaram um procedimento multimalha

adaptativo com decomposição do domínio, mostrando seus resultados grande concordância

com Gartling [57] na componente "u" e alguma discrepância na componente "v".

Truncaram o domínio infinito em diversas posições, mostrando que as soluções tinham

concordância concluindo assim que o efeito de truncamento do canal era de pouca

relevância.

Onur e Baydar [62], fizeram um estudo teórico-experimental do escoamento sobre

um obstáculo quadrado. Para efeitos de validação da análise teórica resolveram o problema

apresentado por Armaly et ai. [54]. Usaram para tal um esquema de diferenças finitas

similar ao do código TEACH [55], considerando as recomendações de Van Doormaal e

Raithby [62] para incrementar a taxa de convergência. Assim obtiveram resultados

numéricos com melhor concordância com os experimentais do que aqueles obtidos por

Armaly et ai. [54].

43

IV.3 - RESULTADOS

Similarmente ao caso de placas paralelas, o erro relativo foi estabelecido em IO"4,

tanto na solução do sistema diferencial ordinário como no esquema de controle global do

erro no campo original. Os casos testados foram computados no VAX.8810 e no CRAY

EL - 94 da COPPE/UFRJ.

Como primeiro caso teste, foi simulado o experimento de Denham e Patrick [51].

Foram escolhidos para o presente estudo,o menor e maior número de Reynolds empregados

naquele estudo (Re = 438 e Re = 1374), em que foram estabelecidas as condições de

entrada.

O comportamento da convergência para o campo de velocidades foi estudado nas

mesmas posições ao longo do duto em que foram feitas as medições experimentais. As

tabelas (4.1) e (4.2) mostram detalhadamente os resultados para ordens crescentes de

truncamento das séries, até atingir-se a tolerância solicitada. Para Re = 438, foram

necessários N = 20 termos na expansão da função-corrente, mas apenas 8 termos são

requeridos para uma convergência em escala gráfica. Para o Re = 1374, foram necessários

mais termos na expansão (N = 22) a fim de satisfazer a tolerância, ao passo que em escala

gráfica 10 termos foram suficientes.

44

TABELA 4.1 - Convergência da componente longitudinal do campo da velocidade

(Re - 438).

::!H§!

4

8

12

16

20

-2.69E-2

-5.13E-2

-5.23E-2

-5.22E-2

-5.23E-2

:. •::x==^333|i|:

-7.76E-2

-8.88E-2

-9.30E-2

-9.27E-2

-9.28E-2

-4.64E-2

-5.09E-2

-5.01E-2

-5.00E-2

-5.00E-2

4.64E-2

5.02E-2

5.30E-2

5.29E-2

5.28E-2

;;r|Iiiiii3.3|i:i|::

0.120

0.125

0.126

0.126

0.127

4

8

12

16

20

1.188

1.221

1.204

1.202

1.203

1.220

1.228

1.221

1.220

1.220

1.226

1.248

1.247

1.246

1.246

1.250

1.226

1.225

1.225

1.225

1.201

1.181

1.180

1.180

1.180

u(x,0.S)

4

8

12

16

20

1.339

1.374

1.362

1.361

1.362

1.263

1.291

1.284

1.284

1.284

1.043

1.069

1.065

1.066

1.066

0.839

0.861

0.858

0.859

0.859

0.736

0.756

0.754

0.755

0.755

45

TABELA 4.2 - Convergência da componente longitudinal do campo da velocidade

(Re =1374).

6

10

14

18

22

6

10

14

18

22

f|ÍÍ:;p.|5»l

-0.103

-6.91 E-2

-7.57 E-2

-7.48 E-2

-7.53 E-2

1.054

1.040

1.018

1.020

1.020

llillill

-0.123

-0.111

-0.113

-0.113

-0.113

-0.144

-0.157

-0.160

-0.161

-0.161

1.046

1.030

1.019

1.020

1.019

1.051

1.032

1.031

1.031

1.031

iiiiiiiOJQOj:!

-0.154

-0.178

-0.181

-0.181

-0.181

lilMíii1.049

1.031

1.032

1.032

1.032

-0.127

-0.126

-0.121

-0.119

-0.120

1.021

1.003

1.001

1.000

1.000

u (x, 0.5)

6

10

14

18

22

1.088

1.067

1.082

1.083

1.081

1.076

1.055

1.060

1.060

1.060

1.013

0.991

0.989

0.988

0.988

0.902

0.860

0.849

0.847

0.847

0.732

0.662

0.642

0.639

0.640

46

Re = 438

(Re = 73 [51])

GITT

O Experimental [51]

90 100 110 130 140 15

1 I

130 140 15

x(mm)

FIGURA 4.2 - Desenvolvimento da componente longitudinal da velocidade ao longo

do canal (Re = 438).

Re = 1374

(Re = 229 [51])

GITT

O Experimental [51]

x(mm)

FIGURA 4.3 - Desenvolvimento da componente longitudinal da velocidade ao longo

do canal (Re = 1374).

47

i.o

Re = 438

(Re = 73 [51])

FIGURA 4.4 - Isolinhas da Função Corrente (Re = 438).

REYNOLDS= 1374

(Re = 229 [51])

FIGURA 4.5 - Isolinhas da Função Corrente (Re = 1374).

48

As figuras (4.2) e (4.3) mostram, comparativamente, o perfil da componente

longitudinal da velocidade em desenvolvimento calculado via GITT e os respectivos

valores experimentais, evidenciando-se que existe uma excelente concordância entre eles.

Para Re = 438, a concordância é mais notória que para Re = 1374, sendo este fato

explicável a partir da condição de contorno assumida na entrada do canal. Denham e

Patrick [ 51 ], mostraram que em Re = 438, a velocidade longitudinal medida no degrau

difere muito pouco do perfil parabólico (escoamento completamente desenvolvido),

demonstrando-se portanto que a componente transversal de velocidade é quase nula, fato

que não ocorre no caso de Re = 1374, onde existe uma diferença mais perceptível entre o

perfil medido e o perfil parabólico, traduzindo-se isto numa componente transversal da

velocidade ainda significativa. Devido a não possuir-se informação da componente

transversal da velocidade, foi assumido como condição de contorno apenas a componente

longitudinal, refletindo-se isto em uma ligeira discrepância entre os resultados teóricos e

experimentais nas vizinhanças da entrada do canal.

Nas figuras (4.4) e (4.5), se apresentam isolinhas da função corrente para ambos os

valores do número de Reynolds. Claramente são mostradas as zonas de recirculação que se

posicionam justamente atrás do degrau. Para o maior número de Reynolds, o comprimento

do vórtice localizado é maior, apresentando além disso o início de uma recirculação

secundária na parede superior.

O segundo passo no presente estudo foi uma covalidação com os resultados

benchmark de Gartling [ 57 ], obtidos via elementos finitos. Para tal efeito, se analisou a

convergência da componente longitudinal da velocidade em duas posições ao longo do

duto, x = 1 4 e x = 3 0 . O processo de convergência até N = 20 termos é mostrado nas

tabelas (4.3) e (4.4), obtendo-se uma excelente concordância com os resultados obtidos por

Gartling, covalidando-se assim os resultados.

49

TABELA 4.3 - Convergência da componente longitudinal do campo de velocidadeemx=14(Re=1600)

-1.0

-0.9

-0.8

-0.7

-0.6

-0.5

-0.4

-0.3

-0.2

-0.1

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

mmm4m0

0.187

0.386

0.591

0.787

0.955

1.076

1.134

1.123

1.045

0.912

0.743

0.559

0.379

0.219

0.089

-0.004

-0.058

-0.075

-0.057

0

0

0.235

0.437

0.621

0.797

0.954

1.068

1.120

1.103

1.020

0.882

0.705

0.516

0.343

0.200

0.091

0.014

-0.034

-0.030

-0.038

0

11111111!0

0.234

0.430

0.615

0.795

0.950

1.063

1.119

1.105

1.023

0.884

0.708

0.520

0.348

0.203

0.092

0.014

-0.032

-0.049

-0.038

0

lillliiill0

0.233

0.429

0.614

0.793

0.949

1.063

1.118

1.105

1.023

0.885

0.709

0.522

0.348

0.203

0.092

0.015

-0.032

-0.049

-0.038

0

iililiiiil0

0.233

0.428

0.614

0.793

0.949

1.063

1.118

1.105

1.024

0.885

0.709

0.522

0.349

0.204

0.093

0.015

-0.032

-0.049

-0.038

0

;!MIÍÍI!:I0

0.232

0.428

0.613

0.792

0.948

1.062

1.118

1.105

1.024

0.885

0.709

0.522

0.349

0.204

0.092

0.015

-0.032

-0.049

-0.038

0

50

TABELA 4.4 - Convergência da componente longitudinal do campo de velocidade

emx=30(Re=1600).

-1.0

-0.9

-0.8

-0.7

-0.6

-0.5

-0.4

-0.3

-0.2

-0.1

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

0

0.105

0.213

0.322

0.432

0.541

0.643

0.733

0.803

0.846

0.857

0.836

0.784

0.709

0.618

0.517

0.412

0.307

0.204

0.101

0

0

0.109

0.218

0.328

0.439

0.547

0.650

0.738

0.805

0.845

0.854

0.831

0.779

0.704

0.612

0.512

0.407

0.304

0.202

0.101

0

0

0.109

0.218

0.328

0.438

0.547

0.649

0.737

0.805

0.844

0.853

0.831

0.779

0.704

0.613

0.512

0.408

0.304

0.202

0.101

0

EÊKmm0

0.109

0.218

0.328

0.438

0.547

0.649

0.737

0.805

0.844

0.853

0.831

0.779

0.704

0.613

0.512

0.408

0.304

0.202

0.101

0

0

0.109

0.218

0.327

0.438

0.547

0.649

0.737

0.805

0.845

0.853

0.831

0.779

0.704

0.613

0.512

0.408

0.304

0.202

0.101

0

0

0.109

0.218

0.328

0.438

0.547

0.649

0.737

0.804

0.844

0.853

0.831

0.779

0.704

0.613

0.512

0.408

0.304

0.202

0.101

0

51

Um outro aspecto estudado foi o esclarecimento das discordâncias entre os

resultados obtidos pelos recentes trabalhos de Gartling [57] Srinivasan e Rubin [61], uma

vez que esses estudos concordaram muito bem para a componente longitudinal da

velocidade, não ocorrendo o mesmo para a componente transversal. Na fig. (4.6)

apresenta-se resultados comparativos para x = 1 4 e x = 30. Claramente é observada uma

excelente concordância entre Gartling e o presente trabalho, já os resultados de Srinivasan

e Rubin apresentam diferenças. Isto deve-se, possivelmente, a uma convergência

insuficiente nos resultados da referência [61], não se constituindo, portanto, em uma

solução benchmark.

As isolinhas da função corrente para este exemplo são também apresentadas na

fig. (4.8). Os vórtices primário (na parede inferior) e secundário ( na parede superior) são

plenamente descritos.

Um terceiro estudo foi feito com base no clássico trabalho de Armaly et ai. [54], no

qual se levanta a questão da validade da formulação bidimensional para valores crescentes

do número de Reynolds. As comparações feitas naquele trabalho, entre os valores medidos

experimentalmente e os simulados através do código TEACH [55] mostram excelente

concordância até Re = 400 (definição de Reynolds de [54], Re = VD/v , V = 2/3 da

velocidade máxima na entrada, D é o diâmetro hidráulico na entrada do canal, equivalente

a duas vezes sua altura, v = viscosidade cinemática), mas mostram-se marcantemente

discordantes para números de Reynolds maiores, como no caso de Re = 1000, apresentado

em [54].

Inicialmente, mostra-se a convergência da componente longitudinal da velocidade

ao longo do duto, requerendo-se 22 termos nas expansões, na pior situação, tal como

mostrado na tabela (4.5). A seguir, compara-se as velocidades assim obtidas pela GITT

com os resultados experimentais da referência [54](fig. 4.7). Uma boa concordância é

alcançada na maioria das posições do duto mostradas, discordando para um segmento em

uma região intermediária ao longo do comprimento do duto. Estas discrepâncias são

52

justificadas pela influência das paredes laterais do canal, as quais não são consideradas

numa formulação bidimensional, e que é mais significativa ao incrementar-se o número de

Reynolds. O efeito tridimensional é marcante, essencialmente nessa região intermediária

do duto, onde as camadas limites das paredes laterais interferem mais fortemente sobre o

escoamento principal que está se desenvolvendo.Sob este aspecto, os resultados obtidos

pela GITT demonstram esta coerência física, a qual não é reproduzida na simulação

numérica feita por Armaly et ai. [ 54 ].

A reprodução das intensas zonas de recirculação primária e secundária, são

apreciadas na fig. (4.9) para este exemplo.

53

TABELA 4.5 - Convergência da componente longitudinal do campo de velocidade

(Re=1942).

S l l l l

4

8

12

16

18

20

22

EÊÍMÈÊÊÊKÍ

-5.98 E-2

-0.105

-0.105

-0.106

-0.106

-0.106

-0.106

ÉIIiiilÉi

-7.28 E-2

-8.78 E-2

-8.70 E-2

-8.77 E-2

-8.79 E-2

-8.78 E-2

-8.78 E-2

liiliiili:liiilliii

0.514

0.484

0.475

0.469

0.466

0.468

0.467

Ililiili:

1.062

1.032

1.030

1.030

1.030

1.030

1.030

IlIIIilI

0.856

0.865

0.862

0.863

0.863

0.863

0.863

liiiiii

0.557

0.638

0.639

0.639

0.639

0.639

0.639

u (x, 0)

4

8

12

16

18

20

22

4

8

12

16

18

20

22

0.444

0.494

0.505

0.505

0.504

0.504

0.504

1.146

1.428

1.429

1.430

1.431

1.431

1.431

0.675

0.757

0.757

0.755

0.754

0.755

0.755

1.320

1.280

1.280

1.282

1.282

1.282

1.282

1.219

1.214

1.217

1.218

1.219

1.218

1.219

!§!l!!l!lf0.376

0.404

0.412

0.419

0.423

0.421

0.421

0.836

0.822

0.824

0.824

0.825

0.824

0.825: , ; • : • • • : . • • : • : • : • ; • : • : - : • ; • : • : • : • : • > : • : • : • ; • : • : • ; • : • : • :

•:-.;;: : •::::-:-:-.;-:-;:x-x-:-:"xo:-:^x':X-

2.88 E-2

5.02 E-2

5.17 E-2

5.12 E-2

5.13 E-2

5.13 E-2

5.13 E-2

0.953

0.937

0.938

0.937

0.937

0.937

0.937

lillllllllll

0.182

0.176

0.177

0.176

0.176

0.176

0.176

1.001

0.980

0.980

0.980

0.980

0.980

0.980

0.439

0.368

0.368

0.367

0.367

0.367

0.367

54

1.00

0.50 —

y 0.00 —

-0.50 —

Re = 1600

GirrO BementosRnitos[57]

0 Elementos Rnitos[ 57]

rj Diferenças Rnitas [ 61 ]

" ' I I I I I I I I I I I I I I I I I I I I I I I I I I I I I T I I I I I I I I I

-0.03 -0.02 -0.01 0.00 0.0v(x,y)

FIGURA 4.6 - Comparação dos Perfis da componente transversal da velocidade ao

longo do duto (x =14 e x = 30 ; Re = 1600).

Re «1942

( Re=1000 [16])

GITT

O Experimental (16]

U (cm/s)

180 180 180 180 180 180 180 180 180 180 180 180 180 180 180

0.00 2.55 4.00 6.00 8.00 10.00 12.00 14.00 15.53 16.12 20.00 22.45 23.88 30.20 49.59x

FIGURA 4.7 - Desenvolvimento da velocidade longitudinal ao longo do duto(Re =1942).

55

•t.o

Re= 1600

(Re = 800 [54])

20 25

FIGURA 4.8 - Isolinhas da Função Corrente (Re = 1600).

4.0

3.0

2.0

1-0,

Re= 1942

(Re = 294 [55])

FIGURA 4.9 - Isolinhas da Função Corrente (Re = 1942).

56

V - ESCOAMENTO EM CANAIS DE GEOMETRIA ARBITRÁRIA

V.l - DEFINIÇÃO DO PROBLHMA

Considera-se agora um duto cujas paredes são definidas por uma geometria

arbitrária, identificada pelas funções que descrevem o contorno, y{ (x) e y2 (x), dentro do

qual escoa um fluido newtoniano em regime laminar permanente. Para maior generalidade

da formulação do problema, assume-se que não existe simetria com relação à linha de

centro do duto como descrito na figura 5.1.

O escoamento é governado pelas equações de continuidade e de Navier-Stokes, que

em sua forma adimensional e em termos de função-corrente somente, são dadas por:

onde o número de Reynolds ( Re ) é definido a partir da semidistância entre as paredes e a

velocidade média, ambas na entrada do duto. V4 é o operador biharmônico definido

anteriormente na equação (2.5).

As hipóteses consideradas são:

- escoamento bidimensional

- escoamento incompressível

- propriedades do fluido constantes

- impermeabilidade e não deslizamento das paredes do duto

- fluido newtoniano

As condições de contorno nas paredes do duto são especificadas ao considerar-se

impermeabilidade e não deslizamento, assim:

58

duto o campo de velocidades é dado pelas componentes nx (y) e Vj (y), que em termos de

função corrente são dados por:

— ; - v , ( y ) = -£- (5.2.a,b)

Integrando-se ambas as equações:

v | / ( O , y ) = k 1 + j u , ( y ) d y ( 5 . 3 . c )

- J v lv | / ( 0 , y ) = k 2 ( y ) - J V l ( y ) d x (5.3.d)o

Dado que estamos interessados em determinar v|/ (0,y), a integral na equação (5.3.d)

é nula e a relação resultante apenas nos informa que a função corrente na entrada depende

de "y". Por outro lado, a equação (5.3.c) descreve o valor da função corrente na entrada do

duto em função de Uj (y).

A partir de (5.3.c) é possível determinar \\i (0,y2):

+ J u , ( y ) d y (5.4)

A integral em (5.4) define a vazão por unidade de comprimento (Q), na entrada do

duto, ou seja:

v}/(0,y2) = k 2 = Q + k, (5.5)

59

V.2 - SOLUÇÃO VIA GITT

Para facilitar a solução do sistema (5.1) fazendo-se uso de transformações integrais,

é conveniente definir uma função "filtro" que permita homogeneizar as condições de

contorno na direção "y". Isto permitirá posteriormente escolher o problema de autovalor

em tal direção. Assim:

\)/(x,y) = \|/ (x,y) + F(x,y) (5.6)

onde y* (x,y) é o potencial desconhecido a ser determinado, e F (x, y) é o filtro, o qual

tem os mesmos valores que i|/(x,y) nas paredes do duto. A função F (x, y) não é

necessariamente uma solução particular de y (x,y).

Logo, o sistema (5.1) pode ser reescrito, considerando-se (5.2), como:

õy õ y õy õ y õy õ y õy õ y

õy Õx3 + õy Õxõy2~ õx Õx2õy~ õx Õy3

a. * ^3-n a. * ÍS3T-' a. * a3 T̂ 5. • 3 3 rõy o r õy õ r õy o r õy o r1 1

õy õx3 õy õxõy2 õx õx2õy õx õy3

'Yn ^3 * a-n o3 * OT^ a3 * ai^ -̂ 3 *

õr õ y or o 11/ ar õ y õr õ yõy õx3 õy õxõy õx ôx2õy õx õy

õy õx3 õy Õxõy2 õx Õx2 õy õx õy3

Õ4F Ô4F2

1 í õ\ Õy Õy ÕF

R e j õx4 +2Õx2ôy2+ Õy4 + õx4+2ôx2õy2 + Õy4}

(5.7.a)

A filtragem das condições de contorno na direção y resulta em:

60

V|/*(x,-y,) = k 1 - F ( x , - y 1 ) (5.7.b)

õy * (x-yl )= 0 (5.7.c)

vj ;*(x,y 2) = k 2 - F ( x , y 2 ) (5.7.d)

= 0 (5.7.e)

A filtragem das condições de contorno na entrada e saída do duto serão analisadas

no capítulo seguinte para um problema específico.

Nota-se que os termos convectivos conservam um núcleo não linear, com a mesma

estrutura da equação original, além de alguns termos com uma estrutura linear devido ao

"filtro". Os termos difusivos resultam em estruturas semelhantes tanto para o filtro F(x,y),

como para a nova variável incógnita \|/ * ( x, y ) .

V.3 - O FILTRO F (x. vi

O filtro F (x, y) deve ser tal que reproduza o valor da função corrente nas paredes

do duto, ao longo do seu comprimento. Esta função pode ser construída ao se considerar

que em cada posição ao longo do duto se tem um perfil de velocidades desenvolvido, o

qual vai se adaptando à irregularidade do canal.

61

FIGURA 5.2 - Geometria que identifica o filtro F (x, y)

A relação entre o sistema de coordenadas original (x, y) e o sistema transformado

(T|, X) é dado, por:

(5.8,9)

onde y3 representa a defasagem entre eixos "y" e "T|" . Assim, em analogia à função

corrente do escoamento desenvolvido, tal como mostrado no anexo A, é possível descrever

a função F como:

F(X,T1) = -ri I 1 TI

(5.10)

ou em termos das coordenadas originais:

62

y-y3 (5.11)

onde Q é a vazão entre A - A'.

A dependência em x do filtro F (x, y) recai sobre y3 (x) e y0 (x) os quais definem as

supostas "placas paralelas" que passam pelas paredes do duto.

Uma verificação do valor de F (x, y) nas paredes fornece:

F ( x , y 2 ) = (5.12.a)

F ( x , - y , ) = (5.12.b)

Fica claro que y pertence ao intervalo [ - y i ( x ) , y 2 ( x ) ] , e r| a

[ - y o ( x ) , y o ( x ) ] , mas é útil definir uma nova variável, £, , a qual permitirá determinar

mais facilmente os coeficientes que aparecem da transformação integral, tal como mostrado

posteriormente.

Assim, define-se:

ri y - y 3 (5.13)

onde o domínio de % e [ —1, 1 ] .

O filtro F (x, y) pode ser emtão representado como:

3Q

4(5.14)

63

V.4 -PROBLEMA PR AUTOVALOR

Dadas as condições homogêneas (5.7.b-c) escolhe-se o problema auxiliar na direção

y. Deve ficar claro que a escolha do problema de autovalor na direção x não é apropriado

pela presença de condições de contorno não homogêneas.

Analogamente ao caso de dutos regulares o problema auxiliar de 4â ordem é dado

por:

, -yi(x)<y<y2(x) (5.15.a)n Í Y , ( x , y

com condições de contorno:

on= 0 (5.15.b,c)

Y i ( x , y 2 ) = 0 ; ^ = 0 (5.15.d,e)on

Considerando a relação (5.13), o sistema (5.15) pode ser reescrito em termos de t,

como:

*4 " ' - (5.16.a)

dY-, (-1)' ^ = 0 (5.16.b,c)

(5.16.d,e)

64

A solução de (5.16) é dada, similarmente a (2.15.a), por:

cosh(ni

cos()j.* ) cosh(fi*

sinh(^j Ç)

sin(u.-) sinh((x-)

i = 1,3,5,.

i = 2,4,6,.

(5.17.a,b)

onde o autovalor |_tj , é encontrado a partir de:

Nas coordenadas originais:

(5.18)

Y i = i

cos[Hj(y-y3 j (y-y3

sinh[Hj ( y - y

i = 1,3,5,.

i = 2,4,6,.

(5.19.a,b)

A norma é calculada a partir de sua definição:

Ns = 1 Yf dy i = 1,2,3,.. (5.20)

Fazendo-se a transformação de variável indicada por (5.13), obtém-se:

i = y 0(5.21)

65

Logo:

N j ( x ) = 2 y o ( x ) (5.22)

Fica também evidente que a norma é independente da ordem i, de forma que daqui

por diante denominaremos a norma apenas por N.

V.4.1 - PROPRIEDADE DE ORTOGONALIDADE

Para determinar a propriedade ortogonalidade de que gozam as autofunções

operamos com J Yj dy ambos os membros da equação (5.15.a):- y i

f yi ô4Y- fJ a 4 J •> ~ i j Y.Yjdy (5.23)

-y,

O lado esquerdo de (5.23) pode ser integrado por partes e, utilizando as condições

de contorno (5.15.b-e), obtém-se:

3 4 Y : f y i d 3 Y . ô Y i Çy> Õ2Y-õ2Y-Yd i i^ i7 d - L ^wL ̂ w (5-24)

Logo:

( 5 - 2 5 )

Intercambiando os sub-índices, i -> j , j -> i :

66

y2 d4 Y : f " f »» a 2 Y-f' Y < d y = JFl

Subtraindo entre si as equações (5.25) e (5.26):

- . J ) íV2

í

ou seja:

~yi 2 y 0 , para i = j

Y,Y j dy = 0 (5.27)

Dado que \x f * jLt j , a integral é nula:

Y iY jdy = 0 ,parai?ij (5.28)

J Y iY jdy = i (5.29.a,b)

V.5 - PAR TRANSFORMADA-INVERSA

Representando a função desconhecida \\i* (x,y) como uma expansão em torno da

autorunção Y^

(5.30)

Operando-se ambos os lados da equação (5.30) com o operador J Yj dy:

67

Çyi Çyi

J Y j M / * d y = J Yj

(

i=i

dy (5.31)

Reordenando o segundo membro da equação (5.31), dado que a integração é em "y"

e o domínio de ijT; é em "x" :

fy2 « j-

J YjM/*dy = 2 J V j J Yi Yj dy (5.32)i=l - y i

Fica evidente reconhecer, pela propriedade de ortogonalidade, que o somatório do

lado direito da equação (5.32) tem apenas um único termo não nulo, para o caso i = j .

Obtém-se desta forma a regra de transformação:

1 fyi= ^7T J Y«<x'y> ̂ W

^\X) -v,(5.33)

V.6 - TRANSFORMAÇÃO INTEGRAL DO PROBLEMA

Tendo estabelecido o par de transformação, transformada e inversa, equações (5.30,

5.33), e havendo sido definidas a autofunção e as propriedades de ortogonalidade do

problema de autovalor escolhido, podemos agora proseguir com a transformação integral

do problema diferencial original. Para tal, opera-se com JYj dy ambos os lados da- yi

equação (5.7). Mostramos a seguir a transformação de cada um dos termos de equação

(5.7), a fim de evidenciar a aplicação da GITT. Em todos os casos, v|/* (x ,y ) foi

representada pela expansão (5.30); as suas derivadas, tanto em x como em y, encontram-se

no apêndice B.

68

fJ

y2 * J r—dy = I1 v

Y,- y i

a2Y k

j=l k=l

j õ 3Yk

1 dydy+

__, f3W :\lf k I

J w

" S Y J-dy +

fT J T K J dy

(5.34.a)

J'- y i

^ a 3 Yk _ a2 Yk _ ,dy =

j=l k=l - y i

1 ay axay:- - P2v

+ VjVk J Y,- y iay ay: -dy

(5.34.b)

69

í r^=í Y=k _

õ2Yk _ ,

*y = 2.2.j=i k=i

- - f y2^7"

1

k

õx2õy Y

VjVk J Y:1 ax Õxôy-dy + - - • f

VjM>k J- y i

i dx

2xj7 jVk J52Yk

Y.Y.-r^-dy- y i

+ ¥jVk J YiYi-rf-dy- y i

"jVk" J• yi

(5.34.C)

í Y:1 ax ay3- fd y = l Y,L k=l

dy =

+VjVk J Y-Y- ^ d y

(5.34.d)

í dy 5 J F1 3y õx 3 -í Y: dy =

-;íy2

ôy

(5.34.e)

70

f -Ay2

Y,t? ay YJ dxõy

dy =

,-i Y;dxõy

•dy

(5.34.Í)

Í- y i 5x

= jy2

- y ij ^ j \õK2õy

U í dY=Y; -dy

- y i

53FYi YJ TTT-dy

yi õx õy(5.34.g)

J Y,6 3 F

dy =

,-A + W J(5.34.h)

71

Ô2Y;j _ _ .

3VJT;32Yj aF

Y, -r^—dy-- y i - yi

1 Ôx

*•; í " V, Yj f dy

(5.3.Í)

í-yi

^"i Y : - — - ^ ^ d y + vjT:" Y:J J_ ' Sxôy dy J J_

(5.34.J)

ôy õxdy

(5.34.k)

í * °o

5y: J JY: (5.34.1)

72

FC, = íY,ÕF Ô 3 F ÕF 5 3 F

i

ox3 3y ôxdy2 dx õx28ydy (5.34.m)

í»

- y i

a4Yj r;J y2

- y i - y i OX'dy

ax-- f

J wY i17d y

Observa-se que nesta última expressão utilizou-se a relação:

(5.34.n)

í- y i

,-•• Í- y i

(5.34.O)

í 1 a,. 4 •íd y = J Ys- y i

dy=nf Nil/, (5.34.p)

FD > J Ys

a4F a 4 F2 a., 2. a x ax ay

dy (5.34.q)

a4F

~õy~'= 0 (5.35)

Fatorando os termos comuns se chega ao seguinte sistema diferencial ordinário:

73

s +\|/jV|/kF10 i jk

j=l k=l

ijk +MTjVk'F12ijk j + X i ^ J 1 7 0 ' ] +vl rj'F1iJ +HTj"F2iJ +vjrj'"F3ij

j"E2 i j +vjTj"E3fl ] f »if NijT, +Nijriiv +FD;

(5.36)

Deixando em evidência o diferencial de 4a ordem e reordenando convenientemente,

tem-se:

j=i k=i

ijk +\jrjVkFllijk +^j NHlii +

(5.37)

onde os coeficientes, que dependem de cada posição em x, são calculados a partir de:

74

F00 : ;

f fJ

» j a 3 Y kd y -

-í- y i - y i

1 ÕX

FO1f

ük = 3 J- y i

• d y -

(5.38.a)

Í2 Y:- y i

1 ox 5xõy-dy

(5.38.b)

F02 i ; t =5 Y J

- y i

j ox ay(5.38.c)

FO3 Ük - J- y i

(5.38.d)

=̂-1 \ V1 53Yk

F10i!lr = - I Y : Y ; — r ^ d y - J " "dy-J Y^i^^dyJ(5.38.e)

F l lrvi Õ 2 Y | . y 2

dy~ J Y.Y ^^ (5.38.Í)

F12 ijk - - JaYk (5.38.g)

75

F0 ; i = íõxõy :

-dy-

fJ Y;

Í Y:-y, ' ax 3

r» a3Yj OF:dy+J Y i T — r —dy-

Í- y i õy

(5.38.h)

Fl 0 =~ J Í 53F

5y;Y i Y j — d y

Í

2 Y;

Í- y i

- y iôxõy dx.

dy

(5.38.Í)

-dy (5.38.J)

F3,,=" J e y

(5.38.k)

Í ax' - y i

(5.38.1)

76

, - íE l i j = 4 j Y i — ^ d y + 4 j Y i - — T d y (5.38.m)ox _y, õxõy

Ç ôY-. Çy> ÕY:

ij=6j Y d 2j YÇ ô2Y-. Ç

E 2 i j = 6 j Y i T T - d y + 2 j Yi — - dy (5.38Jd 5dx

E 3 i i = 4 j Y i — - d y (5.38.0)- y , OX

FCD i = Re (FC i ) - FD i (5.38.p)

HOü=Re(FOu)-EOu (5.38.q)

Hlu = R e ( n u ) - E l u (5.38.r)

H2 i j =Re(F2 i j ) -E2 i j (5.38.s)

-ES,, (5.38.t)

O passo seguinte na aplicação da GITT, é a avaliação dos coeficientes, oriundos da

transformação integral em y.

Uma análise do sistema diferencial ordinário resultante mostra que em boa parte o

custo computacional de obter-se resultados em problemas de domínio irregular via GITT

recai sobre o cálculo dos coeficientes.

Avaliar diretamente os coeficientes em forma numérica traz a desvantagem que eles

podem somente ser calculados dentro do processo de solução do sistema diferencial

ordinário e não previamente. Isto porque os coeficientes têm uma dependência da posição

"x", o que traria um alto custo computacional e como conseqüência uma baixa performance

do código computacional como um todo.

77

Uma outra estratégia para a obtenção dos coeficientes foi estudada, a qual permitiu

uma abordagem analítica.

Em primeiro lugar uma transformação de domínio, em y, foi realizada. Valendo-se

da relação (5.13), foi possível transformar os limites de integração de [—Vj , y2 ] para [—1,

1] e, através da regra da cadeia, foi possível estabelecer as derivadas, em termos de 4 , que

aparecem nos coeficientes.

Em segundo lugar, uma análise da nova estrutura dos coeficientes permitiu

descobrir que eles eram formados por uma combinação de integrais, as quais denominamos

de "núcleos", que apareciam repetidamente nos coeficientes. O interessante destas

integrais núcleo é que elas são independentes da posição "x", permitindo portanto, calcular

os coeficientes em separado, uma única vez, e armazená-los num arquivo. Portanto, o

custo de calcular os coeficientes numa determinada posição foi reduzido substancialmente,

ao ponto de serem necessárias apenas algumas operações aritméticas para a sua obtenção.

Para esta abordagem analítica no cálculo dos coeficientes, se fez uso do software

MATHEMATICA [ 63 ], a fim de facilitar toda a manipulação algébrica necessária.

Por sua vez, para uma simplificação no processo de aplicação da regra da cadeia se

fez uso das relações que são mostradas no apêndice C. No mesmo anexo mostra-se

também a estrutura dos coeficientes assim encontrada.

V.7 - ALGORITMO COMPUTACIONAL

Tal como visto no caso de dutos regulares, é necessário truncar as séries infinitas

em um número de termos suficientemente grande que garanta o erro relativo prefixado para

obtenção dos potenciais originais.

Portanto, o sistema diferencial ordinário pode ser reescrito como:

78

Revví dyk

LL\¥00 +F01 i j k ijT

¥kiF03,ik V- 1 ? - + F101Jk

(5.39)

onde Nt é a ordem de truncamento das séries infinitas.

É também necessário transformar o domínio infinito para um outro finito, usando a

mesma relação que para o caso de dutos regulares, equação (2.32).

A natureza híbrida da GITT é evidenciada ao resolver-se numericamente o sistema

diferencial ordinário não-linear. Este sistema diferencial ordinário possui um

comportamento stiff acentuado, o qual se intensifica ao aumentar-se o número de Reynolds.

O uso de subrotinas matemáticas encontradas em softwares bem estabelecidos, que

possuem controle automático do erro global na solução de um problema diferencial

ordinário é muito conveniente. Este fato, acompanhado da característica do controle do

erro do potencial original pela ordem de truncamento da expansão, fazem da presente

metodologia uma ferramenta muito forte na obtenção de resultados acurados.

Assim, foi selecionada a subrotina DBVPFD do IMSL [ 42 ] para obtenção dos

resultados numéricos, como na seção anterior para dutos regulares.

79

A implementação do programa para uso dessa subrotina torna necessário

transformar o sistema diferencial ordinário de 4â ordem em um outro de lâ ordem. Para

isto é necessário definir-se:

X; =yi (5.40.a)

j dvjj:

~ " Xi+Nt" "17 (5-40-b)

(5-40-d)

dXi+3Nt _ d V idx dx4 (5.40.e)

Devido à transformação ao domínio finito de [0, 1] é conveniente definir-se o

seguinte operador:

dx dr| dx

- ^ = c[l--n] (5.42)dx

Portanto o sistema pode ser reescrito como:

!ÍÍ1L (5.43 .a)dr| (dri/dx)

80

dXi+Nt Xi+2Nt

~A = (A IA ^ (5-43-b)

dr| (dri/dx)

(dT!/dx)

Nt

+X jXk+NtF01 i j k +F02ijk XjXk+2Nt

k=l

F03yk XjX k + 3 N t +F10jjk Xj+ N tXk +F l ly k Xj+ N tXk + N t

ijk Xj+NtXk+2Nt

Nt

1 F C D i a 1 /H25j Xj+2Nt +H3ÍJ Xj+3Nt ] + - — — vt V, /(dTi/dx)N

(5.43.d)

81

VI - ESCOAMENTO NO INTERIOR DE UM DUTO COM

EXPANSÃO GRADUAL

Uma grande variedade de problemas envolvendo dutos irregulares são encontrados

na literatura, os quais têm sido resolvidos por métodos numéricos tradicionais, como

volumes finitos, elementos finitos ou diferenças finitas.

Um caminho para comparar performances de diferentes metodologias de solução na

abordagem de problemas de dutos de geometria complexa foi estabelecido pelo V

Congresso da International Association for Hydraulic Research (IAHR), no qual foi

especificado como caso teste um problema recomendado por Roache [64]. O congresso

foi realizado em Roma em 1982, com participação de 15 grupos, os quais forneceram seus

resultados para tal problema proposto. A compilação dos trabalhos foi feita por Napolitano

e Orlandi [ 65 ], fazendo uma extensa comparação entre eles.

Devido à ampla quantidade de resultados, o problema proposto por Roache foi

também escolhido no presente trabalho como caso teste.

VI. 1 - DEFINIÇÃO DO PROBLEMA

Trata-se de um duto o que se expande gradualmente e mantém simetria ao longo do

eixo x. A geometria das paredes é definida analiticamente, tal como se mostra na fig. 6.1.

Como se observa na fig. 6.1, a geometria do duto é dependente do número de

Reynolds, definido a partir da velocidade média e metade da distância entre as paredes na

entrada do duto.

É evidente que quando o número de Reynolds aumenta o canal fica mais comprido

e estreito. Roache [ 64 ], mostra que para grandes valores do número de Reynolds o

escoamento exibe uma forma similar. Isto deve-se ao fato de que para altos Reynolds , o

duto tende à uma situação de placas paralelas.

82

A condição de contorno, na entrada do duto, é de velocidade uniforme e paralela:

ou em termos de função corrente:

Considera-se que o duto tem um comprimento finito, sendo truncado na posição:

Rexout =

As condições de contorno na saída do duto não foram especificadas, deixando-as

em aberto para a escolha de cada participante do congresso na solução do problema teste.

No presente estudo, o problema teste será resolvido considerando-se tanto o duto

trancado como duto infinito, com as condições de contorno adequadas em cada caso.

No caso de duto finito, as condições de contorno na saída são especificadas

similarmente às usadas por Magi & Napolitano (diferenças finitas), e Quartapelle &

Napolitano (elementos finitos) [ 65 ].

= 0 ; v (x o u t , y ) = 0 (6.4.a,b)

Isto permitirá comparar resultados obtidos por outras metodologias, as quais

fizeram uso das mesmas condições de contorno.

83

u=0 v=0

3̂

2v = 0

(0,1)

(0,-1)

u = 0 v=0

tanh(2)-tanh 2 -30x

V Re

FIGURA. 6.1 - Geometria do Problema teste: Canal com expansão gradual

Para o caso de duto infinito, as condições de contorno especificadas são as do

escoamento completamente desenvolvido, uma vez que a geometria tende para um canal de

placas paralelas para cada caso:

J3_2a - f v(°o,y) = (6.5.a,b)

Assumindo que o duto é infinito, as condições de contorno na saída ficam melhor

definidas, removendo assim a incerteza da posição do truncamento do duto e a correta

especifição das condições de contorno.

Na figura (6.2.a, b) são mostradas as oito condições de contorno para os dois casos

de estudo.

84

^ = 0

t y

&=°ãc = o

(a)

xout = Re/3

5vi/- . 0

(oo, a)

3 y

(oo, -a)

(b)

FIGURA 6.2.a, b - Condições de contorno no duto com expansão gradual.

85

VI.2 - REVISÃO BIBLIOGRÁFICA

Os resultados do problema teste, apresentados no V Congresso do IAHR, foram

resumidos e comparados por Napolitano e Orlandi [ 65 ]. Alguns dos trabalhos desse

encontro podem ser citados:

K.A. Cliffe et ai. [ 66 ] usaram o método de elementos finitos para resolver as

equações de Navier-Stokes, as quais foram formuladas em variáveis primitivas. Utilizaram

elementos isoparamétricos com nove nós e interpolação biquadrática. Os resultados

obtidos por este grupo foram considerados como benchmark.

V. Magi e M. Napolitano [ 67, 68 ] usaram uma formulação transiente de

vorticidade - função corrente. As equações foram discretizadas pelo método de diferenças

finitas, utilizando um esquema ADI.

Quartapelle e M. Napolitano [ 69 ] usaram uma formulação transiente de

vorticidade - função corrente, mas empregaram o método de elementos finitos para a

discretização, usando elementos quadriláteros do tipo bilinear.

A. Wada e K. Ada Chi [ 65 ] apresentaram duas soluções do problema. Empregaram

um método que envolve transformação de coordenadas usando funções analíticas

complexas para a transformação de variáveis e diferenças finitas para sua solução

numérica.

W. Schonauer [ 65 ] usou o método de diferenças finitas com malha autoadaptativa,

para uma formulação de velocidade - vorticidade. Seus resultados mostraram boa

concordância com os de Cliffe et ai. [ 66 ].

I. Demirovic e A. D. Gosman [ 65 ] usaram uma formulação em variáveis primitivas

das equações de Navier-Stokes e o método de diferenças finitas para sua solução. Foi

empregado um ADI iterativo para resolver as equações simultâneas.

86

Existem na literatura outras referências que utilizaram este problema para validação

de códigos gerados para solução de problemas em dutos de domínio arbitrário [ 70 - 74 ]

que mostram similar acuracia nos resultados, com relação aos apresentados no congresso

em questão [ 65 ].

VI.3 - SOLUÇÃO DO PROBLEMA

A solução do problema de escoamento em desenvolvimento no interior de dutos

irregulares foi apresentada no capítulo anterior numa forma geral, deixando-se em aberto a

especificação das condições de contorno na entrada e saída do duto. Portanto, a seguir

mostramos apenas a transformação integral das condições de contorno para os dois casos

analisados.

VI.3.1 - DUTO TRUNCADO

As condições de contorno filtradas são dadas por:

- na entrada (0, y)

v|/* = i|/ — F (6.6.a)

(6.6.b)ôx õx õx

- na saída (xout, y)

dij/ õ\\f õF

õx õx õx

5x3 + õxõy2 ôx3 + õxõy2 õx3 õxõy2

(6.6.c)

87

Operando com JY; dy e ordenando convenientemente:- yi

1 Ç y i

= - J Yi[M/(0,y)-F(0,y)]dy (6.7.a)

1 fy i ra|/(O,y) SF(0,y)l f̂ _ f

Neste ponto, é necessário advertir que na entrada do duto existe uma singularidade,

uma vez que a geometria que define as paredes do duto não tem gradiente nulo,

questionando assim, a validade de ter-se um perfil completamente desenvolvido.

Portanto, admitir escoamento completamente desenvolvido justifica assumir que o

filtro não tem dependência em x, nessa posição, resultando então:

Nota-se que o filtro é exatamente igual ao valor da função corrente em x=0. Desta

forma, as condições de contorno ficam reduzidas a:

ijT, (0) = 0 (6.9.a)

= 0 (6.9.b)dx

A transformação das outras duas condições de contorno fornece as seguintes

relações acopladas:

dx

dx

88

1 fyz OF v ^ _• = — — I Y: — d y - / 11/

ÕY;Y,

j=i1 dx

• = — Y:

fa3õ3F a 3 F V

õxdydy

ax3 axay•

j \ A OUt s I

dx

dx- J- y i

- y i aydy +

(6.10.a)

(6.10-b)

VI.3.2 - PUTO INFINITO

Como no item anterior, o primeiro passo é especificar a variável filtrada:

— na entrada (0, y):

\|/ = \\i - F

d\|/ ÕF

dx 3x ôx

(6.11.a)

(6.1 l.b)

- na saída (co, y):

\\f* - \\i - F (6.12.a)

89

ôy õv/ SF~ã ~~x~~!T~ (6.12.b)õx ôx ôx v '

Tanto na entrada como na saída do duto, o filtro F adota os mesmos valores que a

função corrente em cada uma dessas posições; realizando a transformação integral resulta:

(6.13.a,b)

Y ^ =0 (6.14.a,b)

VI.4 - RESULTADOS

Foi gerado um código em FORTRAN - 77, o qual foi executado no CRAY EL - 94

do Núcleo de Atendimento de Computação de Alto desempenho (NACAD) da

COPPE/UFRJ.

As tarefas numéricas envolvidas no caso de duto irregular, ficaram restritas à

solução do sistema diferencial ordinário não linear de quarta ordem, que surge como

conseqüência da transformação integral do problema. Para tal fim, foi usada a subroutina

DBVPFD do IMSL [ 42 ], na qual um erro relativo global de IO"4 foi sempre solicitado.

Alguns coeficientes núcleo foram calculados numericamente com as subrotinas DQDAG e

DQDAGS, ambas do IMSL. Os outros coeficientes foram calculados analiticamente e em

alguns casos o desenvolvimento analítico foi auxiliado pelo uso do software

MATHEMATICA [ 63 ].

Em primeiro lugar, foi estabelecido o processo de convergência da função corrente,

nos dois tipos de dutos considerados e para os dois valores de número de Reynolds

90

analisados (Re = 10, 100). Os resultados para ambos os casos, ao longo da linha y = 0.5,

são mostrados nas tabelas 6.1 e 6.2, respectivamente.

Foram requeridos somente Nt = 6 termos em todos os casos estudados, para

satisfazer à tolerância solicitada no campo original. Os cálculos foram conduzidos até Nt =

10, para ilustrar como o potencial original se mantém inalterável, verificando-se portanto a

convergência na precisão requerida.

Quando o IAHR definiu o problema teste, tratou que ele fosse de tal forma que

qualquer tratamento numérico dos termos convectivos fosse de segunda importância, dado

que problemas com dominância dos termos convectivos eram extremamente sensíveis ao

processo de discretização [ 75 ] . Este fato explica o baixo número de termos requerido

para a convergência na presente solução, uma vez que a expansão empregada é fortemente

baseada nos termos difusivos.

O efeito do truncamento do duto é também registrado nas tabelas 6.1 e 6.2. Em

posições próximas da entrada, a condição de contorno na saída do duto tem pouca

influência, pois o escoamento, tanto no duto infinito como finito, comporta-se

identicamente, já que a condição de contorno na entrada é a mesma. Para posições um

pouco mais distantes da entrada do duto, os valores da função corrente de ambos os casos

estudados ( duto truncado e infinito ) ficam ligeiramente diferentes, intensificando-se a

diferença para distâncias maiores como conseqüência do truncamento. No caso de Re =

10, a diferença de valores da função corrente, entre o duto truncado e o infinito, é de 1.4%;

já para Re = 100 a diferença se reduz a 0.1%. A melhor concordância para este segundo

caso deve-se ao fato de se ter um canal mais longo e daí uma expansão mais suave.

Nas figuras 6.3 e 6.4, apresenta-se as isolinhas da função corrente. Se observa o

surgimento de vórtices localizados, os quais ficam mais marcantes para o menor valor de

Re.

91

TABELA 6.1 -Convergência da Função Corrente \\> ( x / x out , 0.5)

( R e - 1 0 )

:'/ OUt :

•:'._. ;.;.;;.;.; .; ;.;.;.;..-.;. .-. . i

0

0.05

0.10

0.15

0.20

0.25

0.30

0.35

0.40

0.45

0.50

0.55

0.60

0.65

0.70

0.75

0.80

0.85

0.90

0.95

1.00

1.6875

1.6762

1.6617

1.6446

1.6257

1.6063

1.5869

1.5680

1.5500

1.5330

1.5171

1.5023

1.4886

1.4760

1.4645

1.4539

1.4444

1.4357

1.4280

1.4210

1.4148

Düto Infiiiitb

1.6875

1.6761

1.6618

1.6447

1.6259

1.6063

1.5869

1.5680

1.5500

1.5330

1.5171

1.5023

1.4886

1.4760

1.4645

1.4540

1.4444

1.4358

1.4280

1.4210

1.4148

Éiill

1.6875

1.6762

1.6617

1.6447

1.6259

1.6064

1.5869

1.5680

1.5500

1.5330

1.5171

1.5023

1.4886

1.4760

1.4645

1.4540

1.4444

1.4358

1.4280

1.4210

1.4158

1.6875

1.6761

1.6617

1.6447

1.6259

1.6064

1.5869

1.5680

1.5500

1.5330

1.5171

1.5023

1.4886

1.4761

1.4645

1.4540

1.4444

1.4358

1.4280

1.4210

1.4158

i l l l l1.6875

1.6762

1.6617

1.6446

1.6257

1.6064

1.5870

1.5682

1.5504

1.5336

1.5179

1.5034

1.4902

1.4782

1.4676

1.4582

1.4502

1.4437

1.4388

1.4358

1.4347

um1.6875

1.6762

1.6618

1.6447

1.6260

1.6065

1.5871

1.5683

1.5504

1.5336

1.5179

1.5035

1.4903

1.4783

1.4676

1.4582

1.4503

1.4438

1.4389

1.4358

1.4347

•111!1.6875

1.6762

1.6618

1.6447

1.6259

1.6065

1.5871

1.5683

1.5504

1.5336

1.5179

1.5035

1.4903

1.4783

1.4676

1.4582

1.4503

1.4438

1.4389

1.4358

1.4347

If*-11.6875

1.6762

1.6618

1.6447

1.6259

1.6065

1.5871

1.5683

1.5504

1.5336

1.5179

1.5035

1.4903

1.4783

1.4676

1.4582

1.4503

1.4438

1.4389

1.4358

1.4347

92

TABELA 6.2 - Convergência da Função Corrente \j/ (x / x out , 0.5)

(Re= 100)

W ;í

%/ x p u t ;0

0.05

0.10

0.15

0.20

0.25

0.30

0.35

0.40

0.45

0.50

0.55

0.60

0.65

0.70

0.75

0.80

0.85

0.90

0.95

1.00

1.6875

1.6793

1.6619

1.6340

1.6041

1.5782

1.5559

1.5355

1.5166

1.4991

1.4831

1.4688

1.4561

1.4448

1.4349

1.4261

1.4185

1.4119

1.4061

1.4011

1.3968

1.6875

1.6794

1.6619

1.6340

1.6041

1.5782

1.5559

1.5355

1.5166

1.4991

1.4831

1.4688

1.4561

1.4448

1.4349

1.4261

1.4185

1.4119

1.4061

1.4011

1.3968

wmm.1.6875

1.6794

1.6619

1.6340

1.6041

1.5782

1.5559

1.5355

1.5166

1.4991

1.4831

1.4688

1.4561

1.4448

1.4349

1.4261

1.4185

1.4119

1.4061

1.4011

1.3968

1.6875

1.6794

1.6619

1.6340

1.6041

1.5782

1.5559

1.5355

1.5166

1.4991

1.4831

1.4688

1.4561

1.4448

1.4349

1.4261

1.4185

1.4119

1.4061

1.4011

1.3968

; M : ' : 7 Í V . " : .•••V-^yy-íy-:y>y í f i : f > y y - •: . • . . . ; : : ; ; y i y y y y y y y y . - y •]:•• •• • • : • • . " . - • • • • • • • : • . • : : .

iOíy-: ; í : • í-..'-; y< :ÍV:ÍÍ ;.'i :.;:>.; • i K s f :::::i;:::í:V: : í : i : í i : í f : í í í : íf ™ i •': y lii 'f •-.•'. • y yi\hy : .í:;:í ' í

igmãã

1.6875

1.6793

1.6619

1.6340

1.6041

1.5782

1.5559

1.5356

1.5166

1.4991

1.4832

1.4689

1.4561

1.4448

1.4349

1.4262

1.4186

1.4119

1.4062

1.4012

1.3983

S.1ÜÍI1

1.6875

1.6794

1.6619

1.6340

1.6041

1.5782

1.5559

1.5356

1.5166

1.4991

1.4832

1.4689

1.4561

1.4448

1.4349

1.4262

1.4185

1.4119

1.4062

1.4012

1.3983

isms1.6875

1.6794

1.6619

1.6340

1.6041

1.5782

1.5559

1.5356

1.5166

1.4991

1.4832

1.4689

1.4561

1.4448

1.4349

1.4262

1.4186

1.4119

1.4062

1.4012

1.3983

1.6875

1.6794

1.6619

1.6340

1.6041

1.5782

1.5559

1.5356

1.5166

1.4991

1.4832

1.4689

1.4561

1.4448

1.4349

1.4262

1.4186

1.4119

1.4062

1.4012

1.3983

93

Re= 10

1.0 1.5 2.0 2.5 3.0 X

FIGURA 6.3 - Isolinhas da Função Corrente [ Re = 10 ].

Re = 100

y

1.5

1.0

0.5

n n

^^, -y ' 1~

/ V u ^ ^ ^^-— •

^———' - 1.5-

=—— I I I ^ - - - - ^ ^r • " '__ ' ' 3

, , 1 , , , , 1 . , , . 1 . • . . 1 . . . . 1

•

_ „ -. •

—1.6

—1.4

. i , . , 1 , • ,

IS 20 25 30 X

FIGURA 6.4 - Isolinhas da Função Corrente [ Re =100 ].

94

Em seqüência ao presente estudo, calculou-se a vorticidade nas paredes do duto,

possibilitando comparações com resultados obtidos através de outras metodologias.

O processo de convergência da vorticidade na parede é mostrado nas tabelas (6.3.a,

b) para Re = 10 e (6.4.a, b) para Re = 100. De fato, como pode apreciar-se, um maior

número de termos foi requerido a fim de se respeitar a tolerância imposta em relação à

convergência da função corrente. Em alguns casos, para posições muito próximas da

entrada, é notório que a convergência da vorticidade ainda não foi completamente atingida

com as ordens de truncamento empregadas.

Este comportamento de convergência retardado em relação ao caso da função

corrente era esperado. A segunda derivada da função corrente na direção y, requerida no

cálculo da vorticidade, incide diretamente na autofunção, trazendo ao numerador da

expansão o quadrado do auto-valor, adquirindo maior importância para ordens cada vez

mais crescentes, manifestando-se num retardamento de convergência.

A vorticidade, calculada nos dois tipos de dutos estudados, mostra pouca diferença

entre si para Re = 100 (menos de 1.9%) e uma distância muito mais pronunciada para o

caso de Re = 10 (menos de 19%).

95

TABELA 6.3.a, b - Convergência da vorticidade na parede do duto

(Re =10)

w%/... Xput

0

0.05

0.10

0.15

0.20

0.25

0.30

0.35

0.40

0.45

0.50

0.55

0.60

0.65

0.70

0.75

0.80

0.85

0.90

0.95

1.00

•till

3.000

2.657

1.672

0.478

-0.222

-0.138

-0.087

-0.088

-0.079

-0.046

0.006

0.069

0.135

0.201

0.265

0.232

0.378

0.427

0.470

0.509

0.543

3.000

2.707

1.570

0.362

-0.087

-0.094

-0.076

-0.083

-0.078

-0.047

0.004

0.068

0.136

0.203

0.267

0.327

0.381

0.430

0.473

0.512

0.546

'••yyyy)Xym-y--'''-V.

3.000

2.727

1.536

0.355

-0.060

-0.086

-0.073

-0.082

-0.079

-0.048

0.004

0.068

0.136

0.204

0.268

0.327

0.382

0.430

0.474

0.512

0.546

• • l i3.000

2.736

1.520

0.357

-0.051

-0.084

-0.072

-0.082

-0.079

-0.049

0.003

0.068

0.136

0.204

0.268

0.328

0.382

0.431

0.474

0.513

0.547

lilli!3.000

2.657

1.672

0.477

-0.223

-0.139

-0.087

-0.089

-0.080

-0.047

0.038

0.064

0.129

0.192

0.251

0.303

0.349

0.385

0.412

0.428

0.434

I l l l l l l l

3.000

2.706

1.569

0.362

-0.087

-0.095

-0.076

-0.083

-0.079

-0.049

0.002

0.064

0.130

0.194

0.254

0.307

0.353

0.390

0.417

0.434

0.439

;i;i;i|li||

3.000

2.727

1.535

0.355

-0.601

-0.087

-0.073

-0.082

-0.079

-0.050

0.002

0.064

0.130

0.195

0.255

0.308

0.354

0.391

0.418

0.435

0.440

iiiii3.000

2.736

1.520

0.357

-0.051

-0.084

-0.072

-0.082

-0.080

-0.050

0.001

0.064

0.131

0.195

0.255

0.309

0.355

0.392

0.419

0.435

0.441

96

TABELA 6.4.a, b - Convergência da vorticidade na parede do duto

(Re =100)

: : / Xout :

0

0.05

0.10

0.15

0.20

0.25

0.30

0.35

0.40

0.45

0.50

0.55

0.60

0.65

0.70

0.75

0.80

0.85

0.90

0.95

1.00

3.000

2.541

1.743

0.733

0.102

-0.092

-0.126

-0.095

-0.029

0.054

0.142

0.224

0.299

0.366

0.424

0.474

0.516

0.553

0.584

0.611

0.634

W§Ê§,

3.000

2.530

1.728

0.726

0.102

-0.092

-0.125

-0.095

-0.029

0.055

0.142

0.225

0.300

0.367

0.425

0.475

0.517

0.554

0.585

0.612

0.635

3.000

2.528

1.726

0.725

0.102

-0.092

-0.125

-0.095

-0.029

0.055

0.142

0.225

0.301

0.367

0.425

0.475

0.518

0.554

0.585

0.612

0.635

l l l l l l l lf l

3.000

2.528

1.725

0.724

0.102

-0.092

-0.125

-0.095

-0.029

0.055

0.143

0.225

0.301

0.367

0.425

0.475

0.518

0.554

0.585

0.612

0.635

lillij;

3.000

2.541

1.743

0.733

0.102

-0.092

-0.126

-0.095

-0.029

0.054

0.141

0.224

0.299

0.366

0.423

0.473

0.516

0.553

0.584

0.610

0.623

3.000

2.530

1.728

0.726

0.102

-0.092

-0.125

-0.095

-0.029

0.055

0.142

0.225

0.300

0.367

0.425

0.475

0.517

0.554

0.585

0.611

0.623

iiilii!3.000

2.528

1.726

0.725

0.102

-0.092

-0.125

-0.095

-0.029

0.055

0.142

0.225

0.300

0.367

0.425

0.475

0.517

0.554

0.585

0.611

0.623

3.000

2.528

1.725

0.724

0.102

-0.092

-0.125

-0.095

-0.029

0.055

0.142

0.225

0.301

0.367

0.425

0.475

0.518

0.554

0.585

0.611

0.623

97

Gráficos comparando o valor da vorticidade na parede com aqueles obtidos por

outras metodologias são apresentados nas figuras (6.5) e (6.6). Para o caso de Re = 100, os

resultados do presente trabalho concordam muito bem com os demais, mas para Re = 10

uma ligeira discrepância é observada na região de entrada do duto, possivelmente por uma

diferença na definição da vorticidade na parede do duto. Em [ 65 ], apenas um dos

participantes define que a vorticidade deverá ser avaliada pela expressão:

co = ducos9- õvsinG (6.15)

onde 9 é o ângulo entre a normal à parede e a direção da coordenada x; õ u, õ v são os

valores das normais das componentes da velocidade nas direções x e y.

A implementação desta fórmula forneceu resultados díspares. Uma outra

expressão, definindo a vorticidade como:

õu õvoo = —cosG- — sinG (6.16)

õy ox

onde 9 é o ângulo entre a tangente à curva que define a parede e a direção positiva do eixo

x, forneceu os resultados mostrados nas figuras (6.5) e (6.6).

4.00

-1.00

Re = 10

••oA

OA

GITT (duto infinito)

GITT (duto truncado)

Cliffe et. ai. [ 66 ]

Magi et. ai. [ 67,68]

Quartapelle et. ai. [ 69 J

Demirovic et. ai. [ 65 ]

Wada et. ai. [ 65 ]

Schonaueret. ai. [ 651

0.00 0.20 0.40 0.60X/Xout

0.80 1.00

<B

CO

a.

ICOTO

O

4.00 - n

3.00

2.00

1.00 —

0.00

-100

••oA

O•

Re =100

- GITT (duto infinito)

- GITT (duto truncado)

Cliffe et ai. [66]

Magl et ai. [ 67, 68 ]

Quartapelle et ai. [ 69 ]

Demirovic et ai. [65]

Wada et ai. [65]

Schonaueret ai. [65]

0.00 0.20 0.40 0.60X/Xout

0.80 1.00

FIGURAS 6.5, 6 - Vorticidade ao longo da parede do duto ( Re = 10 , 100)

99

VII - CONCLUSÕES

- A Técnica da Transformada Integral Generalizada, mostrou-se uma ferramenta capaz de

conseguir resultados acurados em problemas de escoamentos dentro de canais de geometria

arbitrária. A característica de controlar automaticamente o erro relativo no campo estudado

faz com que a presente metodologia seja empregada com sucesso como ferramenta para

obter-se resultados benchmark.

Outra vantagem oferecida pela GITT é que pode oferecer resultados bastante

acurados com poucos termos nas expansões que definem os campos originais, traduzindo-

se isto num baixo custo computacional. Esta vantagem é muito interessante para obtenção

de resultados práticos de engenharia.

- A estratégia de se resolver o problema de dutos sem truncamentos, transformando o duto

infinito num domínio finito por meio de uma relação algébrica, mostrou-se efetiva e

prática, evitando-se a incerteza de onde truncar o duto. Além disso, desta forma é possível

obter-se melhores resultados próximos à saída do duto, sem serem perturbados pelo

truncamento.

- A correta especificação das condições de contorno é de crucial importância no processo

de simulação do escoamento. Condições de contorno inadequadas destroem a coerência

física e portanto a acurácia dos resultados.

- No caso de dutos regulares é adequado usar o escoamento completamente desenvolvido

como solução filtro, pois além de homogeneizar o contorno, serve como acelerador de

convergência, mais marcante em zonas distantes da entrada.

- Foi gerado um código capaz de resolver problemas de escoamentos em dutos de

geometria arbitrária para escoamentos bidimensionais em regime laminar e estado

permanente. A flexibilidade do código desenvolvido permite que no "input" se identifique

as condições de contorno na entrada e saída do canal, a geometria do canal e o número de

Reynolds em que se deseja simular o escoamento. Esta experiência de solução pode ser

100

agora estendida para casos mais gerais, que possam ter dependência no tempo, geometria

tridimensional e regime turbulento.

- Em seqüência ao presente estudo, buscar-se-á a implementação do problema térmico, o

qual permitirá estudar problemas de transferência de calor em dutos de geometria irregular.

101

VIII - REFERÊNCIAS

[1] COTTA, R.M., "Integral Transforms in Computational Heat e Fluid Flow". CRC

Press, Boca Raton, FL., 1993.

[2] ÕZISIK, M.N. e MURRAY, R.L., "On the Solution of Linear Diffusion Problems

with Variable Boundary Conditions Parameters", J. Heat Transfer. V.96, pp.48-

51, 1974.

[3] COTTA, R.M. e SANTOS, C.A.C., "Transient Diffusion Problems with Time-

Dependent Boundary Condition Coefficients", J. Eng. Physics. V.61, nfl 5,

pp.1411-1418,1992.

[4] YENER, Y. e ÕZISIK, M.N., "On the Solution of Unsteady Heat Conduction in

Multi-Region Finite Media with Time Dependent Heat Transfer Coefficient",

Anais da 5° Conf. Int. Transf. de Calor. Tóquio, Setembro, 1974.

[5] SANTOS, C.A.C., COTTA, R.M. e ÔZISIK, M.N., "Laminar Forced Convection

Inside Externally Finned Tubes", Anais do 2o Encontro Nacional de Ciências

Térmicas - ENCIT 88. pp.87-90, Rio de Janeiro, 1988.

[6] COTTA, R.M. e ÕZISIK, M.N., "Diffusion Problems with General Time

Dependent Coefficients", Revista Brasileira de Ciências Mecânicas. RBCM, V.9,

na4,pp.269-292, 1987.

[7] COTTA, R.M., LEIROZ, A.J.K. e APARECIDO, J.B., "Steady- State Diffusion

Problems with Variable Equation Coefficients", Int. J.Heat & Technology. V.10,

ns 3-4, pp.104, 1992.

[8] MIKHAILOV, M.D., "On the Solution of the Heat Equation with Time

Dependent Coefficient". Int. J. Heat & Mass Transfer. V.18, pp.334-345,1975.

[9] APARECIDO, J.B., COTTA, R.M. e ÔZISIK, M.N., "Analytical Solutions to

Two-Dimensional Diffusion Type Problems in Irregular Geometries", J. Franklin

Inst. V. 326, pp.421-434, 1989.

[10] APARECIDO, J.B. e COTTA, R.M., "Laminar Flow Inside Hexagonal Ducts",

Computational Mech.. V.6, pp.93-100, 1990.

102

[11] DINIZ, A.J., APARECIDO, J.B. e COTTA, R.M, "Heat Conduction with

Ablation in a Finite Slab", Int. J. Heat & Technology. V.8, pp.30-43, 1990.

[12] APARECIDO, J.B. e COTTA, R.M., "Laminar Thermally Developing Flow

Inside Right Triangular Ducts", Applied Scientific Research. V.49, pp.355-368,

1992.

[13] ÕZISIK, M.N. e GUÇERI, S.I., "A Variable Eigenvalue Approach to the Solution

of Phase-Change Problems", Canadian J. Chem. Eng.. V.55, pp. 145-148, 1977.

[14] LEITE, S.Q.B., ÕZISIK, M.N. e VERGHESE, K., "On the Solution of Linear

Diffusion Problems in Media with Moving Boundaries", Nucl. Science & Eng..

V.76, pp.345-350, 1980.

[15] COTTA, R.M., "Diffusion in Media with Prescribed Moving Boundaries

Application to Metals Oxidation at High Temperatures", Anais do II Congresso

Latino-Americano de Transferência de Calor e Massa. V.I, pp.502-513, São

Paulo, 1986.

[16] COTTA, R.M. e ÕZISIK 5, M.N., "Laminar Forced Convection in Ducts with

Periodic Variation of Inlet Temperature", Int. J. Heat & Mass Transfer. V.29, n°

10, pp.1495-1501, 1986.

[17] APARECIDO, J.B. e COTTA, R.M., "Thermally Developing Laminar Flow

Inside Rectangular Ducts", Int. J. Heat & Mass Transfer. V.33, nQ 2, pp.341-347,

1990.

[18] COTTA, R.M., "On the Solution of Periodic Multidimensional Diffusion

Problems", Int. Comm. Heat & Mass Transfer. V.I6, pp.549-579, 1989.

[19] COTTA, R.M. e ÕZISIK, M.N., "Transient Forced Convection in Laminar

Channel Flow with Stepwise Variations of Wall Temperature", Can. J. Chem.

Eng.. V.64, pp.734-742, 1986.

[20] GUEDES, R.O.C., COTTA, R.M. e BRUM, N.C.L., "Conjugated Heat Transfer

in Laminar Flow Between Parallel-Plates Channels", Anais do X Congresso

Brasileiro de Enegenharia Mecânica, pp.285-288, Rio de Janeiro, Dezembro,

1989.

103

[21] KAKAÇ S., LI W. e COTTA, R.M., "Theoretical e Experimental Study of

Transient Laminar Forced Convection in Duct with Timewise Variation of Inlet

Temperature", ASME Winter Annual Meeting. HTD. V.123, pp.57-63, San

Francisco, 1989.

[22] BAYAZITOGLU, Y. e ÒZISIK, M.N., "On the Solution of Graetz Type

Problems with Axial Conduction", Int. J. Heat & Mass Transfer. V.23, pp. 1399-

1402, 1980.

[23] COTTA, R.M., "Hybrid Numerical Analytical Approach to Nonlinear Diffusion

Problems", Num. Heat Transfer. Part B - Fundamentals. V.I7, pp.217-226, 1990.

[24] SERFATY, R. e COTTA, R.M., "Integral Transform Solutions of Diffusion

Problems with Nonlinear Equation Coefficients", Int. Comm. Heat & Mass

Transfer, V. 17, n~ 6, pp.852-864,1990.

[25] SERFATY, R. e COTTA, R.M., "Hybrid Analysis of Transient Nonlinear

Convection -Diffusion Problems", Int. J. Num. Meth. Heat & Fluid Flow. V.2,

pp.55-62, 1992.

[26] COTTA, R.M. e SERFATY, R., "Integral Transform Algorithm for Parabolic

Diffusion Problems with Nonlinear Boundary e Equation Source Terms", Proc. of

the 7° Int. Conf. on Num. Meth. for Thermal Problems. Part 2. pp.916-926,

Stanford, CA, Julho 1991.

[27] COTTA, R.M. e MIKHAILOV, M.D., "The Integral Transform Method", Appl.

Math. Modeling. V.17, pp.156-161, 1993.

[28] COTTA, R.M. e CARVALHO, T.M.B., "Hybrid Analysis of Boundary Layer

Equations for Internal Flow Problems", Proc. of the 7a Int. Conf. on Num. Meth.

in Laminar & Turbulent Flow. Part 1. pp. 106-115, Stanford, CA, Julho 1991.

[29] CARVALHO, T.M.B., COTTA, R.M. e MIKHAILOV, M.D., "Flow

Development in the Entrance Region of Ducts", Comm. Num. Meth. Eng.. V.9,

pp.503-509, 1993.

104

[30] MACHADO, H.A. e COTTA, R.M., "Integral Transform Method for Boundary

Layer Equations in Simultaneous Heat e Fluid Flow Problems", Int. J. Num.

Meth. Heat & Fluid Flow. V.5, 225-237,1995.

[31] MIKHAILOV, M.D. e COTTA, R.M., "Unified Integral Transform Method", L

Braz. Assoc. Mech, Sciences. RBCM (Invited Paper), V.12, n2 3, pp.301-310,

1990.

[32] LEIROZ, A.J.K. e COTTA, R.M., "On the Solution of Nonlinear Elliptic

Convection -Diffusion Problems Through the Integral Transform Method", Num.

Heat Transfer - Part B: Fundamentals. V.25. 401 -41.1. 1993.

[33] BAOHUA, C. e COTTA, R.M., "Integral Transform Analysis of Natural

Convection in Porous Enclosures", Int. J. Num. Meth. in Fluids. V.I7, 787-801,

1993.

[34] PEREZ GUERRERO, J.S., "Solução das Equações de Navier- Stokes em

Formulação de Função Corrente via Transformação Integral", COPPE/UFRJ,

M.Sc, Tese, Rio de Janeiro, 1991.

[35] PEREZ GUERRERO, J.S. e COTTA, R.M., "Integral Transform Method for

Navier-Stokes Equations in Stream-Function Only Formulation", Int. J. Num.

Methd. in Fluids. V.15, pp.399-409, 1992.

[36] COTTA, R.M., PEREZ GUERRERO, J.S. e SCOFANO NETO, F., "Hybrid

Solution of the Incompressible Navier-Stokes Equations via Integral

Transformation", Proc. of the 2° Int. Conf. Advanced Computational Methods in

Heat Trasnfer. Heat Transfer 92. V.I, pp.735-750, Milão, Itália, Julho 1992.

[37] RIBEIRO, J.W. e COTTA, R.M., "Numerical-Analytical Study of Nonlinear

Drying Problems with Radiative Boundaries", 6 "International Symposium on

Transport Phenomena. V.I, pp.209-214, Seoul, Korea, Maio, 1993.

[38] COTTA, R.M., "Computational Integral Transform Approach in Nonlinear

Diffusion e Convection -Diffusion Problems", Laboratório de Ingegneria

Nucleane di Montecuccolino. Serie Scientifica (Invited Lecture), Universitá degli

Studi di Bologna, Itália, Julho 1992.

105

[39] RAMOS, R., "Análise Bidimensional e Não-Linear de Aletas Longitudinais de

Perfil Variável via Técnica de Transformada Integral Generalizada",

COPPE/UFRJ, M.Sc, Tese, Rio de Janeiro, 1993.

[40] PIMENTEL, L.C.G., "Simulação do Escoamento Turbulento em Dutos via

Modelo Algébrico de Turbulência e Transformação Integral", COPPE/UFRJ,

M.Sc, Tese, Rio de Janeiro, 1993.

[41] PÉREZ GUERRERO, J.S., COTTA, R.M. e SCOFANO NETO, F., "Integral

Transformation of Navier-Stokes Equations for Incompressible Laminar Flow in

Channels", Proc. 8° International Conference on Numerical Methods in Laminar e

Turbulent Flow. Swansea, v.2,pp. 1195-1206, Julho, 1993.

[42] IMSL Library, Math/Lib., Houston, TX 1989.

[43] PEREYRA, V., "PASVA3: An Adaptive Finite-Difference FORTRAN Program

for First Order Non Linear Boundary Value Problems", Lee. Notes in Comp.

Science. 76, Springer-Verlag, Berlim, pp.67-88, 1978.

[44] SHAH, R.K. e LONDON, A.L., "Laminar Flow Forced Convection in Ducts".

Adv. Heat Transfer, Suppl.l, 1978.

[45] WANG, Y.L. e LONGWELL, P.A., "Laminar Flow in the Inlet Section of

Parallel Plates", AICHE J.. V.10, pp.323-329, 1964.

[46] BRET, A. e GILLIS, J., "Magnetohydrodynamic Flow in the Inlet Region of a

Straight Channel". Phys. Fluids. V.9, pp.690-699, 1966.

[47] McDONALD, J.W., DENNY, V.E. e MILLS, A.F. , "Numerical Solutions of the

Navier-Stokes Equations in Inlet Regions", J. Appl. Mech.. V.39, pp.873-878,

1972.

[48] NARANG, B.S. e KRISHNAMOORTHY, G., "Laminar Flow in the Entrance

Region of Parallel Plates". J. Appl. Mech.. V.43, pp.186-188, 1976.

[49] COMINI, G. e DEL GIUDICE, S., "Finite-Element Solution of the

Incompressible Navier-Stokes Equations", Num. Heat Transfer. V.5, pp.463-478,

1982.

106

[50] BODOIA, J.R. e OSTERLE, J.F., "Finite Difference Analysis of Plane Poiseuille

e Couette Flow Developments", Appl. Sci. Res.... V.10, pp.265-276,1961.

[51 ] DENHAM, M.K. e PATRICK, M.A., "Laminar Flow over a Dowmstream-Facing

Step in a Two-Dimensional Flow Channel", Trans. Inst. Chem. Engns.. 52,

pp.361, 1974.

[52] CAMPION -RENSON, A. e CROCHET, M.J., "On the Stream-Function Vorticity

Finite Element Solutions of Navier-Stokes Equations", Int. J. Num. Meth. Eng..

12, pp.1809-1818, 1978.

[53] LESCHZINER, M.A., "Practical Evaluation of Three Finite Difference Schemes

for the Computation of Steady-State Recirculating Flows", Comp. Meth. Appl.

Eng., 23, pp.293-299, 1980.

[54] ARMALY, B.F., DURST F., PEREIRA, J.C.F. e SCHONUNG, B., "Experimental

e Theoretical Investigation of Backward Facing Step Flow", J. Fluid Mech.. 172,

pp.473-496, 1983.

[55] GOSMAN, A.D. e PUN, W.M., Lecture notes for course entilled: "Calculation of

Recirculating Flow", Heat Transfer Rep. HTS/74/2. Imperial College. London,

1974.

[56] THANGAM, S. e KNIGHT-DOYLE, D., "A Computational Scheme in

Generalized Coordinates for Viscous Incompressible Flows", Computers &

Fluids. Vol.18, n24, pp.317-327,1990.

[57] GARTLING, D.K., "A Test Problem for Outflow Boundary Condition: Flow

Over Backward Facing Step", Int. J. Num. Meth. in Fluids. 11, pp.953-967, 1990.

[58] GARTLING, D.K., "NACHOS II - A Finite Element Computer Program for

Incompressible Flow Problems", SE86-1816 e SE86-1817, Seia National

Laboratories, Albuquerque, NM, 1987.

[59] FIDAP Theoretical Manual, Version 4, 5, Fluids Dynamics International,

Evanston, IL, 1988.

107

[60] ROGERS, S. e KWAK, D., "An Upwind Differencing Scheme for the

Incompressible Navier-Stokes Equations", Appl. Num. Math.. V.8, pp.43-64,

1991.

[61 ] SRINIVASAN, K. e RUBIN, S.G.,"Segmented Multigrid Domain Decomposition

Procedure for Incomprssible Viscous Flows", Int. J. Num. Meth. in Fluids. V.I5,

pp.1333-1355, 1992.

[62] ONUR, H.S. e BAYDAR, E., "Laminar Channel Flow Over a Square Step", Int. J.

Eng. Science. V.30. na 9, pp.1109-1116, 1992.

[63] MATHEMATICA Standard Version 2.2, Champaign, Illinois, 1993.

[64] ROACHE, P. ; "Scaling of high Reynolds number weakly separated channel

flows", Symposium on Numerical e Physical Aspects of Aerodynamic Flows,

pp. 1-8, California State Univ., Long Beach, CA, 1981.

[65] NAPOLITANO, M. e ORLANDI, P. ; "Laminar Flow in a Complex Geometry: a

Comparison", Int. J. Num. Meth. Fluids 5, 667-683,1985.

[66] CLIFFE, K.A., JACKSON CP. e GREENFIELD A.C.; "Finite element Solutions

for Flow in a Symmetric Channel with a Smoth Expansion", AERE - R10608.

[67] MAGI, M. e NAPOLITANO, M.; "Calcolo del Flusso Laminare en un Canale de

Geometrea Complessa"; VI Congresso Nazionale AIMETA. Gênova, 7-9,

Octover, 1982.

[68] NAPOLITANO, M., "Efficient ADI e Spline ADI Methods for the Steady state

Navier-Stokes Equations", Int. J. Num. Meth. Fluids. 4,1101-1115, 1984.

[69] QUARTAPELLE, L. e NAPOLITANO, M.; "A method for Solving the

Factorized Vorticity-stream Function Equations by finite elements", Int. J. Num.

Meth. Fluids. 4, 109-125, 1984.

[70] NAPOLITANO, M., "Efficient Solution of Two-dimensional Steady Separated

Flows", Computers & Fluids. 20, 3, 213-222, 1991.

[71] JORDAN, S.A.; "An Interative Scheme for Numerical Solution of Steady

Incompressible Viscous Flows", Computers & Fluids. 21,4, 503-517, 1992.

108

[72] NAPOLITANO, M. e WALTERS, R.W.; "An Incrementai Block-Line-Gauss-

Seidel Method for the Navier-Stokes Equations", AIAA J..V.24. ns 5 , pp. 770-

776, may, 1986.

[73] VRADIS, G., ZALAK, V. e BENTSON, J.; "Simultaneous Variable Solutions of

the Incompressible Steady Navier-Stokes Equations in General Curvilinear

Coordinate Systems., J. Fluids Eng.. 114, 299-305, 1992.

[74] YAO-HSIN, H.; "Arbitrary Domain Velocity Analysis for the Incompressible

Navier-Stokes Equations", J. Comp. Phys.. 110, 134-149, 1994.

[75] SMITH, R.M. e HUTTON, A.G., "The Numerical Treatment of Advection. A

Comparison of Current Methods", Num. Heat Transfer. 5, 4, 439-462, 1982.

[76] LIMA J.A., "Solução das Equações de Navier-Stokes Médias para Escoamento

Turbulento via Transformação Integral", COPPE/UFRJ, M.Sc, Tese, Rio de

Janeiro, 1995.

[77] PEREIRA, L.M.; "Transformação Integral das Equações de Navier-Stokes em

Coordenadas Cilíndricas: Desenvolvimento Hidrodinâmico em Dutos Anulares",

COPPE/UFRJ, M.Sc, Tese, Rio de Janeiro, 1995.

[78] PÉREZ GUERRERO, J.S. e COTTA R.M.; "Benchmark Integral Transform

Results for Flow over a Backward-Facing Step"; Computers & Fluids, in press.

109

APÊNDICE A

Escoamento completamente Desenvolvido entre Placas Paralelas

Considera-se um duto de placas paralelas, tal como o esquematizado na fig A.l.

Quando o escoamento torna-se completamente desenvolvido, a velocidade

transversal é nula, e as equações de Navier-Stokes então se reduzem a:

dy

com as condições de contorno:

u(-b) = 0 ; u(b) = 0

b

h

y u=0 V = 0

u=0 v = 0

FIGURA A.l

Geometria do escoamento completamente desenvolvido

entre Placas Paralelas

110

A solução de (A.l) fornece a componente longitudinal da velocidade no

escoamento desenvolvido:

Kb'u = — • (A.2)

considerando que a vazão por unidade de comprimento no canal é Q:

Q= Judy- b

(A.3)

Integrando-se esta última expressão, resulta que:

- K b 2 3Q

4b(A.4)

Logo o perfil da componente de velocidade u pode ser expresso convenientemente

em função de Q:

u =3Q

4b1- (A.5)

A função corrente é encontrada a partir de:

u =õy

(A.6)

Logo,

J dy = J uk, - b

dy (A.7)

I l l

onde kj é um valor arbitrário do valor de \\i ( - b) .

Integrando-se a equação (A.7), obtém-se:

v|/=-3Q y iy

bJ 3 lb .(A.8)

O valor da função corrente na parede superior, isto é, em y = b, é dado por:

(A.9)

112

APÊNDICE B

Derivadas da Fórmula de Inversão

ôx

„ 5Y; ,„4 ^

(B.l)

(B.3)

Z l -^ÜT; +3 -̂ -\j7V +3—-üTi" + Y:\jT:'" (B.4)

(B, )

ôyn t f Uy11 J (B.6)

õy:

5 2 Y S , 5Y, „ j+ 2 -viTi +—LvjT i (B.8)

113

Õx2õy2

Õ2Y;(B.9)

114

APÊNDICE C

Coeficientes do Sistema Diferencial Ordinário

Tem-se que:

y o ( x )

p - q (C2)

Portanto:

(C5)7 y p q

onde (n), indica a ordem da derivada.

Além disso:

y = 4 y o + y 3 (C6)

Logo:

— - = ^ y 0 p ( n ) +y 3 P - q ( (c-8)ox"

115

ou:

(C.9)

onde:

W n = y o p ( n ) ; S n = y 3 p < n > - q ( n > (C.10,11)

Desta forma, qualquer derivada de autofunção com respeito a x, poderá ser

representada facilmente como uma derivada em Ç, ao se aplicar a regra da cadeia.

A derivada da autofunção com respeito a y é encontrada mais facilmente como:

(C l 2)õya

L = ( p ) n

Õyn

As expressões dadas acima foram fornecidas ao MATHEMATICA, com o qual se

gerou a seguinte representação dos coeficientes:

F00ijk = 2*CeO*p**2*pl - B2eO*p2*sl - 2*CeO*pl*sl**2 + 2*CeO*p*sl*s2 +B2eO*p*s3 -B2el*p2*wl -4*Cel*pl*sl*wl +2*Cel*p*s2*wl -2*Ce2*pl*wl**2 + 2*Cel*p*sl*w2 + 2*Ce2*p*wl*w2 + B2el*p*w3

(Cl 5)

116

F01ijk = CeO*p**3 - 2*B2eO*pl*sl + CeO*p*sl**2 + 3*B2eO*p*s2 - 2*B2el*pl*wl +2*Cel*p*sl*wl + Ce2*p*wl**2 + 3*B2el*p*w2

(C.I 6)

F02ijk = 2*B2eO*p*sl + 2*B2el*p*wl(C.17)

F10ijk = -(AeO*p**3) - BeO*p2 - 2*EeO*pl*sl - AeO*p*sl**2 - EeO*p*s2 -2*Eel*pl*wl - 2*Ael*p*sl*wl - Ae2*p*wl**2 - Eel*p*w2

(C.18)

Fll i j k = 2*Be0*pl -2*EeO*p*sl -2*Eel*p*wl(C.19)

FOy = 2*DeO*fq*p**2*pl - 2*De2*fq*p**2*pl - 4*Dlel*fq*p**2*pl -2*DleO*fq*p2*sl + 2*Dle2*fq*p2*sl - 2*DeO*fq*pl*sl**2 +

2*De2*fq*p*sl*s2 - 4*Dlel*fq*p*sl*s2 + 2*DleO*fq*p*s3 -2*Dle2*fq*p*s3 - 2*Dlel*fq*p2*wl + 2*Dle3*fq*p2*wl -4*Del*fq*pl*sl*wl +4*De3*fq*pl*sl*wl + 8*Dle2*fq*pl*sl*wl +2*Del*fq*p!l:s2*wl - 2*De3*fq*p*s2*wl - 4*Dle2*fq*p*s2*wl -2*De2*fq*pl*wl**2 + 2*De4*fq*pl*wl**2 + 4*Dle3*fq*pli|twl**2 +2*Del*fq*p*sl*w2 - 2*De3*fq*p*sl*w2 - 4*Dle2*fq*p*sl*w2 +2*De2*fq*p*wl*w2 - 2*De4*fq*p*wl*w2 - 4*Dle3*fq*p*wl*w2 +2*Dlel*fq*p*w3 - 2*Dle3*fq*p*w3

(C.20)

Fly = 2*CroeO*fq*p**3 + DeO*fq*p**3 - De2*fq*p**3 - CroeO*fq*p2 +Croe2*fq*p2 + 4*Croel*fq*pl*sl - 2*DleO*fq*pl*sl +2*Dle2*fq*pl*sl + 2*CroeO*fq*p*sl**2 + DeO*fq*p*sl**2 -De2*fq*p*sl**2 + 2*Croel*fq*p*s2 + 3*DleO*fq*p*s2 -3*Dle2*fq*p*s2 + 4*Croe2*fq*pl*wl - 2*Dlel*fq*pl*wl +2*Dle3*fq*pl*wl + 4*Croel*fq*p*sl*wl +2*Del*fq*p*sl*wl -2*De3*fq*p*sl*wl + 2*Croe2*fq*p*wl:|!*2 + De2*fq*p*wl**2 -De4*fq*p*wl**2 + 25|!Croe2*fq*p*w2 + 3*Dlel*fq*p*w2 -3*Dle3*fq*p*w2

(C.21)

F2y = 2*DleO*fq*p*sl - 2*Dle2*fq*p*sl + 2*Dlel'*fq*p*wl " 2*Dle3*fq*p*wl(C.22)

F3jj = CroeO*fq*p - Croe2*fq*p(C.23)

117

j = -4*Autel*fq**2*p**2*pl + 4*Aute3*fq**2*p**2*pl - AuteO*fq**2*p2*sl +2*Aute2*fq**2*p2*sl - Aute4*fq**2*p2*sl +4*Autel*fq**2*pl*sl**2 - 4*Aute3*fq**2*pl*sl**2 -4*Autel*fq**2*p*sl*s2 + 4*Aute3*fq**2*p*sl*s2 +AuteO*fq**2*p*s3 - 2*Aute2*fq**2*p*s3 + Aute4*fq**2*p*s3 -Autel*fq**2*p2*wl + 2*Aute3*fq**2*p2*wl -Aute5*fq**2*p2*wl + 8*Aute2*fq**2*pl*sl*wl -8*Aute4*fq**2*pl*sl*wl - 4*Aute2*fq**2*p*s2*wl +4*Aute4*fq**2*p*s2*wl + 4*Aute3*fq**2*pl*wl**2 -4*Aute5*fq**2*pl*wl**2 - 4*Aute2*fq**2*p*sl*w2 +4*Aute4*fq**2*p*sl*w2 - 4*Aute3*fq**2*p*wl*w2 +4*Aute5*fq**2*p*wl*w2 + Autel*fq**2*p*w3 -2*Aute3*fq**2*p*w3 + Aute5*fq**2*p*w3

(C.24)

= -8*Autel*fq*pl**2 - 8*Autel*fq*p*p2 - 16*AuteO*fq*p*pl*sl -4*AuteO*fq*p**2*s2 - 12*AuteO*fq*sl**2*s2 -6*Autel*fq*s2**2 - 8*Autel*fq*sl*s3 + AuteO*fq*s4 -Aute2*fq*s4- 16*Autel*fq*p*pl*wl - 24*Autel*fq*sl*s2*wl -8*Aute2*fq*s3*wl - 12*Aute2*fq*s2*wl**2-4*Autel*fq*p**2*w2 - 12*Autel*fq*sl**2*w2 -12*Aute2*fq5|;s2*w2 - 24*Aute2*fq*sl*wl*w2 -12*Aute3*fq*wl**2*w2 - 6*Aute3*fq*w2**2 - 8*Aute2*fq*sl*w3 -8*Aute3*fq*wl*w3 + Autel*fq*w4 - Aute3*fq*w4

(C.25)

EOjj = 4*DeO*pl**2 + 4*DeO*p*p2 + 8*D2eO*p*pl*sl + 2*D2eO*p**2*s2 +6*D2eO*sl**2*s2 + 3*DeO*s2**2 + 4*060*81*83 + DleO*s4 +2*CroeO*p:|i*2*sl**2*uj + CroeO*sl**4*uj + 8*D2el*p*pl*wl +12*D2el*sl*s2*wl +4*Del*s3*wl +4*Croel*p**2*sl*uj*wl +4*Croel*sl**3Hsuj*wl + 6*D2e2*s2*wl**2 +2*Croe2*p**2*uj*wl**2 + 6*Croe2*sl**2*uj*wl**2 +4*Croe3*sl*uj*wl**3 + Croe4*uj*wl**4 + 2*D2el*p**2*w2 +6*D2el*sl**2*w2 + 6*Del*s2*w2 + 12*D2e2*sl*wl*w2 +6*D2e3*wl**2*w2 + 3*De2*w2**2 + 4*Del*sl*w3 + 4*De2*wl*w3 +Dlel*w4

(C.26)

Ely - 8*DeO*p*pl + 4*D2eO*p**2*sl + 4*D2eO*sl**3 + 12*DeO*sl*s2 +4*DleO*s3 + 4*D2el*p**2*wl + 12*D2el*sl**2*wl +12*Del*s2*wl + 12*D2e2*sl*wl**2 + 4*D2e3*wl**3 +12*Del*sl*w2 + 12*De2*wl*w2 + 4*Dlel*w3

(C.27)

E2y = 2*DeO*p**2 + 6*060*81**2 + 6*DleO*s2 + 12*Del*sl*wl + 6*De2*wl**2 +6*Dlel*w2

(C.28)

118

E3ÍJ = 4*DleO*sl + 4*Dlel*wl(C.29)

As integrais núcleo são definidas como:

AeO= J YjYjYk'dC (C.30)

= J ^YjYjY^d^ ((Ael-1

-r= J YJYJY;

- i

'= J YiY;Yk

, _ Í ' Y Y Y

- f

A P 9 — I ? Y Y Y d? (C 3?"»4 ~ ' J K ~ v

-1

BeO= I YjYjYkd^ (C.33)

B2eO= J YjYjYkd^ (C.34)

f . .B2el= ^Y:Y:Yk d£ (C.35)

-1

CeO= J YiYjYk dÇ (C.36)- i

Cel= ^Y:Y:Ykd^ (C.37)/ J hv

_ I

Ce2= I l2 Y:Y|Yk ' d£ (C.38)

119

•= J YÍYJY:EeO= I Y i Y j Y ^ (C.39)- i

. .J 1Eel=J ^ Y j Y , ^ (C.40)

-1

•1

DeO= J YjYj" d^ (C.41)

-rDel= J ^Y;Yj d£, (C.42)- i

-fDe2= V Y:Yi d£ (C.43)-1

- fDe3= t; Y:Y: d^ (C.44)-1

•ÍDe4= J ^4 YiYj d§ (C.45)

DleO= J YjY- (C.46)

Dlel= J ^YjYjd^ (C.47)

Dle2= t2 Y:Yi d£ (C.48)-1

Dle3= ^3Y:Y|d4 (C.49)

120

'= Í YiYj"D2eO= I YiYj dÇ (C.50)

D2el= J ÇYjY^dÇ (C.51)- í

, fD2e2= J r YjYj d^ (C.52)- i

-fD2e3= J Ç3 YiYj-'dC (C.53)- i

D2e4= J Ç4 YiYJ"d^ (C.54)

-fv.y,CroeO = J YjYjdÇ (C.55)- i

•fCroel= ^Y:Y=dE (C.56)- i

2

- l

Croe2= J \ l Y ^ d Ç (C.57)

• i

Croe3= J ^ YjYjdÇ (C.58)- i

Croe4= J ^4 YjYjd^ (C.59)

121

AuteO=|Y id^ (C.60)

•A

(C.62)

Aute3= I Ç3 Yj d£ (C.63)- i

Aute4= J Ç4 Y; dÇ (C.64)- i

(C.65)

A família de coeficientes "Aen" e "Een" pode ser calculada ao se fazer integração

por pontos, resultando uma combinação de outros já avaliados.

Aen i j k = J Ç - Y . Y j Y , ; ' ^- i

- J [n^"1 YiYjY,: +^n YJYjY; +^n Y, Y^Y^- i

(C.66)

Portanto:

i j k = - jApO _ _ ??Y Y Y + ? 2 Y Y Y + ? z Y Y Y H? (C 671-i

Ae2 i jk = - ( 2 E e l i j k +Ce2 j i k +Ce2 i j k ) (C.68)

122

Ael i jk =-(EeO i j k +Cel j i k +Cel i j k ) (C.69)

AeOijk =- (CeO j i k +CeO i j k) (C.70)

No outro caso:

Een i j k = J Ç" Y,YjY;d4 =

- J [n^"-1 Y.Yj Yk + Ç" YiYjY, + £" Y, Y^Y ]

(C.71)

Eel i jk =-(BeO i j k +B2e l j i k +B2el i j k ) (C.72)

EeOijk = - (B2eO j i k +B2eO i jk ) (C.73)