Embed Size (px)
Citation preview
Estudo Numérico do Problema
do Golpe de Aríete em Tubulações
Raphael Costa Carvalho
Projeto de Graduação apresentado ao Curso de
Engenharia Mecânica da Escola Politécnica,
Universidade Federal do Rio de Janeiro, como
parte dos requisitos necessários à obtenção do
título de Engenheiro.
Orientador : Helcio Rangel Barreto Orlande, Ph.D.
Rio de Janeiro
Março de 2018
i
Carvalho, Raphael Costa
Estudo Numérico do Problema do Golpe de Aríete em
Tubulações/ Raphael Costa Carvalho. – Rio de Janeiro:
UFRJ/ Escola Politécnica, 2018.
VIII, 47 p.: il.; 29,7 cm.
Orientador: Helcio Rangel Barreto Orlande
Projeto de Graduação – UFRJ/ Escola Politécnica/
Curso de Engenharia Mecânica, 2018.
Referencias Bibliográficas: p. 45-47.
1. Introdução. 2. Revisão Bibliográfica. 3. Modelo
Matemático. 4. Método WAF-TVD 5. Resultados
6.Conclusões e Sugestões para Trabalhos Futuros 7.
Referências I. Orlande, Helcio Rangel Barreto. II.
Universidade Federal do Rio de Janeiro, Escola
Politécnica, Curso de Engenharia Mecânica. III. Estudo
Numérico do Problema do Golpe de Aríete em
Tubulações.
ii
Agradecimentos Agradeço a Deus por ter me dado força e resiliência ao longo da graduação e na realização
deste trabalho.
Agradeço à Nossa Senhora pela intercessão.
Agradeço à minha mãe Tânia e ao meu pai Márcio pela educação, atenção e amor dados a
mim e à minha irmã e também ao apoio e oração durante os cinco anos de graduação.
Agradeço à minha irmã Raphaelle por ser minha melhor amiga durante toda a vida.
Agradeço à minha avó Lourdes e as minhas tias Fátima e Aparecida e ao meu primo Lucas
por todo o carinho.
Agradeço ao professor Helcio Rangel Barreto Orlande pela paciência, atenção e
disponibilidade na orientação deste trabalho e todas as instruções valiosas para realizá-lo.
Agradeço aos meus amigos de faculdade Bruno, Caio, Eduardo, Lucas, Pedro e Rafael
pela convivência e horas de estudo compartilhadas ao longo da graduação.
Agradeço a todos os alunos e funcionários do Laboratório de Tecnologia e Transmissão do
Calor (LTCC) pela amizade e ajuda durante minha iniciação científica no laboratório e em
especial a Diego Estumano que me acompanhou nos dois primeiros anos no laboratório.
Agradeço também a todos os professores do Departamento de Engenharia Mecânica da
UFRJ pelos ensinamentos e experiências passadas durante a graduação.
iii
Resumo do Projeto de Graduação apresentado à Escola Politécnica/UFRJ como parte dos
requisitos necessários para a obtenção do grau de Engenheiro Mecânico
Estudo Numérico do Problema do Golpe de Aríete em Tubulações
Raphael Costa Carvalho
Março/2018
Orientador: Helcio Rangel Barreto Orlande
Curso: Engenharia Mecânica
Escoamentos transientes possuem uma grande importância em problemas de engenharia,
tanto para dimensionar corretamente a tubulação, quanto para a escolha adequada dos
materiais que a compõe. Dessa forma, no presente trabalho, é feito o estudo numérico de
um fenômeno transiente chamado golpe de aríete. Para tanto, é utilizada uma abordagem
numérica baseada no método de volumes finitos usando o esquema WAF-TVD. Esse
esquema é de segunda ordem no espaço e no tempo. Para a verificação e validação do
código, os resultados numéricos obtidos aqui são comparados com resultados presentes na
literatura.
Palavras-chave: Golpe de Aríete, WAF-TVD, Método Numérico.
iv
Abstract of Undergraduate Project presented to POLI/UFRJ as a partial fulfillment of the
requirements for the degree of Mechanical Engineer
Numerical Study of Water Hammer Problem in Pipelines
Raphael Costa Carvalho
Março/2018
Adivisor: Helcio Rangel Barreto Orlande
Course: Mechanical Engineering
Transient flows have a great importance in engineering problems, as for properly
dimensioning the pipelines and for the proper choice of the materials that make it. Thus, in
the present work, the numerical study of a transient phenomenon called water hammer is
done. For this, the numerical approach based on the finite volume method is used through
the WAF-TVD scheme. This is a second-order scheme in space and time. For the
verification and validation of the code, the numerical results obtained here are compared
with results present in the literature.
Keywords: Water Hammer, WAF-TVD, Numerical Method
v
Sumário 1 - Introdução.................................................................................................................... 1
2 – Revisão Bibliográfica ................................................................................................. 3
2.1 Contexto Histórico .................................................................................................. 3
2.2 A Importância da Prevenção do Golpe de Aríete ................................................... 4
2.3 Equipamentos de Proteção...................................................................................... 6
2.4 Métodos Numéricos ................................................................................................ 8
3 - Modelo Matemático .................................................................................................. 10
3.1 Equação da Conservação da Quantidade de Movimento Linear .......................... 10
3.2 Equação da Conservação da Massa ...................................................................... 11
4 - Método WAF-TVD ................................................................................................... 17
4.1 Estudo da Natureza do Problema.......................................................................... 17
4.2 Fracionamento das equações Não Homogêneas ................................................... 18
4.3 Esquema de Fluxo Médio Ponderado (WAF) ...................................................... 19
4.4 Esquemas TVD ..................................................................................................... 21
4.5 Método de Solução ............................................................................................... 21
4.6 Solução Aproximada do Problema de Riemann ................................................... 23
4.7 Volumes Fictícios ................................................................................................. 25
4.8 Critério de Estabilidade ........................................................................................ 26
5 - Resultados ................................................................................................................. 27
5.1 Testes e verificações ............................................................................................. 27
5.2 Caso 1 ................................................................................................................... 27
5.3 Solução para o Caso 1 .......................................................................................... 28
5.4 Caso 2 ................................................................................................................... 32
5.5 Solução para o Caso 2 .......................................................................................... 35
6 - Conclusão e Sugestões para Trabalhos Futuros ........................................................ 43
7 - Referências ................................................................................................................ 45
vi
Nomenclatura
𝑎 celeridade
𝐴 Área da seção transversal do duto
b variável correspondente cada família de onda na solução do problema de
Riemann
𝑐𝑘 número de Courant
D diâmetro interno do duto
e espessura da parede do duto
E módulo de elasticidade do duto
𝑓 fator de atrito
F força
F vetor de fluxo
𝑔 aceleração da gravidade
J matriz Jacobiana
K módulo de elasticidade do líquido
L comprimento do duto
N número de ondas na solução do problema de Riemann
𝑝 pressão
𝑞 vazão volumétrica
r parâmetro de fluxo
𝑅𝑒 número de Reynolds
S operador representando sistema de equações diferenciais ordinárias
S vetor termo fonte
vii
t tempo
TV variação total
u velocidade
U vetor das variáveis conservativas
x posição axial
Z operador representando sistema de equações diferenciais parciais
Letras Gregas
𝜏0 tensão de cisalhamento na parede do duto
Δ variação
𝛾 peso específico
𝛿 diferencial
𝜃 ângulo entre a direção do escoamento e o plano horizontal de referência
𝜆 autovalor
𝜉 função utilizada nas Equações (4.16) e (4.17)
𝜌 massa específica
𝜐 volume específico
𝜙 função limitadora
Subscritos
* região com solução aproximada no problema de Riemann
0 estado inicial
at atrito
D à direita
viii
e entrada
E à esquerda
ext externa
i passo da solução no espaço
k família de onda do problema de Riemann
max máximo
min mínimo
N número de ondas na solução do problema de Riemann
𝑠 saída
Sobescritos
k família de onda do problema de Riemann
n passo da solução no tempo
1
1 - Introdução O escoamento de um líquido dentro de um duto é um fenômeno bastante conhecido e
analisado em diversos problemas de engenharia como pode ser visto, por exemplo, em
FOX et al (2014). Se em algum momento, a velocidade do fluido é alterada bruscamente
em algum ponto na tubulação, surge um fenômeno conhecido como transiente
hidráulico ou golpe de aríete (nome dado devido ao ruído de batida característico do
processo) (LOMBARDI, 2005, MENDES, 2011). Essa alteração de velocidade
normalmente é causada pelo fechamento rápido de válvulas ou devido à partida ou ao
desligamento rápido de bombas.
Segundo LÜDECKE et al. (2006), por causa da inércia do fluido, a velocidade do
líquido não é capaz de se ajustar à nova configuração do escoamento quando ocorre um
fechamento brusco de válvula. Como consequência, parte da energia cinética do fluido é
transformada em energia de pressão, que se converte em trabalho de compressão do
líquido e na deformação das paredes do duto (LOMBARDI, 2005). Devido a esses
fenômenos, surgem ondas de pressão ao longo da tubulação. As altas magnitudes e
velocidades de propagação (em torno de 1000 m/s) dessas ondas de pressão podem
colocar em risco as tubulações e seus acessórios, podendo ocorrer até mesmo
deformações plásticas nos tubos e falhas estruturais.
Assim, se a tubulação não for de material adequado ou não for bem dimensionada
para suportar as mudanças de pressão ocorridas durante o golpe de aríete, podem surgir
sérios problemas para o circuito hidráulico. Portanto, é necessário um melhor
entendimento de como ocorrem e quais as magnitudes e velocidades de propagação das
ondas de pressão.
O presente trabalho possui como objetivo a análise numérica do fenômeno de golpe
de aríete resultante do fechamento repentino de uma válvula. Para tanto, as equações
governantes são resolvidas pelo método de volumes finitos, em um código
computacional desenvolvido nesse projeto. As equações discretas do método de
volumes finitos são resolvidas com uma versão TVD (Total Variation Diminishing) do
esquema de Fluxo Médio Ponderado (WAF), descrito por TORO (1997). Esse é um
método explícito de volumes finitos de segunda ordem no tempo e no espaço. Para a
dedução das equações que governam o fenômeno, é usado o princípio de conservação
da massa e da quantidade de movimento linear.
Outra motivação para o estudo do golpe de aríete é o seu possível uso na detecção
de vazamentos nas tubulações. Existem hoje diversas técnicas que buscam identificar a
presença e a localização desses vazamentos tais como métodos acústicos que utilizam
sensores capazes de medir a energia acústica gerada pelo vazamento assim como
métodos óticos, que identificam as mudanças de temperatura na parede do duto
causadas pelo vazamento (MURVAY e SILEA, 2012).A presença de um vazamento no
duto onde ocorre o golpe de aríete causa uma atenuação no sinal gerado pelas ondas de
pressão, através do escape do fluido pressurizado. Como mostra COLOMBO et al
2
(2009), uma das técnicas para a detecção de vazamentos se baseia justamente na análise
dessas ondas de pressão. No trabalho de COLOMBO et al. (2009), foi feito um
experimento onde foi analisado o sinal obtido pelas ondas de pressão geradas pelo
fechamento de uma válvula com e sem vazamento. Os sinais transientes de pressão ao
longo da tubulação, com e sem vazamento, apresentam comportamentos distintos. Além
disso, o sinal transiente das ondas de pressão no sistema com vazamento decai mais
rapidamente do que aquele sem vazamento, conforme ilustrado na Figura 1
(COLOMBO et al., 2009)
Figura 1: Comparação entre os sinais para os casos com e sem vazamento (adaptado de
COLOMBO et al., 2009)
Dessa forma, a aplicação de métodos numéricos que sejam capazes de resolver as
equações do golpe de aríete pode ser uma importante ferramenta na detecção de
vazamentos, através da estimativa do comportamento das ondas de pressão. Entretanto,
tal análise não é realizada neste trabalho.
O trabalho está dividido da seguinte maneira: no capítulo 2, é feita uma revisão
bibliográfica do fenômeno estudado; no capítulo 3, são deduzidas as equações do
transiente hidráulico através das equações de conservação de massa e quantidade de
movimento linear; no capítulo 4, a natureza do problema é analisada e o método
numérico para a simulação do escoamento é aplicado; no capítulo 5 são apresentados e
analisados casos testes para a verificação e validação da simulação. Finalmente, no
capítulo 6 são apresentadas as conclusões obtidas e sugestões para trabalhos futuros.
3
2 – Revisão Bibliográfica
2.1 Contexto Histórico
No final do século XVIII e começo do Século XIX ocorreu na Europa o processo
conhecido como Revolução Industrial, onde aconteceu a substituição da produção
artesanal pela produção em série realizada por trabalhadores assalariados com o uso
predominantemente de máquinas (COTRIM, 2008). Com o avanço do processo de
industrialização ao longo do século XIX e início do século XX, a demanda por energia
pelas fábricas aumentou, tornando-se necessário o acréscimo do fornecimento de
energia elétrica para o abastecimento dessa necessidade. Além disso, era preciso atender
também a crescente população urbana. A geração hidroelétrica foi uma das responsáveis
pelo o abastecimento energético durante esse período. (GHIDAOUI et al., 2002)
Companhias hidroelétricas contribuíram fortemente para o desenvolvimento de
laboratórios para o estudo de fluidos e turbo máquinas, que analisavam entre outros
fenômenos, o processo de golpe de aríete. Allievi, um dos primeiros pesquisadores a
estudar esse processo, realizou os seus experimentos através do resultado de falhas e
incidentes causados pelo aumento de pressão devido ao rápido fechamento de válvulas
em usinas de potência no norte da Itália (GHIDAOUI et al., 2005)
Como mostra LÜDECKE et al. (2006) e GUIDAOUI et al. (2005), a equação mais
conhecida para modelar o fenômeno do transiente hidráulico foi desenvolvida por
Joukowsky através de uma série de experimentos realizados em sistemas hidráulicos de
Moscou em 1897. Esse modelo é conhecido na literatura como equação fundamental do
golpe de aríete ou simplesmente equação de Joukowsky, sendo mostrada na Equação
(2.1) a seguir
p a u
onde Δ𝑝 se refere a diferença de pressão, 𝜌 é a massa específica do fluido e Δ𝑢 é a
variação de velocidade longitudinal do fluido. O parâmetro a representa a velocidade de
propagação da onda de pressão no fluido, também conhecida como celeridade
(SZYDLOWSKI, 2002).
A teoria geral do golpe de aríete, entretanto, só foi desenvolvida em 1902 por
Allievi, sendo esse autor responsável também pela produção de cartas para os picos de
pressão na seção da válvula devido ao seu fechamento uniforme (GHIDAOUI et al.,
2005). Com o passar do tempo, o fenômeno do transiente hidráulico foi sendo melhor
estudado e as suas equações governantes foram modificadas até resultar nas equações de
conservação de massa e momento mostradas no capitulo 3.
(2.1)
4
2.2 A Importância da Prevenção do Golpe de Aríete
Os riscos que o golpe de aríete pode gerar nas tubulações são consideráveis. A Figura 2
mostra uma tubulação de 12 mm de espessura complemente destruída pelo transiente
hidráulico. Dessa forma, um melhor entendimento das consequências causadas pelo
golpe de aríete é necessário para manter os dutos e todos os componentes do sistema
hidráulico em um bom estado de funcionamento.
Figura 2: Tubulação com paredes de 12 mm de espessura completamente destruída pelo Golpe
de Aríete (LÜDECKE et al, 2006)
A seguir, vemos alguns dos problemas que o golpe de aríete causar nas tubulações:
2.2.1 Altas e Baixas Pressões
Como mostra BOULOS et al. (2005) e LÜDECKE et al. (2006), o fenômeno do
transiente hidráulico gera pontos de pressão máxima e mínima nos dutos. Pressões
máximas durante o regime transiente podem destruir as tubulações, as suas fixações,
válvulas, bombas entre outros componentes do sistema hidráulico. Muitas vezes, o dano
causado pelas ondas de pressão não é visto no mesmo momento em que ocorre o
fenômeno, sendo que as suas consequências não se tornam aparentes até muito tempo
depois do evento. A Figura 3 mostra uma válvula de retenção destruída devido à
sobrepressão.
Figura 3: Válvula de retenção DN800 destruída devido à sobrepresão (LÜDECKE et al, 2006)
5
Baixas pressões também podem ser extremamente danosas ao sistema hidráulico.
Ao se atingir valores menores que a pressão atmosférica, pode ocorrer o colapso de
tubos de parede fina ou de seções reforças de concreto, além da flambagem dos dutos
(BOULOS et al, 2005 e LÜDECKE et al., 2006).
Outro fenômeno relacionado às baixas pressões na tubulação é o processo de
cavitação. Segundo MATTOS et al. (1998), a cavitação ocorre quando a pressão
absoluta em qualquer ponto do sistema hidráulico atinge valor igual ou inferior à
pressão de vapor do líquido, fazendo com que parte do líquido se vaporize. Esse
processo gera uma mistura de líquido e bolhas dentro da tubulação. Quando essa
mistura atingir algum ponto do escoamento onde a pressão absoluta for novamente
superior à pressão de vapor do líquido, haverá o colapso das bolhas com o retorno na
fase líquida. Como o volume específico do líquido é inferior ao volume específico do
vapor, o colapso das bolhas implicará a existência de um vazio, causando o surgimento
de ondas de choque.
Algumas consequências do fenômeno de cavitação para os sistemas hidráulicos são:
vibrações indesejáveis, diminuição do rendimento da bomba, acelerar o processo de
corrosão, além de gerar flexões nos dutos (MATTOS et al., 1998 e ELBASHIR et al.,
2007)
2.2.2 Vibrações
Como mostra BOULOS et al. (2005), o golpe de aríete pode gerar fortes vibrações nos
sistemas hidráulicos, prejudicando os suportes das tubulações, danificando os
equipamentos de controle e medição e também o concreto presentes em dutos
revestidos. Essas vibrações se continuarem por muito tempo, causarão esforços cíclicos
que podem acarretar no surgimento de falhas por fadiga na tubulação.
É preciso estar atento às frequências em que ocorrem as vibrações para que elas não
se aproximem da frequência natural de vibração dos dutos, evitando assim que aconteça
a ressonância e isso seja capaz de destruir todo o sistema hidráulico (BOULOS et al.,
2005).
2.2.3 Poluição da Água
A qualidade da água na tubulação que sofreu o golpe de aríete pode ser também muito
prejudicada. Segundo BOULOS et al. (2005), os eventos transientes podem gerar fortes
cisalhamentos entre o fluido e as paredes da tubulação, causando a raspagem de
partículas presentes nos dutos.
Níveis baixos de pressão ajudam na contaminação da água ao gerar vazamentos nos
dutos ou nas juntas das tubulações, permitindo a intrusão da água do solo ou mesmo de
contaminantes (BOULOS et al., 2005 e GHIDAOUI et al., 2005). Além disso, pontos de
baixa pressão na tubulação induzem o retorno da água no sistema hidráulico, podendo
ocorrer a mistura entre a água doméstica e a água industrial nos sistemas de distribuição
(BOULOS et al., 2005).
6
2.3 Equipamentos de Proteção
Tendo em vista as consequências que o golpe de aríete pode gerar nos sistemas
hidráulicos, alguns equipamentos são usados com o objetivo de minimizar esses efeitos.
2.3.1 Volante de Inércia
Os Volantes de Inércia são instalados na bomba com o objetivo de aumentar o seu
tempo de parada quando ela for desligada, garantindo a sua rotação por mais alguns
instantes. Dessa forma, a variação da velocidade do escoamento será menor, diminuindo
a intensidade do golpe de aríete. (CAMARGO, 1989 e MENDES, 2011)
Uma desvantagem desse método se deve ao aumento de peso da bomba,
necessitando maior potência do motor para acioná-la e consequentemente maior
intensidade de corrente elétrica.
2.3.2 Ventosas
As ventosas são dispositivos que atuam contra os baixos valores de pressão na
tubulação, através da entrada de ar nos dutos, limitando os valores da pressão ao valor
da pressão atmosférica (CAMARGO, 1989) Um tipo comum desse dispositivo é a
ventosa com flutuador esférico conforme mostrado na Figura 4.
Figura 4: Ventosa com flutuador esférico (CAMARGO, 1989).
Conforme explica CAMARGO (1989), através da pressurização da linha, o
flutuador esférico será deslocado para cima na direção do orifício de passagem de ar.
Ao ocorrer o golpe de aríete e o nível da pressão onde está instalada a ventosa cair,
haverá a redução do nível d’água causando a descida do flutuador. Isso abrirá o orifício,
permitindo a entrada de ar na tubulação, evitando regiões de vácuo e impedindo o
colapso do duto.
2.3.3 Chaminés de Equilíbrio
Chaminés de Equilíbrio são dispositivos que permitem atenuar as variações bruscas de
pressão através do fornecimento ou do armazenamento de água. Ao se reduzir a pressão
na tubulação, a água presente na chaminé irá alimentar o sistema, No caso de ocorrer
picos de pressão na tubulação, a água será direcionada da linha para o reservatório,
permitindo o enchimento da chaminé. O topo desse dispositivo é aberto para a
7
atmosfera e normalmente as chaminés são usadas em dutos de alimentação de
hidroelétricas (MENDES, 2011)
A Figura 5 abaixo mostra um esquema de como é uma chaminé de equilíbrio
Figura 5: Chaminé de equilíbrio (adaptado de CAMARGO, 1989).
2.3.4 Reservatórios de Ar comprimido
Reservatórios de ar comprido ou do inglês Air Vessels são reservatórios fechados que
contém água no nível inferior e ar comprido no nível superior, sendo a massa de ar
dentro do reservatório regulada por manômetros e compressores. Esses reservatórios
possuem praticamente a mesma função das chaminés de equilíbrio, atuando nos níveis
altos e baixos de pressão (CAMARGO, 1989 e MENDES, 2011).
2.3.5 Válvulas de Retenção
As válvulas de retenção são dispositivos que impedem que o fluxo se inverta dentro da
tubulação. Dessa forma, atuam contra os altos níveis de pressão causados pelo golpe de
aríete (CAMARGO, 1989). Uma válvula de retenção do tipo portinhola é mostrada na
Figura 6
Figura 6: Válvula de retenção com portinhola (CAMARGO, 1989).
8
Quando o fluxo ocorre da esquerda para a direita, a portinhola se encontra aberta
permitindo assim a passagem normal de fluido. Ao se inverter o sentido do escoamento,
a portinhola irá se fechar, bloqueando fluxo.
2.3.6 Válvula Limitadora de Pressão
As válvulas limitadoras de pressão ou válvula de alívio são válvulas que possuem como
função básica limitar a pressão do sistema hidráulico. Ao se atingir determinado valor
de pressão pré-fixado, essa válvula é aberta, permitindo a saída do fluido até a pressão
retornar ao valor estabelecido (RACINE HIDRÁULICA, 1981 e CAMARGO, 1989).
2.4 Métodos Numéricos
Devido às consequências mostradas que o golpe de aríete pode causar nos sistemas
hidráulicos, há a necessidade de controlar e prever continuamente esse fenômeno. O uso
de métodos numéricos é uma forma de tentar prever como ocorrem as variações de
pressão durante o processo do transiente hidráulico.
Nos trabalhos de STREETER et al. (1968) e (1982), é mostrado alguns métodos
para a solução das equações do golpe de aríete. Os autores recomendam o uso do
método das características ou MOC na abreviatura em inglês, devido a sua simplicidade
na programação, assim como sua boa eficiência e acurácia.
No MOC, as equações diferenciais parciais do golpe de aríete são transformadas em
equações ordinárias ao longo de linhas características como é mostrado em GHIDAOUI
et al (1994). Como explica TORO (1997), as linhas características são curvas y=y(t) no
plano y-t onde as equações diferenciais parciais se transformam em equações
diferenciais ordinárias.
Como mostra GUIDAOUI et al. (1994), soluções acuradas podem ser obtidas pelo
método das características se o número de Courant é bem próximo de 1. Para que essa
condição seja satisfeita, é necessário que a velocidade de propagação da onda seja
constante, o que, em geral, não ocorre no golpe de aríete, devido à compressibilidade do
líquido e da elasticidade da tubulação (SZYDLOWSKI, 2002). Além disso, a condição
de Courant é fortemente influenciada por bifurcações ou junções no sistema, devido à
mudança de vazão nos dutos (WANG et al., 2015).
Além do MOC, outros métodos para resolver as equações do transiente hidráulico
estão presentes na literatura. GUINOT (2000) usa o método de Godunov para resolver
as equações do golpe de aríete. Segundo o autor, a vantagem desse método está na sua
característica conservativa e o fato dele lidar bem com descontinuidades. Nesse
trabalho, é mostrado que as aproximações de segunda ordem são mais rápidas
computacionalmente e fornecem resultados mais precisos que as de primeira ordem para
o método de Godunov.
SZYDLOWSKI (2002) usa o esquema proposto por Roe para a solução do
problema de Riemann como forma de resolver as equações do golpe de aríete utilizando
9
a abordagem de volumes finitos. AFSHAR et al. (2008), aplica o método das
características de forma implícita (IMOC) como forma de contornar as limitações do
MOC convencional.
HWANG et al. (2013) busca resolver as equações do transiente hidráulico pelo que
o autor chama de método de partículas características. Através de casos testes, o autor
compara a solução obtida com as soluções numéricas dadas pelo método de Godunov e
pelo método das características. Já WANG et al.(2015) propõe um método acoplado
explícito-implícito. Nesse trabalho, é mostrado mostra as vantagens e desvantagens de
cada abordagem com justificativa para acoplá-las.
O presente trabalho tem como objetivo a solução das equações do golpe de aríete
pelo esquema WAF-TVD, que é um método explicito de volumes finitos de segunda
ordem no tempo e espaço. O uso desse método foi motivado pela sua eficiência para a
solução de escoamentos compressíveis de gás no interior de dutos, conforme foi
mostrado no trabalho de BARBOSA (2003). Nesse trabalho foi visto que o método
WAF-TVD é bastante robusto para escoamentos com descontinuidades.
O método WAF-TVD foi usado também por MADEIRA (2011) para a simulação
de escoamento de gás em dutos com o objetivo de aplicar a solução obtida na detecção
do fechamento de válvula em gasodutos, através do filtro de partículas SIR. Esse
trabalho foi continuado por LUCUMI (2015) através da comparação entre os resultados
obtidos pelo filtro SIR e ASIR.
10
3 - Modelo Matemático
Para o estudo das equações básicas do golpe de aríete, será usado o princípio da
conservação da massa e da quantidade de movimento linear, bem como a aproximação
de que o escoamento é unidimensional.
As variáveis dependentes são a pressão p e a velocidade média u na seção
transversal. As variáveis independentes são a distância x medida ao longo da tubulação
e o tempo t. Dessa maneira, p=p(x,t) e u=(x,t).
As equações gerais de conservação da massa e da quantidade de movimento linear
são descritas a seguir (STREETER et al.,1982).
3.1 Equação da Conservação da Quantidade de Movimento Linear
Para a aplicação das equações de conservação da quantidade de movimento linear
iremos analisar um elemento de fluido definido por um volume de controle diferencial,
conforme a Figura 7 abaixo.
Figura 7: Volume de controle para a dedução da equação de conservação da quantidade de
movimento linear (adaptado de STREETER et al.,1982)
Aplicando a Segunda Lei de Newton na componente paralela à velocidade do
escoamento para o elemento diferencial da Figura 7 e realizando o balanço de forças,
obtêm-se.
0 sen = A Du
pA pA pA x p x A x D x A xx x dt
(3.1)
11
No lado esquerdo da Equação (3.1), os dois primeiros termos referem-se à diferença
das forças de pressão entre a entrada e a saída do volume de controle, o terceiro termo
refere-se à força de pressão exercida na superfície lateral do volume de controle, o
quarto termo é o componente da força peso paralelo ao escoamento e o último termo é a
força de atrito na parede do duto, sendo 𝜏0 a tensão de cisalhamento. O termo à direita
do sinal de igualdade refere-se à massa do volume de controle (𝜌𝐴𝛿𝑥) multiplicada pela
derivada total da velocidade.
Dividindo todos os membros da Equação (3.1) pela massa e simplificando temos:
041 sen =
p Dug
x D dt
Expandindo o termo da aceleração
Du u uu
dt x t
Substituindo a Equação (3.3) na Equação (3.2) e agrupando os termos no mesmo
lado, obtêm-se:
041 + 0
u u pu g sen
t x x D
3.2 Equação da Conservação da Massa
Para a dedução da equação da conservação da massa, vamos utilizar como modelo a
Figura 8 abaixo
Figura 8: Volume de controle para a dedução da equação da conservação da massa (adaptado de
STREETER et al.,1982)
(3.2)
(3.3)
(3.4)
12
O princípio de conservação da massa é dado por
e s
mm m
t
onde �̇�𝑒 é a vazão mássica que entra no volume de controle e �̇�𝑠 a vazão mássica de
saída. Ao realizar o balanço de massa no elemento diferencial da Figura 8, obtêm-se
A x Au Au Au xt x
Simplificando, temos
Au x A xx t
Após o desenvolvimento da Equação (3.7) e a sua divisão pela massa (𝜌𝐴𝛿𝑥)
obtêm-se:
10
u A A u u
A x A t t x
Podemos reduzir a Equação (3.8) utilizando o conceito de derivada total, onde
1 1DA u A A
A dt A x A t
1 1D u
dt x t
Com isso, a Equação (3.8) pode ser escrita da seguinte forma:
1 10
DA D u
A dt dt x
(3.7)
(3.8)
(3.9)
(3.10)
(3.11)
(3.5)
(3.6)
13
onde o primeiro termo da Equação (3.11) se refere à elasticidade da parede do tubo e à
sua taxa de deformação devido às ondas de pressão, enquanto o segundo termo se refere
à compressibilidade do líquido.
STREETER et al. (1982) nos mostra que a taxa de variação da área pode ser
substituída da seguinte maneira
2 2
DA D Dp DD
dt eE dt
)
onde e se refere à espessura da parede do tubo e E ao módulo de elasticidade do duto.
Dessa forma, o primeiro termo da Equação (3.11) pode ser escrito como
1 DA D Dp
A dt eE dt
Pela definição do módulo de elasticidade volumétrica do fluido (STREETER et al.,
1982)
Dp DpK
D D
)
Podemos, assim, substituir o segundo termo da Equação (3.11) por
1 1D Dp
dt K dt
Aplicando as Equações (3.13) e (3.15) na Equação (3.11) temos
11 0
Dp K D u
K dt E e x
Define-se agora o parâmetro a que representa a velocidade de propagação da onda
de pressão no fluido, também conhecida como celeridade (STREETER et al., 1982 e
LOMBARDI, 2005) A equação desse parâmetro é dada por (SZYDLOWSKI, 2002):
(3.12)
(3.13)
(3.14)
(3.15)
(3.16)
14
1
1a
D
K eE
)
Escrevendo a Equação (3.16) em função da celeridade temos
210
Dp ua
dt x
Expandindo a taxa de variação da pressão, temos
Dp p pu
dt x t
Para aplicações do golpe de aríete, o termo 𝑢𝜕𝑝
𝜕𝑥 é muito menor do que o termo
𝜕𝑝
𝜕𝑥
(STREETER et al., 1982). Dessa forma, a equação da conservação da massa fica
2 0p u
at x
Sabemos que a vazão volumétrica é igual a
q Au 20)
Reescrevendo as equações de conservação da quantidade de movimento linear e da
massa, dadas pelas Equações (3.4) e (3.20), em função da vazão e considerando que a
velocidade do escoamento seja paralela ao plano horizontal de referência, obtêm-se:
2
2
0
0
40
p qa
t x A
Aq qAp
t x A D
)
Existem algumas formas de substituir a tensão de cisalhamento 𝜏0 na parede do
duto. A primeira delas é utilizando a aproximação de Darcy-Weisbach (STREETER et
al., 1982)
(3.17)
(3.18)
(3.19)
(3.20)
(3.21)
(3.22)
15
08
f u u
)
ou, em função da vazão,
0 28( )
f q q
A
onde f é o fato de atrito de Darcy-Weisbach dado pela equação transcendental abaixo
(FOX et al.2014)
1 2.512log
3.7
eD
Df Re f
)
É mostrado na literatura (AXWORTHY et al, 2000, BRUNONE et al.,1995) que as
Equações (3.23) e (3.24) para o cálculo da tensão de cisalhamento não são capazes de
proporcionar dissipação suficiente para escoamentos transientes rápidos em tubulações,
que é o caso do golpe de aríete, ou seja, a tensão de cisalhamento para esse caso é maior
que prevista pelas Equações (3.23) e (3.24), que se referem a escoamentos
hidrodinamicamente desenvolvidos.
Tendo essa questão em vista, será utilizada uma nova fórmula para calcular a força
de atrito na parede do duto. Essa fórmula está presente em WYLIE (1997), baseada no
estudo feito por BRUNONE em 1991 e é igual a
2at u
f u u u uF f a
D t x
onde o parâmetro fu é definido como o fator de atrito transiente. Adicionando essa
fórmula na equação da conservação da quantidade de movimento linear temos
2
02
u
f q qq q u uAp f a
t x A AD t x
Devido à natureza método numérico escolhido para resolver as equações do golpe
de aríete, é conveniente escrevê-las na forma vetorial da seguinte maneira:
(3.23)
(3.24)
(3.25)
(3.26)
(3.27)
16
t x
U F(U)S(U)
onde
p
q
U
)
2
2
qa
A
qAp
A
F(U)
)
Se usarmos a fórmula de Darcy-Weisbach para o atrito, o termo fonte é dado por:
0
2
f q q
AD
S(U)
e se usarmos a fórmula dada por WYLIE (1997), o termo fonte é igual a:
0
2u
f q q u uf a
AD t x
S(U)
Nessas equações, U é o vetor das variáveis de estado, F é o vetor que representa o
fluxo e S é o vetor que representa o termo fonte. Podemos ver que F e S são funções do
vetor de variáveis de estado, ou seja, F=F(U) e S=S(U).
(3.28a)
(3.28b)
(3.28c)
(3.28d)
(3.28e)
17
4 - Método WAF-TVD
4.1 Estudo da Natureza do Problema
Um sistema de equações diferenciais parciais pode ser classificado das seguintes
maneiras (OZISIK et al., 2017):
Sistema Elíptico
Sistema Parabólico
Sistema Hiperbólico
Sistema Misto
Para analisar a natureza do problema do golpe de aríete, escrevemos o seu sistema
de equações na forma quasi-linear (TORO, 1997)
t x
U F US(U)
U
ou
t x
U UJ S(U)
onde J é a matriz Jacobiana do sistema que é igual a
1 1
1 2
2 2
1 2
F F
U U
F F
U U
FJ
U
Para o presente caso formulado pelas Equações (3.28), a matriz jacobiana é dada
por
2
0
2
a
A
qA
A
J
(4.1)
(4.2)
(4.3)
(4.4)
18
cujo autovalores são dados por:
2 2
.1
2 2
.2
u u a
u u a
Podemos ver claramente que os dois autovalores são sempre reais e distintos,
portanto, o sistema pode ser classificado como hiperbólico (TORO, 1997 e OZISIK et
al., 2017).
4.2 Fracionamento das equações Não Homogêneas
A presença do termo fonte na Equação (3.27a) a torna complicada para se aplicar os
métodos numéricos. Uma forma de contornar esse problema é fracionar o sistema não
homogêneo em dois: um sistema homogêneo de equações diferenciais parciais (EDP) e
um sistema não homogêneo de equações diferenciais ordinárias (EDO) (TORO, 1997,
MADEIRA, 2011 e OZISIK et al, 2017)
O sistema é fracionado em cada intervalo de tempo Δ𝑡 com a condição inicial dada
pela solução no tempo anterior, ou seja,
nU U
onde 𝐔 representa o vetor de variáveis de estado e 𝑡𝑛 = 𝑛Δ𝑡. A cada tempo 𝑡𝑛+1, 𝐔n
evolui para 𝐔n+1 com o passo de tempo Δ𝑡 = 𝑡𝑛+1 − 𝑡𝑛
Um fracionamento de primeira ordem é dado por (TORO, 1997 e OZISIK et al.,
2017).
10 n
n
t x
U F(U)
U
U U
1
1
( )n
n
d
dt
US U
U
U U
Podemos ver que a solução do Sistema (4.7) é a condição inicial do sistema (4.8).
Outra forma de reescrever esse fracionamento é
1 t tn nS Z U U
(4.5)
(4.6)
(4.7)
(4.8)
(4.9)
19
onde 𝑍(∆𝑡) e 𝑆(∆𝑡) representam os operadores do sistema de equações diferenciais
parciais (4.7) e o sistema de equações diferenciais ordinárias (4.8) respectivamente. O
sistema de equações representado pelo operador 𝑆(∆𝑡) pode ser resolvido através da
integração da Equação (4.8), obtendo a seguinte expressão:
11 nn t U U S
Para o sistema representado por 𝑍(∆𝑡) , iremos descrever o método utilizado na
próxima seção.
4.3 Esquema de Fluxo Médio Ponderado (WAF)
O método WAF (Weighted Avarage Flux) ou em português, Fluxo Médio Ponderado, é
uma generalização do método de Lax-Wendroff e do método upwind de primeira ordem
de Godunov para sistemas não lineares de leis conservativas. O método WAF tem suas
origens no esquema de fluxo randômico, que é de segunda ordem no tempo e no espaço,
sendo sua abordagem determinística (TORO, 1997)
Será considerada agora a integração da equação homogênea correspondente ao
sistema dado pela Equação (4.7), mostrada na Equação (4.11).
( )0
t x
U F U
Na abordagem de volumes finitos, dividimos o nosso domínio em regiões
chamadas de volumes de controle ou células de acordo com a Figura 9.
Figura 9: Discretização em volumes finitos para o caso unidimensional (adaptado de LUCUMI,
2011)
(4.11)
(4.10)
20
Considerando o domínio [𝑥 −Δ𝑥
2, 𝑥 +
Δ𝑥
2] x [𝑡, 𝑡 + Δ𝑡], sendo Δ𝑥 e Δ𝑡 os intervalos
espaciais e temporais, podemos escrever a forma integral da Equação (4.11) da seguinte
maneira
/2
/2
( , ) ( , ) ( / 2, ) ( / 2, ) 0
x x t t
x x t
x t t x t dx x x t x x t dt
U U F F
1
onde 𝐅(𝑥, 𝑡) = 𝐅(𝐔(𝑥, 𝑡))
O método WAF assume propriedades constantes localmente, ou seja, 𝐔𝑖𝑛 é uma
média integral da solução 𝐔(𝑥, 𝑡𝑛) na célula [𝑥𝑖−
1
2
, 𝑥𝑖+
1
2
] no instante 𝑡 = 𝑡𝑛
/2
/2
1( , )
x x
n n
i
x x
x t dxx
U U
. 12)
Aplicando a Equação (4.13) na Equação (4.12), temos
1 1( / 2, ) ( / 2, )
t t
n n
i i
t
x x t x x t dtx
U U F F
Um esquema numérico explícito para resolver a integral acima pode ser definido
como
1
1 1
2 2
n n n n
i ii i
t
x
U U F F
)
onde 𝐅𝑖+
1
2
𝑛 é uma aproximação numérica para o fluxo 𝐅 calculado na fronteira x=𝑖 +1
2
do volume de controle.
(4.12)
(4.13)
(4.14)
(4.15)
21
4.4 Esquemas TVD
Segundo TORO (1997), os esquemas TVD são uma classe de esquemas de variação
total estáveis e eles são baseados na condição de que a variação total da solução
numérica não aumenta com o tempo. A variação total pode ser definida para uma
solução discreta da seguinte maneira (TORO, 1997):
1( )n n n
i i
i
TV
onde 𝜉 é uma função 𝜉 = 𝜉(𝑥).
Dessa forma, um esquema TVD pode ser definido matematicamente como:
1( ) ( ), n nTV TV n
4.5 Método de Solução
Métodos WAF são de segunda ordem no espaço e no tempo. Nesta seção, aplicamos a
condição TVD no método WAF. Para a aplicação da restrição TVD em sistemas não
lineares, foi observado empiricamente que a solução para sistemas completos pode ser
caracterizada por saltos de uma variável 𝑏, através das ondas do problema de Riemann
(TORO, 1997).
A versão TVD do fluxo WAF é
1 1 1 1
12 2 2
1 1s
2 2
Nk k
i i ki i i
k
ign c
F F F F
)
onde N(=2) é o número de ondas do problema de Riemann para o sistema dado pela
Equação (3.28)
1 1
2 2
k k
i ir
)
é a função limitadora e
kk
k
csign c
c
é a função sinal
(4.16)
(4.17)
(4.18)
(4.19)
(4.20)
22
O parâmetro 𝑐𝑘 representa o número de Courant para as ondas k=1 e k=2 de
velocidade 𝑆𝑘 e é dado pela equação a seguir
kk
tSc
x
com 𝑐0=-1 e 𝑐𝑁+1=1 (TORO, 1997)
O salto do fluxo através da onda k é dado por
1
1 1 11
2 2 2
k k k
i i
F F F
onde 𝐅𝑖+
1
2
(𝑘)= 𝐅(𝐔(𝑘)) e
1
2
*
3
1
n
i E
n
i D
U U U
U U
U U U
onde o cálculo para 𝐔∗ será demonstrado na próxima seção.
O parâmetro de escoamento 𝑟(𝑘) refere-se a onda k da solução 𝑈𝑖+
1
2
(𝑘) do problema de
Riemann e é dado pela razão
12
12
32
12
para 0
para 0
k
i
kk
ik
k
i
kk
i
bc
b
rb
cb
Para as equações de Euler, é recomendado usar como o valor para a variável 𝑏 a
massa específica ou a energia interna específica (TORO, 1997). Para o caso da equação
do golpe de aríete, vamos utilizar a pressão como a variável para 𝑏. Assim,
(i) Δ𝑏𝑖−1/2(𝑘)
representa o salto em 𝑝 através da onda k na solução 𝐔𝑖−
1
2
(𝑥, 𝑡) do
problema de Riemann com dados (𝐔𝑖−1, 𝐔𝑖) em um instante 𝑡𝑛
(4.21)
(4.22)
(4.23)
(4.24)
23
(ii) Δ𝑏𝑖+1/2(𝑘)
representa o salto em 𝑝 através da onda k na solução 𝐔𝑖+
1
2
(𝑥, 𝑡) do
problema de Riemann com dados (𝐔𝑖, 𝐔𝑖+1) em um instante 𝑡𝑛
(iii) Δ𝑏𝑖+3/2(𝑘)
representa o salto em 𝑝 através da onda k na solução 𝐔𝑖+
3
2
(𝑥, 𝑡) do
problema de Riemann com dados (𝐔𝑖+1, 𝐔𝑖+2) em um instante 𝑡𝑛
A função limitadora 𝜙𝑖+1/2(𝑘)
na Equação (4.18) tem como função garantir que o
esquema numérico seja do tipo TVD. TORO (1997) nos sugere algumas funções como:
Superbee, Minbee e Van Leer. Como demonstrado pelos trabalhos de BARBOSA
(2003), MADEIRA (2011) e LUCUMI (2015), vamos utilizar a função Van Leer, que
funciona como um limitador intermediário. Essa função é representada pela equação
abaixo
1 se r 0
, 1 21 se r >0
1
VL r c c r
r
4.6 Solução Aproximada do Problema de Riemann
Como visto na seção anterior, para implementarmos o método WAF-TVD, precisamos
revolver o Problema de Riemann. Originalmente, o problema de Riemann referia-se ao
chamado problema do tubo de choque onde dois gases estacionários com diferentes
pressões são separados por um diafragma. A ruptura desse diafragma causa uma onda
de compressão para o lado de mais baixa pressão e uma onda de expansão (choque) no
sentido oposto (TORO, 1997 e LUCUMI, 2015).
O problema de Riemann foi generalizado para a solução de qualquer problema
unidimensional hiperbólico onde as condições iniciais são dadas por
(0) 0
( ,0) ( ) 0
E
D
se xx x
se x
UU U
U
onde 𝐔𝑑 e 𝐔𝑒 são dois vetores representados pela Figura 10
(3.24)
(4.26)
24
Figura 10: Condições iniciais para o problema de Riemann (LUCUMI, 2015).
Para a solução aproximada do problema de Riemann, TORO (1997) recomenda o
uso do método HLL desenvolvido por Harten, Lax e Van Leer. Nesse método, a
aproximação para o fluxo na fronteira da célula é obtido diretamente. A ideia central do
método assume, para a solução, uma configuração que consiste de duas ondas separando
três estados de acordo com a Figura 11.
Figura 11: Regiões e ondas para o método HLL (TORO, 1997).
Dado um algoritmo que calcule as velocidades de onda 𝑆𝐸 e 𝑆𝐷, podemos calcular o
vetor 𝐔∗ da seguinte maneira:
* hll D D E E E D
D E
S S
S S
U U F FU U
onde 𝐔𝐸 = 𝐔 𝑖𝑛, 𝐔𝐷 = 𝐔 𝑖+1
𝑛 , 𝐅𝐸 = 𝐅(𝐔𝐸) e 𝐅𝐷 = 𝐅(𝐔𝐷)
(4.27)
25
O método que será usado para estimar as velocidades 𝑆𝐸 e 𝑆𝐷 foi proposto por
Davis (TORO,1997)
1 1
2 2
min ,
max ,
E D
E
E D
D
S
S
onde {𝜆1𝐸 , 𝜆1
𝐷} e {𝜆2𝐸 , 𝜆2
𝐷} são os autovalores 𝜆1 = 𝑢 − √𝑢2 + 𝑎2 e 𝜆2 = 𝑢 + √𝑢2 + 𝑎2
da matriz Jacobiana J para as células i e i+1.
4.7 Volumes Fictícios
Em uma solução numérica por volumes finitos na forma explícita, é necessário saber o
valor da solução em todas as células em um instante de tempo anterior. A princípio,
pode-se não saber o valor da solução nas células do contorno. Nesses casos, iremos
utilizar condições de contorno numéricas para volumes fictícios conforme mostrado na
Figura 12.
Figura 12: Volumes finitos e fictícios (LUCUMI, 2015)
LANEY (1998) apresenta uma técnica de extrapolação utilizando dois volumes
fictícios na entrada e na saída da fronteira de acordo com a seguinte equação
0 1
1 2
1
2 1
n n
n n
n n
M M
n n
M M
U U
U U
U U
U U
Para o problema estudado, vamos impor o valor da pressão na entrada e extrapolar
a vazão, enquanto que para a saída, impõe-se o valor da vazão e extrapola-se a pressão.
Essa escolha foi motivada pela configuração dos casos estudados no capítulo 5, onde
temos um tanque na entrada e uma válvula que será fechada na saída.
(4.28)
(4.29)
26
4.8 Critério de Estabilidade
Para que o sistema de solução explícito seja estável, foi escolhido utilizar o critério
proposto por TORO (1997) sendo o incremento de tempo Δ𝑡 dado pela fórmula
max
cfl n
xt C
S
onde Δ𝑥 é comprimento da célula de discretização, 𝐶𝑐𝑓𝑙 é o coeficiente CFL (Courant-
Friedrichs-Lewy) com 𝐶𝑐𝑓𝑙 ∈ (0,1]. 𝑆𝑚𝑎𝑥(𝑛)
é a velocidade de propagação da onda mais
rápida no instante tempo 𝑛 dada por (TORO,1997)
max max
n n n
i iS u a
sendo 𝑖 é a posição do volume finito, incluindo as condições de contorno. 𝑢𝑖𝑛 é a
velocidade do escoamento no local e 𝑎𝑖𝑛 é a velocidade do som local. No caso do
problema do golpe de aríete, 𝑎𝑖𝑛 é a celeridade local.
(4.30)
(4.31)
27
5 - Resultados
5.1 Testes e verificações
Para a verificação do código implementado, vamos analisar dois casos testes estudados
por SZYDLOWSKI (2002). A grande diferença entre esses dois casos é o fato de que o
primeiro não inclui a variação de diâmetro da tubulação e a variação da massa
específica causada pelas ondas de pressão e também considera que o tubo é liso, ou seja,
não há atrito.
Os códigos computacionais desse estudo foram implementados na plataforma
Matlab®, executados em um computador com processador Intel
® Xeon
® de dois núcleos
de 2.40 GHz e memória RAM de 32 GB.
5.2 Caso 1
Para o primeiro caso analisado, temos um fluxo d’água dentro de um duto liso de 500
metros de comprimento e diâmetro igual a 0.1 metros. Inicialmente a pressão em toda
tubulação é igual a 0.5 MPa e a velocidade é igual a 0.4 m/s. Após 0.5 segundos, a
válvula presente no final da tubulação é completamente fechada. Devido a esse
fechamento, surge dentro da tubulação o fenômeno do golpe de Aríete.
Figura 13: Desenho esquemático para o Caso 1
O escoamento foi considerado incompressível com 𝜌=1000 kg/m3 e a parede do
duto rígida, ou seja, a área transversal da tubulação permanece constante durante o
transiente hidráulico. Além disso, a celeridade 𝑎 foi considerada constante e igual a
1000 m/s. As condições do escoamento podem ser resumidas na tabela abaixo
28
Tabela 1: Parâmetros do escoamento para o Caso 1
Pressão Inicial 𝑝0 0.5 𝑀𝑃𝑎
Velocidade Inicial 𝑢0 0.5 𝑚/𝑠
Massa Específica 𝜌 1000 𝑘𝑔/𝑚3
Celeridade 𝑎 1000 𝑚/𝑠
Comprimento do Duto 𝐿 500 𝑚
Diâmetro do duto 𝐷 0.1 𝑚
As condições de contorno para esse problema, segundo Szydlowski (2002), são
constantes durante todo o escoamento. Devido ao fechamento da válvula, a velocidade a
jusante da válvula é igual a 𝑢 = 0 𝑚/𝑠 Para o começo do duto, a pressão é mantida
igual à pressão inicial.
O resultado analítico para os valores de pressão na seção da válvula foi obtido por
SZYDLOWSKI (2002), sendo mostrado na Figura 14
Figura 14: Resultado Analítico (SZYDLOWSKI, 2002)
5.3 Solução para o Caso 1
Para a verificação da solução do Caso 1 utilizando o método WAF-TVD, foi realizada a
simulação do escoamento para 25, 50, 100, 200, 500 e 1000 volumes, com o objetivo de
29
verificar a convergência de malha. Os resultados para esses volumes de controle são
apresentados nas Figuras 15 a 19, respectivamente.
Figura 15: Comparação do resultado analítico (linha vermelha) com o resultado obtido pelo
método WAF TVD (linha azul) para 25 volumes de controle
Figura 16: Comparação do resultado analítico (linha vermelha) com o resultado obtido pelo
método WAF TVD (linha azul) para 50 volumes de controle
30
Figura 17: Comparação do resultado analítico (linha vermelha) com o resultado obtido pelo
método WAF TVD (linha azul) para 100 volumes de controle
Figura 18: Comparação do resultado analítico (linha vermelha) com o resultado obtido pelo
método WAF TVD (linha azul) para 200 volumes de controle
31
Figura 19: Comparação do resultado analítico (linha vermelha) com o resultado obtido pelo
método WAF TVD (linha azul) para 500 volumes de controle
Figura 20: Comparação do resultado analítico (linha vermelha) com o resultado obtido pelo
método WAF TVD (linha azul) para 1000 volumes de controle
32
Analisando a Figuras 15-20, podemos ver que conforme aumentamos o número de
volumes, melhora a concordância entre o resultado analítico e experimental. A solução
com 100 volumes já consegue recuperar bem a magnitude e a forma das ondas de
pressão, entretanto, o resultado numérico encontra-se um pouco defasado em relação ao
resultado analítico.
Na solução com 500 volumes, já conseguimos ter um resultado bastante
satisfatório, com uma defasagem bem pequena entre os resultados analítico e numérico,
mostrando a robustez do método para resolver as equações do golpe de aríete para o
caso sem atrito, incompressível e com a tubulação rígida.
Na Figura 21, comparamos os resultados obtidos pela solução com 500 e 1000
volumes de controle, onde podemos ver qualitativamente a convergência gráfica entre
ambas soluções.
Figura 21: Comparação entre os resultados numéricos obtidos pela solução com 500 volumes de
controle (linha vermelha) e 1000 volumes de controle (linha azul).
5.4 Caso 2
Nesse segundo caso para fim de validação do código desenvolvido, os resultados
numéricos são comparados com os resultados experimentais realizados por
Lewandowski e Makowski (SZYDLOWSKI, 2002). O caso analisado é semelhante ao
anterior.
A principal diferença para esse segundo caso é que, agora, consideramos a
elasticidade do duto e a variação da massa específica do líquido devido ao transiente
hidráulico. O material do tubo é de aço com o módulo de elasticidade igual a 2x1011
Pa,
33
com comprimento de 72 metros, diâmetro inicial igual a 0.042 metros e a espessura da
parede igual a 0.003 metros (SZYDLOWSKI, 2002).
Para quantificarmos a alteração da área transversal do duto, será usada a seguinte
aproximação para o calculo da área presente em SZYDLOWSKI (2002)
00 1
p pDA A
e E
1)
onde 𝐴0 é a área transversal do duto antes de ocorrer o golpe de aríete, 𝑒 refere-se a
espessura da parede do duto, 𝐸 é o módulo de elasticidade do material e 𝑝0 é a pressão
inicial
Para a variação da massa específica do líquido, temos (SZYDLOWSKI, 2002)
00 1
p p
K
.2)
onde K é o módulo de elasticidade da água e 𝜌0 a massa específica inicial. Para esse
problema 𝐾 = 2𝑥109 Pa e 𝜌0 = 1000 𝑘𝑔/𝑚3
Inicialmente, a válvula no final da tubulação está completamente aberta e a água
escoa pela tubulação com velocidade igual a 0.38 m/s por todo o duto. Na região da
válvula a pressão é igual a 0.51 MPa e é aumentada ao longo da tubulação pelo seguinte
acréscimo de pressão resultante do atrito do fluido com a parede do duto
(SZYDLOWSKI, 2002)
0
2
0
0.51
2
p p
uLp f
D
)
onde 𝑢0 é a velocidade inicial e Δ𝐿 é a distância da válvula até uma posição específica
da tubulação.
A partir do instante 𝑡 = 0.38 s, é iniciado o fechamento da válvula até o tempo de
𝑡 = 0.494 𝑠. A partir desse instante, a válvula encontra-se totalmente fechada. É
utilizada uma aproximação linear para a velocidade durante o processo de fechamento
da válvula, que é dado por (SZYDLOWSKI, 2002)
(5.1)
(5.2)
(5.3)
34
0 0
0.38
0.494 0.38
tu u u
Os parâmetros usados na simulação são resumidos na Tabela 2
Tabela 2: Parâmetros para o escoamento para o Caso 2
Pressão Inicial 𝑝0 Eq. (5.3)
Velocidade Inicial 𝑢0 0.4 𝑚/𝑠
Massa Específica inicial 𝜌0 1000 𝑘𝑔/𝑚3
Comprimento do Duto 𝐿 72 𝑚
Diâmetro do duto 𝐷 0.042 𝑚
Espessura da parede 𝑒 0.003 𝑚
Módulo de elasticidade do aço 𝐸 2𝑥1011 𝑃𝑎
Módulo de elasticidade da água 𝐾 2𝑥109 𝑃𝑎
As condições de contorno para as duas extremidades do duto são (SZYDLOWSKI,
2002)
Extremidade oposta à válvula: Pressão igual à pressão inicial
Extremidade da válvula: velocidade igual a zero quando a válvula está
totalmente fechada e igual ao resultado da Equação (5.4) quando a válvula está
no processo de fechamento
O resultado experimental para os valores de pressão na posição da válvula é
mostrado na Figura 22 (SZYDLOWSKI, 2002)
Figura 22: Resultado experimental (SZYDLOWSKI, 2002).
(5.4)
35
5.5 Solução para o Caso 2
Para a solução do problema, vamos considerar as duas abordagens do fator de atrito
como foi descrito no capítulo 3, usando a fórmula de Darcy-Weisbach e a fórmula
presente em WYLIE (1997), dada pelas equações (3.25) e (3.26)
Na fórmula do fator de atrito presente em WYLIE (1997) existem dois termos com
derivadas espaciais e temporais da velocidade. A presença desses termos modifica as
equações do golpe de aríete e seria necessário reescrevê-las. Uma forma de contornar
essa situação é aproximar essas derivadas explicitamente por diferenças finitas usando a
diferença atrasada para o tempo
1n n
i iu uu
t t
e usando a diferença de quarta ordem no espaço. Nos volumes de controle centrais
usamos a diferença centrada conforme a Equação (5.6)
2 1 1 2
1 2 2 1
12 3 3 12
n n n n
i i i iu u u uu
x x
Nos contornos, usamos a diferença avançada na entrada e a diferença atrasada na
saída conforme as Equações (5.7) e (5.8), respectivamente.
1 2 3 4
25 4 14 3
12 3 4
n n n n n
i i i i iu u u u uu
x x
1 2 3 4
25 4 14 3
12 3 4
n n n n n
i i i i iu u u u uu
x x
Como no primeiro caso estudado, era necessário 500 volumes de controle para se
obter a convergência, vamos usar esse número de volumes para resolver o escoamento
utilizando a fórmula de Darcy-Weisbach para o fator de atito. A Figura 23 mostra esse
resultado na posição da válvula
(5.5)
(5.6)
(5.7)
(5.8)
36
Figura 23: Comparação ente os resultados experimentais e numérico usando a fórmula de
Darcy-Weisbach para o atrito.
A Figura 23 mostra que o fator de atrito dado pela fórmula de Darcy-Weisbach não
é capaz de dissipar energia suficiente para escoamentos transientes rápidos.
Sabendo disso, aplicamos agora a fórmula do fator de atrito apresentada em
WYLIE (1997) nas equações do golpe de aríete e usamos a aproximação dada pelas
Equações (5.5), (5.6), (5.7) e (5.8). WYLIE (1997) recomenda alterar o fator de atrito
transiente até se obter uma dissipação adequada para o problema analisado. No nosso
caso, utilizamos 𝑓𝑢 = 0.06.
Para análise de convergência, vamos comparar os resultados obtidos por 500, 1000
e 2000 volumes. O tempo computacional para cada simulação está presente na Tabela 3
abaixo
Tabela 3: Tempos computacionais para a realização das simulações
Número de Volumes Tempo computacional (horas)
500 4.18
1000 14.83
2000 56.46
A comparação gráfica entre as simulações está presente nas Figuras 24 a 26
37
Figura 24: Comparação entre os resultados numéricos para 500 e 1000 volumes de controle
Figura 25: Comparação entre os resultados numéricos para 1000 e 2000 volumes
38
Figura 26: Comparação entre os resultados numéricos para 500, 1000 e 2000 volumes de
controle.
Pela análise das Figuras 24-26, podemos ver qualitativamente que a convergência
gráfica foi obtida entre as simulações. Os valores numéricos entre as soluções estão
presentes nas Tabelas 4 e 5, para os valores de tempo e pressão, respectivamente.
Tabela 4: valores obtidos para o tempo (s) usando 500, 1000 e 2000 volumes de controle.
Número de volumes de controle Posição do pico 500 1000 2000
1º 0,4884 0,4894 0,4900
5º 0,9279 0,9309 0,9326
9º 1,3680 1,3730 1,3760
13º 1,8080 1,8150 1,8190
17º 2,2440 2,2570 2,2620
21º 2,6890 2,7000 2,7050
25º 3,1300 3,1420 3,1428
29º 3,5700 3,5840 3,5920
33º 4,0110 4,0270 4,0350
37º 4,4520 4,4690 4,4780
39
Tabela 5: valores obtidos para a pressão (MPa) usando 500, 1000 e 2000 volumes de controle.
Número de volumes de controle Posição do pico 500 1000 2000
1º 0,9969 1,0010 1,0030
5º 0,8898 0,8988 0,9042
9º 0,8226 0,8313 0,8367
13º 0,7688 0,7769 0,7819
17º 0,7252 0,7325 0,7370
21º 0,6895 0,6959 0,7000
25º 0,6601 0,6658 0,6694
29º 0,6358 0,6408 0,6440
33º 0,6157 0,6201 0,6229
37º 0,5990 0,6027 0,6052
Podemos observar através da análise das Tabelas 4 e 5 que os valores numéricos
obtidos entre as simulações estão bem próximos entre si. Para melhor análise, é
mostrado nas Tabelas 6 e 7 o desvio relativo percentual entre as soluções de 500 e 2000
e 1000 e 2000 volumes de controle para os valores de tempo e pressão respectivamente.
Tabela 6: Desvio relativo percentual para os valores de tempo em relação à solução numérica
para 2000 volumes de controle.
Nº de volumes
Posição do pico 500 1000 1º 0,3265 0,1225
5º 0,5039 0,1823
9º 0,5814 0,2180
13º 0,6047 0,2199
17º 0,7958 0,2210
21º 0,5915 0,1848
25º 0,4072 0,0255
29º 0,6125 0,2227
33º 0,5948 0,1987
37º 0,5806 0,2010
40
Tabela 7: Desvio relativo percentual para os valores de pressão em relação à solução numérica
para 2000 volumes de controle.
Nº de volumes Posição do pico 500 1000
1º 0,6082 0,1994
5º 1,5926 0,5972
9º 1,6852 0,6454
13º 1,6754 0,6395
17º 1,6011 0,6106
21º 1,5000 0,5857
25º 1,3893 0,5378
29º 1,2733 0,4969
33º 1,1559 0,4495
37º 1,0245 0,4131
Observando a Tabela 6, vemos que o desvio relativo entre as soluções para os
valores de tempo são inferiores a 1%, mostrando uma concordância excelente entre os
resultados. Contudo, ao se analisar a Tabela 7 vemos que o desvio relativo entre o
resultado obtido por 500 e 2000 volumes aumenta, ficando em torno de 1.5 %.
Considerando uma convergência para erros relativos inferiores a 1%, podemos
concluir que o resultado para 1000 volumes de controle é o ideal para a solução do caso
estudado, ainda mais considerando o tempo computacional elevado para se obter o
resultado com 2000 volumes de controle. Em determinadas aplicações, contudo, o
tempo computacional de aproximadamente 15 horas para terminarmos a simulação com
1000 volumes pode ser inviável. Para esses casos, é possível usar 500 volumes de
controle, que apresentou resultados com desvios relativos inferiores a 2% e um tempo
computacional de aproximadamente 4 horas.
Dessa forma, vamos comparar os resultados experimentais e os obtidos pelo
método WAF TVD usando 500 e 1000 volumes de controle. Essas comparações estão
presentes nas Figuras 27 e 28.
41
Figura 27: Comparação entre os resultados experimental e numérico para 500 volumes
de controle
Figura 28: Comparação entre os resultados experimental e numérico para 1000 volumes
de controle
42
Podemos observar pelas figuras anteriores que a fórmula presente em WYLIE
(1993) para o fator de atrito produz a dissipação desejada. Comparando com os
resultados numéricos, obtemos a dissipação de energia bem semelhante daquela vista
experimentalmente.
Ao analisarmos os picos de pressão obtidos pelas soluções numéricas, vemos que
logo quando a válvula é fechada, eles são um pouco inferiores em relação aos resultados
experimentais. Uma possível causa para esse fato seja o uso de uma função limitadora
de características intermediárias, como é o caso da função de Van Leer
A partir de 2 segundos aproximadamente, vemos que os picos de pressão medidos
experimentalmente começam a decair mais rápido que aqueles obtidos numericamente.
Isso talvez ocorra pelo modelo escolhido para calcular o atrito entre o fluido e a parede
do duto presente em WYLIE (1997). No tempo de 8 segundos em diante, os resultados
numéricos e experimentais são muito próximos, sendo as possíveis divergências
consequência da dificuldade de se obter os valores experimentais exatos.
43
6 - Conclusão e Sugestões para Trabalhos Futuros O objetivo deste trabalho é a análise e a implementação numérica para a solução das
equações do golpe de aríete. Para isso, foi utilizado a versão TVD do Fluxo Médio
Ponderado (WAF), que é um método explícito de volumes finitos, sendo de segunda
ordem no espaço e no tempo.
As equações do golpe de aríete foram deduzidas a partir das equações de
conservação da massa e da quantidade de movimento linear. Em seguida, foi verificado,
através da análise dos autovalores da matriz jacobiana, a natureza hiperbólica do
problema. Com isso, foi mostrado uma forma de fracionamento das equações
diferenciais parciais com o objetivo de implementar o método WAF TVD para a
solução do sistema homogêneo de equações diferenciais parciais.
Para a verificação e validação do código computacional, a solução numérica foi
comparada com dois casos presentes na literatura. O primeiro caso considerava o
fenômeno do transiente hidráulico em um tubo liso (sem atrito), onde a massa específica
e a área transversal da tubulação não variavam. Foi variado o número de volumes
utilizados e a solução numérica foi comparada com o resultado analítico presente na
literatura, onde pode ser verificado que o uso de 500 volumes produz um resultado
numérico bastante próximo do analítico.
O segundo caso estuda o mesmo fenômeno do caso anterior, com a diferença de
que, agora, é considerada a presença de atrito na parede interna da tubulação, além da
variação da seção transversal do duto e da massa específica do líquido causada pelas
ondas de pressão, características do golpe de aríete. Foi analisada a solução utilizando
duas abordagens para o fator de atrito bem como a convergência da malha.
Para o segundo caso analisado, foi visto que a presença do coeficiente de atito
transiente produz uma dissipação satisfatória de energia, algo que a formulação clássica
de Darcy-Weisbach não foi capaz de obter. A análise da convergência da malha mostrou
que o desvio relativo entre os valores de pressão para as soluções de 500 e 2000
volumes de controle são inferiores a 2% enquanto o desvio relativo entre os valores de
pressão para as soluções de 1000 e 2000 volumes de controle são inferiores a 1%.
Considerando o tempo computacional para realizar as simulações e tendo consciência
do desvio relativo, pode-se usar os resultados numéricos obtidos com 500 ou 1000
volumes de controle.
Comparando os resultados numéricos e experimentais para o segundo caso, foram
vistas algumas discordâncias entre os resultados. A primeira refere-se ao valor máximo
do primeiro pico de pressão. O resultado numérico foi um pouco inferior ao resultado
experimental. A segunda discordância foi o decaimento mais acelerado dos picos de
pressão medidos experimentalmente. Um possível motivo para isso deve-se ao modelo
escolhido para representar o atrito entre o fluido e a parede da tubulação. Entretanto,
devido à complexidade do problema analisado, vemos que o método WAF-TVD foi
muito satisfatório para a solução das equações do golpe de aríete.
44
Como sugestões para trabalhos futuros, podemos citar:
Aplicação do método WAF-TVD para a simulação de mais casos do
transiente hidráulico.
Realização de experimentos para melhor análise do golpe de aríete.
Aplicação do método WAF-TVD na detecção de vazamentos em
tubulações.
Pesquisa e possível desenvolvimento de novas formulações para o fator de
atrito entre a o fluido e a parede da tubulação durante o golpe de aríete.
45
7 - Referências AFSHAR, M. H., ROHANI, M., Water Hammer Simulation by Implicit Method of
Characteristics, International Journal of Pressure Vessels and Piping, Vol. 85, pp. 851-
859, 2008.
AXWORTHY, D. H., GHIDAOUI, M. S., MCINNIS, D. A., Extended
Thermodynamics Derivation of Energy Dissipation in Unsteady Pipe Flow, Journal of
Hydraulic Engineering, vol. 126(4), pp. 276-287, 2000.
BARBOSA, M. V, F. Implementaçaõ do Método WAF TVD para a Simulação de
Escoamentos Compressíveis Quasi-Unidimensionais. Projeto Final de Graduação,
DEM/EE/UFRJ, Rio de Janeiro, RJ, Brasil, 2003.
BOULOS, P. F., KARNEY, B. W., WOOD, D. J., LINGIREDDY, S, Hydraulic
Transient Guideless for Protecting Water Distribution Systems, Journal of American
Water Association, pp.111-124, 2005.
BRUNONE, B., GOLIA, U.M., GRECO, M., Effects of Two-Dimensionality on Pipe
Transients Modeling, Journal of Hydraulic Engineering, vol. 121 (12), pp. 906-912,
1995.
CAMARGO, L.A., O golpe de aríete em Tubulações de Racalque. Análise
Simplificada, XV Encontro de Engenheiros de Assistência Técnica, Joinville, 1989.
COLOMBO, A. F., LEE, P., KARNEY, B. W., A selective review of transient-based
leak detection methods, Journal of Hydro-environment, vol. 2, pp. 212-227, 2009.
COTRIM, G., Historia Global: Brasil e Geral: volume único, 9ª ed., Saraiva, 2008.
ELBASHIR, M. A. M., AMOAH, S. O. K., Hydraulic Transient in a Pipeline Using
Computer Model to Calculate and Simulate Transient, Master Tesis, Lund University,
Sweden, 2007.
FOX, R. W., MCDONALD, A. T., LEYLEGIAN, J.C., Introdução à Mecânica dos
Fluidos, 8ª ed. LTC, 2014.
GHIDAOUI, M. S., KARNEY, B. W., Equivalent Differential Equations in Fixed-Grid
Characteristics Method, Journal of Hydraulic Engineering, Vol. 120, 1994.
GHIDAOUI, M.S., ZHAO, M., MCINNIS, D. A., AXWORTHY, D. H., A Review of
Water Hammer Theory and Practice, Applied Mechanics Reviews, Vol. 58, pp.49-76,
2005.
GUINOT, V., Riemmann Solvers for Water Hammer Simulations by Gudonov Method,
International Journal for Numerical Methods in Engineering, Vol. 49, pp. 851-870,
2000
46
HWANG, Y., Development of a Characteristic Particle Method for Water Hammer
Simulation, Journal of Hydraulic Engineering, Vol. 139, pp. 1175-1192, 2013.
LANEY, C. B., Computational Gasdynamics, Cambridge, Cambridge University Press,
1998.
LOMBARDI, J.C., Análise de distribuição de pressão em válvulas de diafragma
poroso, 125 p., Tese (Doutorado Engenharia Mecânica), INPE, São José dos Campos,
2005.
LUCUMI, M. A.R. Comparação dos Algoritmos de Filtros de Partículas SIR e ASIR na
Detecção de Fechamento de Válvulas em Gasodutos. 96 p. Tese (Doutorado em
Engenheira Mecânica) – Instituto COPPE de Engenheira Mecânica, Universidade
Federal do Rio de Janeiro, Rio de Janeiro, 2015.
LÜDECKE, H. J., KOTHE, B., KSB Know- how Water Hammer, Volume 1, KSB,
2006.
MADEIRA, I.M. Detecção do fechamento de válvula em gasodutos pela utilização de
filtros de partículas. 150 p. Tese (Doutorado em Engenheira Mecânica) – Instituto
COPPE de Engenheira Mecânica, Universidade Federal do Rio de Janeiro, Rio de
Janeiro, 2011.
MATTOS, E. E., DE FALCO, R., Bombas Industriais, 2ª Ed., Interciência, 1998.
MENDES, L. F. M., Métodos Clássicos de Protecção de Sistemas Elevatórios Contra o
Golpe de Aríete, p. 95, Tese (Mestrado Engenharia Ambiental), Universidade Nova de
Lisboa, Lisboa, 2011.
MURVAY, P., SILEA, I., A survey on gas leak detection and localization techniques,
Journal of Loss Prevention in the Process Industries, Vol. 25, pp. 966-973, 2012.
OZISIK, N., ORLANDE, H. R. B., COLAÇO, M. J., COTTA, R. M., Finite Difference
Methods in Heat Transfer, 2ª Ed., CRC Press, 2017.
RACINE HIDRÁULICA, Manual de Hidráulica Básica, 3ª Ed., 1981.
STREETER, V. L. and WYLIE, E. B., Hydraulic Transients, McGraw-Hill, New York,
1967.
STREETER, V. L., WYLIE, E. B., Mecânica dos Fluidos, 7ª Ed., Editora McGraw-Hill
do Brasil, 1982.
SZYDLOWSKI, M., Finite Volume Method for Water Hammer Simulation, I
International Scientific and Technical Conference on Technology, Automation and
Control of Wastewater and Drinking Water Systems, TiASWiK'02, pp. 159-165, 2002.
TORO, E. F., Riemann Solvers and Numerical Methods for Fluid Dynamics, 2ª ed.
Berlin, Springer, 1997.
47
WANG, C., YANG, J., Water Hammer Simulation Using Explicit-Implicit Coupling
Methods, Journal of Hydraulic Engineering, Vol. 141, 2015
WYLIE, E. B., Frictional effects in unsteady turbulent pipe flows, Applied Mechanics
Reviews, Vol. 50, nº 11, pp. 241-244, 1997.
