Embed Size (px)
Citation preview
Departamento de Mecânica
1
MODELO DE DOIS FLUIDOS DE ALTA ORDEM PARA PREVER A
FORMAÇÃO DE GOLFADAS
Aluno: Fabricio Ciotti Chamma
Orientadora: Angela Ourivio Nieckele
Introdução
O transporte de um ou mais fluidos envolvendo diferentes fases em tubos está presente
em diversas áreas da engenharia, em especial na indústria petrolífera (Figura 1). A
complexidade do tema está associada à existência de diferentes arranjos das fases presentes
causando a necessidade de modelagens diferenciadas para os fenômenos de transferência
interfacial dependendo do padrão de escoamento [1]. Para a compreensão dos fenômenos
envolvidos, assim como nas etapas de projeto e de operação dos equipamentos é fundamental
a previsão acurada do escoamento.
Figura 1: Sistema de produção de petróleo em águas profundas. Fonte: Petrobras
A Figura 2 apresenta os principais tipos de padrões para escoamentos gás-líquido para
tubulações na horizontal e vertical: bolhas dispersas, estratificado, bolhas alongadas, golfada,
anular. Escoamentos horizontais apresentam adicionalmente, o padrão estratificado. A
transição de um padrão para o outro é um tema ainda em discussão teórica e experimental.
Diversos fatores podem influenciar os padrões como: a inclinação da tubulação, o diâmetro,
as velocidades das fases e suas frações mássicas.
Departamento de Mecânica
2
Figura 2: Desenho esquemático de alguns dos principais padrões de escoamento gás-líquido na configuração horizontal e vertical [2].
Dentre os diversos regimes de escoamento bifásico, o padrão de golfadas se destaca
entre os padrões observados por sua frequência nas indústrias, principalmente a petrolífera
(Fig. 3). Este padrão é o que requer maior esforço em sua caracterização e modelagem, devido
às características marcantes da distribuição espacial entre fases, que gera intermitência ao
escoamento. Este escoamento ocorre em larga faixa de vazões de gás e líquido em tubulações
de diâmetro médio e pequeno, com variação periódica da densidade, fração de vazio e
pressões na seção transversal da tubulação.
Figura 3: Padrão de escoamento de golfada em uma linha submarina
A correta previsão deste tipo de escoamento é fundamental para o projeto adequado de
tubulações, e dos reservatórios de separação das fases. O desenvolvimento de uma ferramenta
numérica também é muito útil para auxiliar na operação de dutos.
Este projeto se insere na linha de pesquisa em Escoamentos no Padrão de Golfadas que
vem sendo desenvolvida pelo Grupo de Dinâmica dos Fluidos Computacional.
Departamento de Mecânica
3
Objetivos
O presente projeto consiste em implementar um algoritmo de discretização temporal de
segunda ordem para prever numericamente o escoamento bifásico tanto no padrão estratificado
como em golfadas, utilizando o Modelo de Dois Fluidos [3]. É feita uma comparação entre o
desempenho da integração de primeira ordem (Euler implícito) e de segunda ordem (Crank-
Nicholson) [4].
Os resultados das simulações serão comparados com dados experimentais disponíveis
na literatura e com dados obtidos no Departamento de Engenharia Mecânica da PUC-Rio.
Metodologia
Para determinar os escoamentos no padrão estratificado e no padrão de golfadas,
utilizou-se o Modelo de Dois Fluidos [3]. Foi realizada a simulação numérica de diferentes
casos de escoamento água-ar. No primeiro introduziu-se uma perturbação controlada na
entrada e investigou-se o efeito do passo de tempo na taxa de amplificação das perturbações.
No segundo tipo, investigou-se o escoamento no padrão de golfadas. Neste caso, após a
obtenção do regime estatisticamente permanente, as grandezas características das golfadas são
determinadas. Em ambos os casos, utilizou-se água e ar.
Modelo Matemático
O Modelo de Dois Fluidos [3] consiste na solução de um conjunto de equações de
conservação de massa e de quantidade de movimento linear para cada uma das fases
envolvidas, assim como equações de transporte através das interfaces.
Considerou-se o escoamento unidimensional, isotérmico, líquido incompressível, e gás
ideal, ao longo de uma tubulação inclinada de β com a horizontal, conforme ilustrado na Fig. 4.
Figura 4: Escoamento estratificado
As equações de conservação para o Modelo de Dois Fluidos são obtidas através de uma
média nas fases. Adicionalmente, na formulação unidimensional, médias na seção transversal
precisam ser obtidas.
Na presente modelagem, baseada nos trabalhos [1 e 4], as equações de conservação são
obtidas, considerando o escoamento base como estratificado (Fig.4).
Conservação de Massa para o Gás e Líquido
( )
( )
(1)
( )
( )
(2)
Departamento de Mecânica
4
Conservação de Quantidade de Movimento Linear para o Gás e Líquido:
( )
(
)
(3)
( )
(
)
(4)
onde L, G e i referem-se à fase líquida, gasosa e interface. representa a velocidade média
na seção transversal, a fração volumétrica e é a massa específica. representa a pressão,
é a área da seção transversal do duto ( A = πD2/4 ), onde D é o diâmetro da tubulação, é a
inclinação do mesmo e e correspondem às coordenadas axial e temporal, respectivamente.
representa a tensão de cisalhamento com a parede, a tensão de cisalhamento da
interface gás-líquido. é o perímetro molhado da fase k e é o perímetro da interface gás-
líquido, e hL é o nível de líquido na seção transversal. A aceleração da gravidade é
representada por .
Tal modelo matemático necessita de correlações empíricas para avaliar a transferência
de quantidade de movimento linear τi através das interfaces. Como o modelo utilizado aqui é
uni-dimensional, necessita também de correlações empíricas para avaliar as tensões
cisalhantes τwk de cada fase com a parede da tubulação. As tensões são avaliadas assumindo
escoamento localmente desenvolvido, em função do fator de atrito, o qual depende do número
de Reynolds e regime de escoamento.
| | ;
| |( ) (5)
Existe uma série de correlações existentes na literatura para a determinação do fator de
atrito. Como exemplo pode-se citar: Taitel e Dukler (1976), Issa e Kempf (2003) e Issa e
Bonizzi (2003). As correlações apresentadas na Tabela 1 se mostraram adequadas por
Nieckele et al, 2013.
Tabela 1 – Fórmulas para o cálculo do fator de atrito.
ReG, ReL, Rei ≤ 2100
(Laminar)
ReG, ReL, Rei > 2100
(Turbulento)
fL
sLRe
24
(Hand, 1991)
139,0)Re(0262,0 sLl
(Spedding e Hand, 1997)
fG
GRe
16
(Hagen–Poiseuille)
2500460 ,)Re(,
G
(Taitel e Dukler, 1976)
fi
GRe
16
(Hagen–Poiseuille)
2500460 ,)Re(,
i
(Taitel e Dukler, 1976)
Departamento de Mecânica
5
Nas equações apresentadas na Tabela 1, os números de Reynolds sLRe , ReG e Rei são
definidos, de acordo com Taitel e Dukler (1976) como:
L
sLLsL
DU
Re
(6)
GiG
GGGG
SS
UA
)(
4Re
(7)
GiG
GLGG
iSS
UUA
)(
4Re
(8)
As definições apresentadas são baseadas na hipótese de que o gás escoa num canal fechado,
uma vez que viaja a uma velocidade muito maior que a do líquido, o qual, por sua vez, escoa
como se estivesse num canal aberto. Os números de Reynolds do gás e da interface são
baseados no diâmetro hidráulico do gás, calculado através de sua área de escoamento, AG, e
perímetro molhado, SG. Além disso, é a viscosidade dinâmica da fase, D é o diâmetro da
tubulação e UsL é a velocidade superficial do líquido, definida como:
LLsL UU (9)
As frações volumétricas de líquido e gás possuem a seguinte relação:
, (10)
A pressão na interface do líquido é relacionada com a pressão na interface do lado do
gás, pela tensai superficial s e raio de curvatura, o qual pode ser aproximado por
=
, (11)
O líquido é incompressível, e o gás segue a lei dos gases ideais
, (12)
onde R é a constante do gás e Tref é a temperatura de referência. A pressão do gás na interface
é considera como constante em, toda secao transversal, logo, = .
Os parâmetros geométricos necessários para o modelo, são obtidos considerando o
escoamento estratificado, pois é a partir deste que o padrão de golfadas evolui.
A seção transversal de um escoamento estratificado é representado na Fig. 5, que ilustra
os parâmetros geométricos associados a cada fase (perímetros e áreas da seção transversal de
cada fase).
As relações geométricas das Fig. 5 são apresentadas a seguir:
[ ( ) √ ] (13)
( ) ; ; √ (14)
Departamento de Mecânica
6
[ ] ;
(15)
(
) ;
(16)
Figura 5:–Geometria para escoamentos estratificados.
Método Numérico
As equações de conservação (1) a (4) são discretizadas pelo método dos Volumes
Finitos [5]. Este método tem como base a divisão do domínio computacional em volumes de
controle, e, em cada um desses volumes, são aplicadas as leis de conservação, as quais são
integradas no espaço e no tempo, garantindo assim a conservação global de todas as
grandezas que estão sendo avaliadas.
O volume de controle principal é onde todas as grandezas escalares são armazenadas
Fig. 6(a). Os pontos nodais são representados por letras maiúsculas, sendo P o nó principal, E
o vizinho leste e W o vizinho oeste. As faces são representadas por letras minúsculas, sendo e
a face leste e w a face oeste. Seguindo a recomendação de Patankar [5], as velocidades são
armazenadas em volumes de controle deslocados em relação ao volume de controle principal
Fig. 6(b), de forma que soluções oscilatórias e irrealistas sejam evitadas.
Figura 6(a) – Volume de controle principal ; Figura 6(b) – Volume de controle das velocidades
Uma vez que as equações de conservação não apresentam termos difusivos, aplica-se a
aproximação “upwind” de primeira ordem [5] nos fluxos obtidos após a integração espacial.
Para realizar a integração temporal de uma variável genérica é preciso supor a
variação de no intervalo de tempo. A Figura 7 ilustra os três caminhos mais usuais: Euler
𝛾
𝑆𝐿
𝐴𝐺
𝐴𝐿
𝑆𝐺
𝑆𝑖
ℎ𝐿
𝐷
Departamento de Mecânica
7
explícito, Euler implícito e Crank-Nicholson. A integração de uma grandeza genérica no
tempo pode ser representada por
∫ ( ( ) ) (17)
onde o sobrescrito o indica o valor conhecido no instante de tempo anterior.
Figura 7 - Caminho de integração no Tempo
Os esquemas de Euler implícito e explícito consideram uma variação em degrau, e são
esquemas de primeira ordem, enquanto que o método de Crank-Nicholson considera um perfil
linear, sendo de segunda ordem de precisão. A Tabela 2 ilustra o coeficiente f correspondente
a cada método.
Tabela 2 – Esquemas de Integração Temporal
Esquema: f =
Totalmente Implícito 1
Crank-Nicholson 0,5
Totalmente Explícito 0
Para ilustrar os coeficientes das equações de discretização, a integração da equação de
conservação de massa do gás, Eq. (1) é apresentada.
Integrando a Eq.(1) no volume de controle escalar (de volume xAd )ao longo do
intervalo de tempo t, de modo que:
0
e
w
GGGtt
t
tt
t
GGe
w
dtdxx
UAdxdt
tA
)()(
(18)
A integração da derivada temporal se dá primeiramente no tempo, depois no espaço;
considerando que as variáveis armazenadas no ponto nodal do volume de controle escalar
prevalecem em todo domínio de integração. O termo da derivada espacial é integrado
primeiro no espaço depois no tempo. Um esquema upwind de interpolação é utilizado para
avaliar o valor da fração volumétrica nas faces do volume de controle. A Eq. (11) é aplicada
para avaliar a integração temporal. Assim, a equação discretizada assume a forma:
Departamento de Mecânica
8
01
owGwGGeGeGG
wGwGGeGeGG
oPGGPGG
AUAUf
AUAUfxAt
,,
,,
ˆ)(ˆ)()(
ˆ)(ˆ)()()(
(19)
onde o subscrito o indica que as variáveis são referentes ao instante de tempo anterior. A
fração volumétrica de gás avaliada nas faces do volume de controle, de acordo com o
esquema upwind, é dada por:
EGeGPGeGeG UU ,,,,, ),(),(ˆ 0sinal0sinal (20)
PGwGWGwGwG UU ,,,,, ),(),(ˆ 0sinal0sinal (21)
Nas Eqs. (14) e (15), o símbolo ba, denota o máximo valor entre a e b. Para avaliar o valor
da massa específica do gás nas faces do volume de controle, um esquema upwind também foi
utilizado, fornecendo:
EGGPGGG eee UU ,, ),(),(ˆ,,, 0sinal0sinal (22)
PGGWGGG www UU ,, ),(),(ˆ,,, 0sinal0sinal (23)
Pode-se definir os “pseudo” fluxos convectivos como:
AUFAUF wGwGweGeGe ,,,,ˆ
~;ˆ
~ (24)
Com as definições acima, a Eq. (13) é reescrita da seguinte maneira:
0001
00100
00
o
PGwWGw
o
EGePGePGwWGw
EGePGe
oPGGPGG
FFf
FFfFFf
FFfxAt
,,
,,,,
,,
,~
,~
)(
,~
,~
)(,~
,~
,~
,~)()(
(25)
Assim, o sistema de equações algébricas resultante para a fração volumétrica de gás possui a
seguinte forma:
baaa WGWEGEPGP ,,, (26)
onde os coeficientes são dados pelas seguintes expressões:
Departamento de Mecânica
9
00 ,~
;,~
wWeE FfaFfa ; (27)
0101 ,~
)(;,~
)(o
woW
oe
oE FfaFfa ; (28)
)~~
)((,
ow
oe
oW
oE
oPG
oP FFfaa
t
xAa 1
(29)
)~~
(, weEWPGP FFfaa
t
xAa
(30)
oWG
oW
oEG
oE
oPG
oP aaab
,,, (31)
A Fig. 7 mostra as dependências dos vizinhos ao ponto nodal principal para cada um dos
métodos. Os métodos de Euler Explícito é mais barato, porém pode apresentar instabilidade
numérica, se o passo de tempo não for pequeno o suficiente. O método de Euler Implícito é
sempre estável, e necessita da solução do sistema algébrico para determinar as variáveis de
interesse. O método de Crack-Nicholson é de segunda ordem, e envolve mais vizinhos na
integração. Também requer solução do sistema de equações. É considerado um método
estável pelo ponto de vista matemático, implicando que eventuais oscilações numéricas,
provindas de passos de tempo não pequenos o suficiente, decairão com o tempo.
(a) – Implícito (b) – Cranck-Nicholson (c) - Explícito
Figura 8 – Dependência dos vizinhos
O código desenvolvido por [2,6] utiliza o método de Euler implícito de primeira ordem,
contudo, tendo em vista obter uma maior precisão nos resultados bem como uma diminuição
no tempo de simulação computacional, foi implementado o método de Crank-Nicholson,
correspondente à integração de segunda ordem.
A cada passo de tempo simulado, o acoplamento velocidade-pressão é resolvido através
de um algoritmo baseado no método de PRIME [6]. A fração de vazio é obtida da equação de
conservação de massa de gás, enquanto a de velocidade do gás e do líquido são obtidas a
partir da solução das respectivas equações de quantidade de movimento. A pressão é obtida
indiretamente a partir da solução da equação de conservação de massa da mistura (soma das
equações de conservação de massa de cada fase).
O sistema algébrico resultante é resolvido pelo algoritmo TDMA [5].
O modelo computacional utilizado consiste de um programa em linguagem Fortran e foi
Departamento de Mecânica
10
desenvolvido pelo núcleo de Dinâmica dos Fluidos Computacional do Departamento de
Engenharia Mecânica da PUC-Rio. Foi utilizado também um programa de pós processamento
em MATLAB que caracteriza as taxas de amplificação e as golfadas de acordo com diferentes
parâmetros já citados.
Resultados
Dois tipos de escoamento foram avaliados. No primeiro, estratificado, introduziu-se
uma perturbação controlada na entrada e investigou-se o efeito do passo de tempo na taxa de
amplificação das perturbações.
No segundo tipo, investigou-se o escoamento no padrão de golfadas. Neste caso, após a
obtenção do regime estatisticamente permanente, as grandezas características das golfadas são
determinadas.
Em ambos os casos, utilizou-se água e ar. A pressão atmosférica foi imposta na saída
para todos os casos e um par de velocidades superficiais (Usk=Qk/A, Q=vazão volumétrica,
A=área transversal da tubulação) foi prescrito na entrada, juntamente com a fração
volumétrica do gás, αg. A temperatura de referência do caso estratificado foi igual a 298.15 K
e do escoamento no padrão de golfadas foi 281,8 K. A Tabela 3 apresenta os parâmetros
geométricos do caso estratificado e de golfada.
Tabela 3 – Configurações da Tubulação
Parâmetros da tubulação
Estratificado Golfada
Diâmetro 50,8 mm
24 mm
Comprimento 5 m 10 m
O caso de escoamento estratificado foi analisado para o seguinte par de velocidades
superficiais.: UsL=0,2m/s e UsG=1 m/s. A fração de gás de equilíbrio G = 0,22174 foi imposta
na entrada.
Para determinar a fração de gás de equilíbrio, considera-se que as variações temporais e
axial do nível de líquido são nulas. As Eq. (3) e (4) se simplificam para
0sen
A
S
A
Sg
x
p iiGwG
GG
iG
G
(32)
0sen
A
S
A
Sβg
x
p iiLwLLL
iGL
(33)
Combinando as Eqs. (26) e (27), e eliminando o gradiente de pressão obtém-se a seguinte
expressão:
011
sen
GL
ii
G
GwG
L
LwLGLL
A
S
A
S
A
Sg
(34)
Departamento de Mecânica
11
Os valores do hold-up de equilíbrio podem ser obtidos através da Eq. (28), dados os valores
para as velocidades superficiais de líquido e gás, por meio de um método iterativo qualquer
(uma vez que se trata de uma equação não-linear).
Investigou-se a taxa de amplificação de perturbações introduzidas na entrada da
tubulação, a qual foi definida
)sin( thheqLL 1 (35)
onde é a amplitude da perturbação e é a frequência. A amplitude imposta foi de 1% e a
frequência igual a =2 Hz.
A Figura 1 ilustra a taxa de amplificação para o Caso 1, utilizando os esquemas de
primeira e segunda ordem. Em ambos os casos utilizou-se Δx/D =0,05
O passao de tempo está representado em função do adimentsional número de Courant,
Co= Umax Δt/Δx (36)
onde Umax é a velocidade máxima, e Δx é o espaçamento da malha.
Note que o esquema de 2ª ordem praticamente não apresenta dependência com o passo
de tempo e que só com um passo de tempo 100 vezes menor foi possível obter a mesma
solução com o método de 1ª ordem.
Figura 9 – Caso 1 – Estratificado
O caso do escoamento no regime de golfada foi obtido considerando dois casos,
ilustrado na Tabela 4. Utilizou-se uma malha com espaçamento x/D=1, número de Courant
Co=0,05. A fração de gás de equilíbrio foi prescrita na entrada, mas este valor não é relevante
para o escoamento estatisticamente permanente. Este casos foram medidos
experimentalmente por Fonseca [6].
Os resultados para as grandezas estatisticamente permanente da golfada,
correspondentes ao comprimento da bolha de Taylor Lb, velocidade de translação da golfada
Us, e frequência obtidos para o Caso 2, escoamento no padrão de golfadas, são apresentados
na Tabela 5 Na mesma tabela, os dados experimentais de Fonseca [6]. Observa-se que melhor
concordância com os dados experimentais foi obtida com o esquema de 2a. ordem
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
0 0.2 0.4 0.6 0.8 1 1.2
Taxa
de
Am
plif
icaç
ão
Courant
USL=0,2m/s ; USG=1 m/s ; Δx/D=0,05 ; Freq=2 Hz
Primeira Ordem
Segunda Ordem
Departamento de Mecânica
12
Tabela 4 – Velocidades superficiais. Caso de Golfadas
Casos de Golfada
Casos UsG (m/s) UsL (m/s)
2.1 0,788 0,295
2.2 0,788 0,393
Tabela 5 – Caso 2 – Golfadas
Caso 2.1 Caso 2.2
1a Ordem 2a ordem Experimental 1a Ordem 2a ordem Experimental
LB/D 150,38 117,335 109.69 79,38 54,15 68,01
Us 1,21 1,27 1,31 1,34 1,42 1,35
0,23 0,35 0,60 0,57 0,87 0,95
Conclusão
Os resultados obtidos mostraram um melhora significativa no emprego do esquema de
integração de segunda ordem. Mostraram ainda que é possível utilizar um passo de tempo
maior de forma a obter a solução em um tempo menor.
Como as comparações com os dados experimentais e com as correlações empíricas
forneceram resultados satisfatórios, conclui-se que o Modelo de Dois Fluidos é um boa
ferramenta para a previsão do padrão de golfadas.
Referências
1. Nieckele, A.O., Carneiro, Chucuya, R.C., , Azevedo, J.H.P., 2013, Initiation and
Statistical evolution of horizontal slug flow with a Two-Fluid Model; ASME J. Fluids
Engineering., Vol. 135, pp. 121302-1,121302-11.
2. Carneiro, J.N.E., 2006. Simulação Numérica de Escoamentos Bifásicos no Regime de
Golfadas em Tubulações Horizontais e Levemente Inclinadas. Dissertação de
Mestrado, Dept. Engenharia Mecânica, PUC-RJ.
3. Ishii, M., Hibiki, T. Thermo-fluid Dynamics of Two-Phse Flow, Springer, 2006.
4. Issa, R.I.; Kempf, M.H.W. Simulation of slug flow in horizontal and nearly horizontal
pipes with the two-fluid model. International Journal of Multiphase Flow, Volume 29,
Issue 1, 69-95, 2003.
5. Taitel, Y.; Dukler, A. E. A model for prediction of flow regime transition in horizontal
and near horizontal gas-liquid flow. Aiche Journal 22, 547-555, 1976.
6. Issa, R.I.; Bonizzi, M. On the simulation of three-phase slug flow in nearly horizontal
pipes using the multi-fluid model. International Journal of Multiphase Flow, Volume
29, Issue 11, November 2003, 1719-1747, 2003
7. Patankar, S.V. Numerical Heat Transfer and Fluid Flow. Taylon & Francis, 1980.
8. Ortega, A. J. ; Nieckele, A. O. Simulation of Horizontal Two-Phase Slug Flows Using
the Two-Fluid Model with a Conservative and Non-Conservative Formulation.
Proceedings of Congresso Brasileiro de Engenharia Mecânica. Ouro Preto, MG, 2005.
9. Fonseca Jr, R. Medição do Campo Instantâneo de Velocidade do Líquido no
Escoamento Bifásico Intermitente em Tubos Horizontais e Inclinados; Dissertação de
Mestrado, Dept. Eng. Mecânica, 2009.
