UNIVERSIDADE TECNOLÓGICA FEDERAL DO PARANÁ
CAMPUS CURITIBA
PROGRAMA DE PÓS-GRADUAÇÃO EM ENGENHARIA MECÂNICA E
DE MATERIAIS - PPGEM
WILLIAN SEGALA
SIMULAÇÃO NUMÉRICA DO ESCOAMENTO MONOFÁSICO
NO PRIMEIRO ESTÁGIO DE UMA BOMBA CENTRÍFUGA DE
DUPLO ESTÁGIO
CURITIBA
MARÇO – 2010
WILLIAN SEGALA
SIMULAÇÃO NUMÉRICA DO ESCOAMENTO MONOFÁSICO
NO PRIMEIRO ESTÁGIO DE UMA BOMBA CENTRÍFUGA DE
DUPLO ESTÁGIO
Dissertação apresentada como requisito parcial à
obtenção do título de Mestre em Engenharia, do
Programa de Pós-Graduação em
Engenharia
Mecânica e de Materiais, Área de Ciências
Térmicas, do Departamento de Pesquisa e PósGraduação, do Campus de Curitiba, da UTFPR.
Orientador: Prof. Rigoberto E. M. Morales, Dr. Eng.
Co-Orientador: Prof. Cezar O. R. Negrão, PhD.
CURITIBA
MARÇO – 2010
TERMO DE APROVAÇÃO
WILLIAN SEGALA
SIMULAÇÃO NUMÉRICA DO ESCOAMENTO MONOFÁSICO NO PRIMEIRO
ESTÁGIO DE UMA BOMBA CENTRÍFUGA DE DUPLO ESTÁGIO
Esta Dissertação foi julgada para a obtenção do título de mestre em engenharia,
área de ciências Térmicas, e aprovada em sua forma final pelo Programa de Pósgraduação em Engenharia Mecânica e de Materiais.
__________________________________________
Prof. Giuseppe Pintaúde, D.Sc.
Coordenador do Programa
Banca Examinadora
_____________________________
Prof. Rigoberto E. M. Morales, Dr. Eng.
PPGEM/UTFPR-Orientador
___________________________
Valdir Estevam, Dr, Eng.
E&P / PETROBRAS
_____________________________
Prof. Antônio Carlos Bannwart, Dr. Eng.
DEP/FEM/UNICAMP
___________________________
Admilson Teixeira Franco, Dr. Eng.
PPGEM / UTFPR
Curitiba, 26 de Março de 2010
iv
AGRADECIMENTOS
Agradeço em primeiro lugar a Deus, pelos desafios que Ele me proporciona e
também pelas ferramentas oferecidas para a solução dos mesmos.
Agradeço à Universidade Tecnológica Federal do Paraná, ao Laboratório de
Ciências Térmicas (LACIT) e ao Programa de Pós-Graduação em Engenharia de
Mecânica e de Materiais (PPGEM) pela possibilidade de realização desse trabalho.
Agradeço à Petrobrás e ANP pelo suporte técnico e financeiro para o
desenvolvimento do tema.
Agradeço aos professores Rigoberto E. M. Morales e Cezar O. R. Negrão,
pela oportunidade, apoio, incentivo, troca de conhecimentos no enriquecimento do
trabalho e amizade durante todo o período de convívio.
Agradeço aos amigos, de maneira especial aos colegas de mestrado Hendy,
Henrique, Leandro que de alguma maneira contribuíram com esse trabalho.
Dedico esse trabalho à minha mãe, Lilia Segala, pelo incentivo e atenção
incondicionais, à minha noiva Tatiane, pelo apoio, carinho e compreensão durante
esse período, a minha irmã Debora e ao meu pai, Osias Segala (In Memoriam).
v
SEGALA, Willian. Simulação Numérica do Escoamento Monofásico no Primeiro
Estágio de uma Bomba Centrífuga de Duplo Estágio, 2010. Dissertação (Mestrado
em Engenharia) - Programa de Pós-graduação em Engenharia Mecânica e de
Materiais, Universidade Tecnológica Federal do Paraná, Curitiba, 104p.
RESUMO
Bombas centrífugas são equipamentos muito utilizados em plantas industriais.
Porém, os esforços para a elaboração de curvas de desempenho confiáveis (altura
de elevação e eficiência versus vazão) ainda são requisitados, principalmente
quando o fluido apresenta variação de suas características, como o petróleo, por
exemplo. Atualmente é comum utilizar fatores de correção para prever o
comportamento de uma bomba centrífuga quando submetida a fluidos distintos do
fluido padrão (água). Existem ainda controvérsias do uso desses fatores de
correção, por exemplo, para condições operacionais diferentes das de projeto. Por
isso a utilização de dinâmica dos fluidos computacional pode ser uma ferramenta
muito útil para melhor entender o escoamento no interior de bombas centrífugas.
Imaginando esse cenário, o presente trabalho apresenta resultados de uma
simulação numérica do escoamento monofásico no interior do primeiro estágio de
uma bomba centrífuga comercial de duplo estágio operando com água. A bomba
centrífuga estudada é similar a uma bomba centrífuga radial utilizada em operações
de extração de petróleo. Foram escolhidas quatro rotações para avaliação do
escoamento: 1150rpm (rotação nominal), 1000rpm, 806rpm e 612 rpm. Foi utilizado
um modelo de turbulência e o escoamento foi considerado em regime transiente. Os
resultados
numéricos obtidos foram
comparados a
valores
experimentais,
apresentando boa concordância.
Palvras-chave: Bomba centrífuga, Dinâmica dos fluidos Computacional, Testes
Experimentais.
vi
SEGALA, Willian. Simulação Numérica do Escoamento Monofásico no Primeiro
Estágio em uma Bomba Centrífuga de Duplo Estágio, 2010. Dissertação (Mestrado
em Engenharia) - Programa de Pós-graduação em Engenharia Mecânica e de
Materiais, Universidade Tecnológica Federal do Paraná, Curitiba, 104p.
ABSTRACT
Centrifugal pumps are among the most used equipment in industrial plants.
Besides that, efforts to infer the so-called pump performance curves are still required,
mainly when the pump is used to transfer fluids that have a history of properties
change, as the viscosity or density and chemical composition. Nowadays, it is
common to use correction factor to prevent the centrifugal pump behavior when
delivering different fluids. There are controversies between the uses of theses
correction factors, for example, to the pump off-design conditions. Thus,
computational fluid dynamics (CFD) modeling may become a reasonable alternative
to understand the flow behavior inside centrifugal pumps. Focusing on these aspects
this work presents numerical results of the flow inside the first stage of a commercial
double stage centrifugal pump delivering water. This pump is similar to a radial
Electrical submersible pump (ESP) commonly used in oil wells. It was select four
different impeller angular velocities: 1150rpm (nominal), 1000rpm, 806rpm and
612rpm. It was used a turbulence model and the flow was considered in transient
regimen. The numerical results was compared to experimental data and shown good
agreement.
Keywords: Centrifugal pump, Computicional fluid dynamics, experimental tests.
vii
LISTA DE FIGURAS
Figura 1.1 – Exemplo de ábacos fornecido por fabricantes de bombas centrífugas.
(fonte: KSB Bombas)............................................................................................2
Figura 1.2 – Curvas de desempenho de uma bomba centrífuga convencional. O
ponto sobre a curva de eficiência é denominado ponto de eficiência máxima ou
ponto de projeto ...................................................................................................2
Figura 1.3 – Bomba centrífuga IMBIL ITAP 65-330/2. (b) Fotografia do primeiro
estágio da bomba centrífuga estudada. ...............................................................4
Figura 1.4 – Esquema do modelo numérico implementado no programa
computacional Ansys CFX 11.0 para a simulação do escoamento no interior do
primeiro estágio da bomba centrífuga Imbil ITAP 65-330/2. ................................5
Figura 2.1 – Desempenho de uma bomba centrífuga ao bombear fluidos de
diferentes viscosidades (Li, 1999)........................................................................9
Figura 2.2 – Interação entre difusor aletado e rotor, no escoamento de uma bomba
centrífuga. (Feng et al. (2007))...........................................................................12
Figura 2.3 – Estudo da fluidodinâmica do escoamento em uma bomba centrífuga
apresentado por Cheah et al. (2007). (a) Vazão na condição de projeto e (b)
vazão acima da condição de projeto..................................................................13
Figura 2.4 – Volume de controle e componentes de velocidades do escoamento na
entrada e saída de um rotor genérico ................................................................13
Figura 2.5 – Geometria e notação usadas para desenvolver diagramas de
velocidade em bombas centrífugas. (a) Distribuição dos perfis de velocidade ao
longo da pá do rotor e (b) os triângulos de velocidades na entrada e saída do
rotor....................................................................................................................15
Figura 2.6 – Componentes (tangencial e radial) da velocidade absoluta na secção de
saída do rotor da bomba ....................................................................................17
Figura 2.7 – Corte axial do rotor de uma bomba centrífuga ......................................18
Figura 2.8 – Curva característica idealizada de uma bomba centrífuga....................19
viii
Figura 2.9 – Rotores estudados numericamente por Benra (2001). (a) Pás com bordo
de entrada horizontal e (b) pás com bordo de entrada inclinadas em relação à
entrada do rotor..................................................................................................20
Figura 3.1 – Vários sub-domínios co-existentes em um programa comercial de
simulação numérica de escoamentos envolvendo máquinas de fluxo...............23
Figura 3.2 – Sistema de coordenadas rotativo aplicado a uma bomba centrífuga ....24
Figura 3.3 – Propriedade genérica de um escoamento em regime turbulento em
função do tempo, t .............................................................................................26
Figura 3.4 – Esquema das condições de contorno impostas aos sub-domínios do
primeiro estágio da bomba estudada. ................................................................31
Figura 4.1 – Rotor do primeiro estágio da bomba IMBIL ITAP 65-330/2. (a) Rotor na
sua forma original com coroa e (b) rotor sem a coroa superior para facilitar a
visualização de partes internas. .........................................................................33
Figura 4.2 – Rotor real sem a tampa inferior (cubo). Rotor utilizado na determinação
da curvatura das pás e confecção do modelo virtual. ........................................34
Figura 4.3 – Domínio fluido de cálculo do rotor .........................................................34
Figura 4.4 – Detalhe das pás do rotor e aletas do difusor. Tanto as pás do rotor
quanto as aletas do difusor não tocam a superfície de contato entre os domínios
...........................................................................................................................35
Figura 4.5 – Difusor do primeiro estágio da bomba centrífuga Imbil ITAP 65-330/2 .36
Figura 4.6 – Modelo virtual do difusor da bomba centrifuga Imbil ITAP 65 330/2 .....36
Figura 4.7 – Domínio fluido do difusor da bomba centrífuga Imbil ITAP 65-330/2 ....37
Figura 4.8 – Extensão do domínio à saída do difusor. ..............................................38
Figura 4.9 – Tipos de malhas computacionais aplicadas em mecânica dos fluidos e
transferência de calor e massa. (a) Malhas Cartesianas, (b) Malhas CilíndricoPolares, (c) Malhas Ajustadas ao Corpo e (d) Malha não-estruturada...............40
ix
Figura 4.10 – Tipos de elementos finitos utilizados na confecção de uma malha
computacional não-estruturada. (a) Tetraedro, (b) Hexaedro, (c) prisma
triangular e (d) pirâmide. ....................................................................................41
Figura 4.11 – Volume de controle de uma malha não-estruturada............................41
Figura 4.12 – Elemento de malha .............................................................................43
Figura 4.13 – Exemplo de funções de forma utilizadas em elementos tetraédricos
para ponderação de valores internos.................................................................45
Figura 4.14 – Fluxograma simplificado de resolução do sistema de equações do
programa Ansys CFX 11.0 .................................................................................47
Figura 5.1 – Configuração 01 (“Radial”) de escoamento entre rotor e difusor da
bomba centrífuga IMBIL, modelo ITAP 65-330/2 ...............................................49
Figura 5.2 – Configuração 02 (“Axial para baixo”) de escoamento entre rotor e difusor
da bomba centrífuga IMBIL, modelo ITAP 65-330/2 ..........................................49
Figura 5.3 – Configuração 03 (“Extendido”) de escoamento entre rotor e difusor da
bomba centrífuga IMBIL, modelo ITAP 65-330/2 ...............................................50
Figura 5.4 – Configuração 4 (“Indutor sem aletas”) de escoamento entre rotor e
difusor da bomba centrífuga IMBIL, modelo ITAP 65-330/2...............................51
Figura 5.5 – Modelo numérico da bomba centrífuga Imbil ITAP 65-330/2 ................54
Figura 5.6 – Perfil de velocidade relativa em um dos canais da bomba centrífuga
estudada. Comparativo de diferentes critérios de parada ..................................56
Figura 5.7 – Variação de pressão na direção tangencial, no raio de saída do rotor.
Comparativo de diferentes critérios de parada...................................................57
Figura 5.8 – Corte transversal da malha computacional utilizada no rotor e difusor.
Detalhe ampliado mostra a distribuição dos prismas ao longo das pás do rotor e
aletas do difusor.................................................................................................58
Figura 5.9 – (a) Malha no indutor. (b) Malha no tubo de entrada ..............................59
Figura 5.10 – Perfil de velocidade relativa versus direção tangencial, tomado no raio
R=40,0mm e z=6,0mm (metade da largura do canal do rotor)...........................60
x
Figura 5.11 – Perfil de pressão versus direção tangencial, tomado no raio R=40,0mm
e z=6,0mm (metade da altura do canal do rotor) ...............................................61
Figura 5.12 – Perfil de velocidade relativa do escoamento no rotor versus a direção
axial, em R=64,41mm e entre duas pás consecutivas do rotor..........................62
Figura 5.13 – Perfil de pressão no difusor versus direção tangencial, tomado no raio
R=110,0mm e z=6,0mm.....................................................................................63
Figura 5.14 - Teste de malha temporal (Parte 1).......................................................66
Figura 5.15 – Teste de malha Temporal (Parte 2).....................................................67
Figura 5.16 – Pressão versus tempo na entrada do difusor do primeiro estágio da
bomba centrífuga Imbil ITAP 65-330/2, na vazão de 36m3/h e rotação de
1150rpm. ............................................................................................................69
Figura 6.1 – Ganho de pressão no primeiro rotor da bomba centrífuga Imbil ITAP 65330/2. Dados experimentais obtidos por Amaral (2007) ....................................71
Figura 6.2 – Ganho de pressão no difusor da bomba centrífuga Imbil ITAP 65-330/2.
Dados experimentais obtidos por Amaral (2007) ...............................................72
Figura 6.3 – Altura de elevação do primeiro estágio da bomba centrífuga Imbil ITAP
65-330/2. Dados experimentais obtidos por Amaral (2007) ...............................73
Figura 6.4 – Altura de elevação total versus vazão volumétrica da bomba centrífuga
IMBIL ITAP 65-330/2, na rotação de 1150 rpm. Dados fornecidos no catálogo do
fabricante (Fonte: site Imbil)...............................................................................74
Figura 6.5 – Aplicação das equações de similaridade às curvas de 1000, 806 e 612
rpm transpondo-as para a rotação de 1150 rpm. ...............................................75
Figura 6.6 – Torque numérico versus vazão no primeiro rotor da bomba centrífuga
Imbil ITAP 65-330/2, para uma rotação de 1150 rpm ........................................79
Figura 6.7 – Potência de entrada numérica versus vazão no primeiro rotor da bomba
Imbil ITAP 65-330/2, para um rotação de 1150 rpm ..........................................80
Figura 6.8 – Potência de saída numérica versus vazão no primeiro estágio da bomba
centrífuga Imbil ITAP 65-330/2, para rotação de 1150 rpm................................81
xi
Figura 6.9 – Eficiência hidráulica numérica do primeiro estágio da bomba centrífuga
Imbil ITAP 65-330/2 e Eficiência global fornecida pelo fabricante......................82
Figura 6.10 – Instante de tempo no qual a pá do rotor e a aleta do difusor estão
alinhadas. Definição das Faces de Pressão (FP) e Faces de Sucção (FS) .......83
Figura 6.11 – Pressão versus direção tangencial no raio de saída do rotor para o
caso de pá do rotor alinhada com a aleta do difusor..........................................84
Figura 6.12 – Periodicidade da pressão com a direção tangencial na saída do
primeiro rotor da bomba centrífuga Imbil ITAP 65-330/2. Vazão 35m3/h e rotação
de 1150rpm. .......................................................................................................84
Figura 6.13 – Pressão versus tempo em um ponto localizado no bordo de entrada da
aleta do difusor da bomba centrífuga Imbil ITAP 65-330/2, para a rotação de
1150 rpm e vazão de 35m3/h. ............................................................................86
Figura 6.14 – Localização do bordo de entrada e saída do rotor da bomba. ............87
Figura 6.15 – Distribuição de pressão na pá do rotor desde a entrada até a saída na
altura da raiz da pá (ou cubo) ............................................................................87
Figura 6.16 - Distribuição de pressão na pá do rotor desde a entrada até a saída na
altura da coroa, para o rotor girando a 1150 rpm e na vazão de 35m3/h ...........88
Figura 6.17 – Campo de pressão no interior do primeiro estágio da bomba centrífuga
Imbil ITAP 65-330/2, vazão de 20m3/h e rotação de 1150 rpm..........................89
Figura 6.18 – Campo de pressão no interior do primeiro estágio da bomba centrífuga
Imbil ITAP 65-330/2, vazão de 35m3/h e rotação de 1150 rpm..........................90
Figura 6.19 – Campo de pressão no interior do primeiro estágio da bomba centrífuga
Imbil ITAP 65-330/2, vazão de 50m3/h e rotação de 1150 rpm..........................90
Figura 6.20 – Planos de corte da visualização do escoamento, obtido
numericamente, no rotor e difusor do primeiro estágio da bomba centrífuga Imbil
ITAP 65-330/2 ....................................................................................................92
Figura 6.21 - Esquema da saída do rotor e entrada do difusor que mostra a diferença
de áreas entre esses dois componentes............................................................94
xii
Figura 6.22 – Linhas de corrente do escoamento no interior da bomba centrífuga
Imbil ITAP 65 330/2, na rotação nominal e vazão de 50m3/h. (a) escoamento no
plano T1, (b) escoamento no plano T2 e (c) escoamento no plano T3. .............95
Figura 6.23 - Linhas de corrente do escoamento no interior da bomba centrífuga
Imbil ITAP 65 330/2, na rotação nominal e vazão de 44m3/h. (a) escoamento no
plano T1, (b) escoamento no plano T2 e (c) escoamento no plano T3. .............95
Figura 6.24 - Linhas de corrente do escoamento no interior da bomba centrífuga
Imbil ITAP 65 330/2, na rotação nominal e vazão de 35m3/h. (a) escoamento no
plano T1, (b) escoamento no plano T2 e (c) escoamento no plano T3. .............96
Figura 6.25 – Linhas de corrente do escoamento no interior da bomba centrífuga
Imbil ITAP 65 330/2, na rotação nominal e vazão de 20m3/h. (a) escoamento no
plano T1, (b) escoamento no plano T2 e (c) escoamento no plano T3. .............96
Figura 6.26 - Linhas de corrente do escoamento no interior do difusor bomba
centrífuga Imbil ITAP 65 330/2, na rotação nominal e vazão de 20m3/h. (a)
Vazão 20 m3/h, (b) Vazão 35 m3/h (c) Vazão 44 m3/h e (d) Vazão 50 m3/h...97
xiii
LISTA DE SÍMBOLOS
Símb.
Descrição
Unidade (S.I.)
α
Ângulo entre vetor de velocidade absoluta do escoamento
e a direção tangencial
β
rad
Ângulo entre o vetor de velocidade relativa do escoamento
no rotor e a direção tangencial (medido no sentido contrário
ao giro do rotor)
ε
rad
Componente isotrópica da taxa de dissipação de energia
turbulenta
m 2 ⋅ s −3
ηh
Eficiência hidráulica da bomba centrífuga
-
φ
Propriedade genérica
-
φ
Média temporal de uma propriedade genérica
-
φ'
Flutuação temporal de uma propriedade genérica
-
ϕ
Variável angular
κ
Energia cinética turbulenta
m 2 ⋅ s −2
κv
Constante de Von Karman
-
µ
Viscosidade dinâmica
kg ⋅ m-1 ⋅ s−1
µt
Viscosidade turbulenta
kg ⋅ m-1 ⋅ s−1
ν
Viscosidade cinemática
m2 ⋅ s-1
νt
Viscosidade cinemática turbulenta
m2 ⋅ s-1
θ
Variável angular
ρ
massa específica do fluido
rad
rad
kg ⋅ m −3
xiv
σε
Número de Prandtl difusivo para a taxa de dissipação
-
σk
Número de Prandtl difusivo para a energia cinética turbulenta
-
τw
Tensão de cisalhamento sobre a parede
ω
Velocidade angular do rotor
kg ⋅ m-1 ⋅ s −2
rad ⋅ s −1
a, b, c, d, e Constantes oriundas da determinação do ajuste de curvas
referente aos resultados numéricos obtidos
-
Cε1,Cε 2
Constantes de fechamento do modelo de turbulência κ − ε
-
Cµ
Constante de fechamento associada à viscosidade turbulenta
-
Dε
Termo difusivo da equação de transporte de ε
m 2 ⋅ s −4
Dκ
Termo difusivo da equação de transporte de κ
m 2 ⋅ s −3
FS
Face de sucção da pá do Rotor ou da aleta do Difusor
FP
Face de pressão da pá do Rotor ou da aleta do Difusor
g
Vetor aceleração da gravidade
H
Altura de elevação
m
p
Pressão termodinâmica
Pa
p
Pressão com aplicação das médias de Reynolds
Pk
Produção de energia cinética turbulenta ( κ )
m 2 ⋅ s −3
Q
Vazão volumétrica
m 3 ⋅ s −1
r
Vetor posição de uma partícula fluida em relação ao
m ⋅ s −2
kg ⋅ m-1 ⋅ s −2
Sistema de coordenadas não-inercial
m
R2
Quadrado do coeficiente de correlação
-
S
Termo fonte
-
t
Variável tempo
s
xv
T
Torque
N.m
u'
Flutuação de velocidade na direção x
m ⋅ s-1
u
Velocidade média de Reynolds na direção x
m ⋅ s-1
V
Vetor de velocidade
m ⋅ s −1
v'
Flutuação de velocidade y
m ⋅ s-1
v
Velocidade média de Reynolds na direção y
m ⋅ s-1
w'
Flutuação de velocidade na direção z
m ⋅ s-1
w
Velocidade média de Reynolds na direção z
m ⋅ s-1
Wɺ
Potência
Watt
RESRSM
Valor do resíduo numérico médio de uma variável do escoamento, calculada pelo programa computacional
RESMAX
-
Valor do resíduo numérico máximo de uma variável do escoamento, calculada pelo programa computacional
-
x,y,z
Direções coordenadas de um sistema não-inercial
m
X,Y,Z
Direções coordenadas de um sistema inercial
m
ξ1, ξ2 e ξ3
Direções coordenadas de um sistema genérico
-
Operadores
∆
Variação de uma propriedade do escoamento
∇
Operador nabla
∂
Operador diferencial parcial
∫
Operador integral
xvi
SUMÁRIO
Agradecimentos ......................................................................................................... iv
Resumo ....................................................................................................................... v
Abstract ...................................................................................................................... vi
Lista de Figuras......................................................................................................... vii
Lista de Símbolos......................................................................................................xiii
Sumário .................................................................................................................... xvi
1
2
Introdução ............................................................................................................1
1.1
Objetivos .......................................................................................................4
1.2
Justificativas..................................................................................................5
1.3
Estrutura do trabalho.....................................................................................6
Revisão Bibliográfica............................................................................................8
2.1
3
Teoria de máquinas de fluxo aplicado a bombas centrífugas .....................13
Modelagem Matemática .....................................................................................23
3.1
Equações de conservação da massa e da quantidade de movimento
escritas, para um sistema de coordenadas rotativo...............................................24
3.2
4
Modelagem da turbulência..........................................................................25
3.2.1
Modelo de Turbulência a duas equações κ − ε padrão .......................29
3.2.2
Condições de contorno do problema ...................................................30
Modelagem Geométrica e Numérica..................................................................32
4.1
Definição das Geometrias do primeiro estágio da bomba centrífuga..........32
4.1.1
Geometria do rotor do primeiro estágio da bomba centrífuga..............32
4.1.2
Geometria do difusor do primeiro estágio da bomba centrífuga ..........35
4.1.3
Geometrias auxiliares utilizadas nos estudos numéricos.....................37
4.2
Discretização das equações de conservação .............................................38
xvii
5
Avaliação de Parâmetros e Testes de Malha.....................................................48
5.1
A escolha da configuração geométrica adequada para as simulações
numéricas ..............................................................................................................48
6
5.2
Critério de parada .......................................................................................55
5.3
Teste de Malha ...........................................................................................57
5.4
Teste de Passo de Tempo ..........................................................................64
5.5
Influência da condição inicial ......................................................................68
Resultados .........................................................................................................70
6.1
Altura de elevação versus vazão volumétrica .............................................71
6.2
Obtenção de uma curva unificada de altura de elevação para qualquer
vazão e rotação .....................................................................................................73
6.3
Potência e Eficiência Hidráulica ..................................................................78
6.4
Propriedades do escoamento no interior da bomba centrífuga estudada ...83
7
conclusões e sugestões de trabalhos futuros ....................................................99
8
Referências Bibliográficas................................................................................102
1
1 INTRODUÇÃO
A bomba centrífuga é um dos equipamentos mais utilizados no mundo e é
usada nas mais diversas áreas como, por exemplo, na irrigação, no abastecimento
de água, na indústria petroquímica, no setor automobilístico, entre outros. O princípio
de funcionamento de uma bomba centrífuga consiste em transferir a energia cinética
proveniente de um rotor ao líquido bombeado, fazendo com que o fluido escoe no
interior de um duto.
Devido a sua importância, segundo a literatura, pesquisas envolvendo
bombas centrífugas foram intensificadas a partir da primeira metade do século XX.
Os primeiros estudos foram desenvolvidos utilizando técnicas experimentais e
modelos teóricos simplificados, com destaque aos desenvolvidos por Alexey Joakim
Stepanoff (Amaral, 2007).
As principais variáveis utilizadas para o dimensionamento e seleção de
ɺ , e
bombas centrífugas são a altura de elevação, H, potência consumida, W
ent
eficiência, η , sendo essas variáveis mensuradas para cada vazão de trabalho
imposta. A altura de elevação quantifica a capacidade da bomba centrífuga de
elevar a pressão do fluido bombeado, desde a entrada até a saída da bomba. A
potência e a eficiência contabilizam a quantidade de energia necessária para o
funcionamento da bomba centrífuga e a parcela dessa energia efetivamente
transferida ao líquido bombeado, respectivamente.
Catálogos de fabricantes de bombas centrífugas apresentam ábacos
semelhantes aos mostrados na Figura 1.1. No eixo das abscissas tem-se a vazão
volumétrica e no eixo das ordenadas, a altura de elevação. As regiões delimitadas
do ábaco correspondem ao campo de aplicação de cada modelo de bomba
produzido pelo fabricante para 1750 rpm. Os números no interior de cada região
delimitada do gráfico, por exemplo 40-315, indicam o diâmetro de entrada da bomba
e o diâmetro externo máximo do rotor. Em aplicações, onde há, por exemplo, uma
necessidade especifica de vazão e altura de elevação, a Figura 1.1 pode indicar se
2
esse fabricante irá ou não possuir uma bomba centrífuga que atenda a essa
necessidade.
Figura 1.1 – Exemplo de ábacos fornecido por fabricantes de bombas centrífugas.
(fonte: KSB Bombas)
Uma vez selecionada a bomba centrífuga, o fabricante fornece uma carta
específica, semelhante à mostrada na Figura 1.2. Nessa carta, constam informações
sobre a eficiência, potência consumida e altura de elevação.
Eficiência
Eficiência
Altura Elevação
Potência
Altura Elevação
Potência
Vazão
Figura 1.2 – Curvas de desempenho de uma bomba centrífuga convencional. O
ponto sobre a curva de eficiência é denominado ponto de eficiência máxima ou
ponto de projeto
3
De um modo geral as curvas de operação de bombas centrífugas são
determinadas pelos fabricantes utilizando água como fluido padrão. Quando se
deseja transportar um fluido diferente do padrão, fatores de correção (Hydraulic
Institute, por exemplo) são utilizados para ajustar a curva de desempenho na nova
condição de operação, entretanto, estudos mostram que nem sempre esses fatores
de correção representam o comportamento real desses equipamentos (Amaral,
2007).
Quando bombas centrífugas são utilizas em aplicações específicas, como por
exemplo, no bombeio de fluidos de viscosidade elevada ou ainda no bombeio de
misturas
bifásicas
(líquido-gás),
pode
haver
alterações
significativas
no
comportamento das curvas de desempenho desses equipamentos. Na área
petrolífera, por exemplo, com a descoberta de novas jazidas, na camada pré-sal,
possivelmente bombas centrífugas sejam submetidas a altas pressões e altas
vazões envolvendo misturas multifásicas, o que também caracteriza um grande
desafio de aplicação.
Atualmente, com o desenvolvimento tecnológico na área da informática, a
utilização de simulações numéricas, através de dinâmica dos fluidos computacional,
tem sido uma ferramenta importante no estudo e compreensão de detalhes do
escoamento no interior de bombas centrífugas (Cheah et al., 2007 e Feng et. al.,
2007).
Nesse cenário, propõe-se no presente trabalho estudar numericamente o
escoamento no rotor e no difusor do primeiro estágio de uma bomba centrífuga
comercial de duplo estágio da marca Imbil, modelo ITAP 65-330/2 mostrado na
Figura 1.3. Esse modelo possui dois rotores de oito pás curvadas para trás (em
relação ao sentido de giro do rotor), um difusor de doze aletas no primeiro estágio e
uma voluta na saída do segundo estágio (Amaral, 2007).
4
(a)
(b)
Figura 1.3 – Bomba centrífuga IMBIL ITAP 65-330/2. (b) Fotografia do primeiro
estágio da bomba centrífuga estudada.
1.1 Objetivos
O objetivo do presente trabalho é simular numericamente o escoamento no
interior do rotor e difusor do primeiro estágio de uma bomba centrífuga de duplo
estágio, cujo esquema é mostrado na Figura 1.4.
Para atingir o objetivo proposto, foram resolvidas numericamente as
equações de conservação da massa e da quantidade de movimento, considerando o
escoamento como transiente, newtoniano, monofásico, incompressível e de
propriedades constantes. A turbulência do escoamento foi modelada utilizando o
modelo de turbulência de duas equações κ − ε padrão.
A geometria do primeiro estágio da bomba centrífuga foi construída
virtualmente utilizando o programa SolidWorks 2007. O modelo geométrico obtido foi
implementado no programa comercial de dinâmica dos fluidos computacional
ANSYS CFX 11.0, plataforma que foi utilizada para a solução numérica do
escoamento no interior da bomba centrífuga.
5
A partir dos resultados numéricos obtidos, foram elaboradas curvas de altura
de elevação, potência e eficiência em função da vazão volumétrica, no primeiro
estágio da bomba. Os resultados foram comparados com medidas experimentais
realizadas no LABPETRO/UNICAMP (Amaral, 2007).
Tubo
Rotor
Difusor
Extensão Difusor
Figura 1.4 – Esquema do modelo numérico implementado no programa
computacional Ansys CFX 11.0 para a simulação do escoamento no interior do
primeiro estágio da bomba centrífuga Imbil ITAP 65-330/2.
1.2 Justificativas
O conhecimento da dinâmica do escoamento que ocorre em bombas é de
fundamental importância para o projeto e manutenção desses equipamentos. Com
informações detalhadas sobre o escoamento, ou seja, com o perfil de velocidades,
campo de pressão e tensão de cisalhamento nas paredes, é possível identificar
fatores que possam aperfeiçoar o projeto e escolha do material da bomba.
Do ponto de vista acadêmico, estudar o comportamento do escoamento em
uma bomba centrífuga é um desafio, já que existe uma combinação de efeitos no
escoamento, como a curvatura e rotação, alem da complexidade da geometria, que
dificulta o desenvolvimento da modelagem matemática e numérica do escoamento.
Muitas vezes, são necessários ensaios experimentais para a validação dos
6
resultados numéricos ou teóricos. Os resultados numéricos obtidos no presente
trabalho foram comparados com dados experimentais obtidos por Amaral (2007).
Nesse cenário, a contribuição do presente trabalho está no desenvolvimento
de uma metodologia numérica para o estudo detalhado do escoamento monofásico
em uma bomba centrífuga. Metodologia que pode ser estendida para o estudo do
escoamento em bombas centrífugas sob diferentes condições de operação, como
por exemplo, presença de fluidos viscosos ou de escoamentos bifásicos líquido-gás,
situações que são comumente encontradas na produção de petróleo e gás em
águas profundas.
O presente trabalho é pioneiro de outros que estão sendo desenvolvidos no
Laboratório de Ciências Térmicas (LACIT) da Universidade Tecnológica Federal do
Paraná (UTFPR). Pesquisas estão sendo desenvolvidas, considerando como
referência a metodologia desenvolvida no presente trabalho, com a finalidade de
reproduzir detalhes do escoamento, monofásico ou bifásico, no interior de Bombas
Centrífugas Submersas (BCS).
1.3 Estrutura do trabalho
O presente trabalho está organizado da seguinte maneira: No capítulo 1 uma
introdução e contextualização do tema são feitas, bem como apresentação do foco
dos estudos na bomba centrífuga. O capítulo 2 apresenta uma revisão de trabalhos
reportados da literatura que abordam o estudo do escoamento em bombas
centrífuga, sob o enfoque experimental e numérico. No capítulo 3 são apresentadas
as equações matemáticas que regem o escoamento no interior da bomba centrífuga,
assim como o modelo de turbulência e hipóteses utilizadas. No capítulo 4 é descrita
a metodologia numérica de modelagem do domínio, foco de estudo como, por
exemplo, a confecção e delimitação da geometria utilizada nas simulações
numéricas. No capítulo 5 são mostrados os testes numéricos de parâmetros
importantes nas simulações como a malha espacial adequada, critério de parada,
passo de tempo adequado para as simulações transiente, etc. No capítulo 6 estão
mostrados os resultados obtidos das simulações numéricas do escoamento no
interior do primeiro estágio da bomba centrifuga estudada. Os capítulos 7 e 8
7
reportam, respectivamente, às conclusões, juntamente com sugestões de trabalhos
futuros e referências bibliográficas utilizadas.
8
2 REVISÃO BIBLIOGRÁFICA
Nesse capítulo é mostrada uma revisão de alguns estudos, experimentais e
numéricos, envolvendo bombas centrífugas. Como mencionado, até meados do
século XX, os estudos envolvendo bombas centrífugas eram essencialmente
experimentais e ocorreu um certo esgotamento de possibilidades de se otimizar a
eficiência com melhorias de projeto. Uma revisão sobre fundamentos e estudos
experimentais desenvolvidos com a finalidade de levantar ábacos e curvas de
desempenho global em bombas centrífugas pode ser encontrada em Stepanoff
(1957), Pfleiderer (1960), Hydraulic Institute (1983), Macintyre (1997) e Amaral
(2007).
Junto com desenvolvimento tecnológico, as técnicas experimentais foram
sendo aprimoradas. Técnicas avançadas de visualização com a finalidade de
observar detalhes do escoamento no interior das bombas centrífugas, como o LDV
(Laser-Doppler Velocimetry, um anacrônico em língua inglesa que significa
Velocimetria por Laser-Doppler), PIV (Particle Image Velocimetry, ou Velocimetria
por imagens de partículas), câmeras de alta velocidade de obturação, entre outros,
foram desenvolvidas. Essas informações fornecem ao projetista da bomba centrífuga
subsídios para confecção de equipamentos mais eficientes.
Neste contexto, Li (1999) fez um estudo experimental sobre a influência da
viscosidade do fluido no desempenho de bombas centrífugas. Através de ensaios,
utilizando medição de velocidade por LDV ele avaliou o desempenho de uma bomba
centrífuga operando com dois fluidos distintos: primeiramente com água ( ρ =1000
kg/m3 e ν =1mm2/s) como fluido de trabalho e, posteriormente com óleo ( ρ =851
kg/m3 e ν =48mm2/s). Ele mostrou que o aumento da viscosidade do fluido, aumenta
a altura de elevação da bomba centrífuga, diminui a eficiência e aumenta a potência
consumida (ver Figura 2.1). Segundo Li (1999), o desempenho da bomba diminui,
quando opera com fluidos de altas viscosidades, devido ao aumento nas perdas por
atrito, tanto na sucção como na saída do rotor.
2
1.5
1
0.5
0
14
Óleo
Água
Altura
12
10
8
Eficiência
6
4
Potência
2
0
0
2
4
6
Vazão [L/s]
8
100
90
80
70
60
50
40
30
20
10
0
Eficiência [%]
Potência [ kW ]
2.5
Altura de Elevação [m]
9
10
Figura 2.1 – Desempenho de uma bomba centrífuga ao bombear fluidos de
diferentes viscosidades (Li, 1999)
Wuibaut et al. (2005) buscaram desenvolver técnicas experimentais de
medição em bombas centrífugas submersas (BCS). Os experimentos foram feitos
em um rotor de BCS conhecido da literatura (rotores SHF), sobre o qual muitos
dados experimentais estão disponíveis, principalmente para validação de simulações
numéricas. Os experimentos de Wuibaut et al. (2005) foram realizados utilizando as
técnicas de visualização PIV-2D para a BCS operando com ar e LDV-2D para BCS
operando com água. Devido à diferença das propriedades dos fluidos utilizados,
bombas de tamanhos diferentes tiveram que ser construídas, porém, considerando
que os escoamentos fossem semelhantes. Para a bomba operando apenas com ar,
Wuibaut et al. (2005) não utilizaram difusor, nem voluta, diferentemente do que foi
feito para o caso da bomba operando com água, onde um difusor e uma voluta
industrial foram acoplados ao rotor. Apesar das pequenas diferenças de medição
encontradas entre os dois testes, os campos de velocidades obtidos por eles foram
similares. Eles concluíram que ambas as técnicas de medição utilizadas são
satisfatórias para avaliação de campos de velocidade no interior de bombas
centrífugas, por exemplo, para validação de técnicas numéricas.
Güilich (1999) estudou o comportamento de bombas centrífugas operando
com óleos pesados, utilizando um enfoque analítico-experimental. Em geral, o
desempenho de bombas operando com óleos viscosos é estimado por comparação
ao desempenho da bomba operando com água, através de fatores de correção da
10
vazão, altura da bomba e eficiência. Segundo ele, esses fatores de correção nem
sempre correspondem à realidade. Ele fez um tratamento aprofundado das perdas
existentes nas bombas, e as separou em quatro parcelas: perdas volumétricas,
perdas hidráulicas, perdas por atrito e perdas mecânicas. Segundo Güilich (1999),
as parcelas mais influenciadas pelo aumento da viscosidade foram as perdas
hidráulicas e por atrito. Por meio de relações semi-empíricas, foram estimados os
coeficientes de correção considerando somente essas perdas. Tais correlações
foram validadas com sucesso. Ele concluiu que seus resultados foram bons e
enfatiza que fazer um tratamento das perdas é vantajoso por poder ser aplicado a
qualquer modelo de bomba.
Amaral (2007), estudou experimentalmente o comportamento de uma bomba
centrifuga comercial (Imbil ITAP 65 330/2) e dois modelos de bombas centrífugas
submersas (BCS). O estudo realizado foi concentrado na determinação da faixa
operacional
de
uma
BCS
operando
com
fluidos
viscosos.
Ele
levantou
experimentalmente as curvas de desempenho para diversas condições de operação
e apresentou uma metodologia para a correção das curvas de desempenho de
bombas centrífugas sob a presença de fluidos viscosos.
Com o aparecimento da mecânica dos fluidos computacional e, também, com
o desenvolvimento da informática, a pesquisa e o desenvolvimento de bombas
centrífugas ganhou força. Surgiram diversos programas computacionais comerciais
com enfoque específico na simulação de máquinas de fluxo (Ansys CFX, Flow3D,
Fluent, entre outros) (Amaral, 2007).
Os programas computacionais de dinâmica dos fluidos computacional
possuem uma diversidade de modelos matemáticos que buscam descrever os
comportamentos físicos presentes no escoamento real. Porém, sob certas
condições, como por exemplo, em escoamentos turbulentos se faz necessária a
utilização de modelos que busquem representar adequadamente esse fenômeno.
Sempre que se utiliza um programa de mecânica dos fluidos computacional para
simular escoamentos em bombas centrífugas, com fins de projeto ou de
desenvolvimento, é necessário, posteriormente, validar os resultados obtidos com
11
um protótipo real. Isso é feito porque podem ocorrer pequenas divergências entre o
modelo numérico e o fenômeno físico real.
Nesse contexto, Asuaje et al. (2005) estudaram numericamente uma bomba
centrífuga comercial com o objetivo de avaliar a influência da voluta nos campos de
pressão e velocidade da bomba. Eles testaram três modelos de turbulência
diferentes, k − ε , k − ω e SST, e concluíram que não houve diferenças significativas
entre os resultados dos campos médios de velocidade e pressão. Para a bomba
operando na vazão de projeto, eles mostraram que o escoamento é praticamente
isento de recirculações desde a entrada até a saída da voluta. Para vazões menores
que a condição de projeto, verificaram que ao sair do rotor o escoamento está sujeito
a um gradiente adverso pronunciado, o que gera instabilidades e recirculações e
conseqüentemente, resulta em perdas de rendimento. Por outro lado seus estudos
mostram que quando a bomba centrífuga está operando em vazões maiores que a
de projeto, não há um gradiente de pressão intenso, porém, devido à inércia do
escoamento ser elevada ocorrem recirculações no interior da voluta, resultado de
perdas de rendimento da bomba. Essas recirculações, aliadas à forma da voluta,
geram uma força resultante que o fluido irá transmitir ao mancal onde a bomba
centrífuga está apoiada. Asuaje et al. (2005) mostraram ainda que, para a vazão de
projeto, a magnitude dessa força é menor quando comparada às demais.
Muitos pesquisadores atentaram para a possível interação existente entre
rotor e difusor de bombas centrífugas. Essa interação foi avaliada de diversas
maneiras como, por exemplo, observando as variações das forças atuantes no rotor
e difusor quando da alteração do número de aletas do difusor, das pás do rotor, da
distância relativa entre esses componentes, etc.
Feng et al. (2007) estudaram o comportamento transitório de uma bomba
centrífuga, juntamente com a interação existente entre rotor e um difusor aletado.
Eles fizeram um estudo de independência de malha e uma validação experimental
da bomba estudada, em termos de altura de elevação versus vazão, onde obtiveram
boa concordância. De maneira análoga aos estudos de Asuaje et al. (2005), Feng et
al. (2007) simularam numericamente uma bomba centrífuga na condição de projeto e
em duas vazões distintas (uma maior e outra menor). Seus resultados também
12
mostram que, para vazões inferiores ou superiores à de projeto, ocorrem
recirculações no rotor e no difusor. Nos resultados apresentados por Feng et al.
(2007), um estudo da distância relativa entre rotor e difusor foi feito. Para isso, eles
simularam numericamente o escoamento em uma bomba centrífuga com frestas
entre rotor e difusor de 3% e 10% do raio externo do rotor, respectivamente (Ver
Figura 2.2 (a) . Os resultados obtidos por eles são mostrados na Figura 2.2 (b), onde
as curvas se referem à flutuação da altura de elevação em função da direção
tangencial. Da Figura 2.2 (b), Feng et al. (2007) mostraram que quanto mais distante
o difusor se encontra do rotor (maior fresta), menores são as amplitudes de
oscilação das variáveis calculadas (altura de elevação, torque, etc) com o tempo, ou
H [m]
seja, menor é a interação entre rotor e difusor.
6.4
6.3
6.2
6.1
6
5.9
5.8
5.7
5.6
Fresta 3%
Fresta 10%
0
20
(a)
40
60
80
ϕ [graus]
(b)
100 120
Figura 2.2 – Interação entre difusor aletado e rotor, no escoamento de uma bomba
centrífuga. (Feng et al. (2007)).
Posteriormente,
Cheah
et
al.
(2007)
simularam
numericamente
a
fluidodinâmica do escoamento em uma bomba centrífuga convencional de um
estágio com rotor aberto e voluta em espiral. Na Figura 2.3, é evidenciada a
diferença de escoamento entre a condição de projeto (Figura 2.3(a)) e uma condição
de vazão maior (Figura 2.3(b)). Na condição de projeto, eles mostram que o
escoamento é praticamente isento de recirculações, que se reflete em poucas
perdas e, conseqüentemente, em maior eficiência. Para o caso onde a vazão é
13
maior, eles mostram que há recirculação na saída da voluta, reflexo de perdas de
energia no interior da bomba centrífuga e conseqüente redução de eficiência.
Velocidade Absoluta
(a)
Velocidade Absoluta
(b)
Figura 2.3 – Estudo da fluidodinâmica do escoamento em uma bomba centrífuga
apresentado por Cheah et al. (2007). (a) Vazão na condição de projeto e (b) vazão
acima da condição de projeto.
2.1 Teoria de máquinas de fluxo aplicado a bombas
centrífugas
Devido à geometria complexa das partes internas de uma bomba centrífuga,
alguns trabalhos focam no estudo ou modelagem de parâmetros específicos desse
equipamento. Através da teoria de turbomáquinas são definidos diversos ângulos e
termos que serão utilizados ao longo do presente trabalho.
Figura 2.4 – Volume de controle e componentes de velocidades do escoamento na
entrada e saída de um rotor genérico
14
A Figura 2.4 mostra o desenho de um rotor, onde estão definidos o raio da pá,
as componentes da velocidade do escoamento tanto na secção de entrada
(denotada pelo índice 1) como na secção de saída do rotor (índice 2). Através do
princípio da quantidade de movimento angular é possível obter uma equação que
expressa o torque de eixo, Teixo, em função das componentes da velocidade do
escoamento à entrada e saída do rotor, como mostra de forma escalar a equação
(2.1).
ɺ
Teixo = ( r2 Vt2 − r1Vt1 ) m
(2.1)
A equação (2.1), denominada de equação de Euler para Turbomáquinas,
pode ser utilizada para determinar a potência fornecida ao rotor, através da relação
ɺ = ωT :
escalar W
eixo
ɺ = ωT = ω (r V − r V ) m
ɺ = (u V − u V ) m
ɺ ⇒W
ɺ
W
eixo
2 t2
1 t1
2 t2
1 t1
(2.2)
ɺ , se obtém a magnitude da energia
Ao dividir-se a equação (2.2) por mg
transferida ao rotor em unidades de comprimento:
H=
ɺ
W
1
= (U2 Vt2 − U1Vt1 )
ɺ
mg g
(2.3)
A equação (2.3) apresentada a altura de elevação como função das
componentes de velocidade na entrada e saída do rotor. Obviamente, as equações
(2.1) a (2.3) exprimem de forma muito simplificada o comportamento de uma bomba
centrífuga, uma vez que uma série de fenômenos não estão sendo considerados
nessa modelagem. Porém, dessa modelagem é possível notar a importância das
componentes de velocidade na obtenção dos valores de torque, potência e altura de
elevação. Isso sugere que é necessário definir a forma com que o escoamento entra
e sai do rotor, através do formato das pás.
A Figura 2.5 mostra que na condição de projeto de uma bomba centrífuga, o
escoamento é admitido como entrando e saindo tangencialmente ao perfil da pá do
15
rotor. O ângulo β é medido a partir da direção tangencial em direção à pá do rotor,
no sentido oposto ao de giro do rotor. Esse ângulo é construtivo, ou seja, depende
somente da geometria da pá do rotor. O ângulo α, por sua vez, é ditado pela
condição de operação da bomba, ou seja, depende da magnitude da rotação
imposta ao rotor. O ângulo α também é dependente da vazão que flui pela bomba,
vazão esta, que é controlada por válvulas no sistema de bombeamento a que a
bomba está conectada.
(a)
(b)
Figura 2.5 – Geometria e notação usadas para desenvolver diagramas de
velocidade em bombas centrífugas. (a) Distribuição dos perfis de velocidade ao
longo da pá do rotor e (b) os triângulos de velocidades na entrada e saída do rotor.
Com a definição das velocidades do escoamento, e os ângulos que elas
formam, pode-se então avaliar o valor das componentes tangenciais, Vt1 e Vt2, em
função das velocidades absolutas, V1 e V2, pelas seguintes relações:
Vt 1 = V1 cos (α1 )
Vt 2 = V2 cos (α 2 )
(2.4)
Substituindo as relações dadas pela equação (2.4) na equação (2.1), do
torque no eixo do rotor da bomba, tem-se:
ɺ ( r2 V2 cos (α 2 ) − r1V1 cos (α1 ) )
Teixo = m
(2.5)
Com a equação (2.5) pode-se, assim, determinar a altura de elevação da
bomba, através de parâmetros do escoamento e dimensões características do rotor.
Portanto, substituindo a equação (2.5) na equação (2.3):
16
H=
ω
ɺg
m
Teixo =
ω
g
(r V
2
2
cos (α 2 ) − r1V1 cos (α1 ) )
(2.6)
Por a equação (2.6) se tratar de uma idealização do escoamento em uma
bomba, é conveniente alterar a nomenclatura da altura de elevação para sinalizar
que o resultado obtido será para um caso Ideal. Com isso, define-se Ht∞ como
sendo a altura de elevação teórica infinita. O índice t indica um processo sem perdas
(teórico) e o índice ∞ sinaliza que a bomba possui infinitas pás. Assim:
Ht∞ =
ω
ɺ
mg
Teixo =
ω
g
(r
2
V2 cos (α 2 ) − r1 V1 cos (α1 ) )
(2.7)
Da análise da equação (2.7), pode-se notar que a energia transferida ao fluido
varia proporcionalmente com a velocidade angular (quanto maior a magnitude da
rotação do rotor, maior é a quantidade de energia transferida ao fluido). Os dois
termos entre parênteses têm sinal invertido, por isso, suas contribuições à
quantidade de energia transferida são opostas. Portanto, a quantidade de energia
específica, Ht∞ , transferida ao fluido será máxima quando o termo negativo for nulo.
Isto ocorre quando o ângulo α1=90°. Isso faz com que a velocidade absoluta do
fluido, na entrada do rotor seja perpendicular à direção tangencial. Na prática, devido
às perdas envolvidas no processo de transferência de energia α1 não é exatamente
igual a 90°, porém, é próximo desse valor, fazendo com que o termo negativo da
equação (2.7), de fluxo de quantidade de movimento angular na entrada do rotor
seja pequeno quando comparado com o fluxo de quantidade de movimento angular
na saída do rotor, (V2 r2 cosα2). Dessa forma, tem-se:
Ht∞ =
ω
ɺ
mg
Teixo =
ω
g
r2 V2 cos (α 2 ) =
u2
Vt 2
g
(2.8)
onde Vt2 é a componente tangencial da velocidade absoluta do fluido (na saída do
rotor) e u2 é a velocidade periférica do rotor, devido à rotação. A Figura 2.6 mostra o
triangulo de velocidades na saída do rotor, com as componentes, tangencial e
normal (ou radial), do vetor velocidade absoluta do fluido.
17
Figura 2.6 – Componentes (tangencial e radial) da velocidade absoluta na secção de
saída do rotor da bomba
Observando a Figura 2.6 e a equação (2.8), pode-se notar que quanto maior
a rotação ou o tamanho do rotor, maior será u2 e quanto maior for a componente
tangencial da velocidade absoluta da bomba, na saída do rotor, maior a altura de
elevação, Ht∞ .
A formulação apresentada na equação (2.8), ainda não leva em conta
características operacionais como a vazão nem mesmo construtivas do rotor como o
valor do ângulo β, responsável pela geometria da pá do rotor. Para explicitar esses
parâmetros na equação (2.8), com o auxílio da Figura 2.6, pode-se escrever a
componente tangencial da velocidade absoluta em função componente normal da
velocidade absoluta, ou seja:
Vn2 = w 2 sen( β 2 ) ⇒ w 2 =
Vn2
sen( β 2 )
(2.9)
Vt 2 = u2 − w 2 cos( β 2 ) ⇒ Vt2 = u2 − Vn2 cot g( β 2 )
Substituindo a equação (2.9) na equação da altura de elevação ideal
(equação (2.8)):
Ht∞ =
ω
ɺ
mg
Teixo =
ω
g
r2 V2 cos (α 2 ) =
u2
(u2 − Vn2 cotg ( β2 ) )
g
(2.10)
A componente normal da velocidade absoluta na saída do rotor, Vn2, está
relacionada com a vazão mássica que flui pela bomba. A Figura 2.7 apresenta um
desenho esquemático da geometria do rotor de uma bomba centrífuga. Utilizando a
18
equação da conservação da massa e os dados fornecidos na Figura 2.7, pode-se
escrever que:
Figura 2.7 – Corte axial do rotor de uma bomba centrífuga
ɺ ⇒ Vn2 =
ρ Vn2 A 2 = m
ɺ
ɺ
Q
Q
⇒ Vn2 =
A2
2π r2b2
(2.11)
Substituindo a equação (2.11) na equação (2.10), tem-se:
Ht∞ =
ɺ

u2 
Q
u
−
cotg ( β 2 ) 
 2
g 
2π r2b2

(2.12)
A equação (2.12) é mais interessante pois apresenta, explicitamente, as
principais variáveis operacionais da bomba tais como vazão, geometria e rotação
(apresentada através da velocidade tangencial). Analisando a equação (2.12), podese notar que a altura de elevação recebe uma influência quadrática da variação da
rotação. Nessa formulação simplificada, a altura de elevação varia linearmente com
o aumento da vazão. A forma da pá (ditada pelo ângulo β2) determina se a vazão
influencia negativamente (β2<90°) ou positivamente ( β2>90°) na altura de elevação.
Se β2=90° a vazão não influencia na magnitude do valor d a altura de elevação.
A curva característica de uma máquina de fluxo é, por definição, a curva que
representa a dependência que existe entre a quantidade de energia transferida pela
máquina (real ou idealizada) e a vazão do fluido. Ela se aplica a uma bomba de
características geométricas conhecidas (valores especificados de r2 e b2, dimensões
geométricas que aparecem na equação (2.12)), operando com rotação ω , também
19
especificada. A Figura 2.8 mostra a relação entre a altura de elevação e a vazão
volumétrica da bomba, para diferentes ângulos β2 .
Figura 2.8 – Curva característica idealizada de uma bomba centrífuga
A curva característica real de uma bomba centrífuga difere substancialmente
destas curvas idealizadas, qualquer que seja o valor do ângulo β2. A curva para
β2<90º apresenta a mesma tendência da curva real: redução de Ht∞ à medida em
que a vazão aumenta. As curvas para β2=90º e β2>90º, entretanto, têm
comportamento inverossímil. Na medida em que a vazão aumenta, é de se esperar
que, nos escoamentos reais (viscosos), a energia dissipada (em perdas hidráulicas,
por exemplo) aumente com o quadrado da vazão. Assim, parcela substancial da
potência disponível no eixo é irreversivelmente dissipada, e a energia específica
transferida não pode, indefinidamente, aumentar, ou mesmo se manter constante,
com o aumento da vazão.
A influência da magnitude do ângulo β2 sobre a curva característica da
bomba, e sobre as formas construtivas das pás dos rotores, entretanto, deve ser
objeto de análise. As bombas centrífugas quase sempre apresentam rotores com
pás curvadas para trás em relação ao sentido de rotação do rotor, isto é, β2<90º. Os
valores usuais para β2 estão por volta dos 30º. Em bombas mais antigas, de
pequena potência e baixa responsabilidade, pode-se ainda encontrar rotores com
20
pás inteiramente radiais, com β1=β2=90º. Atualmente, este tipo de rotor está em
desuso.
Benra (2001) apresenta uma metodologia numérica para o desenvolvimento
de rotores eficientes a serem utilizados em bombas centrífugas. Seus estudos se
basearam na modelagem e simulação numérica do escoamento em rotores radiais
isolados, ou seja, sem contabilizar volutas ou difusores. Seu principal objetivo foi
avaliar a eficiência hidráulica, através da comparação do escoamento em dois
rotores de mesmas dimensões, porém, com diferenças na configuração das pás. O
primeiro rotor estudado por ele possuía pás com bordo de entrada alinhados com a
entrada do rotor, como mostra a Figura 2.9 (a). Já o segundo rotor possuía pás com
bordo de entrada inclinados em relação à entrada do rotor (Figura 2.9 (b). Essa
diferença na geometria da entrada dos rotores reflete na forma com que a pá evolui
da entrada até a saída do rotor.
Bordo de entrada
das pás
(a)
(b)
Figura 2.9 – Rotores estudados numericamente por Benra (2001). (a) Pás com bordo
de entrada horizontal e (b) pás com bordo de entrada inclinadas em relação à
entrada do rotor.
No primeiro rotor, Benra (2001) observou que a evolução do ângulo β
(ângulo construtivo do rotor) desde a entrada até a saída do rotor, no cubo, possuía
uma grande oscilação. Já o segundo rotor não continha tais variações, sendo a
transição de β suave desde a entrada até a saída. Como conseqüência dessas
diferenças de geometria, o escoamento no segundo rotor apresentou menores
instabilidades e recirculações para vazões fora da condição de projeto. Isso levou à
conclusão de que a forma da entrada do rotor influencia no desempenho da bomba.
Bacharoudis et al. (2008) estudaram o comportamento do escoamento em
uma bomba centrífuga (composta de um rotor e uma voluta). O foco do trabalho foi
avaliar as mudanças no comportamento da altura de elevação e eficiência da
21
bomba, quando da alteração do ângulo de saída das pás do rotor, β2 . Eles
mostraram que quando há um aumento do ângulo de saída, há um aumento na
altura de elevação da bomba, porém com redução na eficiência. Entretanto, esse
comportamento não é observado em todo o campo operacional da bomba. Eles
mostraram que para a bomba operando em vazões maiores que a de projeto, ao
utilizar um ângulo de saída da pá maior, ocorre um aumento na eficiência em relação
ao rotor original. Já para vazões inferiores à de projeto, utilizando um ângulo de
saída da pá menor, observaram uma redução na eficiência, também em relação ao
rotor original. Shoajaee e Boyaghchi (2007) também mostraram que aumentando o
ângulo de saída de um rotor centrífugo, há um aumento na altura de elevação.
Esses resultados estão em acordo com a equação (2.12), uma vez que quanto maior
o ângulo β 2 , menor será a contribuição do termo de sinal negativo no lado direito
dessa equação.
É importante salientar que a análise de resultados oriundos de simulações
numéricas do escoamento em componentes isolados pode levar a conclusões que
não correspondem efetivamente ao que acontece em uma bomba centrífuga real.
Isso porque há interações entre os componentes desse tipo de equipamento, o que
foi comprovado experimentalmente por Amaral (2007) e numericamente por Feng et
al. (2007).
Recentemente, Segala et al. (2008) apresentam um estudo preliminar do
escoamento no rotor e difusor de uma bomba centrífuga em regime permanente.
Nesse estudo, foi discutido que, para resolver o acoplamento entre rotor e difusor,
uma
metodologia
numérica
adequada
deve
simular
ambos
os
domínios
simultaneamente. Eles também apontam a necessidade de dados experimentais
para o desenvolvimento de uma metodologia numérica confiável.
Com o contínuo desenvolvimento da tecnologia, e a crescente capacidade de
processamento dos computadores, tem sido possível simular numericamente
escoamentos em equipamentos reais de interesse de engenharia, como bombas
centrífugas, por exemplo. Os primeiros estudos numéricos envolvendo bombas
centrífugas eram bastante simplificados devido à limitação computacional existente.
22
Atualmente é possível simular numericamente escoamentos em rotores com
geometrias mais próximas do modelo real, o que torna esse tipo de abordagem
versátil e interessante, pois, permite testar diversos parâmetros e configurações de
rotores, por exemplo, de maneira simples e relativamente rápida.
O presente trabalho visa desenvolver uma metodologia de simulação
numérica de escoamentos em bombas centrífugas. Para isso, será utilizada a
geometria de uma bomba centrífuga comercial, de onde foram extraídos dados
experimentais por Amaral (2007). Esses dados experimentais servem de base de
comparação e validação para a metodologia numérica desenvolvida no presente
trabalho.
23
3 MODELAGEM MATEMÁTICA
Para simular numericamente a fluidodinâmica do escoamento em bombas
centrífugas, as equações de conservação da massa e da quantidade de movimento
são utilizadas. Nesse tipo de equipamento, têm-se simultaneamente domínios
rotativos (rotores) e estacionários (tubo de admissão da bomba centrífuga, carcaça,
difusor, voluta, entre outros). Os programas comerciais de simulação numérica de
escoamentos quando aplicados ao estudo de bombas centrífugas oferecem
soluções de sistemas de múltiplos domínios, isto é, há um sub-domínio específico
para cada parte do equipamento em estudo (seja rotativo ou estacionário), como
mostra a Figura 3.1. Nos domínios rotativos, os efeitos da rotação são impostos por
meio de termos fontes, aplicados às respectivas componentes da equação de
quantidade de movimento.
Figura 3.1 – Vários sub-domínios co-existentes em um programa comercial de
simulação numérica de escoamentos envolvendo máquinas de fluxo.
Para acoplar os domínios rotativos aos estacionários existem modelos de
interface que transferem as informações de um domínio ao outro. Nesse capítulo
será apresentada a modelagem matemática das equações de conservação da
massa e da quantidade de movimento utilizadas pelo programa computacional na
simulação numérica do escoamento de uma bomba centrífuga, bem como as
equações do modelo de turbulência escolhido.
24
3.1 Equações de conservação da massa e da quantidade de
movimento escritas, para um sistema de coordenadas
rotativo
Para estudar o escoamento no interior de domínios providos de rotação é
conveniente utilizar um sistema de coordenadas que acompanhe o giro desse
domínio. Isso facilita a implementação de diversos parâmetros do programa
numérico, como por exemplo, construção da malha, aplicação das condições de
contorno e processamento dos resultados. Para contabilizar os efeitos de rotação
em domínios desse tipo, necessita-se de termos adicionais nas componentes
equações da quantidade de movimento, como está mostrado a seguir.
Considerando o rotor mostrado na Figura 3.2 onde em seu eixo de rotação é
acoplado um sistema de coordenadas não-inerciais (x,y,z), ou seja, que gira na
mesma rotação do rotor. Um sistema de coordenadas inerciais (X,Y,Z) é
estabelecido conforme mostrado na Figura 3.2.
Z
z
y
Y
x
X
Figura 3.2 – Sistema de coordenadas rotativo aplicado a uma bomba centrífuga
Da Figura 3.2 e utilizando conceito de velocidade e aceleração relativa é
possível obter as equações de conservação da massa e quantidade de movimento
25
escritas utilizando o sistema de coordenadas não-inercial como referencial, o que
resulta em:
∇⋅Vxyz = 0
(3.1)
DVxyz
1
2
− ∇p + ν∇ Vxyz + g = 2 ω× Vxyz + ω× ω× r +
ρ
Dt
(
)
(3.2)
onde ρ é a massa específica do fluido, p é a pressão hidrostática, ν é a viscosidade
cinemática, Vxyz é a velocidade do fluido no sistema de coordenadas não-inercial, g
é a aceleração da gravidade, ω é a velocidade angular do rotor e r é a posição de
uma partícula fluida em relação à origem do sistema de coordenadas não-inercial.
Os termos do lado esquerdo da Equação (3.2) são, respectivamente, o gradiente de
pressão, a dissipação viscosa e o termo gravitacional. O último termo do lado direito
da Equação (3.2) representa a representa a aceleração temporal e convectiva do
fluido. O termo 2 ω× Vxyz é denominado de aceleração de Coriollis e surge devido à
mudança do sistema de coordenadas inercial para o não-inercial. O termo
ω× ω× r , chamado de aceleração centrípeta também surge devido à mudança de
(
)
referencial.
O escoamento a ser estudado é considerado incompressível, com
propriedades constantes, isotérmico e o rotor gira com velocidade angular constante
em torno de um eixo fixo (sem movimento de translação). A equação da quantidade
de movimento para os domínios estacionários (tubo de entrada, difusor) pode ser
obtida da Equação (3.2) bastando atribuir ω = 0 .
3.2 Modelagem da turbulência
Da literatura sabe-se que a fluidodinâmica do escoamento no interior de
máquinas de fluxo ocorre em regime turbulento. Portanto, é necessário incorporar
tais efeitos na modelagem do problema estudado. De maneira resumida, os
escoamentos em regime turbulento são intrinsecamente transientes, ou seja, os
26
valores de propriedades desse tipo de escoamento flutuam com o tempo. Porém,
utilizando métodos estatísticos é possível observar uma certa coerência nesse tipo
de escoamento.
φ
φ'
φ
t
Figura 3.3 – Propriedade genérica de um escoamento em regime turbulento em
função do tempo, t
Osbore Reynolds em 1895 observou que era possível escrever o valor de
uma propriedade genérica do escoamento turbulento através de duas parcelas: uma
média e uma flutuação, como está representado na Figura 3.3 e na Equação (3.3)
(Wilcox, 1993).
φ(t ) = φ + φ '(t )
(3.3)
onde φ(t ) representa uma propriedade genérica instantânea do escoamento em um
ponto do espaço, φ é o valor médio no tempo dessa propriedade e φ '(t ) é o valor da
flutuação temporal dessa propriedade. Com essa forma de abordagem, Reynolds
propôs um estudo da turbulência aplicando médias temporais às equações da
conservação da massa e da quantidade de movimento resultado da substituição da
Equação (3.3) nas variáveis velocidade e pressão, resultando em:
ρ
onde
∇⋅Vxyz = 0
(3.4)
DV
= −∇ p + ∇τij
Dt
(3.5)
27
 ∂u ∂u 
τij = µ  i + j  − ρui′u ′j
(3.6)
 ∂x j ∂xi 


onde o símbolo u i representa a componente do vetor de velocidade na direção
coordenada “i” (x, y ou z). Analisando a Equação (3.5), pode-se dizer que os termos
turbulentos de inércia se comportam como se a tensão total, τij , no sistema fosse
composta de uma tensão viscosa newtoniana mais um tensor tensão turbulenta
aparente ρui′u ′j , denominado Tensor Tensão de Reynolds. Esse produto das
flutuações médias é uma razão média temporal da transferência de quantidade de
movimento devido à turbulência, sendo a prescrição desse fluxo de quantidade de
movimento o grande desafio na modelagem da turbulência. A modelagem desses
termos é necessária para resolver o sistema de equações (fechamento do problema
da turbulência).
Para resolver o problema de fechamento da turbulência Boussinesq propôs
um modelo onde se assume que o produto médio das flutuações de velocidade
( ρui' u 'j ) é proporcional à taxa de deformação média do fluido (Wilcox,1993), ou seja:
 ∂u ∂u j 
(3.7)
−ρui'u 'j = µt  i +

 ∂x
 j ∂xi 
onde i e j são índices que representam as direções coordenadas (x, y ou z). O
coeficiente µt é denominado viscosidade dinâmica turbulenta e é uma característica
do escoamento estudado.
Para o presente trabalho, assume-se como hipótese escoamento turbulento
de fluido newtoniano governado pelas equações de conservação da massa e da
quantidade de movimento. Utilizam-se as equações médias de Reynolds e a
hipótese de Boussinesq para a modelagem da turbulência. Sendo assim, as
equações da conservação da massa e da quantidade de movimento em x, y e z
assumem, respectivamente a forma:
∂u ∂v ∂w
+
+
=0
∂x ∂y ∂z
(3.8)
28
∂u ∂uu ∂uv ∂uw
1 ∂p ∂ 
 ∂u  
+
+
+
=−
+
( ν + ν t )  2
 +
∂t
∂x
∂y
∂z
ρ ∂x ∂x 
 ∂x  
 ∂u ∂v   ∂ 
∂ 
∂u ∂w 
+
+
( ν + νt )  +   + Sx
( ν + ν t ) 
 +

∂y 
 ∂z ∂x  
 ∂y ∂x   ∂z 
 ∂v  
∂v ∂vu ∂vv ∂vw
1 ∂p ∂ 
+
+
+
=−
+
( ν + ν t )  2
 +
∂t
∂x
∂y
∂z
ρ ∂y ∂y 
 ∂y  
 ∂u ∂v   ∂ 
 ∂v ∂w  
∂ 
+ ( ν + ν t ) 
+
+
+
( ν + ν t ) 
  + Sy

∂x 
 ∂y ∂x   ∂z 
 ∂z ∂y  
(3.9)
(3.10)
∂w ∂wu ∂wv ∂ww
1 ∂p
∂ 
 ∂w ∂u  
+
+
+
=−
+ + ( ν + νt ) 
+
 +
∂t
∂x
∂y
∂z
ρ ∂z
∂x 
 ∂x ∂z  
(3.11)
 ∂w ∂v   ∂ 
∂ 
 ∂w  
+
+
( ν + ν t ) 
 +
( ν + νt )  2 ∂z  
∂y 


 ∂y ∂z   ∂z 
onde u , v e w são as componentes médias de velocidade nas direções x, y e z, ν
é a viscosidade cinemática do fluido, νt é a viscosidade cinemática turbulenta dada
por νt = µt / ρ , Sx e Sy são termos fonte. Para o caso de um sistema de
coordenadas inerciais, esses termos fontes são iguais a zero. Considerando um
sistema de coordenadas rotativo, com velocidade angular constante, tem-se:
Sx = −2 ωz .v − ωz 2 x
(3.12)
Sy = 2 ωz .u − ωz 2 y
(3.13)
As barras sobre as componentes de velocidade e pressão são mantidas para
denotar a aplicação das médias temporais de Reynolds. Para modelar o valor de
µt = ρ.ν t em cada ponto do domínio, surgiram diversos modelos ao longo dos anos.
Dentre esses, os denominados modelos de turbulência a duas equações são
bastante utilizados. No presente trabalho, é utilizado o modelo de turbulência a duas
equações κ − ε padrão de Launder e Spalding (1974).
29
3.2.1 Modelo de Turbulência a duas equações κ − ε padrão
O modelo de turbulência a duas equações κ − ε padrão é, talvez, o mais
conhecido dos modelos de turbulência. A idéia principal desse modelo é estabelecer
uma relação entre a variável ν t e duas outras propriedades turbulentas do
escoamento, que para esse modelo são as variáveis κ e ε :
νt =
Cµ κ 2
(3.14)
ε
onde κ é a energia cinética turbulenta, ε é a componente isotrópica da taxa de
dissipação de energia, Cµ é uma constante do modelo. Como já mencionado, duas
equações diferenciais adicionais são utilizadas nesse modelo, uma para calcular κ e
outra para calcular ε . As equações para κ e ε em coordenadas cartesianas são,
respectivamente,
Dκ
= Dκ + Pκ − ε
Dt
(3.15)
Dε
ε
= Dε + (Cε1Pκ − Cε 2ε )
(3.16)
Dt
κ
onde o lado esquerdo das equações (3.15) e (3.16) estão contemplados os termos
de variação temporal e transporte convectivo das propriedades turbulentas, Dκ e Dε
são termos difusivos de κ e de ε , respectivamente e Pκ é o termo de produção de
κ . Cε1 e Cε 2 são coeficientes do modelo. Os termos Dk , Dε e Pk são dados
respectivamente por:
Dκ =
νt
∂ 
 ν +
∂x 
σκ
 ∂κ  ∂ 
νt
 ν +
 +
σκ
 ∂x  ∂y 
 ∂κ  ∂ 
ν t  ∂κ 
 ν +
 +
 
σ κ  ∂z 
 ∂y  ∂z 
(3.17)
Dε =
ν t  ∂ε  ∂ 
ν t  ∂ε  ∂ 
ν t  ∂ε 
∂ 
 ν +   +
 ν +   +
 ν +  
∂x 
σε  ∂x  ∂y 
σε  ∂y  ∂z 
σε  ∂z 
(3.18)
30
 ∂u ∂v  2  ∂u ∂w 2  ∂v ∂w  2
Pκ = ν t 
+
+
+
 +
 +
 +
 ∂y ∂x   ∂z ∂x   ∂z ∂y 
(3.19)
 ∂u  2  ∂v  2  ∂w  2  
+2  
 +
 +
 
 ∂x   ∂y   ∂z   
onde σk e σε são os números de Prandtl turbulento para κ e ε , respectivamente.
Os coeficientes adotados pelo modelo κ − ε padrão são Cµ = 0,09 , Cε1 = 1,44 ,
Cε 2 = 1,92 , e os números de Prandtl turbulentos são σk = 1,0 e σε = 1,30 . Uma
restrição a esse modelo é que ele não pode ser utilizado em regiões muito próximas
à parede. Uma discussão sobre as conseqüências dessa limitação é feita com
maiores detalhes no capítulo 6.
3.2.2 Condições de contorno do problema
As condições de contorno do problema, mostradas esquematicamente na
Figura 3.4, são: condição de não-deslizamento e de impermeabilidade sobre todas
as superfícies sólidas do domínio ( u = v = w = 0 )parede. O valor da energia cinética
turbulenta κ , deve também ser nulo sobre todas as superfícies sólidas. A condição
de contorno para a dissipação da energia cinética turbulenta, ε , em superfícies
sólidas é calculada através da relação empírica mostrada por Versteeg e
Malalasekera (1995) dada por ε = Cµ3 / 4 .κ3p / 2
( ∆y .κv ) . À entrada do domínio da bomba
centrífuga, é imposta uma pressão de referência constante e igual a zero (Pref.= 0
Pa) e à saída é imposta vazão constante, sob forma de uma velocidade média
constante em toda área de saída do domínio. As distribuições de k e ε devem ser
conhecidas na entrada, porém, quando não se conhece essa distribuição (maioria
dos casos), faz-se uma estimativa. No programa Ansys CFX, há diversas maneiras
de como isso pode ser feito. Uma dessas maneiras usa o conceito de intensidade
turbulenta (I). O valor de I na entrada é normalmente assumido entre 0 e 0,10. O
31
valor de κ é calculado através da relação κ = ( 3 2 ) I2 V
ε = ρCµ κ 2 (1000 I µ ) .
Nas
interfaces
dos
domínios
2
de
e ε é calculado por
cálculo
(tubo/rotor,
rotor/difusor, etc) são impostos modelos de interface para possibilitar a conexão de
um domínio ao outro (detalhes desses modelos são mostrados no capítulo 5).
Entrada Fluido
Pressão Ref.
(P = 0 Pa)
Superfícies
Sólidas
(u = v = w = κ = 0 )
Interface
dos Domínios
(Modelos Interface)
Rotor
Difusor
ω=cte
Saída de Fluido
Vazão Constante
(Velocidade Média)
Extensão Difusor
Figura 3.4 – Esquema das condições de contorno impostas aos sub-domínios do
primeiro estágio da bomba estudada.
O programa computacional utilizado nas simulações numéricas possui uma
interface gráfica onde são inseridos os dados de entrada tais como: geometria do
problema, a malha computacional, parâmetros do escoamento, modelo de
turbulência utilizado, etc. O capítulo seguinte descreve como foram confeccionados
os domínios de interesse e atenta para algumas simplificações feitas nas geometrias
da bomba estudada.
32
4 MODELAGEM GEOMÉTRICA E NUMÉRICA
Uma vez estabelecidas as hipóteses e as equações governantes que regem o
problema estudado, os domínios geométrico e numérico, que no presente trabalho
compreende o primeiro rotor e o difusor da bomba centrífuga Imbil ITAP 65-330/2,
são definidos. Contudo, apesar do foco dos estudos estar no rotor e no difusor, dois
domínios auxiliares são utilizados nas simulações numéricas, um à entrada do rotor
e outro à saída do difusor, com o intuito de minimizar efeitos das condições de
contorno nos domínios de interesse.
Na primeira parte deste capítulo é descrita a construção dos domínios de
cálculo do problema. Em seguida apresenta-se a metodologia empregada pelo
programa computacional Ansys CFX 11.0, utilizado na solução numérica do
problema.
4.1 Definição das Geometrias do primeiro estágio da bomba
centrífuga
Antes de se implementar a geometria desejada no programa computacional, é
necessário proceder alguns passos quanto à modelagem mecânica das partes
envolvidas (rotor, difusor) e também das extensões do domínio à entrada e saída.
Uma descrição detalha é feita a seguir.
4.1.1 Geometria do rotor do primeiro estágio da bomba centrífuga
O rotor estudado é radial do tipo fechado, isto é, possui tampas que envolvem
ambos os lados das pás. A essas tampas dão-se os nomes de cubo e coroa. O cubo
é a região inferior do rotor e que está em contato direto com o eixo que transmite
33
movimento. A coroa localiza-se na região oposta ao cubo, como mostra a Figura 4.1.
As pás do rotor estudado são voltadas para trás em relação ao sentido de giro do
rotor e são em número de oito. O primeiro rotor da bomba possui diâmetro de
entrada de 80,1 mm, diâmetro de saída de 206,6 mm e altura do canal na saída de
12,0 mm.
A Figura 4.1 mostra o primeiro rotor da bomba centrífuga Imbil ITAP 65-330/2
modelado em realidade virtual com o auxílio do programa computacional de
SolidWorks 2007. Para fins ilustrativos, a Figura 4.1(a) mostra o rotor com as duas
tampas (cubo e coroa) e a Figura 4.1(b) apresenta o mesmo rotor sem a tampa
superior (coroa), onde é possível visualizar as pás e o cubo.
Pás
Coroa
Cubo
(a)
(b)
Figura 4.1 – Rotor do primeiro estágio da bomba IMBIL ITAP 65-330/2. (a) Rotor na
sua forma original com coroa e (b) rotor sem a coroa superior para facilitar a
visualização de partes internas.
Uma réplica do rotor utilizado nos testes experimentais foi obtida junto à
Universidade Estadual de Campinas (UNICAMP). Nesta réplica, foram efetuadas
diversas medições com intuito de determinar a geometria das partes internas do
rotor. A Figura 4.2 mostra o rotor do primeiro estágio da bomba, utilizado nos
estudos numéricos.
34
Figura 4.2 – Rotor real sem a tampa inferior (cubo). Rotor utilizado na determinação
da curvatura das pás e confecção do modelo virtual.
Após uma modelagem geométrica preliminar do rotor, diversos ajustes foram
efetuados para que o rotor fosse o mais fidedigno possível ao modelo real. Isso foi
feito devido à grande complexidade da geometria das pás do rotor. A geometria das
partes sólidas do rotor (pás, cubo e coroa) foi utilizada para obtenção do domínio
fluido de cálculo, mostrado na Figura 4.3. Esse domínio foi utilizado na geração da
malha computacional que finalmente foi fornecida ao programa de dinâmica dos
fluidos computacional Ansys CFX 11.0 para a simulação numérica do escoamento.
Figura 4.3 – Domínio fluido de cálculo do rotor
Uma aproximação feita no domínio numérico do rotor foi considerar as pás
menores que o diâmetro externo do rotor, como mostra a Figura 4.4. Isso foi feito
com objetivo de melhorar o acoplamento entre o domínio de rotor e difusor. O
manual do programa Ansys CFX sugere que o valor numérico da área na secção de
35
saída do rotor seja igual ao da área de secção de entrada do difusor. O programa
tolera pequenas variações, porém, campos distorcidos podem ser gerados caso haja
diferenças muito grandes na área de contato em rotor e difusor. Para evitar esse
problema, o rotor foi confeccionado com as pontas das pás 2mm menores, em
diâmetro, do que o modelo real. Procedimento análogo foi feito ao difusor para a
área de contato entre os dois domínios fosse igual.
Figura 4.4 – Detalhe das pás do rotor e aletas do difusor. Tanto as pás do rotor
quanto as aletas do difusor não tocam a superfície de contato entre os domínios
4.1.2 Geometria do difusor do primeiro estágio da bomba centrífuga
Para a modelagem do difusor, localizado à saída do primeiro rotor da bomba,
são utilizadas técnicas de medição semelhantes às aplicadas na modelagem do
rotor. O difusor é composto por doze aletas, que tem a função principal de organizar
o escoamento na saída do rotor. As aletas possuem 16,1mm de altura por 3,0mm de
espessura. A Figura 4.5 mostra uma réplica do difusor utilizado nos testes
experimentais.
36
Figura 4.5 – Difusor do primeiro estágio da bomba centrífuga Imbil ITAP 65-330/2
Em comparação ao rotor, o difusor não apresenta grandes dificuldades para a
determinação de sua geometria. O difusor mostrado na Figura 4.5 é uma peça única
confeccionada por fundição. Por isso as aletas do difusor possuem variações
dimensionais entre si. Isso não foi levado em consideração no modelo virtual, ou
seja, foram tomados valores médios das medições. Por exemplo, a espessura e o
comprimento das aletas foram escolhidos de modo que todos os canais do modelo
virtual fossem exatamente iguais. A Figura 4.6 mostra o resultado dessa
modelagem.
Figura 4.6 – Modelo virtual do difusor da bomba centrifuga Imbil ITAP 65 330/2
Analogamente ao que foi realizado para determinação do domínio fluido do
rotor, um domínio fluido do difusor foi confeccionado, para também ser utilizado
como dado de entrada no programa Ansys CFX 11.0. O resultado final desse
tratamento é mostrado na Figura 4.7.
37
Figura 4.7 – Domínio fluido do difusor da bomba centrífuga Imbil ITAP 65-330/2
O primeiro rotor e o difusor da bomba centrífuga Imbil ITAP 65 330/2
constituem o foco de estudos do presente trabalho. Porém, para eliminar efeitos de
condições de contorno nesses componentes, pequenas extensões nos domínios à
entrada e à saída foram implementadas. Detalhes dessa técnica estão mostrados na
secção 4.1.3, a seguir.
4.1.3 Geometrias auxiliares utilizadas nos estudos numéricos
Como foi visto no Capítulo 2, quando a bomba centrífuga está operando fora
da condição de projeto, surgem recirculações no interior do rotor e difusor. No caso
de haver recirculações em regiões próximas às fronteiras do domínio (entradas ou
saídas), pequenas extensões do domínio devem ser aplicadas na modelagem.
García e Vicent (2007) realizaram estudos numéricos preliminares na bomba
centrífuga Imbil ITAP 65 330/2 e mostraram que os resultados são melhorados com
a aplicação de um tubo de pequeno comprimento na região de entrada da bomba.
Essa técnica também foi utilizada por Mañez e Vicent (2009) que além de um tubo
de entrada, utilizou uma pequena extensão à saída do difusor.
Em concordância aos trabalhos realizados por García e Vicent (2007) e
Mañez e Vicent (2009), dois domínios foram implementados. O primeiro trata-se de
um tubo com diâmetro idêntico a entrada do primeiro rotor e um comprimento de
100,0 mm e o segundo é uma extensão da saída do difusor.
38
No caso da extensão à saída do difusor, foram necessários alguns estudos
para se determinar a extensão adequada do domínio. Na saída do difusor real, há
uma peça denominada de Indutor. O indutor é composto de seis aletas e tem a
função de preparar o escoamento para o segundo rotor. Como já mencionado, o
objetivo do estudo do presente trabalho é avaliar o escoamento no primeiro rotor e
no difusor. Por isso, foi confeccionada uma peça simplificada que se assemelha à
forma do indutor, porém que não possui aletas, como mostra a Figura 4.8. Uma
análise da escolha da peça que melhor se adequasse à realidade do escoamento à
saída do difusor está mostrada na secção 5.1.
Saída do
Indutor
Entrada do
Indutor
Figura 4.8 – Extensão do domínio à saída do difusor.
4.2 Discretização das equações de conservação
Nessa secção é mostrado como as equações de conservação da massa,
equação da quantidade de movimento e as provenientes dos modelos de turbulência
são tratadas numericamente.
A metodologia numérica utilizada consiste na divisão dos domínios de
interesse em diversos pequenos volumes de controle, gerando uma malha
computacional. Em cada um desses volumes de controle, são aplicadas as
39
equações que regem o escoamento na bomba (conservação da massa, quantidade
de movimento e modelo de turbulência). Essas equações são linearizadas e
integradas ao longo de cada um dos volumes de controle do domínio, no tempo e no
espaço, resultando em um sistema algébrico de equações.
Há, basicamente, quatro tipos de malhas computacionais que podem ser
utilizadas para mapear problemas numéricos de mecânica dos fluidos: malha
cartesiana, malha cilíndrico-polar, malhas ajustadas ao corpo e malhas nãoestruturadas, como mostrado na Figura 4.9. As malhas cartesianas, cilíndricopolares e ajustadas ao corpo são classificadas como malhas-estruturadas, pois
todos os volumes de controle internos ao domínio possuem a mesma quantidade de
vizinhos. Nas malhas cartesianas (Figura 4.9(a)), os volumes de controle são
paralelepípedos, sendo as variáveis espaciais dadas por x, y e z, orientadas nas
direções das três arestas desses volumes de controle. Nas malhas cilíndrico-polares
(Figura 4.9(b)), os volumes de controle são prismas retos de base formada por um
setor circular, sendo as variáveis espaciais dadas por r (raio), θ (direção angular) e
z (direção axial). Malhas ajustadas ao corpo (Figura 4.9(c)) possuem volumes de
controle com formatos análogos a prismas, porém, suas arestas não são formadas
por retas, mas sim curvas que se ajustam ao formato da geometria estudada. As
variáveis espaciais ξ 1 , ξ 2 e ξ 3 mapeiam o domínio com malhas ajustadas ao
corpo e são obtidas pela teoria de transformação de coordenadas. As malhas nãoestruturadas, mostrada na Figura 4.9(d) possuem volumes de controles com
orientações variadas e não há uma lei de formação definida para esses volumes,
como ocorre para os demais tipos de malhas.
40
Figura 4.9 – Tipos de malhas computacionais aplicadas em mecânica dos fluidos e
transferência de calor e massa. (a) Malhas Cartesianas, (b) Malhas CilíndricoPolares, (c) Malhas Ajustadas ao Corpo e (d) Malha não-estruturada.
No caso apresentado pela Figura 4.9(d), devido as particularidades da
geometria, não há garantia de que os volumes de controle possuam uma quantidade
fixa de vizinhos, e também não há o conceito de orientação do volume de controle
em determinada direção coordenada. A escolha do tipo de sistema de coordenadas
mais adequado depende da forma geométrica do problema estudado. O programa
computacional Ansys CFD possui um gerador de malha denominado Ansys ICEM
CFD que é capaz de gerar tanto malhas estruturadas como não-estruturadas.
O programa computacional Ansys CFX 11.0 utiliza o método de discretização
dos Volumes Finitos baseado em Elementos Finitos. Esse método consiste em
mapear a geometria com elementos tetraédricos, hexagédricos, prismáticos ou
piramidais (vide Figura 4.10).
A Figura 4.11 mostra de maneira esquemática um desses volumes de
controle da malha gerada. Nos vértices ou nós de cada elemento gerado são
armazenadas todas as variáveis do problema e propriedades do fluido. No centróide
das faces de cada elemento é discriminado um ponto, denominado Elemento Face-
41
Centrado. Ao redor de cada nó é construído um volume de controle unindo os
Elementos Face-centrados ao ponto médio da aresta de cada elemento
circunvizinho ao nó considerado.
(a)
(b)
(c)
(d)
Figura 4.10 – Tipos de elementos finitos utilizados na confecção de uma malha
computacional não-estruturada. (a) Tetraedro, (b) Hexaedro, (c) prisma triangular e
(d) pirâmide.
Elemento
Superfície do Volume de Controle
Nó
Figura 4.11 – Volume de controle de uma malha não-estruturada
42
As equações de conservação, apresentadas no capítulo 3, são integradas ao
longo do volume de controle “hachurado” da Figura 4.11 e o teorema de divergência
de Gauss é aplicado para converter algumas integrais de volume em integrais de
superfície. Assim, as equações resultantes assumem a seguinte forma:
∫ (V dn ) = 0
j
j
(4.1)
SC
∂
∂t
 ∂V ∂V j 
= − ∫ Pdn j + ∫ µeff  i +
 dn j + ∫ S dV
 ∂x
∂
x
j
i
VC
SC
SC
SC
VC


 ∂φ 
∂
ρφdV + ∫ ρV j φ ⋅ dn j = − ∫ Pdn j + ∫ Γeff 
dn + S dV
∫
 ∂x  j ∫ φ
∂t VC
j
SC
SC
SC
VC


∫ ρV dV + ∫ ρV V ⋅ dn
j
j
i
j
(4.2)
(4.3)
onde VC e SC denota, respectivamente, integração no volume de controle e na
superfície de controle, dn j são diferenciais das componentes cartesianas do vetor
normal de área que aponta para a fora da superfície de controle.
O primeiro passo na solução numérica dessas equações diferenciais é criar
um sistema acoplado de equações lineares algébricas. Isso é feito convertendo cada
termo das equações (4.1) a (4.3) na forma discreta. A Figura 4.12 mostra uma das
faces de um elemento da malha e evidencia os denominados Pontos de Integração,
Pin, que estão localizados na fronteira entre dois volumes de controle, no centro de
cada segmento que compõem a face que circunda o nó.
43
Nó2
Nó1
Ponto Integração
Pi1
Elemento face-centrado
Setores
Pi3
Pi2
Nó3
Figura 4.12 – Elemento de malha
Da Figura 4.12, para integrar os termos volumétricos da Equação (4.2) e (4.3)
de um nó são contabilizadas as contribuições de cada setor a que esse nó está
circunscrito. Os termos de fluxo são convertidos na sua forma discretizada
primeiramente aproximando seus fluxos através dos Pontos de Integração e
posteriormente contabilizados ao nó a que estão circunscritos. A forma discretizada
das equações de conservação é:
∑ (V ∆n )
j
(V − V ) +
ρV
0
i
i
∆t
Pi
j Pi
=0

ɺ
m
V
=
P
∆
n
+
 µeff
(
)
(
)
∑
∑
∑
Pi
i Pi
i Pi

Pi
Pi
Pi

(φ
ρV
i
− φi 0
∆t
)+
∑ mɺ ( φ )
Pi
Pi
i Pi
(4.4)
 ∂Vi ∂V j
+

∂
x
∂xi
j



 ∆n j  + SVi V (4.5)


Pi


∂φ
= ∑  Γeff i ∆n j  + Sφ V


∂x j
Pi 
Pi
(4.6)
onde V é o volume de controle, e ∆t é o passo de tempo, ∆n j é componente
discretizada do vetor de área da superfície de controle, o sob-escrito “ Pi ” denota a
avaliação da variável no ponto de integração. O termo temporal é discretizado com
um esquema de interpolação de primeira ordem que utiliza o resultado da variável no
passo de tempo anterior (sobre-escrito 0). A vazão mássica discretizada através de
44
ɺ Pi nas equações (4.5) e (4.6), é dada
uma superfície de controle, simbolizada por m
por:
ɺ Pi = ( ρV j ∆n j )
m
Pi
(4.7)
Para resolver as equações (4.4) a (4.6) é necessário determinar o valor das
propriedades do escoamento fora dos nós, onde estão armazenados os valores das
variáveis. Para isto, são utilizadas as funções de forma do Método dos Elementos
Finitos para calcular o valor de uma variável φ no interior de um elemento da
seguinte forma:
φ=
Nnós
∑N φ
i
i =1
i
(4.8)
onde Ni é uma função de forma para o nó “i” e φi é o valor da variável φ no nó “i”. O
comportamento da função de forma Ni é tal que:
N=
Nnós
∑N
i =1
i
=1
(4.9)
De acordo com a Equação (4.9), as funções de forma “ N” foram concebidas
para que nos vértices dos elementos o valor da variável φ seja exatamente o valor
φi daquele nó. Pode-se concluir que a equação (4.8) representa uma interpolação
de todos os vértices do elemento em relação ao ponto interno onde se está
determinando a propriedade desejada. Essas funções de forma são escritas em
termos de variáveis paramétricas s, t e u que assumem valores reais entre 0 e 1.
Cada tipo de elemento (tetraédrico, hexaédrico, prismático ou triangular) possui um
conjunto de funções de forma específico para a interpolação de φ no interior
daquele elemento. Um exemplo dessas funções de forma para o elemento
tetraédrico é mostrado na Figura 4.13.
45
Nó4
Nó1
Funções de forma
(tetraédros):
Nó2
Nó3
N1(s,t,u)=1-s-t-u
N2(s,t,u)=s
N3(s,t,u)=t
N4(s,t,u)=u
Figura 4.13 – Exemplo de funções de forma utilizadas em elementos tetraédricos
para ponderação de valores internos.
Com auxílio dessas funções é possível calcular o valor de qualquer variável
dentro de uma posição arbitrária do elemento considerado, inclusive em termos de
gradientes. Dessa maneira, os termos de gradiente de pressão e os gradientes
difusivos da Equação (4.5) e (4.6) são determinados.
Para completar a discretização das equações de conservação, é necessário
avaliar o valor assumido pelos termos advectivos. O programa Ansys CFX possui um
esquema de interpolação de Alta Ordem, que pondera escolha do valor de φ para
cada ponto de integração, utilizando princípio análogo ao esquema híbrido de
interpolação descrito por Patankar (1980), porém, com a adição de termos e funções
interpoladoras de alta ordem. Esse esquema foi escolhido para as simulações
numéricas do presente trabalho. O programa computacional resolve o sistema de
algébrico de equações através de um método de solução proposto por Rhie e Chow
(1983). A equação da conservação da massa é modificada, estabelecendo uma
equação para a pressão, que é resolvida de maneira implícita com a equação da
quantidade de movimento. Cada nó tem um sistema de equações do tipo: (Maliska,
2009)
∑a
i
viz
onde
viz
φi viz = bi
(4.10)
46
aiviz
 auu
a
vu
=
awu

 apu
φviz
i
auv
avv
awv
apv
aup 
avp 
awp 

app 
auw
avw
aww
apw
u 
v 
= 
w 
 
p
(4.11),
viz
(4.12)
e
 bu 
b 
v
bi =  
(4.13)
 bw 
 
 bp  i
φ é a variável calculada, os “a” são os coeficientes oriundos das equações de
conservação discretizadas, i representa o número do nó considerado (cada nó é
identificado por um número inteiro que o discrimina dos demais nós da malha
computacional), viz representa a contribuição dos vizinhos ao nó considerado e
também a contribuição do próprio nó. Todas as matrizes dos nós são resolvidas de
forma acoplada e simultânea através de um sistema da forma:
[ai ]

[ ]






[ ]




[] []
[] [] []
[]
[]
[]
.
.
.
.
[] [] []
.
.
[ φ ] [ bi ]

 i  
 [ ]  [ ] 

 [ ] 

 [ ]  


. 


 .  

 . 

. .  = 




.

 .  


 


 

.
 


 

.

 


(4.14)
47
INÍCIO
AVANÇ AR PASSO
DE TEMPO
RESOLVER AVANÇO DA MALHA
RESOLVER SISTEMA HIDRODINÂMICO
RESOLVER EQUAÇÕES DE TURBULÊNCIA
SIM
NÃO
TRANSIENTE
NÃO
NÃO
TEMPO FÍSICO
MÁXIMO ATINGIDO
SIM
NÃO
CRITÉRIO DE PARADA /
N° ITERAÇÕES ATINGIDO
CRIT. PARADA / N° ITERAÇÕES
ATINGIDOS
SIM
FIM
SIM
Figura 4.14 – Fluxograma simplificado de resolução do sistema de equações do
programa Ansys CFX 11.0
Os coeficientes “a” das equações de algébricas (4.10) a (4.14) são calculados
através de parâmetros da malha (tamanho dos volumes de controle, etc) e do valor
variáveis do escoamento (pressão, velocidade, temperatura, etc). Quando o
programa computacional resolve a matriz dada pela equação (4.14), novos valores
dessas variáveis são obtidos e então, os coeficientes “a” são atualizados e o
processo iterativo se reinicia até que a convergência seja alcançada.
A Figura 4.14 mostra de maneira esquemática e simplificada o fluxograma de
solução das equações no programa Ansys CFX. Tendo definido o problema a ser
estudado numericamente (geometria, malha computacional, etc) e com os dados de
entrada, inseridos no programa, a simulação numérica é iniciada. O programa lê um
campo inicial das variáveis a serem calculadas (estimativa fornecida pelo usuário ou
gerada pelo próprio programa) e segue uma série de testes lógicos até que a
convergência do problema seja alcançada, como mostra a Figura 4.14. A seguir é
mostrado como foi definida geometria a ser simulada numericamente no presente
trabalho, bem como são mostrados os diversos parâmetros numéricos que foram
determinados como dados de entrada nas simulações (critério de parada, número de
nós da malha computacional, etc).
48
5 AVALIAÇÃO DE PARÂMETROS E TESTES DE MALHA
Neste capítulo, é descrita a metodologia utilizada para determinação dos
parâmetros utilizados nas simulações numéricas. Mostra-se como foi definida a
geometria à saída do difusor, o teste de critério de convergência, o teste de malha, o
teste do passo de tempo e o teste da influência do arquivo iniciador nas simulações
em regime transiente.
5.1 A escolha da configuração geométrica adequada para as
simulações numéricas
Para a escolha da geometria do domínio a ser anexada a saída do difusor,
foram propostas quatro configurações diferentes e alguns testes numéricos foram
realizados para observar a influência de cada uma delas nos resultados. A primeira
configuração, mostrada na Figura 5.1, simula o escoamento no rotor e no difusor da
bomba, sem a presença de extensão de domínio à saída do difusor, impondo
pressão de referência igual a zero na entrada do tubo de admissão e vazão na
carcaça do difusor (sob forma de um perfil uniforme de velocidades ao longo da área
de saída), como se a direção do escoamento à saída do difusor estivesse em um
plano transversal ao eixo de rotação do rotor.
49
Pref = 0Pa
Tubo de entrada
Rotor
Difusor
Vazão
Constante
Velocidade
Constante
Figura 5.1 – Configuração 01 (“Radial”) de escoamento entre rotor e difusor da
bomba centrífuga IMBIL, modelo ITAP 65-330/2
A segunda configuração é semelhante à primeira, porém a saída do difusor
tem direção axial, como mostrada na Figura 5.2. Como conseqüência dessa saída
axial, no raio externo do difusor (carcaça), foi considerado o escoamento junto a uma
parede sólida com condição de não-deslizamento.
Pref = 0Pa
Tubo de entrada
Rotor
Difusor
Parede com condição
de não-deslisamento
Vazão
Constante
Velocidade
Constante
Figura 5.2 – Configuração 02 (“Axial para baixo”) de escoamento entre rotor e difusor
da bomba centrífuga IMBIL, modelo ITAP 65-330/2
A terceira configuração, mostrada na Figura 5.3, é uma extensão do domínio
do difusor com objetivo semelhante ao tubo de admissão do rotor, ou seja, afastar a
condição de contorno da região de interesse. Como a bomba é composta por dois
estágios, um dimensionamento preliminar do indutor (região localizada entre o
50
difusor e o segundo rotor) foi feito e as medidas principais (como altura do indutor e
diâmetro externo) foram utilizadas para confecção dessa extensão de domínio.
A quarta e última configuração é uma variante da configuração 3, como
mostrada na Figura 5.4. A diferença entre essas duas configurações está no
diâmetro de saída do escoamento, ou seja, na configuração 4 foi utilizado um
diâmetro de saída menor do que na configuração 3. No indutor da bomba centrífuga
estudada, há uma redução de área na região superior do domínio (ver Figura 5.4)
que também foi determinada e implementada no modelo numérico. No indutor real,
há também um conjunto de 6 aletas igualmente espaçadas ao longo da direção
tangencial, entretanto, essas aletas não estão contabilizadas na configuração 4. A
função dessa extensão continua sendo a de afastar a condição de contorno de saída
da região de interesse.
Pref = 0Pa
Tubo de entrada
Rotor
Difusor
Parede com condição
de não-deslisamento
Velocidade
Constante
Vazão Constante
Extensão Difusor
Figura 5.3 – Configuração 03 (“Extendido”) de escoamento entre rotor e difusor da
bomba centrífuga IMBIL, modelo ITAP 65-330/2
51
Pref = 0Pa
Tubo de entrada
Rotor
Difusor
Parede com condição
de não-deslisamento
Velocidade
Vazão
Constante
Extensão Difusor
Figura 5.4 – Configuração 4 (“Indutor sem aletas”) de escoamento entre rotor e
difusor da bomba centrífuga IMBIL, modelo ITAP 65-330/2
Nas Figuras 5.1 a 5.4, pode-se observar a presença de vários sub-domínios
distintos, que juntos compõem o domínio completo a ser simulado numericamente.
Além das condições de contorno comumente utilizadas como vazão, pressão de
referência e condição de não deslizamento em superfícies sólidas, é necessário
informar ao programa computacional o acoplamento entre os domínios (tubo de
entrada, rotor, difusor).
No problema estudado existem dois tipos de interfaces a serem analisados. A
primeira trata-se de um domínio fixo em contato com outro domínio também fixo,
como é a interface da saída do difusor com a entrada do indutor. Para esse caso, o
programa Ansys CFX oferece um modelo simples que comunica ambos domínios
sem maiores dificuldades. O segundo caso diz respeito a um domínio em movimento
de rotação em relação a outro.
Para interfaces de domínios onde há movimento relativo de rotação, o
programa computacional Ansys CFX 11.0 possui três modelos distintos: O modelo
Estágio ou Segmento (do inglês, “Stage”), O modelo Rotor-Congelado (do Inglês,
“Frozen-Rotor”) e o modelo Transiente.
O modelo Estágio e o modelo Rotor-Congelado são modelos de regime
permanente, isto é, o rotor assume uma posição fixa em relação ao difusor. Nesses
52
dois modelos de interface, não há variação da posição do rotor em relação ao
estator, ou seja, o domínio numérico que compreende o rotor não altera sua posição
angular com o tempo, gerando um campo médio de pressão e velocidade. A
diferença básica entre esses dois modelos de interface está na maneira com que
ambos tratam a variável pressão na interface dos domínios. O modelo Estágio aplica
uma média circunferencial da pressão na interface dos domínios de modo que não
haja variação da pressão na direção tangencial da interface dos domínios, o que não
é feito no modelo Rotor-Congelado. Como conseqüência, o modelo Estágio não é
indicado para situações onde ocorram recirculações no interior do domínio.
O modelo Transiente leva em conta o efeito da variação da posição angular
do rotor em relação ao difusor. A desvantagem desse modelo frente aos dois
primeiros modelos citados está no alto custo computacional e também no grande
volume de dados gerados nas simulações numéricas.
Para verificar o efeito das quatro configurações mostradas nas Figuras 5.1 a
5.4 nos resultados foi utilizado o modelo de interface Rotor-Congelado. Uma malha
computacional grosseira foi construída para realização dos testes. Apesar disso,
espera-se uma certa concordância com os resultados experimentais, uma vez que
as geometrias de rotor e difusor, bem como a malha gerada (ainda que não seja
refinada), são próximas da geometria do modelo real. A tabela 5-1 mostra os
resultados de diferença de pressão no rotor e difusor para as várias configurações,
os resultados experimentais e os desvios em relação aos valores experimentais,
para a vazão de projeto, 36m3/h, e rotação nominal, 1150rpm.
53
Tabela 5-1 – Tabela comparativa da diferença de pressão no rotor e difusor da
bomba centrífuga estudada para 4 configurações de saída do escoamento. vazão de
36m3/h de água
Config.
01
Geometria
Experimental
∆PRotor
Desvio (%)
∆PDifusor
Desvio (%)
[Pa]
∆PRotor
[Pa]
∆PDifusor
62549,6
-
17647,3
-
0,85
9809,6
44,41
(Amaral, 2007)
Saída Radial 62018,4
02
Axial p/ baixo
59352,0
5,11
9964,1
43,54
03
Extendido
61404,1
1,83
15534,1
11,97
Indutor s/ aletas 60785,2
2,82
14898,9
15,57
04
O valor dos desvios relativos foram calculados utilizando a Equação (5.1)
 (Valor Numérico ) − (Valor Experimental ) 
Desvio (%) = Absoluto
 ⋅ 100%
(Valor Experimental )


(5.1)
Observa-se na Tabela 5-1 que as variações dos resultados para a diferença
de pressão no rotor são inferiores a 6,0 %. Já para a diferença de pressão no difusor
houve maiores discrepâncias. As configurações 01 e 02 tiveram variações muito
pronunciadas em relação ao valor experimental além de apresentarem alto custo
computacional. Várias intervenções por meio de alteração de parâmetros das
simulações foram necessárias para obtenção da convergência. Por essas razões,
estas configurações foram desconsideradas. Para escolha da geometria mais
indicada, foi feita uma simulação numérica adicional com uma vazão diferente da de
projeto, 46m3/h, também com rotação de 1150 rpm. Para essa vazão os resultados
estão apresentados na Tabela 5-2, para as configurações 03 e 04.
54
Tabela 5-2 – Tabela comparativa da diferença de pressão no rotor e no difusor para
as configurações 3 e 4. Vazão de 46m3/h de água
Config.
-
∆PRotor
Desvio
∆PDifusor
Desvio (%)
[Pa]
(%) ∆PRotor
[Pa]
∆PDifusor
57538,4
-
14898,6
-
57606,1
0,12
12523,1
15,94
Indutor s/ aletas 57290,8
0,43
12936,9
13,17
Geometria
Experimental
(Amaral, 2007)
03
04
Extendido
Analisando as tabelas 5-1 e 5-2, nota-se que os resultados entre as
configurações 3 e 4 são muito semelhantes. No presente trabalho, foi escolhida a
configuração 4, pois a condição de saída está mais longe da região de interesse
(rotor e difusor) em relação a configuração 3.
Uma vez defina a geometria a ser utilizada nas simulações numéricas (ver
Figura 5.5) constituída por tubo de admissão, rotor do primeiro estágio, difusor e a
extensão à saída do difusor, diversos outros parâmetros foram avaliados de maneira
a verificar a consistência dos resultados. Nas secções seguintes, avalia-se o critério
de parada, a malha numérica utilizada e o passo de tempo.
Figura 5.5 – Modelo numérico da bomba centrífuga Imbil ITAP 65-330/2
55
5.2 Critério de parada
O critério de parada é um parâmetro utilizado para avaliar a convergência
alcançada nas simulações numéricas. O critério de parada indica a magnitude dos
resíduos obtidos na solução das equações discretizadas. A literatura sugere utilizar
como critério de parada o valor máximo dos resíduos (RESMAX) das variáveis
calculadas (u, v, w e P) como sendo menor que 5.10-4 (Asuaje et al, 2005).
Para avaliar o critério de parada adequado, foram feitos testes da bomba
centrífuga operando na vazão de projeto fornecida pelo fabricante (vazão de 36m3/h
na rotação de 1150 rpm). A Tabela 5-3 mostra os resultados para diferentes critérios
de parada.
Tabela 5-3 – Comparativo de resultados da diferença de pressão no rotor e difusor
versus critério de parada do programa Ansys CFX 11.0
∆PRotor
∆PDifusor
Desvio %
Desvio %
-1
[Pa]
46944,6
[Pa]
13381,4
(Rotor)
-
(Difusor)
-
1x10-2
60951,3
15177,7
22,98
11,84
1x10-3
60768,4
15026,1
0,30
1,01
1x10-4
60774,1
15029,9
0,01
0,03
RESMAX
1x10
Os desvios percentuais que aparecem na Tabela 5-3 são calculados de
acordo com a Equação (5.2). Essa equação compara dois critérios de parada
consecutivos, por exemplo, 1x10-1 (i) e 1x10-2 (i+1).
(
) (
)
 ∆P i − ∆P i +1 
 ⋅ 100%
(5.2)
Desvio (%) = Absoluto
∆P i +1


Da tabela 5-3, pode-se notar que a variação entre o valor da diferença de
(
)
pressão tanto no rotor, quanto no difusor não sofreu alteração significativa a partir da
utilização de um RESMAX de 1x10-3. Tal valor foi adotado como sendo o critério de
parada para as simulações do presente trabalho.
Para avaliar se esse valor de critério de parada é realmente aceitável, foram
levantados perfis de velocidade e pressão em algumas regiões do domínio da
56
bomba. A Figura 5.6 mostra o perfil de velocidade relativa no interior de um dos
canais do rotor em uma posição arbitrária. No eixo das ordenadas, a cota z
corresponde à direção axial (paralela ao eixo da bomba) e está adimensionalizada
pela largura do canal na saída do rotor, b. Da Figura 5.6 pode-se observar que
praticamente não há distinção dos perfis de velocidade para RESMAX de 1x10-3 e de
1x10-4, o que se conclui que utilizar um critério de parada em 1x10-3 não trás
prejuízos aos resultados.
1
z/b
0.8
0.6
Res Max 0,1
Res Max 0,01
Res Max 0,001
Res Max 0,0001
0.4
0.2
0
0
1
2
3
Velocidade Relativa [m/s]
4
Figura 5.6 – Perfil de velocidade relativa em um dos canais da bomba centrífuga
estudada. Comparativo de diferentes critérios de parada
A Figura 5.7 apresenta a variação de pressão com a direção tangencial, no
raio de saída do rotor. Nessa figura é possível notar que a partir da utilização de um
critério de parada inferior à 1x10-2 não há mudanças significativas nos resultados.
57
80000
Pressão [ Pa ]
70000
60000
Res Max 0,1
Res Max 0,01
Res Max 0,001
Res Max 0,0001
50000
40000
30000
20000
0
30
60
90
Direção Tangencial [ ° ]
120
Figura 5.7 – Variação de pressão na direção tangencial, no raio de saída do rotor.
Comparativo de diferentes critérios de parada
5.3 Teste de Malha
Esta é uma etapa importante na análise do problema para evitar que os
resultados fiquem dependentes da malha utilizada. Primeiramente, foi necessário
adquirir experiência com a utilização de malhas não-estruturadas, uma vez que
algumas particularidades são impostas na confecção dos volumes de controle desse
tipo de malha. A grande influência que o parâmetro de tamanho da aresta dos
volumes tetraédricos exerce no número final de pontos de cálculo do domínio é um
exemplo da sensibilidade adquirida na confecção de malhas não-estruturadas, no
presente trabalho. Pequenas variações nesse parâmetro implicavam na geração de
um número exagerado de volumes no domínio.
O programa utilizado para a geração da malha computacional foi o Ansys
ICEM CFD. Com esse programa é possível gerar malhas tanto estruturadas como
não-estruturadas. No caso das malhas estruturadas, a complexidade da geometria é
determinante para a utilização dessa opção.
58
Os programas comerciais de geração de malhas não-estruturadas trabalham
basicamente com quatro tipos de elementos ou volumes de controle: os elementos
tetraédricos e hexaédricos que têm função de preenchimento e são versáteis a
praticamente qualquer tipo de geometria ou complexidade (cantos vivos, frestas,
furos, etc); os elementos prismáticos que são utilizados como forma de refinar a
malha em regiões próximas a superfícies sólidas e; os elementos piramidais que
formam uma transição entre os elementos prismáticos e tetraédricos em geometrias
muito complexas.
De um modo geral, a quantidade de volumes de controle envolvida em uma
simulação numérica quando se utilizam malhas não-estruturadas é muito maior em
comparação ao que se utiliza para malhas estruturadas. Feng et al. (2007), por
exemplo, em seus estudos utilizaram cerca de 1,5 milhões de elementos em suas
simulações numéricas. A geometria utilizada por eles é semelhante à utilizada no
presente trabalho. Portanto, esse valor foi tomado como referência na quantidade de
volumes de controle do domínio. Também se observa que em termos de altura de
elevação ou diferença de pressão entre entrada e saída da bomba, malhas
grosseiras produzem resultados bem semelhantes a malhas mais refinadas.
Figura 5.8 – Corte transversal da malha computacional utilizada no rotor e difusor.
Detalhe ampliado mostra a distribuição dos prismas ao longo das pás do rotor e
aletas do difusor.
59
A Figura 5.8 mostra um plano transversal ao eixo de rotação da bomba,
evidenciando o formato dos volumes de controle e a disposição da malha ao longo
do domínio do rotor e do difusor. É possível notar a presença de camadas de
prismas próximo às pás do rotor e, também, das aletas do difusor.
No tubo de admissão da bomba centrífuga e no domínio localizado a saída do
difusor, por serem geometrias simples, foi possível gerar malhas com relativa
facilidade. A malha gerada nesses dois domínios (ver Figura 5.9), apesar de
aparentar ser estruturada, é também não-estruturada devido à arquitetura numérica
utilizada pelo programa computacional Ansys CFX 11.0.
(a)
(b)
Figura 5.9 – (a) Malha no indutor. (b) Malha no tubo de entrada
O teste de malha foi efetuado de maneira que a quantidade de elementos de
cada refinamento fosse aproximadamente o dobro da malha anterior. Na tabela 5-4,
está mostrada a quantidade de elementos de cada malha testada, bem como as
quantidades de pontos em cada um dos seus componentes (tubo de admissão,
rotor, difusor e indutor):
60
Tabela 5-4 – Quantidade de elementos utilizada nos testes de malha
Elementos (Tetraedros, Prismas, Pirâmides e Hexaedros)
Domínio
Malha Inicial
Malha 01
Malha 02
Malha 03
Tubo
14.535
26.061
55.836
114.240
Rotor
85.616
209.740
461.098
743.831
Difusor
53.104
114.419
203.847
465.353
Indutor
36.100
51.832
51.832
215.232
Total
189.355
402.052
772.613
1.538.656
Para avaliar a independência do problema frente à malha computacional,
foram levantados perfis de velocidade e pressão em algumas posições do rotor e do
difusor e com isso foi possível comparar se a cada refino de malha ocorria alguma
mudança significativa nos resultados obtidos. As Figuras 5.10 e 5.11 mostram,
respectivamente, um perfil de velocidade relativa no rotor e pressão, tomados no raio
de 40,0 mm e na metade da largura axial do canal do rotor. Pode se notar que os
resultados entre as malhas 02 e 03 praticamente não sofrem alteração.
Velocidade Relativa [ m/s ]
5
4
Malha Inicial
Malha 01
Malha 02
Malha 03
3
2
1
0
0
45
90
θ [ °]
135
180
Figura 5.10 – Perfil de velocidade relativa versus direção tangencial, tomado no raio
R=40,0mm e z=6,0mm (metade da largura do canal do rotor)
61
Ainda na Figura 5.11, o valor médio de pressão entre a malha inicial e a
malha 01 sofreu uma alteração de 28%. Essa variação caiu para 3% comparando os
valores médios entre as malhas 01 e 02 e ficou em 0,47% comparando as malhas 02
e 03. Nas Figuras 5.10 e 5.11, é possível notar que os perfis tanto de velocidade
quanto de pressão não possuem o mesmo comportamento de um canal do rotor
para o outro. Isso acontece devido a presença do difusor à saída do rotor e também
da diferença na quantidade de aletas do difusor em relação às pás do rotor, 12
aletas versus 8 pás. Essa configuração sugere uma periodicidade do escoamento a
cada 90°, como pode ser observado nas Figuras 5.10 e 5.11.
Em uma outra região do rotor, foi traçado um perfil de velocidade relativa do
escoamento versus a posição axial, em um raio de 64,41mm e uma posição angular
aproximadamente eqüidistante de duas pás consecutivas do rotor e os resultados
estão mostrados na Figura 5.12. Nessa Figura, o eixo das ordenadas aparece
adimensionalizado pela altura do canal do rotor, b=12,0mm.
15000
Pressão [Pa]
10000
Malha Inicial
Malha 01
Malha 02
Malha 03
5000
0
-5000
0
60
θ[°]
120
180
Figura 5.11 – Perfil de pressão versus direção tangencial, tomado no raio R=40,0mm
e z=6,0mm (metade da altura do canal do rotor)
62
1
z/b[]
0.8
Malha 03
Malha 02
Malha 01
Malha Inicial
0.6
0.4
0.2
0
0
1
2
3
4
5
Velocidade Relativa [m/s]
6
Figura 5.12 – Perfil de velocidade relativa do escoamento no rotor versus a direção
axial, em R=64,41mm e entre duas pás consecutivas do rotor.
O difusor apresentou comportamento semelhante em relação à sensibilidade
da malha. De maneira análoga ao que foi feito no rotor, um perfil de pressão versus
direção tangencial foi avaliado em um raio de 110,0mm (o que corresponde a
aproximadamente metade do comprimento do canal do difusor). Discrepâncias em
termos de pressão média foram ainda menores quando comparadas as diversas
malhas simuladas, mas mostram novamente que os resultados das malhas 02 e 03
estão próximos, (ver Figura 5.13). A pressão média da malha inicial quando
comparada com a malha 01 sofreu uma variação 1,75%. Já na comparação relativa
entre a malha 01 com a malha 02 essa variação é de 1,09% e para as malhas 02 e
03, de 0,12%. Com isso, a malha 02 é satisfatória para avaliações tanto de perfis de
velocidade quanto de pressão, sem que os resultados obtidos sejam dependentes
da malha.
63
76000
Pressão [Pa]
74000
Malha Inicial
Malha 01
Malha 02
Malha 03
72000
70000
68000
66000
0
30
60
90
θ [ °]
120
150
180
Figura 5.13 – Perfil de pressão no difusor versus direção tangencial, tomado no raio
R=110,0mm e z=6,0mm
A independência dos resultados em relação à malha foi testada tomando o
valor da pressão em um ponto qualquer do domínio e comparados os resultados de
todas as malhas. Os resultados estão mostrados na tabela 5-5.
Tabela 5-5 – Valor de pressão em um ponto específico do difusor
Malha
Pressão
Desvio
(%)
Inicial
[Pa]
68142,5
01
69429,1
1,85
02
68544,7
1,29
03
68852,2
0,45
Outra maneira de avaliar essa independência de malha é calculando uma
variável secundária do escoamento no rotor da bomba, como por exemplo, o torque
que o eixo da bomba estará submetido na direção axial.
64
Tabela 5-6 – Valor calculado do torque que o escoamento aplica ao eixo da bomba,,
na vazão de projeto e para as várias malhas propostas.
Torque
Desvio
Inicial
[N.m]
-7,867
relativo
(%)
01
-8,219
4,29
02
-8,120
1,22
03
-8,107
0,17
Malha
Das tabelas 5-5 e 5-6 pode se inferir que a malha 02 é suficiente para os
cálculos de engenharia tais como torque no eixo da bomba, altura de elevação, etc.
Outro ponto a ser salientado é a economia computacional ao se utilizar a malha 02
para as simulações numéricas, principalmente quando a simulação for transiente.
Todas as simulações numéricas realizadas foram feitas em um computador pessoal
com processador Intel Core 2 Quad de 2,4GHz e 3,23GB de memória RAM. Para a
malha Inicial, o tempo de simulação foi de 11 minutos, partindo de um campo inicial
de pressão e velocidade igual a zero. Para a malha 01, o tempo de simulação foi de
17 minutos utilizando como ponto de partida os resultados obtidos na simulação da
Malha Inicial. O tempo de simulação da malha 02 foi de 56 minutos utilizando como
campo inicial o campo da malha 01 e o tempo de simulação para a malha 03 foi de
96 minutos, utilizando o campo da malha 02 com ponto de partida.
5.4 Teste de Passo de Tempo
Em todas as simulações numéricas realizadas nos testes de malha da seção
anterior o modelo de interface Rotor-Congelado do programa computacional Ansys
CFX 11.0 foi utilizado para conectar o rotor ao difusor e ao tubo de entrada. Porém,
uma limitação desse tipo de modelo de interface é não permitir contabilizar em uma
mesma simulação numérica os efeitos da posição relativa entre rotor e difusor
65
(devido à presença de aletas no difusor). Por isso, no presente trabalho optou-se por
utilizar o modelo de interface Transiente, com o objetivo de capturar essas
oscilações. Apesar da operar em regime estacionário, ou seja, com uma vazão e
rotação constantes, espera-se que, devido à presença de um número finito de pás
do rotor e aletas do difusor, a pressão em um ponto qualquer na bomba oscile com o
passar do tempo e que essa oscilação seja periódica, isto é, se repetirá após um
determinado período.
No programa Ansys CFX 11.0, há uma opção para simulação da bomba
centrífuga em regime transiente. Entenda-se aqui regime transiente como sendo a
variação da posição angular do rotor em relação ao difusor em cada instante de
tempo. Portanto, é necessário informar quanto tempo físico se deseja calcular e qual
é o intervalo de tempo mínimo entre dois instantes consecutivos.
Como ponto de partida para a escolha desses parâmetros utilizou-se o valor
da rotação nominal da bomba (1150rpm) e com isso foi calculado o tempo físico
necessário para que o rotor completasse duas voltas completas. A sub-divisão desse
tempo para contabilizar os instantes intermediários foi feita de modo que, em uma
simulação numérica preliminar, o rotor deslocasse 30° a cada in stante de tempo.
Nessa configuração, para cada volta do rotor foram contabilizados 12 instantes
iguais de tempo, armazenados pelo programa para posterior avaliação. O manual do
programa Ansys CFX 11.0 informa que para capturar os efeitos transientes ou
periódicos necessita-se de no mínimo 20 passos de tempo dentro do período.
Para quantificar a magnitude dessa periodicidade do escoamento, a
geometria e a configuração rotor-difusor foi analisada. O rotor possui oito pás e o
difusor possui doze aletas ao longo dos 360°. Sabendo q ue o máximo divisor comum
entre oito e doze é quatro, imagina-se que o escoamento seja periódico a cada 90°,
como foi comentado anteriormente. Portanto é de se esperar que essa primeira
simulação numérica utilizando passos de tempo equivalentes ao giro de 30° não seja
suficiente para representar a periodicidade temporal do escoamento. Nas
simulações seguintes, foram feitos refinos no passo de tempo, de modo que uma
nova simulação tivesse metade do passo de tempo da anterior.
66
As Figuras 5.14 e 5.15 mostram os resultados da diferença de pressão entre
entrada e saída do difusor para os vários casos simulados. Foram utilizados como
condições iniciais das simulações numéricas transientes os resultados obtidos para
regime permanente (com o modelo de interface Rotor-Congelado). Por isso, há uma
influência oriunda dessa condição inicial ao longo da primeira metade de volta do
rotor. Isso ocorre por o escoamento utilizando o modelo Rotor Congelado não
contabilizar os efeitos temporais da simulação transiente.
Os resultados foram separados em duas figuras e mostrados a partir da
primeira volta do rotor para facilitar a visualização. Note que os resultados para os
três primeiros casos (30°, 15° e 7,5°) ainda não são consis tentes, ou seja, há ainda
muita diferença de comportamento das curvas obtidas. Na Figura 5.15 pode-se notar
uma tendência de estabilização dos resultados à medida que passos de tempos
menores foram utilizados.
Diferença Pressão Difusor [Pa]
30000
20000
7,5°
15°
30°
10000
0
1
1.2
1.4
1.6
N° Voltas do Rotor
1.8
2
Figura 5.14 - Teste de malha temporal (Parte 1)
67
Diferença Pressão Difusor [Pa]
30000
7,5°
3,75°
1,88°
0,94°
20000
10000
0
1
1.2
1.4
1.6
N° Voltas do Rotor
1.8
2
Figura 5.15 – Teste de malha Temporal (Parte 2)
Os resultados médios da diferença de pressão no difusor estão mostrados na
Tabela 5-7. Na primeira coluna dessa tabela é mostrada a magnitude do intervalo de
tempo utilizado nas simulações numéricas transientes. Na segunda coluna, tem-se o
ângulo percorrido pelo rotor a cada passo de tempo. Na terceira coluna, pode-se
observar a quantidade de intervalos de tempo necessária para o rotor percorrer duas
voltas. A quarta coluna mostra a quantidade de intervalos de tempo presente a cada
90° (Deslocamento angular onde se imagina uma periodicidade do escoamento com
o tempo). Pode-se observar que os desvios relativos da diferença de pressão média
no difusor para duas linhas consecutivas são menores de 2% a partir de intervalos
de tempo menores do que 1,087x10-3s (giro de 7,5° do rotor). No presente trabalho
escolheu-se o passo de tempo de 0,5435x10-3 (giro de 3,75° do rotor). Nos estudos
de Feng et al. (2007), por exemplo, foram utilizados 290 passos de tempo para cada
volta do rotor o que equivale a aproximadamente 0,8° por ins tante de tempo. Porém,
foi observado, no presente trabalho, que a partir de passos de tempo menores, a
diferença dos resultados é muito pequena e o custo computacional se torna muito
alto. Na prática, para o caso do passo de tempo escolhido, são necessários
aproximadamente 1 dia e 15 horas de cálculo computacional, já para o caso em que
o rotor percorre 0,94º a cada passo de tempo, o tempo de cálculo computacional
68
sobe para 5 dias. Assim, têm-se definido o critério de parada, a malha espacial e o
passo de tempo adequados.
Tabela 5-7 – Comparativo dos resultados médios de diferença de pressão no difusor
da bomba centrífuga estuda para diferentes intervalos de tempo.
∆t x 10-3
N° total
∆t’s a ∆Pdifusor Desvio
de
cada 90° médio relativo
Intervalos
(%)
[Pa]
24
3
16096
-
[s]
Giro do
rotor a cada
∆t [°]
4,3478
30
2,1739
15
48
6
11039
45,81
1,0870
7,5
96
12
13160
16,12
0,5435
3,75
192
24
13405
1,82
0,2717
1,88
384
48
13583
1,32
0,1359
0,94
768
96
13843
1,88
5.5 Influência da condição inicial
Em todas as simulações transientes, foram utilizados os resultados de
campos de velocidade e de pressão em regime permanente obtidos com o modelo
de interface Frozen-Rotor (que é um modelo de regime permanente). Foi observado
que nas simulações utilizando o modelo de interface Transient-Rotor-Stator há uma
certa influência dessa condição inicial nos resultados para os primeiros passos de
tempo das simulações transientes, como mostra a Figura 5.16.
69
90
80
Pressão [kPa]
70
60
50
40
30
20
10
0
0.00
0.50
1.00
1.50
2.00
2.50
3.00
3.50
4.00
t / Período de uma volta
Figura 5.16 – Pressão versus tempo na entrada do difusor do primeiro estágio da
bomba centrífuga Imbil ITAP 65-330/2, na vazão de 36m3/h e rotação de 1150rpm.
Na simulação numérica correspondente à Figura 5.16, foram impostas quatro
voltas completas no rotor para observar a variação dos resultados com relação ao
tempo. O ponto monitorado está localizado na interface de rotor e difusor. Foi
observado que após um giro de cerca de meia volta do rotor essa influência já não é
mais sentida nos resultados. Com isso, pode-se garantir que ao simular duas voltas
do rotor tem-se uma solução periódica e, conseqüentemente, independente da
condição inicial.
As etapas descritas no capítulo 5 são fundamentais para se observar e
controlar as diversas influências existentes em uma simulação numérica. Dessa
análise preliminar são extraídos os principais dados numéricos de entrada que
posteriormente serão utilizados nas simulações numéricas do problema, com
posterior economia de tempo. No capítulo seguinte, de posse dos parâmetros
adequados de malha, passo de tempo e critério de convergência, são mostrados os
resultados obtidos nas simulações numéricas efetuadas no primeiro estágio da
bomba, em 4 rotações distintas e diversas vazões.
70
6 RESULTADOS
Uma vez definida a malha utilizada nas simulações numéricas, bem como o
passo de tempo necessário para capturar os efeitos transientes do escoamento na
bomba, os valores de vazão e rotação do rotor foram estabelecidos para as
simulações numéricas do escoamento no rotor e difusor. A bomba centrífuga
estudada no presente trabalho foi testada experimentalmente por Amaral (2007). Por
isso, para facilitar comparações entre os resultados do presente trabalho e os
obtidos por Amaral (2007), foram utilizadas as mesmas rotações e a mesma faixa de
vazão.
Nesse capítulo estão mostrados os resultados obtidos das simulações
numéricas realizadas tendo como foco o primeiro estágio da bomba centrífuga Imbil
ITAP 65-330/2. São mostrados gráficos de altura de elevação, potência e eficiência
em função da vazão volumétrica para o primeiro estágio. Os campos de pressão e
de velocidade do escoamento no interior do rotor e do difusor para várias vazões
volumétricas
simuladas
são
também
apresentados
juntamente
com
o
comportamento transiente do escoamento no interior da bomba.
Nos trabalhos realizados por Amaral (2007), a bomba centrífuga foi testada
em quatro rotações distintas: 1150rpm (rotação nominal), 1000rpm, 806rpm e
612rpm. Ele levantou experimentalmente as curvas de ganho de pressão versus
vazão volumétrica em cada componente da bomba, isto é, diferença de pressão no
primeiro rotor, no difusor, no indutor, no segundo rotor e na voluta. A faixa de vazão
simulada vai de 10,0 a 50,0m3/h. Para todos os casos mostrados a seguir foi
realizada simulação em regime transiente com o modelo de interface transiente para
conectar o rotor ao difusor e o rotor ao tubo de entrada.
71
6.1 Altura de elevação versus vazão volumétrica
O objetivo principal do presente trabalho é reproduzir o comportamento da
bomba centrífuga operando com água ( ρ = 998kg / m3 e µ = 0,001Pa.s ). A Figura 6.1
mostra o ganho de pressão no primeiro rotor, para as quatro rotações escolhidas e
diversas vazões volumétricas. Pode-se observar na Figura 6.1 que os resultados
numéricos seguem a mesma tendência dos respectivos resultados experimentais. O
desvio médio relativo entre esses resultados numéricos e os valores experimentais é
da ordem de 7,0%.
A Figura 6.2 mostra a diferença de pressão no difusor em função da vazão
volumétrica, para as quatro rotações simuladas. Analisando os resultados mostrados
na Figura 6.2, pode-se notar que a tendência dos resultados numéricos é
semelhante a dos resultados experimentais, entretanto, o desvio médio relativo entre
esses resultados numéricos e os valores experimentais é da ordem de 20%.
∆PRotor [ kPa ]
80
Presente Trabalho - 1150rpm
Amaral (2007) - 1150rpm
Presente Trabalho - 1000rpm
Amaral (2007) - 1000rpm
Presente Trabalho - 806rpm
Amaral (2007) - 806rpm
Presente Trabalho - 612rpm
Amaral (2007) - 612rpm
60
40
20
0
10
20
30
40
Q [ m3 / h ]
50
60
Figura 6.1 – Ganho de pressão no primeiro rotor da bomba centrífuga Imbil ITAP 65330/2. Dados experimentais obtidos por Amaral (2007)
72
20
Presente Trabalho - 1150rpm
Amaral (2007) - 1150rpm
Presente Trabalho - 1000rpm
Amaral (2007) - 1000rpm
Presente Trabalho - 806rpm
Amaral (2007) - 806rpm
Presente Trabalho - 612rpm
Amaral (2007) - 612rpm
∆PDifusor [ kPa ]
16
12
8
4
0
10
20
30
40
Q [ m3 / h ]
50
60
Figura 6.2 – Ganho de pressão no difusor da bomba centrífuga Imbil ITAP 65-330/2.
Dados experimentais obtidos por Amaral (2007)
Para se obter a curva da altura de elevação correspondente ao primeiro
estágio da bomba, basta somar o resultado das curvas das Figuras 6.1 e 6.2 e
utilizar a Equação (6.1) para converter o valor de pressão em termos de altura de
elevação.
Para fins práticos de engenharia, expressar o ganho de pressão sob a forma
de altura de elevação facilita na visualização da capacidade de bombeamento do
equipamento estudado, devido a maior sensibilidade sobre a variável de altura em
relação à variável de pressão.
∆Protor + ∆Pdifusor
(6.1)
ρg
onde g é aceleração da gravidade. A Figura 6.3 mostra o resultado os resultados de
H=
altura de elevação da bomba, tanto numérico quanto experimental. Foi observada
boa concordância entre ambos os resultados.
73
10
Presente Trabalho - 1150rpm
Amaral (2007) - 1150rpm
Presente Trabalho - 1000rpm
Amaral (2007) - 1000rpm
Presente Trabalho - 806rpm
Amaral (2007) - 806rpm
Presente Trabalho - 612rpm
Amaral (2007) - 612rpm
H[m]
8
6
4
2
0
10
20
30
40
Q [ m 3/ h ]
50
60
Figura 6.3 – Altura de elevação do primeiro estágio da bomba centrífuga Imbil ITAP
65-330/2. Dados experimentais obtidos por Amaral (2007)
6.2 Obtenção de uma curva unificada de altura de elevação
para qualquer vazão e rotação
De uma maneira geral, em aplicações industriais onde haja a necessidade de
um sistema de bombeamento, nem sempre se tem uma bomba centrífuga disponível
no mercado com as características exatas de rotação e altura de elevação
desejadas. Por isso, os fabricantes de bombas centrífugas costumam disponibilizar
como opções, rotores de tamanhos diferentes para cada modelo de bomba, de
maneira a aumentar o campo de aplicação de cada equipamento sem perda
significativa de desempenho. A Figura 6.4 mostra a curva de altura de elevação total
em função da vazão volumétrica fornecida pelo fabricante da bomba centrífuga Imbil
ITAP 65-330/2.
74
H[m]
25
20
240
15
190
10
220
260
160
5
0
0
10
20
30
Q [ m3 / h ]
40
50
60
Figura 6.4 – Altura de elevação total versus vazão volumétrica da bomba centrífuga
IMBIL ITAP 65-330/2, na rotação de 1150 rpm. Dados fornecidos no catálogo do
fabricante (Fonte: site Imbil)
Na Figura 6.4, os números que acompanham as linhas tracejadas são as
opções oferecidas pelo fabricante do equipamento para o diâmetro externo do
segundo rotor com unidades em milímetros. Como o segundo rotor está dentro de
uma voluta, há flexibilidade de tamanhos que, conseqüentemente, irão influenciar na
altura de elevação total disponibilizada pelo equipamento, como está mostrado na
Figura 6.4.
Ainda assim, outra maneira de modificar o campo de aplicação de uma
bomba centrífuga está em alterar a rotação de trabalho do equipamento. Isso pode
ser usado, quando se deseja reutilizar uma bomba centrífuga em uma outra
aplicação, evitando assim a necessidade de compra de uma nova bomba.
Quando se utiliza uma mesma bomba centrífuga em rotações diferentes da
nominal devem-se utilizar relações de similaridade para se obter a curva de altura de
elevação corrigida para a nova rotação. Essas relações de similaridade são
apresentadas nas equações (6.2) a (6.4) (Fox e MacDonald, 2001):
Q2 ω2
=
Q1 ω1
H2  ω2 
= 
H1  ω1 
(6.2),
2
(6.3)
75
3
Wɺ 2  ω2 
= 
Wɺ 1  ω1 
(6.4),
onde os sub-escritos “1” representam os valores de vazão, altura de elevação e
potência, conhecidos em uma determinada rotação (por exemplo, na rotação
nominal fornecida pelo fabricante do equipamento) e “2”, a nova condição de
ɺ representa a potência necessária para o funcionamento da
trabalho. O símbolo W
bomba em uma determinada vazão e rotação.
Com o objetivo de verificar se os resultados numéricos apresentados na
Figura 6.3 são similares entre si, as equações (6.2) e (6.3) são aplicadas à curva de
altura de elevação para a rotação de 1000 rpm de modo a obter a transposição
desses resultados para a rotação de 1150rpm. O mesmo é realizado para as curvas
de 806 e 612rpm, gerando como resultado a Figura 6.5. Vale ressaltar que a
aplicação das equações (6.2) a (6.4) está baseada na hipótese de isoeficiência do
escoamento nas duas condições.
10
Hsimilar [ m ]
8
Presente Trabalho - 1150rpm
Presente Trabalho - 1000rpm
Presente Trabalho - 806rpm
Presente Trabalho - 612rpm
6
4
2
0
10
20
30
40
Q [ m3 / h ]
50
60
Figura 6.5 – Aplicação das equações de similaridade às curvas de 1000, 806 e 612
rpm transpondo-as para a rotação de 1150 rpm.
Da Figura 6.5 é possível notar que os resultados para as curvas de 1000, 806
e 612 rpm, são muito próximos àqueles apresentados para a curva de 1150 rpm.
Isso indica que os resultados são similares. Note que as curvas de 612, 806 e
1000rpm quando transpostas por similaridade à rotação de 1150rpm apresentam
76
uma faixa de vazão diferente das suas curvas originais. Isso acontece devido à
aplicação da Equação (6.2) que tem a função de garantir que o ponto transposto à
nova rotação seja semelhante ao original, ou seja, o comportamento do escoamento
no interior da bomba centrífuga à 10m3/h e 612rpm é semelhante ao escoamento
18,8m3/h e 1150rpm.
Para avaliar se essa premissa estava sendo observada no modelo numérico,
primeiramente foi considerado o par ordenado (Q ; H
)612rpm = ( 29,0[m3 / h ]; 1,78[m]) ,
extraído da curva de 612rpm que, quando transposto à rotação de 1150rpm gerou o
par ordenado (Q ; H
)1150rpm = ( 54,49[m3 / h ]; 6,28 [m] ) .
Posteriormente, utilizando a
vazão de 54,49m3/h, foi realizada uma simulação numérica adicional para a rotação
de 1150rpm de modo a conferir se o valor obtido para altura de elevação seria
próximo de 6,28m. O valor obtido foi de 6,34m, o que corresponde a um desvio de
1,0% em relação ao valor obtido anteriormente por similaridade.
Uma vez que os pontos são similares é possível obter uma curva de altura de
elevação que seja função da vazão e também da rotação de modo a contemplar, por
exemplo, rotações diferentes das originalmente simuladas. Com os resultados
numéricos obtidos através das simulações da bomba operando a 1150 rpm, foi
realizado um ajuste utilizando um polinômio de segundo grau, obtendo-se:
H = aQ 2 + bQ + c
(6.5)
onde ‘a’, ‘b’ e ‘c’ são valores constantes que foram determinadas e estão mostradas
na tabela 6-1, a altura de elevação é dada em metros e a vazão, Q, em m3/h. A
Equação (6.5) foi determinada para a rotação de 1150rpm, porém, para que seja
possível obter valores de altura de elevação para outras rotações são necessárias
algumas manipulações matemáticas que envolvem as Equações (6.2), (6.3) e (6.5),
de onde obtém-se a seguinte relação:
H (Q , ω ) = aQ 2 + d ωQ + e ω2
(6.6)
77
sendo ‘d’ e ‘e’ valores constantes que foram determinados e também mostrados na
tabela 6-1.
Tabela 6-1 – Valores das constantes das equações (6.5) e (6.6)
a
-0,001675
b
0,0626
c
7,966
d
5,441x 10
e
6,023 x 10
-5
-6
A Equação (6.6) mostra a altura de elevação como função da vazão
volumétrica e, também, da rotação. Observando as constantes ‘a’ e ‘e’ das equações
(6.5) e (6.6) pode-se avaliar o comportamento da altura de elevação com as
variáveis Q e ω . Como a < 0 , espera-se que a altura de elevação diminua com o
aumento da vazão. Em contrapartida, com o aumento da rotação, valores mais
pronunciados de altura de elevação são esperados, pois e > 0 . A faixa de validade
da Equação (6.6) é dada para valores de vazão e rotação de modo que:
10
Q
54
< <
1150 ω 1150
ou
(6.7)
0,00870. ω < Q <0,0470. ω
onde a vazão é dada em m3/h e a rotação em rpm. Como pode-se notar da Figura
6.4, a faixa de vazão contemplada pelo fabricante vai de 10m3/h a 54m3/h, para o
caso do segundo rotor com 260,0 mm (rotor utilizado por Amaral, 2007). A Equação
(6.7) tem o objetivo de transpor essa faixa de vazão para qualquer rotação, de modo
a manter as mesmas características do escoamento na nova condição de operação.
O ajuste utilizando do polinômio de segunda ordem da Equação (6.6)
apresentou um desvio médio relativo de 0,8% em relação aos valores de altura de
elevação obtidos pelas simulações numéricas e um coeficiente de determinação
R2=0,995. A Equação (6.6) ainda está sob a hipótese de que o rendimento hidráulico
é o mesmo para qualquer rotação utilizada. Para que essa hipótese seja sempre
78
atendida, não se devem utilizar valores de rotação muito superiores ou muito
inferiores à 1150rpm.
6.3 Potência e Eficiência Hidráulica
Outros parâmetros de engenharia importantes no contexto de máquinas de
fluxo são a potência necessária para o acionamento da bomba centrífuga e sua
eficiência. O cálculo dessas grandezas depende de diversos parâmetros do sistema,
como altura de elevação, vazão, massa específica do fluido, torque sofrido pelo eixo
do rotor e rotação. Para o cálculo da eficiência global da bomba, são contabilizadas
todas as perdas envolvidas no processo de conversão da energia fornecida pelo
motor da bomba em energia de pressão, porém, no presente trabalho não há como
mensurar perdas por atritos nos mancais da bomba ou perdas por fugas ou
vazamento do equipamento, mas somente a perda hidráulica da bomba, como está
mostrado na Equação (6.8):
Wɺ
ρg Q H
(6.8),
ηh [%] = ɺ saída =
⋅ 100%
Tω
WEntrada
onde T é o torque no motor da bomba centrífuga. Na Equação (6.8) é mostrada a
relação entre potência de saída e de entrada fornecida ao escoamento através da
rotação do rotor. No programa Ansys CFX é possível extrair o torque que o eixo do
rotor está sentindo devido à rotação e a passagem do escoamento por esse
componente. Com isso foram contabilizados todos os torques que as superfícies
sólidas que compõem o rotor estavam sentindo, em relação ao eixo de giro do rotor,
para a rotação de 1150 rpm e o resultado está mostrado na Figura 6.6.
79
12
Tnumérico [ N.m ]
10
8
6
4
2
0
10
20
30
Q [ m3 / h ]
40
50
Figura 6.6 – Torque numérico versus vazão no primeiro rotor da bomba centrífuga
Imbil ITAP 65-330/2, para uma rotação de 1150 rpm
Com os torques obtidos é possível calcular a potência de entrada numérica
através da relação Wɺ Entrada = T . ω . Para aplicação da relação anterior, uma hipótese,
assumida nesse trabalho, é de que toda a potência fornecida pelo motor é
disponibilizada ao rotor. A Figura 6.7 mostra um gráfico de potência de entrada
versus vazão volumétrica para a rotação de 1150rpm. Não é possível estabelecer
uma comparação direta entre resultados numéricos e experimentais para a potência
de entrada, uma vez que não se tem a informação de qual percentual do torque
global medido experimentalmente foi proveniente do primeiro estágio da bomba.
80
1400
Wentrada [ Watt ]
1200
1000
800
600
400
200
0
10
20
30
Q [ m3 / h ]
40
50
Figura 6.7 – Potência de entrada numérica versus vazão no primeiro rotor da bomba
Imbil ITAP 65-330/2, para um rotação de 1150 rpm
Com os resultados de altura de elevação versus vazão, obtidos das
simulações numéricas, é possível construir a curva da potência de saída do primeiro
estágio da bomba. A potência de saída ajuda a identificar qual é a parcela de
energia entregue ao rotor que efetivamente é transferida ao fluido sob forma de
ɺ
energia de pressão, através da relação W
Saída = ρ gQH . A Figura 6.8 mostra uma
comparação entre as potências de saída experimental e numérica. Como esperado
há uma forte semelhança entre as curvas, já que a potência está intrinsecamente
ligada aos valores de altura de elevação, cuja comparação foi muito boa, conforme
observado na Figura 6.3.
81
1400
Wsaída [ Watt ]
1200
1000
800
600
400
Presente Trabalho
Amaral (2007)
200
0
10
20
30
Q [ m3 / h ]
40
50
Figura 6.8 – Potência de saída numérica versus vazão no primeiro estágio da bomba
centrífuga Imbil ITAP 65-330/2, para rotação de 1150 rpm.
Uma vez determinada a curva de potência de entrada e de saída do primeiro
estágio da bomba, é possível determinar a eficiência hidráulica desse estágio, com
auxílio da Equação (6.8). A curva de eficiência hidráulica do primeiro estágio da
bomba está mostrada na Figura 6.9, onde também está mostrada a eficiência global
da bomba fornecida pelo fabricante. Novamente, não é possível uma comparação
direta entre os resultados de eficiência obtidos numericamente e os fornecidos pelo
fabricante, pois, estes últimos incluem efeitos globais da bomba, bem como
irreversibilidades mecânicas, como atritos dos mancais, perdas por vazamentos, etc.
Como foi mencionando, o segundo rotor possui diâmetro externo maior do
que o primeiro rotor, porém, ambos rotores possuem mesmo diâmetro de entrada.
Isso faz com que os rotores não sejam geometricamente idênticos. Além disso, há
componentes distintos à entrada e saída dos dois rotores, fazendo com que os
estágios sejam diferentes entre si. Observe na Figura 6.9 que a vazão onde a
eficiência hidráulica é máxima, obtida como resultado das simulações numéricas no
primeiro estágio da bomba, está próxima de 44m3/h, diferente da vazão de projeto
apresentada pelo fabricante, em 36m3/h. Uma das razões possíveis dessa
discrepância está no fato de se contabilizar apenas o primeiro estágio da bomba no
estudo numérico.
82
90
80
η[%]
70
60
50
Presente Trabalho
Fabricante
40
30
0
10
20
30
40
Q [ m3 / h ]
50
60
Figura 6.9 – Eficiência hidráulica numérica do primeiro estágio da bomba centrífuga
Imbil ITAP 65-330/2 e Eficiência global fornecida pelo fabricante.
Analisando os resultados experimentais apresentados por Amaral (2007),
para vazões altas, o indutor (elemento localizado à saída do difusor) apresenta
perdas de pressão grandes. Para vazões acima de 36m3/h, essas perdas possuem a
mesma ordem de grandeza da variação de pressão no difusor, isto é, para vazões
elevadas uma grande parcela do acréscimo de pressão adquirida pelo fluido no
difusor é perdida no indutor. Com essas observações, espera-se que ao simular
numericamente o escoamento na bomba centrífuga Imbil ITAP 65-330/2,
contabilizando todos os componentes (os dois rotores, difusor, indutor e voluta) os
resultados de eficiência hidráulica tenham comportamento mais próximos aos
apresentados pelo fabricante.
83
6.4 Propriedades do escoamento no interior da bomba
centrífuga estudada
Um dos objetivos do presente trabalho é capturar a periodicidade do
escoamento que, como mencionado no capítulo 5, foi assumida estar a cada 90°,
devido a características geométricas do rotor e do difusor.
A Figura 6.10 mostra de maneira esquemática o rotor e o difusor da bomba
centrífuga estudada, evidenciando o instante de tempo no qual a pá do rotor está
alinhada com a aleta do difusor. São definidas as faces de pressão (FP) e sucção
(FS), tanto da pá do rotor quanto da aleta do difusor. Essas faces têm essa
nomenclatura devido à distribuição da pressão ao longo da pá ou da aleta, onde,
para o caso de rotores com pás curvadas para trás a face convexa apresenta
pressões maiores em relação à face côncava.
Figura 6.10 – Instante de tempo no qual a pá do rotor e a aleta do difusor estão
alinhadas. Definição das Faces de Pressão (FP) e Faces de Sucção (FS)
A Figura 6.11 mostra a variação da pressão na interface entre rotor e difusor,
no centro do canal do rotor, no instante de tempo descrito pela Figura 6.10. Na
região inferior do gráfico, os quadrados indicam a posição das aletas do difusor, e na
região superior, os losangos sinalizam a posição das pás do rotor. Para analisar a
periodicidade do escoamento, os valores correspondentes às faixas 90° < θ < 180° ,
84
180° < θ < 270° e 270° < θ < 360° foram transportados para a faixa 0° < θ < 90° ,
como está mostrado na Figura 6.12.
PInterface Rotor/Difusor [ kPa]
90
80
70
60
50
40
30
20
0
30
60
90 120 150 180 210 240 270 300 330 360
θ [ °]
Figura 6.11 – Pressão versus direção tangencial no raio de saída do rotor para o
caso de pá do rotor alinhada com a aleta do difusor
PInterface Rotor/Difusor [ kPa]
90
ω
80
(FP)
(FS)
70
0° a 90°
90° a 180°
180° a 270°
270° a 360°
60
50
40
30
(FS)
(FP)
20
0
10
20
30
40 50
θ [ °]
60
70
80
90
Figura 6.12 – Periodicidade da pressão com a direção tangencial na saída do
primeiro rotor da bomba centrífuga Imbil ITAP 65-330/2. Vazão 35m3/h e rotação de
1150rpm.
85
Apesar das pequenas variações observadas na Figura 6.12, é nítida a
periodicidade angular do escoamento no primeiro estágio da bomba. Nota-se que
em regiões próximas à face de pressão da aleta do difusor, a pressão atinge valores
mais elevados, pois o fluido ao abandonar o rotor se choca com as aletas do difusor,
sendo desacelerado nessa região e conseqüentemente recuperando pressão. Como
o movimento das pás do rotor se processa da esquerda para a direita (na Figura
6.12) é de se esperar que os picos de pressão se localizem na face convexa da
aleta do difusor (face de pressão), como pode ser observado para os ângulos de 30
e 60°. Para os ângulos de 0 e 90° esse pico não é a centuado pois, a presença da pá
do rotor nessas posições contribui para que haja maior interação com a aleta do
difusor.
A Figura 6.12 comprova que o escoamento possui simetria espacial. Com
relação à variável tempo, a Figura 6.13 mostra o resultado da pressão em um ponto
localizado no bordo de entrada da aleta como função do tempo em uma simulação
numérica com vazão de 36m3/h e rotação de 1150rpm. Pode-se notar que após a
primeira volta do rotor a flutuação da pressão é periódica com relação ao tempo,
como era esperado. A partir da segunda metade da primeira volta os resultados
numéricos mostram que a bomba centrífuga atingiu o regime de operação, ou seja,
apesar da variação de pressão não ser a de regime, essa variação é periódica. Não
significa dizer que após a partida de uma bomba real, depois do deslocamento
angular de apenas meia volta do rotor, a bomba já estará operando em regime, pois,
nessa situação há parâmetros que não foram considerados no presente trabalho,
como a aceleração angular (no intervalo de tempo da partida do rotor da bomba), ou
ainda, uma mudança repentina de vazão. Os valores de pressão como função do
tempo da primeira metade da primeira volta do rotor não são reais e sim
representam o efeito da condição inicial imposta nas simulações numéricas.
86
PBordo entrada aleta difusor [ kPa]
90
80
70
60
50
40
30
20
10
0
0
0.2
0.4
0.6 0.8
1
1.2 1.4
t / período de uma volta [ ]
1.6
1.8
2
Figura 6.13 – Pressão versus tempo em um ponto localizado no bordo de entrada da
aleta do difusor da bomba centrífuga Imbil ITAP 65-330/2, para a rotação de 1150
rpm e vazão de 35m3/h.
É possível notar também, que o período de oscilação da pressão com o
tempo é determinado pela passagem das pás pela aresta da aleta na entrada do
difusor, uma vez que no intervalo da primeira para a segunda volta são identificados
oito ciclos de oscilação da pressão, sinalizando a passagem de todas as oito pás do
rotor por esse ponto. Vale ressaltar que o valor da pressão absoluta apresentada na
Figura 6.13 tem pouca relevância quando analisada isoladamente no presente
trabalho, pois como hipótese das simulações numéricas foi considerada escoamento
incompressível e sem mudança de fase. Isso torna o valor pontual da pressão pouco
expressivo, pois o interesse maior se concentra em diferenças de pressão entre
entrada e saída da bomba. No entanto, em estudos numéricos onde o efeito de
cavitação em máquinas de fluxo é considerado, por exemplo, a condição de
contorno de pressão, bem como os valores absolutos de pressão são relevantes.
As Figuras 6.15 e 6.16 mostram a distribuição de pressão ao longo da pá,
desde a entrada até a saída do rotor. Na Figura 6.15 foram tomados os valores de
pressão na interseção da pá com o cubo do rotor, e na Figura 6.16 esses valores
foram tomados na intersecção da pá com a coroa (parte superior do rotor). O eixo
das abscissas nas Figuras 6.15 e 6.16 são definidos, de acordo com o esquema
mostrado na Figura 6.14, de modo que em zero tem-se a entrada da pá do rotor e
87
em 1 tem-se a saída da pá do rotor, tanto para a região do cubo (Figura 6.15) quanto
para a região da coroa (Figura 6.16).
Bordo de entrada
0
0
Coroa
1
1
Bordo de saída
Cubo
Figura 6.14 – Localização do bordo de entrada e saída do rotor da bomba.
Nas Figuras 6.15 e 6.16 é possível observar a interação existente entre rotor
e difusor, principalmente na região próxima à saída do rotor e quando a pá do rotor
está praticamente alinhada com a aleta do difusor. Nota-se quando a pá do rotor se
afasta da aleta do difusor a pressão praticamente não é alterada.
70
Alinhado
3,75 °
7,50 °
11,25 °
15,00 °
60
P [ kPa ]
50
40
Superfície de Pressão
30
20
10
Superfície de Sucção
0
-10
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
( r - RPá entrada cubo ) / (RSai rotor - RPá entrada cubo ) [ ]
0.8
0.9
1
Figura 6.15 – Distribuição de pressão na pá do rotor desde a entrada até a saída na
altura da raiz da pá (ou cubo)
Os resultados de distribuição de pressão ao longo da pá do rotor, no cubo
(Figura 6.15) e na coroa (Figura 6.16) mostram poucas diferenças entre si. A
diferença mais significativa está região de entrada da pá do rotor. Na entrada, a pá é
88
fixada no cubo, em um diâmetro de aproximadamente 36,0 mm, e também na coroa
do rotor, em um diâmetro de 80,1 mm. Por isso, a diferença de pressão entre a face
de pressão e a face de sucção da pá do rotor é menos significativa na região de
entrada do cubo, por haver um campo centrífugo menor quando comparado a região
da coroa.
Na saída do rotor, nota-se que a pressão sofre alterações a cada instante de
tempo. Ao longo da pá do rotor, essas mudanças na magnitude da pressão são
pouco influenciadas pela posição relativa entre rotor e difusor.
70
Alinhado
3,75°
7,50°
11,25°
15,00°
60
P [ kPa ]
50
40
30
Superfície de Pressão
20
Superfície de Sucção
10
0
-10
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
( r - RPá entrada coroa ) / ( RSai rotor - RPá entrada coroa )
0.9
1
Figura 6.16 - Distribuição de pressão na pá do rotor desde a entrada até a saída na
altura da coroa, para o rotor girando a 1150 rpm e na vazão de 35m3/h
As Figuras 6.17 a 6.19 mostram o campo de pressão no interior do primeiro
estágio da bomba, para três vazões distintas na rotação de 1150 rpm: 20, 35 e
50m3/h. Comparando o campo de pressão gerado para essas três vazões pode-se
visualizar qualitativamente que na vazão de 20m3/h o ganho de pressão no primeiro
estágio da bomba é maior quando comparado com as demais vazões (35 e 50m3/h).
Pode-se observar, também a interação existente entre rotor e difusor, por exemplo,
pela deformação das linhas de iso-pressão na saída do rotor.
É possível notar, também, para todas as vazões que, na superfície côncava
do difusor, no bordo de entrada da aleta há uma região onde a pressão sofre uma
89
leve queda e que, como foi visto, altera o perfil de pressão na ponta da pá do rotor
ao passar por essa região. O motivo dessa queda localizada de pressão pode estar
relacionado a um ajuste no ângulo da velocidade do escoamento ao sair do rotor,
sendo o fluido obrigado a contornar o bordo de entrada da aleta do difusor. As
partículas fluidas que são defletidas para a parte côncava da aleta ficam sujeitas a
um ressalto hidráulico (que pode ser de grande ou de pequena intensidade,
dependendo da magnitude da vazão). Para vazões menores que a de projeto, como
20m3/h, devido à inclinação insuficiente do escoamento em relação as aletas do
difusor, esse ressalto hidráulico pode ser uma fonte de instabilidades e recirculações
no interior do difusor.
Figura 6.17 – Campo de pressão no interior do primeiro estágio da bomba centrífuga
Imbil ITAP 65-330/2, vazão de 20m3/h e rotação de 1150 rpm.
90
Figura 6.18 – Campo de pressão no interior do primeiro estágio da bomba centrífuga
Imbil ITAP 65-330/2, vazão de 35m3/h e rotação de 1150 rpm.
Figura 6.19 – Campo de pressão no interior do primeiro estágio da bomba centrífuga
Imbil ITAP 65-330/2, vazão de 50m3/h e rotação de 1150 rpm.
91
Observando as Figuras 6.17 a 6.19 que mostram que a medida que a vazão
aumenta e o ganho de pressão diminui, observa-se a possibilidade de cavitação em
vazões mais elevadas. Obviamente que os estudos do presente trabalho não foram
feitos levando em consideração tal efeito, pois, se considerou escoamento
incompressível, isotérmico e sem mudança de fase. Em todas as simulações
numéricas realizadas, foi imposta condição de contorno de pressão constante de
referência ( Pref = 0 [Pa] ) na entrada do tubo de admissão da bomba. Entretanto,
como é conhecido na prática, há uma redução da pressão na entrada da bomba que
é função de parâmetros geométricos do sistema como a diferença de nível entre o
reservatório e a entrada da bomba, pressão do reservatório e principalmente em
função características dinâmicas como perdas por atrito estabelecidas ao longo da
tubulação de entrada devido à presença de válvulas, acessórios e conexões. Essas
perdas por atrito são tanto maiores quanto maior é a vazão estabelecida.
Todas essas perdas podem, eventualmente fazer com que haja condições de
operação indesejáveis como a cavitação, principalmente para vazões muito acima
da condição de projeto onde a pressão na entrada tende a ser reduzida a valores
críticos. Nos experimentos apresentados por Amaral (2007), por exemplo, uma
bomba centrífuga adicional foi instalada à entrada da bomba centrífuga Imbil ITAP
65-330/2 (utilizada nos testes experimentais) com o objetivo único de garantir que a
pressão absoluta na entrada da bomba testada não atingisse valores críticos de
cavitação.
Na busca de maior compreensão dos fenômenos presentes no escoamento
monofásico existente no interior da bomba, os campos de velocidades no difusor e
no rotor foram analisados, nas três vazões mencionadas anteriormente, 20, 35 e
50m3/h e também para a vazão onde se obteve a maior eficiência hidráulica
numérica, 44m3/h. Devido à complexidade da geometria foram estabelecidos
diversos planos de corte ao longo dos canais do rotor e do difusor para melhor
visualização do comportamento do escoamento. A Figura 6.20 mostra a forma como
esses planos foram distribuídos ao longo dos canais do rotor e do difusor. Os planos
N1, N2, N3 e N4 estão dispostos de maneira a capturar o comportamento do campo
de velocidades em planos aproximadamente normais as aletas do difusor. Os planos
92
T1, T2 e T3 estão dispostos, respectivamente, próximo à superfície inferior de rotor e
difusor, no centro do canal e próximo à superfície superior.
Figura 6.20 – Planos de corte da visualização do escoamento, obtido
numericamente, no rotor e difusor do primeiro estágio da bomba centrífuga Imbil
ITAP 65-330/2
As Figuras 6.22 a 6.26 mostram o comportamento do escoamento no interior
do rotor e do difusor sob a forma de linhas de corrente dispostas nos planos de corte
mostrados na Figura 6.20. As Figuras 6.22 e 6.26(d) dizem respeito à vazão de
50m3/h. Nessa vazão, devido à alta inércia do escoamento, nota-se que
praticamente não há recirculações ou perturbações no escoamento tanto no rotor
como no difusor. O mesmo pode ser observado para a vazão de melhor eficiência
hidráulica numérica, 44m3/h, através das Figuras 6.23 e 6.26(c).
Para vazões mais baixas, o aparecimento de instabilidades e perturbações no
escoamento é observado, principalmente no interior do difusor. Há várias regiões e
fenômenos envolvidos na nucleação dessas instabilidades.
A face convexa da aleta pode ser um desses nucleadores, pois ao sair do
rotor, o escoamento não recebe energia de nenhuma fonte, escoando apenas por
inércia da energia que recebeu do rotor. Em vazões mais baixas essa inércia
recebida pelo fluido é menor e a face convexa da aleta tende a ser um obstáculo ao
escoamento, podendo haver descolamento do escoamento da face da aleta, com
conseqüente aparecimento de recirculação nessa região, como é mostrado nas
Figuras 6.24 e 6.25 no plano T2.
Outro fator está relacionado com a configuração da interface entre rotor e
difusor (ver Figura 6.21). Na bomba centrífuga estudada, há uma mudança repentina
93
de área entre a saída do primeiro rotor e entrada do difusor. Em todas as vazões, é
observada, no plano T3, uma perturbação do escoamento, em maior ou menor
intensidade. Essa perturbação do escoamento nessa região pode estar associada a
duas razões. A primeira delas, pela mudança de áreas na interface de rotor e difusor
e a segunda pela configuração da geometria do difusor, que tem sua saída de fluido
na região inferior, fazendo com que a região mostrada por T3 seja mais afetada pela
presença de instabilidades no escoamento. A presença da saída do difusor na
região inferior do mesmo, no caso do plano T1, age como um inibidor de
instabilidades, uma vez que o escoamento é forçado a defletir para essa região, e
assim minimizando os efeitos de mudança de áreas nessa região do escoamento.
Outro efeito importante é de que quando uma bomba centrífuga opera fora de
sua condição ideal de trabalho, o escoamento não entra no difusor tangenciando as
aletas, isto é, o vetor velocidade absoluta do escoamento não é tangente à entrada
das aletas do difusor. Portanto, se a vazão volumétrica assumida for maior que a
vazão de projeto, o escoamento irá tender a se chocar com a face côncava da aleta
do difusor. Para vazões de trabalho menores que a vazão de projeto, a angulação
do escoamento será insuficiente e o escoamento tenderá a colidir com a face
convexa da aleta do difusor.
As conseqüências dessas duas situações são distintas. Para vazões maiores,
o escoamento tende a ser defletido, em maior intensidade, pela superfície côncava
da aleta e, também por esse motivo, não são observadas recirculações ou
instabilidades, pois, superfícies côncavas tendem a inibir perturbações no
escoamento. Entretanto, para vazões mais baixas, a angulação insuficiente do
escoamento em relação à entrada da aleta do difusor força o escoamento a colidir
com a face convexa da aleta. Ao ser defletido por essa face (que possui uma
tendência de nucleação de instabilidades) e aliado ao fato de o escoamento possuir
baixa inércia, recirculações tendem a ser mais pronunciadas.
94
Figura 6.21 - Esquema da saída do rotor e entrada do difusor que mostra a diferença
de áreas entre esses dois componentes.
95
Figura 6.22 – Linhas de corrente do escoamento no interior da bomba centrífuga Imbil ITAP 65 330/2, na rotação nominal e vazão de 50m3/h. (a) escoamento no plano T1, (b) escoamento no
plano T2 e (c) escoamento no plano T3.
Figura 6.23 - Linhas de corrente do escoamento no interior da bomba centrífuga Imbil ITAP 65 330/2, na rotação nominal e vazão de 44m3/h. (a) escoamento no plano T1, (b) escoamento no
plano T2 e (c) escoamento no plano T3.
96
Figura 6.24 - Linhas de corrente do escoamento no interior da bomba centrífuga Imbil ITAP 65 330/2, na rotação nominal e vazão de 35m3/h. (a) escoamento no plano T1, (b) escoamento no
plano T2 e (c) escoamento no plano T3.
Figura 6.25 – Linhas de corrente do escoamento no interior da bomba centrífuga Imbil ITAP 65 330/2, na rotação nominal e vazão de 20m3/h. (a) escoamento no plano T1, (b) escoamento no
plano T2 e (c) escoamento no plano T3.
97
(b)
N4
N3
N2
N1
( c)
(d)
Figura 6.26 - Linhas de corrente do escoamento no interior do difusor bomba centrífuga Imbil ITAP 65 330/2, na rotação nominal e vazão de 20m3/h. (a) Vazão 20 m3/h, (b) Vazão 35 m3/h (c)
Vazão 44 m3/h e (d) Vazão 50 m3/h
98
Da Figura 6.26 pode-se notar a evolução do escoamento ao passar pelo
difusor. No plano N1, para as vazões de 44 e 50 m3/h, percebe-se que o
escoamento possui componentes de velocidade
essencialmente
horizontal,
possivelmente devido à geometria do difusor nessa região e também a grande
inércia do escoamento ao sair do rotor. À medida que o escoamento atinge os
planos N2 e N3, já há uma tendência do escoamento em defletir em direção à região
inferior do difusor. Por essa razão, ao observar as vazões de 35 e 20 m3/h nos
planos N2 e N3, nota-se que a magnitude da velocidade na região superior dos
planos é menor, quando comparada à região inferior, tornando assim a região
evidenciada pelo plano T3 mais suscetível às instabilidades do escoamento. Essa
sensibilidade existente na região mostrada pelo plano T3 poderia ser corrigida,
talvez se a parede superior do difusor apresentasse uma curvatura suave em
direção à saída do difusor ao invés da sua forma original plana. O plano N4 mostra a
região onde há comunicação direta entre difusor e indutor. É possível notar nessa
região uma forte tendência do escoamento em defletir para baixo. Para a vazão de
20 m3/h há ainda um turbilhonamento do escoamento, possivelmente conseqüência
das instabilidades geradas ao longo de todo o canal do difusor.
Sabe-se que com o modelo de turbulência utilizado nos estudos numéricos do
presente trabalho ( κ − ε padrão) não é possível localizar, nem tampouco mensurar,
de maneira precisa as recirculações formadas no escoamento, pois tal modelo não
foi desenvolvido para esse fim. No entanto, é possível concluir, ainda que de
maneira qualitativa, que ocorrerão perturbações do escoamento para baixas vazões
e ocorreram descolamentos e recirculações ao longo do difusor para vazões
menores que a vazão de projeto. Vale ressaltar que mesmo na vazão de projeto
fornecida pelo fabricante (36m3/h) também foram observadas recirculações (vide
Figura 6.24 e 6.26(b)) nas superfícies de pressão das aletas, próximas à saída do
difusor. Isso dá indícios de que mesmo o valor de eficiência máxima global da
bomba fornecido pelo fabricante ser dado para vazão em torno de 36m3/h, isso não
implica que o ponto de eficiência hidráulica máximo do primeiro estágio da bomba
centrífuga estudada estará também dado nessa vazão.
99
7 CONCLUSÕES E SUGESTÕES DE TRABALHOS FUTUROS
Neste trabalho, foi simulada numericamente a fluidodinâmica do escoamento
no interior do primeiro estágio de uma bomba centrífuga de duplo estágio,
considerando um escoamento monofásico turbulento em regime transiente e água
com
propriedades
constantes como fluido
de trabalho.
Foram
simuladas
numericamente quatro velocidades angulares diferentes do rotor: 1150, 1000, 806 e
612rpm.
A partir dos resultados numéricos, foram elaboradas as curvas da diferença
de pressão no rotor e no difusor em função da vazão, apresentando boa
concordância com os resultados experimentais obtidos por Amaral (2007). As
discrepâncias entre os resultados de diferença de pressão numéricos e
experimentais para o primeiro rotor foram pequenas, da ordem de 7%, porém, para o
caso do difusor os desvios chegaram a cerca de 20%. Um dos motivos para as
discrepâncias maiores pode estar associado ao modelo de turbulência utilizado. É
sabido que o modelo a duas equações κ − ε padrão não é aconselhável para
escoamentos sobre superfícies curvas e na presença de recirculações. Como no
difusor apresentou maiores instabilidades e recirculações quando comparado ao
rotor, é possível que isso tenha contribuído para o aumento dos desvios de diferença
de pressão nesse componente. Por outro lado, a ordem de grandeza dos valores de
diferença de pressão no rotor é maior do que no difusor fazendo que globalmente os
resultados apresentados tenham boa concordância com os dados experimentais em
termos de altura de elevação versus vazão.
Os resultados das curvas de altura de elevação para 1000, 806 e 612rpm
foram transportados utilizando as equações de similaridade à rotação de 1150rpm e
mostraram-se muito próximos aos obtidos das simulações numéricas à 1150rpm.
Com isso, um ajuste de curvas foi feito e uma equação para altura de elevação como
função da vazão e da rotação foi determinada para o primeiro estágio da bomba.
Com os resultados numéricos obtidos foram levantadas curvas de potência de
entrada e de saída do primeiro estágio da bomba. Para a potência de entrada não foi
100
possível estabelecer uma comparação direta com os resultados experimentais. Isso
porque nos testes experimentais, o torque de entrada foi medido de maneira global
(que reflete na potência de entrada da bomba inteira), não sendo possível mensurar
a parcela de torque consumida em cada estágio. Para a potência de saída, os
resultados numéricos estão em concordância com os experimentais medidos por
Amaral (2007).
Com os resultados de potência de entrada e de saída do primeiro estágio da
bomba, a curva de eficiência hidráulica foi determinada. Novamente, uma
comparação direta com os resultados experimentais não pôde ser estabelecida
devido ao estudo do presente trabalho estar focado apenas no primeiro estágio da
bomba e os dados fornecidos pelo fabricante serem globais. Apesar disso, os
resultados numéricos mostram que a eficiência hidráulica máxima do primeiro
estágio ocorre em uma vazão maior do que a eficiência global máxima fornecida
pelo fabricante.
A partir dos resultados numéricos foi analisado o comportamento do campo
de pressão e velocidade no rotor e difusor da bomba. Foi mostrado que o
escoamento no interior da bomba centrífuga estudada possui periodicidade tanto no
espaço quanto no tempo. A periodicidade espacial está ligada ao número de pás do
rotor (oito) em relação ao número de aletas do difusor (doze). Por essa configuração,
observou-se uma repetitividade do escoamento a cada 90º.
Considerando os resultados alcançados no presente trabalho, sugere-se para
trabalhos futuros:
•
Simular numericamente a bomba centrífuga completa, os dois estágios,
para determinar a eficiência hidráulica de todo o conjunto;
•
Simular o escoamento para outros fluidos, de viscosidades mais
elevadas e observar o comportamento das curvas de desempenho da
bomba;
•
Simular o escoamento na bomba utilizando modelos de turbulência
avançados, como por exemplo os modelos SST e LES (simulação de
grandes escalas).
101
•
Desenvolvimento da metodologia numérica para a simulação do
escoamento bifásico líquido-gás.
102
8 REFERÊNCIAS BIBLIOGRÁFICAS
AMARAL, Gilmar D. L.: “Modelagem do Escoamento Monofásico em Bomba
Centrífuga Submersa Operando com Fluidos Viscosos”. Campinas: Faculdade
de Engenharia Mecânica, Universidade Estadual de Campinas, 2007. 233p.
Dissertação (Mestrado).
Ansys CFX-Solver Theory Guide, “Manual do programa computacional Ansys CFX”
ASUAJE, M., BAKIR, F., KOUIDRI, S., KENYERY, F., REY, R.. “Numerical
Modelization of the Flow in Centrifugal Pump: Volute Influence in Velocity and
Pressure Fields”. International Journal of Rotating Machinery, p. 244-255, 2005.
BACHAROUDIS, E. C., FILIOS, A. E., MENTZOS, M. D., MARGARIS, D. P.,
“Parametric Study of a Centrifugal Pump Impeller by Varying the Outlet Blade
Angle”, The Open Mechanical Engineering Journal, Volume 2, p.75-83, 2008.
BENRA, I. F.-K., “Economic Development of Efficient Centrifugal Pump
Impellers by Numerical Methods”, World Pumps, p.48-53, 2001.
CHEAH, K. W., LEE, T. S., WINOTO, S. H. and ZHAO, Z. M.. “Numerical Flow
Simulation in a Centrifugal Pump at Design and Off-Design Conditions”,
International Journal of Rotating Machinery, 8 p., Singapure, 2007.
ESTEVAM, V. “Uma Análise Fenomenológica da Operação de Bomba
Centrífuga com Escoamento Bifásico”. Campinas: Faculdade de Engenharia
Mecânica, Universidade Estadual de Campinas, 2002. 265 p. Tese (Doutorado).
FENG, J., BENRA, F.K., DOHMEN H. J.. “Numerical Investigation on Pressure
Fluctuations
for
Different
Configurations
of
Vaned
Diffuser
Pumps”.
International Journal of Rotating Machinery, Volume 2007, 10 p. Hindawi Publishing
Corporation, 2007.
FOX, R. W., McDONALD, A. T.. “Introdução à Mecânica dos Fluidos”, 5ª edição,
Editora LTC, 2001.
103
GARCÍA, N. R., VICENT, S. C.. “Simulación, Análisis y Rediseño de Bombas
Centrífugas”, Universidad Jaume I, Castellon de la Plana, (Titulação: Engenharia
Industrial), 2007, Espanha.
GÜILICH, J. F., “Pumping Highly Viscous Fluids with Centrifugal Pumps - Part
1”, World Pumps, pp 30 – 34, 1999a.
GÜILICH, J. F., “Pumping Highly Viscous Fluids with Centrifugal Pumps - Part
2”, World Pumps, pp 39 – 42, 1999b.
Hydraulic Institute “Hydraulic Institute Standards for Centrifugal, Rotary &
Reciprocating Pumps” 14th Edition, 1983.
IMBIL, site www.imbil.com.br, acesso em 16/05/2009
KBS bombas, site www.ksb.com.br, acesso em 10/02/2010
LAUNDER, B.E. e SPALDING, D.B., “The Numerical Computation of Turbulent
Flows”, Comp Meth Appl Mech Eng, 3:269-289, 1974.
LI, W. G.. “Effects of Viscosity of Fluids on Centrifugal Pump Performance and
Flow Pattern in the Impeller”, International Journal of Heat and Fluid Flow, V.21, pp
207-212, 2000.
MALISKA, C. R. “Transferência de Calor e Mecânica dos Fluidos” Editora LTC,
472p, 2009.
MACINTYRE, A. J.. “Bombas e Instalações de Bombeamento”, Editora LTC,
782pg., 1997.
MAÑES, J. P., VICENT, S. C.. “Imbil Ita 65 330/2: Ensayo Experimental y
Simulación Mediante Técnicas de Mecánica de Fluidos Computacional”,
Universidad Jaume I, Castellon de la Plana, (Titulação: Engenharia Industrial), 2009,
Espanha.
PATANKAR, S.V. “Numerical Heat Transfer and Fluid Flow”. Hemisphere
Publishing Corp, 1980.
PFLEIDERER, C., “Bombas Centrífugas y turbocompressores”. Barcelona:
Editorial Labor, 1960, 631 p.
104
PETROBRAS, site www.petrobras.com.br, acesso em 15/05/2009
RHIE, C. M. e CHOW, W. L., “A Numerical Study of Turbulent Flow Past an
Isolated Airfoil with Trailing Edge Separation”, AIAA, 1982.
SEGALA, W., “Centrifugal Pump Performance: Numerical Simulation and
Experimental Data Comparisons”. 12th Brazilian Congress of Thermal Engineering
and Sciences, Belo Horizonte – Brazil, 2008.
SHOAJAEE F., M. H., BOYAGHCHI, F. A., “Studies on the Influence of Various
Blade Outlet Angles in a Centrifugal Pump when Handling Viscous Fluids”.
American Journal of Applied Sciences, p. 718-724, 2007.
STEPANOFF, A. J. “Centrifugal and Axial Flow Pumps – Theory, Design and
Application”. Second edition. John Wiley & Sons, New York, 1957.
VERSTEEG, H.K., MALALASEKERA W. “Introduction to computational fluid
dynamics. The finite volume method”, Longman Scientific & Technical, 1995,
Malásia.
WILCOX, D.C., 1993, “Turbulence modeling for CFD”, 1ª edição, DCW Industries,
La Cañada, Califórnia.
WUIBAUT, G., BOIS, G., EL HAJEM, M., AKHRAS, A., CHAMPAGNE, J. Y.,
“Optical PIV and LDV Comparisons of Internal Flow Investigations in SHF
Impeller”, International Journal of Rotating Machinery, Artigo ID 69521, pp 1-9, 2006.