MODELAGEM DINÂMICA DO ESCOAMENTO BIFÁSICO EM … · Bol. Téc. Petrobras, Rio de Janeiro, 47 (1): 49-87, jan./mar., 2004 49 MODELAGEM DINÂMICA DO ESCOAMENTO BIFÁSICO EM RISERS

Bol. Téc. Petrobras, Rio de Janeiro, 47 (1): 49-87, jan./mar., 2004 49

MODELAGEM DINÂMICA DO ESCOAMENTO BIFÁSICO EM RISERS DE EXPLORAÇÃO DE PETRÓLEO EM ÁGUAS PROFUNDAS

DYNAMIC MODELING OF TWO PHASE FLOW IN

DEEP WATER OIL EXPLORATIONS RISERS

MODELADO DINÁMICO DEL ESCURRIMIENTO BIFÁSICO EN RISERS DE EXPLORACIÓN DE PETRÓLEO EN AGUAS PROFUNDAS

Jaime Neiva Miranda de Souza1 José Luiz de Medeiros1

André Luiz Hemerly Costa1

RESUMO Uma etapa fundamental na produção de petróleo em águas profundas é a elevação do óleo de profundidades superiores

a 1 000 metros até a superfície. O principal componente de todas as tubulações envolvidas nesse transporte é um duto vertical denominado riser, de natureza rígida ou flexível, que conecta um trecho de tubulação submarina à plataforma.

A complexidade nessa operação é bastante elevada já que compreende não-apenas o escoamento do petróleo, mas também de gás e, em algumas situações, água e areia. O conhecimento do comportamento dinâmico do escoamento em

um riser de produção de petróleo é de fundamental importância no controle das operações de uma plataforma. Esse trabalho propõe o desenvolvimento de uma ferramenta computacional para simulação dinâmica de redes de escoamento

bifásico correspondente a uma fase líquida (incompressível) e uma fase gasosa (compressível) que possibilita a determinação do comportamento temporal da pressão, da fração de ocupação das fases (hold up) e das vazões mássicas

de cada fase ao longo do comprimento dos dutos.

ABSTRACT A fundamental stage in the deep water production of oil is the elevation of the oil from depths over 1 000 meters up to

the surface. The main component of all the piping involved in this transport is a vertical duct called riser, of rigid or flexible nature, which connects a section of subsea piping to the rig. The complexity of this operation is very high, since

it involver not only the flow of oil, but also the flow of gas and, in some situations, water and sand. The knowledge of the dynamic behavior of the flow in an oil production riser is fundamentally important in the control of the operations

of a rig. Herein we propose the development of a computer tool for dynamic simulation of two phase flow networks corresponding to one liquid phase (non compressible) and one gas phase (compressible), allowing the determination of the timely behavior of the pressure, of the occupation fraction of the phases (hold up) and the mass flows of each phase

along the length of the ducts.

RESUMEN Una etapa fundamental en la producción de petróleo en aguas profundas es la elevación del petróleo de profundidades

1 Laboratório SDV – Dutos, Escola de Química, Universidade do Federal do Rio de Janeiro. e-mail:[email protected]


superiores a 1 000 metros hasta la superficie. El principal componente de todas las tuberías involucradas en ese transporte es un conducto vertical denominado riser, de naturaleza rígida o flexible, que conecta un tramo de tubería

submarina a la plataforma. La complejidad en esa operación es bastante elevada, ya que comprende no solo el escurrimiento del petróleo sino también de gas y, en algunas situaciones, agua y arena. El conocimiento del

comportamiento dinámico del derrame en un riser de producción de petróleo es de fundamental importancia en el control de las operaciones de una plataforma. Se propone, aquí, el desarrollo de una herramienta computacional para

simulación dinámica de redes de derrame bifásico correspondiente a una fase líquida (incompresible) y una fase gaseosa (compresible), que posibilita la determinación del comportamiento temporal de la presión, de la fracción de

ocupación de las fases (hold up) y de los caudales másicos de cada fase a lo largo de los conductos.

OBJETIVO O estudo aqui apresentado possui os seguintes objetivos: (i) desenvolvimento de um modelo dinâmico de rede escoamento bifásico; (ii) desenvolvimento do software MF-SIM para a simulação de redes de escoamento multifásico

estratificado ou anular; (iii) aplicação desse modelo na simulação de um sistema de exploração de petróleo em águas profundas. 1. INTRODUÇÃO A Petrobras é a líder mundial em exploração de petróleo em águas profundas e ultraprofundas. Essa hegemonia fez-se necessária devido às características das reservas nacionais que, ao final de 1999, compreendiam 17,3 bilhões de barris onde: 14% em terra firme, 11% em águas rasas, 25% em águas profundas e 50% em águas ultraprofundas. Informação essa reforçada pela descoberta, em dezembro de 2002, no Espírito Santo, do campo batizado como Jubarte, perfurado a uma profundidade de 1 500 metros, com reservas avaliadas em 300 milhões de barris de petróleo. Diante desse contexto, cada vez mais se torna necessário o desenvolvimento de tecnologias de produção de petróleo em águas profundas. Uma etapa fundamental na produção de petróleo em águas profundas é a elevação do óleo de profundidades superiores a 1 000 metros até a superfície. O principal componente de todas as tubulações envolvidas nesse transporte é um duto vertical denominado riser, de natureza rígida ou flexível, que conecta um trecho de tubulação submarina à plataforma. A complexidade nessa operação é bastante elevada já que compreende não-apenas o escoamento do petróleo, mas também de gás e, em algumas situações, água e areia. O conhecimento do comportamento dinâmico do escoamento em um riser de produção de petróleo é de fundamental importância no controle das operações de uma plataforma. Com essa motivação, propõe-se, aqui, o desenvolvimento de uma ferramenta computacional para simulação dinâmica de redes de escoamento bifásico correspondente a uma fase líquida (incompressível) e uma fase gasosa (compressível) que possibilita a determinação do comportamento temporal da pressão, da fração de ocupação das fases (hold up) e das vazões mássicas de cada fase ao longo do comprimento dos dutos. O modelo está baseado nas equações de continuidade e de momento de cada fase, compondo um sistema rígido de quatro equações diferenciais parciais ordinárias não-lineares não-homogêneas de primeira ordem. A solução desse sistema é feita via métodos Runge-Kutta, com ajuste de passo após discretização dos trechos de tubulação via Método dos Elementos Finitos, aplicando o formalismo de Galerkin no espaço com funções do tipo pulso triangular. Essa metodologia também fornece a solução estacionária para o escoamento multifásico através da solução do sistema de equações não-lineares obtidas após a discretização no espaço. Simulações dinâmicas do escoamento de gás e líquido simultaneamente em um duto envolvem a elaboração de códigos extremamente complexos que se refletem nos preços dos simuladores comerciais usualmente utilizados na indústria de gás e óleo. Como citado por Masella et al. (1) e Taitel et al. (2), os simuladores OLGA, PLAC, TACITE e TRAFLOW utilizam aproximações do tipo quasi-estacionária ou do tipo drift flux. É desconhecida a existência de simuladores comerciais que utilizem uma abordagem rigorosa semelhante à apresentada neste trabalho.


2. REVISÃO BIBLIOGRÁFICA O escoamento estacionário bifásico foi amplamente estudado nas últimas décadas. Os primeiros trabalhos correspondem a análises experimentais para a determinação de padrões de escoamento bifásico (3). Muitos outros pesquisadores restringiram-se a sistemas ar-água desenvolvendo modelos empíricos de escoamento (4 - 6). O desenvolvimento de simuladores dinâmicos bifásicos iniciou-se na indústria nuclear, porém com um enfoque totalmente diferente do buscado nas indústrias de petróleo e gás. No contexto da indústria nuclear, o maior interesse reside nos transientes de elevada velocidade, permitindo a consideração de quasi-equilíbrio nas equações de balanço de momento, simplificando em muito o modelo. Taitel et al. (7) propôs um modelo em que apenas a equação de continuidade do liquido é rigorosamente tratada, enquanto a equação de continuidade do gás e as equações de momento são assumidas como quasi-estacionárias. Essa abordagem possibilita a obtenção de boa estimativa em escoamentos, onde a fase gás escoa em altas velocidades, porém impede a verificação de acúmulo de gás dentro na linha. Seguindo o processo evolutivo, Taitel et al. (2) propuseram um novo modelo, agora considerando equações de continuidade rigorosas para ambas as fases, mas mantendo a hipótese quasi-estacionária para as equações de momento. Outra linha de pesquisa buscou a unificação das duas equações de momento e a inclusão de relações algébricas entre variáveis do modelo. Estes modelos são denominados drift-flux (1). Pires Neto et al. (8) apresentam o desenvolvimento de um modelo dinâmico de escoamento compressível, adiabático, unidimensional, de um fluido homogêneo com temperatura e viscosidade constantes. Esse trabalho utiliza metodologias propostas por Pires Neto tanto na descrição do modelo termodinâmico de escoamento quanto na resolução numérica do modelo. 3. MODELAGEM DINÂMICA DE REDES ESCOAMENTO BIFÁSICO O modelo matemático utilizado para a descrição do comportamento transiente do escoamento bifásico em um duto é composto por equações de continuidade e momento para cada fase. Devido à complexidade deste sistema físico, algumas simplificações foram utilizadas de modo a reduzir o esforço computacional, são elas: (i) Fluxo unidimensional: as variáveis de estado assumem um valor único em cada posição axial ao longo do duto; (ii) Relação termodinâmica de escoamento para a fase gás de natureza simplificada; (iii) viscosidades e diâmetro do duto constantes ao longo do escoamento; (iv) fluidos imiscíveis; (v) líquido incompressível; (vi) padrão de escoamento anular ou estratificado ao longo de todo o duto. A aplicação dos dois conjuntos de equações, considerando-se as premissas anteriores, gera um sistema de quatro equações diferenciais parciais não-lineares, não-homogêneas de ordem 1, referentes à fração de área de seção transversal ocupada pelo gás (α), a pressão (P), a vazão mássica da fase líquida (qL) e vazão mássica da fase gás (qG). 3.1. Equações da Continuidade e da Quantidade de Movimento das Fases As equações da continuidade e da quantidade de movimento são obtidas aplicando-se o princípio da conservação de massa e momento, respectivamente, em um volume de controle.


Fig. 1 - Volume de controle. Fig. 1 - Control volume. Considerando a situação descrita na figura 1, chega-se à apresentação do modelo (Apêndice 1) contendo, respectivamente, as expressões finais da continuidade da fase líquida e da fase gás e as expressões finais da quantidade de movimento da fase líquida e da fase gás:

(1)

∂

∂+

∂

∂−

−=∂

∂

xLq

LxGq

G

TP

ATP

tTP

ρρ

γ

α

γ 1

0

1

(2)

)1)((

)1(

2)1(02)1(

2

0

αθρψψ

αρα

α

αρ

−−−−

+∂

∂

−−

∂

∂−−

∂

∂

−−=

∂

∂

gsenLAILIKLLK

xLq

LALq

xTP

APxLA

LqTPAP

tLq

(3)

(4)

Definindo-se um vetor de estados [ ]GLTT qqPtxy α=),( , o modelo pode ser expresso em termos

residuais:

(5)

=

1098

654

312

1

00

00000

)(

FFFFFF

FFF

yA

,

=

11

7

00

)(

FF

yB

x x+∆x

θ qG

qL

AG

AL

γ

γγγ

αθρψψ

αρα

αργαρα

1

111

)(

2

0

02

0

2

010

2

0

TGIGIGG

G

TG

G

TG

GT

T

TG

GG

PgsenAKK

xq

PAq

xPAqPAP

xP

PAqAP

tq

−−−

+∂

∂−

∂∂

−−

∂∂

−−=

∂∂

+

xq

AtL

L ∂∂

=∂∂

ρα 1

)()( yBxy

yAty

+

∂

∂+

∂

∂=ℜ


AF

Lρ1

1 −=; αρ

γAPF

L

T=12 , αρ

γ γγ

0

1

3G

T

APF

−

=,

−

−−= 2

2

04 )1( αρL

LT A

qPAPF,

05 )1( PAF α−=

)1(2

6 αρ −=

L

L

AqF

, AgKKF LILILL )1)(sen(7 αθρψψ −++= ,

−=

γαρ12

0

2

08TG

GT

PAqPAPF

,

−= + γαγρ

α 110

2

09TG

G

PAqAPF

, γαρ

1

010 2

TG

G

PAqF =

,

AgPKKF TGIGIGG αθρψψ γ )sen(1

011 ++=

(6)

(7) Estas equações constituem generalizações do modelo de Pires Neto et al. (8) para escoamento dinâmico compressível. 3.2. Resolução Numérica do Modelo de Escoamento O modelo representado pelas equações (1-4) ou (5-7) foi resolvido no espaço através do método dos elementos finitos usando o formalismo Galerkin (Cook, 1981). Esse procedimento permite a redução das equações diferenciais parciais do modelo a um sistema de equações diferenciais ordinárias no tempo. O método é aplicado ao duto k, através da discretização em Nk elementos de tamanho kkk NLh /= delimitados pelo conjunto de nós kNj ,...,1,0= , sendo o elemento j delimitado pelos nós j-1 e j. Define-se um vetor de estados nodais empilhados:

[ ]TTN

TTk k

yyyY L10=

O aproximador para o vetor de estados ao longo do duto é definido como um truncamento em séries de Fourier a partir da base )(xjφ correspondente a uma função pulso triangular como mostrado a seguir:

∑=j

jjxtytxy ))(()(),(~ θφ , onde kk hjhxx /)()( −=θ (8)

≥<≤−<≤−+

−<

=

.1)(,0,1)(0),(1,0)(1),(1

,1)(,0

))((

xxxxx

x

xj

θθθθθ

θ

θφ (9)

onde yj(t) representa o vetor de estado para o nó j da malha de discretização. O procedimento Galerkin prevê a substituição de (8) em (5), gerando-se o vetor de resíduos em cada elemento j:

(10)

)~(

~)~(

~)(

)()(

)()( j

jj

jj yB

xy

yAty

+

∂

∂+

∂

∂=ℜ


Aplicando as propriedades das funções base φj(θ), explicitadas no Anexo 2 deste trabalho, a (10) chega-se a:

(11)

Impondo-se a ortogonalidade entre cada função pulso triangular com ℜ :

(12)

00}1{

)1( =ℜ∫ dxφ

Definindo-se as seguintes expressões para cálculo:

dxH jj

jj φ∫ℜ=

}{

)( (13)

Substituindo-se (11) em (13):

dxYBydxYAydxYAdxydxyH jj

j

jj

jjj

jj

jjj

jjjj

jjjjj φφφφφφφφφ ∫∫∫∫∫ +

+

++=

−−−−}{

)(

}{

)(

1}{

1)(

}{}{11

)()()( &&&&

dxYBydxYA

ydxYAdxydxyJ

jj

j

jj

jjj

jj

jjj

jjjj

jjjjj

φφφ

φφφφφφ

∫∫

∫∫∫

+

+

+

+

++

++

++++

+

+

+

++=

}1{

)1(

}1{

)1(

1}1{

1)1(

}1{}1{11

)()(

)(

&

&&&

(14)

Aplicando-se as restrições das funções-base contidas no Anexo 2:

(15) Dessa forma, a equação (12) fica:

0=+ jj JH para (16)

)()()( )()(11

)(11

)( jjj

jjj

jjjjj

j YByYAyYAyy ++++=ℜ −−−−φφφφ &&&&

11 −≤≤ KNj0}1{

)1(

}{

)( =ℜ+ℜ ∫∫+

+ dxdx jj

jj

j

j φφ

0}{

)( =ℜ∫ dxk

K

kN

N

N φ

kNj ..1=

dxJ jj

jj φ∫

+

+ℜ=}1{

)1(

1..0 −= kNj

( ) θφθφφ dYBhyydYAyhyhJ jj

jkjj

jjj

j

jk

jk

j ∫∫+

+

++

++

++−

++=

}1{

)1(

1}1{

1)1(

1)()(

36&&&

1..1 −= kNj


0=NkH

00 =J para

para kNj = Admitindo o formato matricial:

[ ] [ ]{ } kkkkkkk MYMY Ξ−Λ−Λ= −− 11 00& (17)

Onde:

=

IIIII

III

IIIIII

II

hM k

k

24

4

44

2

6OOO

Com dimensionamento (Nk+1)x(Nk+1) em blocos de 4x4, e I corresponde à matriz identidade 4x4. O vetor

Yk é composto pelos vetores de variáveis de estado nodais empilhados com dimensão 4(Nk+1)x1. A matriz 0 tem dimensionamento 4(Nk+1)x4 e a matriz seguinte tem dimensionamento (Nk+1)xNk em blocos 4x4:

=Λ

−−−

kk

kkkk

NN

NNNN

k

AAA

AAAA

AAA

111

3433

2322

1211

01

OO

E:

++

=Ξ

kk NN

k

B

BBBB

B

hM

2322

1211

01

com dimensão (Nk+1)x1 em blocos 4x1. A obtenção das matrizes

jjA e

1+jjA e dos vetores e se faz pela integração

numérica através do método de quadratura Gauss-Lagrange de ordem 3, 5, 10 ou 12.

0=j

jjB 1+jjB


A resolução da equação (16) fornece equações diferenciais no tempo do vetor de valores nodais para cada duto k com as seguintes condições de contorno e condições iniciais:

kk

kT

kk

kT

k

ZYttEPtE

Lx

tEPtE

x

=⇒=

==

⇒=

==

⇒=

)0(0)()(

)()(

0

4

3

2

1

α

α

(17)

Como verificado em (17), é necessário o conhecimento de α e PT nas extremidades do duto, reduzindo-se de 4 o número de incógnitas do vetor Yk(t). Define-se, então, um vetor de estado reduzido Wk(t) com tamanho 4Nkx1. O sistema de equações diferenciais ordinárias e suas condições iniciais em termos do vetor de estado reduzido adquirem a seguinte forma:

[ ] [ ]{ }{ }

)0()0(

)()()(

00 11

kk

Wk

kk

Ekk

Sk

kkkkkkkWk

YSW

tEStWStY

MYMSW

=

+=

Ξ−Λ−Λ= −−&

(18)

Onde as matrizes seletoras

kSS ,

kES e

kWS são definidas da seguinte forma:

(i) 1,

=

nCkSS , se o n-ésimo estado reduzido corresponde à c-ésima variável nodal do tubo;

(ii) 1,

=

CnkWS , se a c-ésima variável nodal do tubo é o n-ésimo estado reduzido;

(iii) 1,

=

mCkES , se a m-ésima condição de contorno corresponde a c-ésima variável nodal do tubo; e

iguais a zero, caso contrário. 3.3. Tratamento dos Vértices da Rede Os vértices, correspondentes aos trechos de junção entre dois ou mais dutos, são tratados com base em dois princípios: (i) balanços transientes de massa para cada fase; e (ii) álgebra de grafos orientados. Dessa forma, as seguintes equações de vértices são obtidas:

( ) ( )LLESLSEL QqMqMVdiagdtd

+−−= −1)(ρα

[ ] ( )[ ] ( )GGESGSET

G

LLESLSETL

T

QqMqMPdiagdiagVdiag

QqMqMdiagVdiagPdiagdtPd

+−+

++−=

+−

−

)()()(

)()()(

1*1

0

1

1γα

ργ

αργ

(19)

A formulação detalhada da equação (19) encontra-se no Apêndice 3.


Classificando-se os vértices em duas categorias: A e B. A classe A correspondendo aos vértices com que a vazão mássica ao longo do tempo é conhecida e a classe B em que a fração de gás e a pressão ao longo do tempo são conhecidas. Além disso, os vértices devem ser caracterizados por seus volumes ( AV , BV ) que são essenciais para as suas modelagens dinâmicas. Seguem-se as expressões para as classes de vértices A e B, como conseqüência do desdobramento de (19):

onde Ai e Bi são, respectivamente, vetores de índices de vértices da classe A e B; Eq e Sq representam

respectivamente as vazões de entrada e saída de tubos; e as funções ),,(ATAAL PtQ α , ),,(

ATAAG PtQ α , e

)(tPBT são as especificações da rede.

3.4. Sistema em Rede As redes de escoamento são constituídas de um conjunto de elementos como tubulações, vértices, válvulas e bombas interconectados, com o objetivo de transportar fluidos entre pontos de fornecimento e pontos de consumo. No contexto desse trabalho serão considerados apenas dutos e vértices de Classe A e B, representados graficamente por uma matriz de incidência com orientações atribuídas aos dutos de modo que valores negativos de vazão indicam um escoamento no sentido contrário à orientação adotada. A estrutura na forma de redes implica na utilização de especificações de vértices para representar todo o comportamento da rede, no caso desse trabalho, corresponde a duas especificações por vértices, de acordo com a classe de cada vértice. A modelagem de uma rede de escoamento multifásico corresponde a quatro equações diferenciais parciais por trecho de tubulação (equações 1-4) e duas equações diferenciais ordinárias por vértice do tipo A (equação 20). Através do Método dos Elementos Finitos e usando o formalismo de Galerkin, descrito em seções anteriores, é possível transformar esse sistema de equações diferenciais parciais em um sistema de quatro equações diferenciais ordinárias por elemento finito por duto e duas equações diferenciais ordinárias por vértice da Classe A. Compõe-se, então, o grande vetor de estados reduzidos:

[ ]TT

ATT

AT

mT PWWX αL1=

Desta forma, o grande sistema de equações diferenciais e seus respectivos contornos é:

)(tBα

( ) ( )[ ] ( )

[ ] ( )

[ ] (20) )()(

)()()(:),(:),(),,(

)(:),(:),(),(

)0()0(0

),,(:),(:),()()()(

),,(:),(:),()()()(

),,(:),(:),()(

1

1

1

*0

11*0

1*10

1

1

dtdPdiagVdiag

dtPd

PdiagdiagVdiagqiMqiMPtQ

dtdVdiagqiMqiMtQ

PPt

PtQqiMqiMPdiagdiagVdiag

PtQqiMqiMdiagVdiagPdiagdt

Pd

PtQqiMqiMVdiagdt

d

BATBG

BT

ATBBGLSBELEBSBTBBG

BBLLSBELEBSBBL

ATATAA

ATAAGGEASGSAEATAAG

ATAALLEASLSAEAALATAT

ATAALLEASLSAEALA

αγρ

γαρα

αρα

αα

αγαρ

ααργ

αρα

γ

γ

γ

+

++−=

−−=

===

+−+

++−=

+−−=

−−

−−

−

−


[ ]TTP

TTm

TX ΠΠΠΠ= αL&1 (21)

[ ]TTT

TA

Tm

T PWWX )0()0()0()0()0( 1 αL=

A integração no tempo dessas equações diferenciais é feita através de métodos Runge-Kutta adaptativos de ordem 15/23 (solvers ODE-15s e ODE-23s do MATLAB 6.1). O sistema de equações diferenciais ordinárias descritos em (21) produz como subproduto um modelo estacionário para o escoamento bifásico em uma rede de dutos:

[ ] 01 =ΠΠΠΠTT

PTT

mT

αL (22) Esse sistema de equações algébricas foi resolvido utilizando-se o método de Newton-Raphson, em que a estimativa inicial foi gerada através da simulação estacionária do escoamento incompressível de cada uma das fases isoladamente. O modelo de escoamento estacionário incompressível está descrito no Anexo 4 deste trabalho. 4. SIMULAÇÕES E RESULTADOS Foi modelada uma rede de dutos correspondente a um sistema de produção de petróleo em águas profundas. Essa rede é composta por um grande duto (Tubo 1) de 20 cm de diâmetro e 50 km de comprimento, a uma profundidade de 1 000 metros. Uma das extremidades do duto (Vértice 3) corresponde ao poço de petróleo modelado como um vértice de Classe B, ou seja, em que se conhece a pressão da linha e a fração ocupação de gás. A outra extremidade (Vértice 1) corresponde a um vértice de Classe A, ou seja, são conhecidas as vazões mássicas das fases que nesse caso são nulas. O que implica que esse vértice tem a única função de conectar o grande duto submarino ao riser de produção. O riser (Tubo 1) também possui um diâmetro de 20 cm e possui sua extremidade superior (Vértice 2) fixada à plataforma e caracterizada por um vértice de Classe B. Na figura 2 mostra-se a tela principal do software desenvolvido neste trabalho e demonstra a estrutura desse sistema de transporte.

Fig. 2 - Estrutura da rede analisada. Fig. 2 - Analyzed network structure. Os fluidos em escoamento anular correspondem a: (i) gás de viscosidade 0,02 cP, densidade de referência 280 Kg/m3 a uma pressão de referência de 100 bar e expoente politrópico 1,2; e (ii) óleo de viscosidade 8 cP e densidade 800 Kg/m3. A rugosidade dos dutos foi suposta 250 µm. O comportamento dinâmico do escoamento é verificado pela alteração da pressão do poço ao longo do tempo que se inicia a 150 bar, atingindo um valor máximo de 170 bar e retornando ao valor inicial. Além disso, fez-se uma alteração da fração de gás na plataforma ao longo do tempo partindo de 0,25, atingindo 0,28 e, em seguida, retornando a


0,26. Foram utilizados 20 elementos finitos em cada trecho de tubulação, o que corresponde a um sistema de 162 equações diferenciais ordinárias. É importante ressaltar que a metodologia desenvolvida é capaz de representar redes com múltiplas conexões (e.g. multi-well trunkline). Nas figuras 3, 4, 5 (a e b) e 6 descreve-se o comportamento no tempo e no espaço da fração de ocupação, da pressão, da vazão mássica de gás e da vazão mássica de líquido, respectivamente. Nas figuras 7, 8 e 9 descreve-se o comportamento dos vértices. As demais variáveis descrevem o comportamento do escoamento como tensões de cisalhamento, número de Reynolds e velocidade das fases reunidas no Apêndice 5 desse trabalho.

Fig. 3 - Fração de ocupação ao longo dos dutos. Fig. 3 - Occupation fraction along the ducts.


Fig. 4 - Pressão ao longo dos dutos. Fig. 4 - Pressure along the ducts.

Fig. 5a - Vazão mássica da fase líquida ao longo dos dutos. Fig. 5a - Mass flow of the liquid phase along the ducts.


Fig. 5b - Vazão mássica da fase líquida ao longo dos dutos. Fig. 5b - Mass flow of the liquid phase along the ducts.

Fig. 6 - Vazão mássica da fase gasosa ao longo dos dutos. Fig. 6 - Mass flow of the gas phase along the ducts.


Fig. 7 - Comportamentos do vértice 1. Fig. 7 - Behavior of vertex 1.

Fig. 8 - Comportamentos do vértice 2. Fig. 8 - Behavior of vertex 2.


Fig. 9 - Comportamentos do vértice 3. Fig. 9 - Behavior of vertex 3. Nas figuras 10 11 e 12 ilustram-se telas que compõem a interface gráfica do software desenvolvido ao longo desse trabalho.

Fig. 10 - Gráficos de acompanhamento da simulação dinâmica. Fig. 10 - Charts for follow up of the dynamic simulation.


Fig. 11 - Interface para plotagem dos resultados. Fig. 11 - Interface for plotting the results.


Fig. 12 - Telas de configuração do simulador. Fig. 12 - Configuration screens of the simulator. REFERÊNCIAS BIBLIOGRÁFICAS (1) MARSELLA, J. M.; TRAN, Q. H.; FERRE, D.; PAUCHON, C. Transient Simulation of Two-Phase

Flows in Pipes, Int. J. Multiphase Flow, 24, 739-755, 1998. (2) TAITEL, Y.; BARNEA, D. Simplified Transient Simulation of Two-Phase Flow using Quasi-

Equilibrium Momentum Balances, Int. J. Multiphase Flow, 23, 3, 493-501,1997. (3) GOLAN, L. P.; STENNING, A. H. Two-Phase Vertical Flow Map, Proc. Instn. Mech. Engrs., 184, 3C,

108-114, 1969. (4) TAITEL, Y.; Dukler, A. E. A Model for Predicting Flow Regime Transitions in Horizontal and Near

Horizontal Gas-Liquid Flow, AIChE Journal, 22, 1, 47-54, 1976. (5) TAITEL, Y.; BORNEA, D.; DUKLER, A. E. Modelling Flow Pattern Transitions for Vertical Steady

Upward Gas-Liquid Flow in Vertical Tubes, AIChe Journal, 26, 3, 345-354, 1980. (6) BARNEA, D.; SHOHAM, O.; TAITEL, Y.; DUKLER, A. E. Gas Liquid in Inclinjed Tubes: Flow

Pattern Transitions for Upward Flow, Chem. Eng. Sci., 40, 1, 131-136, 1985. (7) TAITEL, Y.; SHOHAM, O.; BRILL, J. P. Simplified Transient Solution and Simulation of Two Phase

Flow in Pipes, Chem. Eng. Sci., 44, 1353-1359, 1989. (8) PIRES NETO, J. P. Modelagem Dinâmica em Redes de Escoamento Compressível para Aplicações à

Detecção de Vazamentos em Tempo Real. Rio de Janeiro: UFRJ/EQ, 2001. (9) WONGWISES, S.; KONGKIATWANITCH, W. Interfacial Friction Factor in Vertical Upward Gas-

Liquid Annular Two-Phase Flow, Int. Comm. Heat Mass Transfer, 28, 3, 323-336, 2001.


APÊNDICE 1 MODELO DINÂMICO DE ESCOAMENTO MULTIFÁSICO EM DUTOS Seja um elemento de volume de fluido V∆ , mostrado na figura , composto por uma fase líquida (L) e uma fase gás (G). Esse elemento de volume pode ser decomposto em:

GL VVV ∆+∆=∆ (A1-1) Que, em termos de área da seção transversal, corresponde a:

xAxAxA GL ∆+∆=∆ (A1-2) As áreas preenchidas pelas diferentes fases relacionam-se através da fração ocupada por gás (α) pelas seguintes expressões:

AAG α= (A1-3)

( )AAL α−= 1 (A1-4) Fazendo-se balanços de massa nesse volume de controle para cada fase em escoamento, considerando-se a inexistência de escoamentos secundários (ao longo da seção transversal do duto), ou seja, considerando-se que as variáveis do modelo sejam função apenas do tempo e da posição ao longo do comprimento do duto, tem-se:

(A1-5) Substituindo-se (A1-2) em (A1-5):

( )

( ))()(

)()(

xxqxqtA

x

xxqxqtAx

GGGG

LLLL

∆+−=∂

∂∆

∆+−=∂

∂∆

ρ

ρ

(A1-6)

Fazendo-se 0→∆x na expressão (A1-6), tem-se:

( )

( )x

qtA

xq

tA

GGG

LLL

∂∂

−=∂

∂∂

∂−=

∂∂

ρ

ρ

(A1-7)

Substituindo-se (A1-3) e (A1-4) em (A1-7):

(A1-8)

( )

( ) )()(

)()(

xxqxqt

V

xxqxqt

V

GGGG

LLLL

∆+−=∂∆∂

∆+−=∂∆∂

ρ

ρ

( )

( )x

qAt

xq

AtGG

LL

∂∂

−=∂

∂∂

∂−=

∂−∂

1

1)1(

αρ

αρ


Utilizando-se a hipótese de incompressibilidade do líquido:

(A1-9) Fazendo-se balanços de quantidade de movimento no volume de controle para cada fase em escoamento, considerando, novamente, a inexistência de escoamentos secundários ao longo da seção transversal do duto, tem-se:

( )

( ) ( )xxAgxKxK

xxvxxxxAxvxxA

xxPxxAxPxAt

vV

LLILILL

LLLLLL

LLLLL

∆−∆−∆−+∆+∆+∆+−+

+∆+∆+−=∂∆∂

)()sen()()()()()()(

)()()()(

22

θρψψρρ

ρ

(A1-10)

( )

( ) ( )xxAgxKxK

xxvxxxxAxvxxA

xxPxxAxPxAt

vV

GGIGIGG

GGGGGG

GGGGG

∆−∆−∆−+∆+∆+∆+−+

+∆+∆+−=∂∆∂

)()sen()()()()()()(

)()()()(

22

θρψψρρ

ρ

(A1-11)

Sendo KL, KG e KI os perímetros correspondentes à fase líquida, à fase gás e à interface, respectivamente. Como esses perímetros dependem do padrão de escoamento, já que: (i) no caso de escoamento anular o líquido está em contado com toda a parede e com o gás, enquanto este está em contato apenas com o líquido; e (ii) no caso de escoamento estratificado, tanto o líquido quanto o gás estão em contato entre si e com a parede. Logo, a utilização deste modelo requer a adoção de regime de escoamento estratificado ou anular. Como a formulação proposta é aplicável a ambos regimes, esta seguirá utilizando os perímetros de contato. A distinção será feita ao final de todo o desenvolvimento. Sendo iiii qxvV ∆=∆ρ para { }VL,∈i e fazendo-se 0→∆x :

( )LLILILLL

LLLLL AgKKx

vAxPA

tq )sen(

2

θρψψρ−−−

∂∂

−∂

∂−=

∂∂

(A1-12)

( )

GGIGIGGGGGGGG AgKK

xvA

xPA

tq )sen(

2

θρψψρ−−−

∂∂

−∂

∂−=

∂∂

(A1-13)

Sendo:

ii

iiii A

qvAρ

ρ2

2 = (A1-14)

Substituindo-se (A1-14) em (A1-12) e (A1-13), tem-se:

)1)(sen()1(

1)1( 2

αθρψψρα

α−−−−

−∂

∂−

∂−∂

−=∂

∂ gAKKqxAx

PAt

qLILILLL

L

LL

(A1-15)

xq

Att

xq

AtG

GG

LL

∂∂

−=∂∂

+∂

∂∂

∂=

∂∂

1

1

αρρα

αρ


αθρψψαρ

α )sen(1 2

gAKKqxAx

PAt

qGIGIGGG

G

GG −−−

∂∂

−∂

∂−=

∂∂

(A1-16)

Utilizando-se a hipótese de incompressibilidade do líquido em (A1-15):

)1)(sen(1

1)1( 2

αθρψψαρ

α−−−−

−∂

∂−

∂−∂

−=∂

∂ gAKKqxAx

PAt

qLILILLL

L

L

L

(A1-17)

Derivando-se os termos da equação (A1-17):

)1)(sen()1(

2)1(

)1( 2

2

αθρψψαρ

ααρ

α

−−−−

+∂∂

−−

∂∂

−

−+∂∂

−−=∂

∂

gAKKxq

Aq

xAqAP

xPA

tq

LILILLL

L

L

L

L

LL

(A1-19)

Derivando-se os termos da equação (A1-16), tem-se:

αθρψψ

ραραρ

αρα

α

)sen(

2 2

2

2

2

gAKKxA

qx

qA

qxA

qAPxPA

tq

GIGIGGG

G

G

GG

G

G

G

GG

−−−

+∂

∂+

∂∂

−∂∂

+−+

∂∂

−=∂

∂

(A1-20)

A partir deste momento, torna-se necessária a utilização de uma expressão para relacionar a densidade do gás às demais variáveis de estado do modelo, já que a utilização de equações de estado traria grande aumento na complexidade do sistema de equações diferenciais pela inclusão de balanços de energia para ambas as fases. Com o objetivo de evitar o uso de equações de balanço de energia, foi adotada uma relação termodinâmica de escoamento que se traduz em uma forma semi empírica do tipo CPVG =γ , que estabelece uma relação entre a pressão e a densidade do gás:

CPMPVG

== γ

γγ

ρ (A1-21)

Diferenciando-se com respeito a uma variável qualquer η, tem-se:

ηργρ

η γ

γ

∂∂

=

∂∂ −

GG

MCP 1

(A1-22)

Substituindo-se (A1-21) em (A1-22), tem-se:

ηρ

ργ

η ∂∂

=

∂∂ G

G

PP

(A1-23)

Obtendo-se como subproduto as seguintes relações:


=

∂∂

GG

PPργ

ρ (A1-24)

xP

PxGG

∂∂

=

∂∂

γρρ

(A1-25)

tP

PtGG

∂∂

=

∂∂

γρρ

(A1-26)

Com o objetivo de simplificar as equações e escalonar as variáveis de estado, será utilizado um estado de referência P0 e ρG0, e variáveis adimensionais PT e ρGT obtidas através das relações a seguir:

TPPP 0= (A1-28)

GTGG ρρρ 0= (A1-29)

γρGTTP = (A1-30)

Substituindo-se (A1-26), (A1-28), (A1-29) e (A1-30) em (A1-9B), tem-se:

xq

AP

tP

tP G

G

TTT

∂∂

−=∂∂

+∂

∂−

αργα

αγ γ

γ

0

1

(A1-31)

Substituindo-se (A1-26), (A1-28), (A1-29) e (A1-30) em (A1-19), tem-se:

)1)(sen()1(

2)1(

)1( 2

2

00

αθρψψαρ

ααρ

α

−−−−

+∂∂

−−

∂∂

−

−+∂∂

−−=∂

∂

gAKKxq

Aq

xAqPAP

xPPA

tq

LILILLL

L

L

L

L

LT

TL

(A1-32)

Substituindo-se (A1-26), (A1-28), (A1-29) e (A1-30) em (A1-20), tem-se:

αθρψψ

αρα

ραα

γαργ

γγγ

)sen(

2

1

111

0

002

2

0010

2

gPAKK

xq

PAq

xPAqPAP

xPAP

PAq

tq

TGIGIGGG

G

TG

G

TG

GT

T

TG

GG

−−−

+∂

∂−

∂∂

+−+

∂∂

−=

∂∂

+

(A1-33) Após as modificações descritas acima, obtém-se o seguinte conjunto de equações diferenciais:

xq

AtL

L ∂∂

=∂∂

ρα 1

(A1-9A)

xq

AP

tP

tP G

G

TTT

∂∂

−=∂∂

+∂

∂−

αργα

αγ γ

γ

0

1

(A1-31)


)1)(sen()1(

2)1(

)1( 2

2

00

αθρψψαρ

ααρ

α

−−−−

+∂∂

−−

∂∂

−

−+∂∂

−−=∂

∂

gAKKxq

Aq

xAqPAP

xPPA

tq

LILILLL

L

L

L

L

LT

TL

(A1-32)

αθρψψ

αρα

ραα

γαργ

γγγ

)sen(

2

1

111

0

002

2

0010

2

gPAKK

xq

PAq

xPAqPAP

xPAP

PAq

tq

TGIGIGGG

G

TG

G

TG

GT

T

TG

GG

−−−

+∂

∂−

∂∂

+−+

∂∂

−=

∂∂

+

(A1-33)

Representando-se o modelo de forma simplificada:

xqF

tL

∂∂

−=∂∂

1α

(A1-34)

xqPF

tPF

tP G

TTT

∂∂

−=∂∂

+∂

∂ ),(),( 32 ααα (A1-35)

),,(),()(),,( 7654 LTL

LT

LTL qPF

xqqF

xPF

xqPF

tq ααααα −

∂∂

−∂

∂−

∂∂

−=∂

∂

(A1-36)

),,(),,(),,(),,( 111098 GTG

GTT

GTGTG qPF

xqqPF

xPqPF

xqPF

tq

ααααα −∂

∂−

∂∂

−∂∂

−=∂

∂

(A1-37) Onde:

AF

Lρ1

1 −= (A1-38)

αγα T

TPPF =),(2

(A1-39)

αργα

γγ

0

1

3 ),(G

TT A

PPF−

= (A1-40)

( )

−−−= 2

2

04 1),,(

αρα

L

LTLT A

qPAPqPF (A1-41)

05 )1()( PAF αα −= (A1-42)


)1(2),(6 αρ

α−

=L

LL A

qqF (A1-43)

AgKKqPF LILILLLLT )1)(sen(),,(7 αθρψψα −++= (A1-44)

−=

γραα 1

02

2

08 ),,(TG

GTGT

PAqPAPqPF

(A1-45)

−= + γγαρ

αα 110

2

09 ),,(TG

GGT

PAqAPqPF

(A1-46)

γαρα 1

010 2),,(

TG

GGT

PAqqPF =

(A1-47)

AgPKKqPF TGIGIGGGGT αθρψψα γ )sen(),,(

1

011 ++= (A1-48)

Substituindo-se (A1-34) em (A1-35), obtém-se:

xqPF

xqPFF

tP G

TL

TT

∂∂

−=∂

∂−

∂∂ ),(),( 321 αα

(A1-49)

Que corresponde a:

xqPF

xqPF

tP G

TL

TT

∂∂

−∂

∂−=

∂∂ ),(),( 312 αα

(A1-50)

),(),( 2112 TT PFFPF αα −= (A1-51)

O modelo modificado passa a ser descrito pelas seguintes equações:

xqF

tL

∂∂

−=∂∂

1α

(A1-34)

xqPF

xqPF

tP G

TL

TT

∂∂

−∂

∂−=

∂∂ ),(),( 312 αα

(A1-52)

),,(),()(),,( 7654 LTL

LT

LTL qPF

xqqF

xPF

xqPF

tq ααααα −

∂∂

−∂

∂−

∂∂

−=∂

∂

(A1-36)

),,(),,(),,(),,( 111098 GTG

GTT

GTGTG qPF

xqqPF

xPqPF

xqPF

tq

ααααα −∂

∂−

∂∂

−∂∂

−=∂

∂

(A1-37)

Escrevendo o modelo segundo um vetor de resíduos


+∂

∂+

∂∂

+∂∂

+∂

∂

+∂

∂+

∂∂

+∂∂

+∂

∂∂

∂+

∂∂

+∂

∂∂

∂+

∂∂

=ℜ

),,(),,(),,(),,(

),,(),()(),,(

),(),(

111098

7654

312

1

GTG

GTT

GTGTG

LTL

LT

LTL

GT

LT

T

L

qPFx

qqPFx

PqPFx

qPFt

q

qPFx

qqFx

PFx

qPFt

qx

qPFx

qPFt

Px

qFt

ααααα

ααααα

αα

α

(A1-53) Define-se o vetor de estado y contendo a fração ocupada pelo gás α , pressão adimensional PT , a vazão mássica da fase líquida qL , a vazão mássica da fase líquida qG .

==

G

L

T

qqP

ytxy

α

),(

(A1-54)

Representando-se o modelo na forma matricial, tem-se:

+

∂

∂

+

∂

∂=ℜ

),,,(),,,(

00

),,(0),,(),,(0),()(),,(

),(),(00000

11

7

1098

654

312

1

GLT

GLT

GTGTGT

LLT

TT

qqPFqqPFx

y

qPFqPFqPFqFFqPF

PFPFF

ty

αα

αααααα

αα

(A1-55) Introduzindo-se as seguintes definições:

=

),,(0),,(),,(0),()(),,(

),(),(00000

)(

1098

654

312

1

GTGTGT

LLT

TT

qPFqPFqPFqFFqPF

PFPFF

yA

αααααα

αα

(A1-56)

=

),,,(),,,(

00

)(

11

7

GLT

GLT

qqPFqqPF

yB

αα

(A1-57)

)()( yBxy

yAty

+

∂

∂+

∂

∂=ℜ

(A1-58)


Modelagem das tensões de cisalhamento As tensões de cisalhamento da fase gás, da fase líquida e da interface são definidas por:

8LLL

LL

vvf

ρψ =

(A1-59)

8GGG

GG

vvf

ρψ =

(A1-60)

( )GGAGIIG φφφφφψψ ↑↓↑↑ +⋅= (A1-61)

( )LLAGIIL φφφφφψψ ↑↓↑↑ +−⋅=

(A1-62) Onde:

( ) ( )↑↓↑↑

++

−= φ

ρφ

ρψ

22

22LGG

ILGG

II

vvf

vvf

(A1-63)

L

LL A

qvρα )1( −

= (A1-64)

γρα1

0 TG

GG

PAqv =

(A1-65)

A intensidade da tensão cisalhante da interface depende dos sentidos dos escoamentos de cada fase, por isso

assume uma forma que depende de funções seletoras φ, como mostrado a seguir:

( ) 1)exp(1 −

↑↑ −+= LGqqζφ (A1-66)

( ) 1)exp(1 −

↑↓ += LGqqζφ (A1-67)

( ) 1)exp(12 1 −−+= −

GG qζφ (A1-68)

( ) 1)exp(12 1 −−+= −

LL qζφ (A1-69)

( )[ ]{ } 1exp12 1 −−−+= −

LGAG vvζφ (A1-70)


Onde ζ é um fator de escala que deve ser muito maior que as vazões mássicas e que as velocidades das fases. Com essa formatação, os sinais das tensões cisalhantes de interface das fases estão mostrados na tabela a seguir:

TABELA A1 VALORES DAS FUNÇÕES SELETORAS E SINAIS DAS

TENSÕES DE INTERFACE TABLE A1

VALUES OF THE SELECTION FUNCTIONS AND SIGNALS OF THE INTERFACE STRESS

qG>0 qG<0 qL>0 qL<0 vG>vL vG<vL

qL<0 qL>0 vG>vL vG<vL

φ↑↑ 1 1 0 0 1 1 φ↑↓ 0 0 1 1 0 0 φAG 1 -1 X X 1 -1 φL 1 1 -1 1 -1 -1 φG 1 1 1 -1 -1 -1 ΨIL (-) (+) (+) (+) (-) (+) ΨIG (+) (-) (-) (-) (-) (+)

Para o cálculo dos fatores de atrito das fases líquido e gás utiliza-se o modelo proposto por Churchill (1977), que corresponde a:

( )

21

23

1Re88

12

++

=

iiii

BAf , para { }L G,∈i (A1-80)

16

9.0

27.0Re7

1ln457.2

+

⋅=

D

A

i

iε

(A1-81)

16

Re37530

=

iiB (A1-82)

Sendo:

Dq

L

LL πµα 5.0)1(

4Re−

= (A1-83)

Dq

G

GG πµα 5.0

4Re = (A1-84)


Para o cálculo dos fatores de atrito na interface utiliza-se o modelo proposto por Wongwises e Kongkiatwanitch (9) que corresponde a:

crb

SGrI Daf

= }{

}{ Reδ

(A1-85)

para: a = 17.172, b = – 0.768 e c = – 0.253. Onde:

Dq

G

GSG πµ

α 5.04Re = (A1-86)

Perímetros de Contato para Regime Anular e Regime Estratificado A determinação dos perímetros de contato e da espessura do filme líquido é feita de acordo com o regime de escoamento, diferenciando-se caso seja anular ou estratificado: A) Para regime estratificado: Fig. A1.1 - Representação da seção transversal de um duto em escoamento estratificado. Fig. A1.1 – Representation of the cross section of a duct under stratified flow.

−=

2cos1

2}{βδ D

est (A1-87)

{ } 2DK estL

β= (A1-88)

{ } ( )DK estG βπ −= 2 (A1-89)

{ }

=

2sen βDK estI (A1-90)

Onde:

)sen()1(2 ββαπ −=− (A1-91) B) Para regime anular:

AG

AL δ{est}

D

β

A

AL

D


Fig. A1.2 - Representação da seção transversal de um duto em escoamento anular. Fig. A1.2 – Representation of the cross section of a duct under annular flow.

( )αδ −= 12}{D

an (A1-92)

{ } DK anL π= (A1-93)

{ } 0=anGK (A1-94)

{ } απDK anI = (A1-95)


APÊNDICE 2 PROPRIEDADES DAS FUNÇÕES BASE )(θφ j Usando-se a notação {j} para denotar elemento j ⇒ [ ]0,1−∈θ , como mostrado na figura A2.1:

Fig. A2.1 – Funções bases pulso triangulares Fig. A2.1 – Triangular pulse base functions. • Propriedades de )(θφ j :

{ } 2)( hdx

ii =∫ θφ (A2-1)

{ } 21)( =∫

i

ii dx

dxdφ

θφ (A2-2)

{ } 3)()( hdx

iii =∫ θφθφ (A2-3)

{ }1ou para ,0)( +><=∫ jijidx

ij θφ (A2-4)

{ },...,ljijidx

i

lj 21 e 1ou para ,0)( =+><=∫ θφ (A2-5)

{ },...,ljijidx

dxd

i

jlj 21 e 1ou para ,0)( =+><=∫

φθφ (A2-6)

{ },...,,...,ljijidxm

ilj 21m e 21 , 1ou para ,0)()( ==+><=∫

∀

θφθφ (A2-7)

j

x

{j} {j+1}j j+1 j-1

... ...

0 1 NkNk-1 {0} {Nk}

0 hk Lk jhk (j+1)hk (j-1)hk Lk-hk

)(xjφ

0

1


{ } 2)(

1

hdxj

j =∫+

θφ (A2-8)

{ } 21)(

1

−=∫+j

jj dx

dxdφ

θφ (A2-9)

{ } 3)()(

1

hdxj

jj =∫+

θφθφ (A2-10)

{ } 6)()( 1

hdxj

jj =∫ − θφθφ (A2-11)

{ }0)()(

11 =∫

+−

jjj dxθφθφ (A2-12)

{ }0)()( 1 =∫ +

jjj dxθφθφ (A2-13)

{ } 6)()(

11

hdxj

jj =∫+

+ θφθφ (A2-14)

{ } 21)( 1

1 −=∫ −−

j

jj dx

dxdφ

θφ (A2-15)

{ } 21)(1 =∫ −

j

jj dx

dxdφ

θφ (A2-16)

{ } 2)(1

hdxj

j =∫ − θφ (A2-17)

{ } 21)(

1

1 =∫+

+

j

jj dx

dxdφ

θφ (A2-18)

{ } 21)(

11 −=∫

++

j

jj dx

dxdφ

θφ (A2-19)

hdxL

j =∫0

)(θφ (A2-20)

2)(

00

hdxL

=∫ θφ (A2-21)

2)(

0

hdxL

Nk=∫ θφ (A2-22)


2j-i onde ,0)()(0

≥=∫L

ji dxθφθφ (A2-23)

k0

1 Njou 0 onde ,6

)()( <≥=∫ + jhdxL

jj θφθφ (A2-24)

6)()(

01

hdxL

jj =∫ − θφθφ (A2-25)

21)(

)(0

1 −=∫ −L

jj dx

dxd θφ

θφ (A2-26)

21)(

)(0

1 =∫ −

Lj

j dxdx

d θφθφ (A2-27)

32)()(

0

hdxL

jj =∫ θφθφ (A2-28)

3)(

0

20

hdxL

=∫ θφ (A2-29)

• Comportamento de funções-base por elementos:

Em {j} : θφθφ

−=−=

−1

1

j

j (A2-30)

Em {j+1} : θφθφ

=−=

−1

1

j

j (A2-31)


LGDuto 1 Duto 3

Duto 2

Duto 4

Vértice 1 QG

QL

APÊNDICE 3 MODELO DINÂMICO DE VÉRTICES BIFÁSICOS O modelo dinâmico dos vértices das redes corresponde a balanços de massa em cada uma das fases em escoamento. No caso do escoamento bifásico de fluidos imiscíveis, o modelo será constituído de um sistema de duas equações diferencias ordinárias. Na figura A3.1 ilustra-se a representação física de um vértice bifásico: Fig. A.3.1 - Representação de um vértice bifásico. Fig. A.3.1 - Representation of a two phase vertex. As seguintes hipóteses foram utilizadas para a elaboração desse modelo de vértice: (i) a fração de ocupação do gás (α) nos tubos conectados ao vértice são iguais à fração volumétrica de gás existente dentro desse vértice; e (ii) a pressão é a mesma em todo o vértice, o que implica que a pressão nas extremidades dos dutos conectados ao vértice é igual à pressão do vértice. As equações de balanço de massa para cada fase no vértice ilustrado na figura A.3.1 ficam:

GGGGGG

LLLLLL

Qqqqqdt

dM

Qqqqqdt

dM

+−−+=

+−−+=

4321

4321

(A3-1)

Que em termos de volume corresponde a:

GGGGGGG

LLLLLLL

Qqqqqdt

Vd

Qqqqqdt

Vd

+−−+=

+−−+=

4321

4321

ρ

ρ

(A3-2)

Usando a relação termodinâmica apresentada no Anexo 1, tem-se:

GGGGGGTG

LLLLLLL

Qqqqqdt

VPd

Qqqqqdt

Vd

+−−+=

+−−+=

4321

/10

4321

γρ

ρ

(A3-3)

Que em termo de fração de ocupação de gás, fica:


GGGGGTG

LLLLLL

Qqqqqdt

VPd

Qqqqqdt

Vd

+−−+=

+−−+=−

4321

/10

4321)1(

αρ

αρ

γ (A3-4)

Derivando-se os termos de (A3-4), tem-se:

( )

GGGGGTGTTG

LLLLLL

QqqqqdtdVP

dtdPPV

QqqqqVdt

d

+−−+=+

+−−+−=

+−

4321/1

0

/110

43211

αργ

αρ

ρα

γγ

(A3-5)

Que, finalmente, corresponde a:

( )

( ) ( )LLLLLL

TGGGGG

G

TT

LLLLLL

QqqqqV

PQqqqqV

Pdt

dP

QqqqqVdt

d

+−−+++−−+=

+−−+−=

−

432143210

/11

43211

αργ

αργ

ρα

γ (A3-6)

Devido à complexidade das redes de escoamento, uma notação matricial deve ser utilizada para representar matematicamente a sua conectividade. Para isso, são necessárias matrizes de incidência de entrada e de saída. A matriz de incidência de entrada (em relação ao vértice) EM , de tamanho NV x m, possui a

seguinte lei de formação: [ ] 1

,=

knEM , caso o k-ésimo duto esteja chegando ao vértice n-ésimo; e igual a zero caso contrário.

A matriz de incidência de saída (em relação ao vértice) SM , de tamanho NV x m, possui a seguinte lei de

formação: [ ] 1

,=

knSM , caso o k-ésimo duto esteja saindo ao vértice n-ésimo; e igual a zero caso contrário.

O modelo de vértice no formato matricial fica:

[ ] ( )( ) ( ) ( )( )

( ) ( ) ( )[ ] ( )LELSSLELT

GEGSSGEGTT

LELSSLEL

QqMqMdiagVdiagPdiag

QqMqMdiagVdiagPdiagdtPd

QqMqMVdiagdtd

+−+

++−=

+−−=

−

−

−

1

0/11

1

)(

αργ

αργ

ρα

γ (A3-7)

onde qSL e qSG são vetores de vazões mássicas de saída em relação ao duto de, respectivamente, líquido e gás; e qEL e qEG são vetores de vazões mássicas de entrada em relação ao duto de, respectivamente, líquido e gás.


APÊNDICE 4 MODELO ESTACIONÁRIO DE ESCOAMENTO INCOMPRESSÍVEL O modelo estacionário de escoamento incompressível em redes de dutos tem grande importância neste trabalho, já que possibilita a obtenção de uma estimativa inicial na determinação do estado estacionário no escoamento multifásico. Este modelo baseia-se na resolução de um sistema de equações algébricas não-lineares através de um algoritmo Newton-Raphson com matriz Jacobiana analítica. O modelo matemático adotado para a simulação estacionária envolve equações de balanço material nos vértices da rede e equações de conservação da energia ao longo dos elementos (Mah, 1990). Seja uma rede composta por N vértices ( ),...,1 Nt = interligados através de S tubulações ( Sk ,...,1= ). As variáveis do sistema são as vazões nas tubulações, representadas pelas variáveis qk (para cada duto é atribuída uma certa direção, se o fluido escoa ao longo desta direção previamente arbitrada, então 0>kq , se escoa na direção contrária, 0<kq ); as pressões nos vértices, representadas pelas variáveis Pt, onde as pressões nos vértices terminais de um elemento k também podem ser representadas por Pk1 e Pk2; e, as vazões externas nos vértices que indicam a transferência de fluido entre a rede e o exterior, representadas pela variável wt ( 0>tw se há entrada de fluido e 0<tw se há saída de fluido). No total, o problema de simulação estacionária de uma rede envolve 2N S+ variáveis. As equações do modelo da rede estão apresentadas a seguir: A) Balanço material nos vértices:

∑ ∑∈ ∈

=+−int

outtQk

tQk

kk wqq 0

para Nt ,...,1= (A2-1)

Onde in

tQ e outtQ são, respectivamente, os conjuntos de elementos que entram e saem do vértice t de acordo

com a direção arbitrada. B) Equações de conservação da energia ao longo dos elementos:

P P F q g z zk k k k k k1 2 2 1 0− − − − =( ) ( )ρ para k S= 1,..., (A2-2) onde zk1 e zk2 são as elevações nos vértices terminais. Para um trecho de tubulação k, a função Fk corresponde ao termo de perda de carga. Este termo pode ser calculado pela equação de Darcy-Weisbach, representada aqui em função da vazão de fluido:

F q fLD

q qk k kk

kk k( ) =

82 5

ρπ

Onde fk é o fator de atrito de Darcy, ρ é a densidade do fluido, Lk e Dk são o comprimento e diâmetro da tubulação. Finalmente, para efetuar a simulação é necessário especificar o valor de N variáveis. Usualmente, em cada vértice é fixado ou o valor da vazão externa ou o valor da pressão.


APÊNDICE 5 GRÁFICOS

Fig. 5.1 - Gráficos da tensão de cisalhamento na fase líquida. Fig. 5.1 – Shear stress charts in the liquid phase.

Fig. 5.2 - Gráficos da tensão de cisalhamento na interface. Fig. 5.2 – Shear stress charts in the interface.

Fig. 5.3 - Gráficos de Número de Reynolds da fase líquida. Fig. 5.3 –Reynolds number charts of the liquid phase.


Fig. 5.4 - Gráficos de Número de Reynolds da fase gás. Fig. 5.4 – Reynolds number charts of the gas phase.

Fig. 5.5 - Gráficos de Número de Reynolds na interface. Fig. 5.5 – Reynolds number charts in the interface.

Fig. 5.6 - Gráficos de Número de Reynolds na interface. Fig. 5.6 – Reynolds number charts in the interface


Fig. 5.7 - Gráficos de Número de Reynolds na interface. Fig. 5.7 – Reynolds number charts in the interface.


NOTAÇÃO LETRAS LATINAS SÍMBOLO DESCRIÇÃO A Área da seção transversal (m x 1) A Coeficiente de Churchill A Matriz de coeficientes do modelo (4x4) B Coeficiente de Churchill B Vetor de coeficientes do modelo (4x1) C Constante termodinâmica do gás D Diâmetro da tubulações (mx1) E Especificações do modelo (4x1) F Fator de atrito, coeficiente do modelo de escoamento h Comprimento dos elementos finitos (mx1) Hk vetor de resíduos do duto k (4Nkx1) Jk vetor de resíduos do duto k (4Nkx1) K Perímetro de contato V Volume M Massa M Matriz de incidência (NVxm)

KM Matriz de ponderação do MEF para o duto k [4(Nk+1)x4(Nk+1)]

P Pressão q Vazão mássica Q Vazão mássica dos vértices (NVx1)

WS Matriz seletora (Y→W) [4Nkx4(Nk+1)]

SS Matriz seletora (W→Y) [4(Nk+1) x4Nk]

ES Matriz seletora (E→Y) [4(Nk+1) x4]

t Tempo v Velocidade W Vetor de estados reduzidos (4x1) x Coordenada espacial y Variáveis de estado (4x1) Y Vetor de estados nodais (4Nkx1) LETRAS GREGAS SÍMBOLO DESCRIÇÃO α Fração de ocupação do gás β Arco em contato com a fase líquida γ Expoente politrópico δ Espessura da camada de líquido θ Inclinação do trecho de tubulação φ Função pulso triangular; função seletora de tensão ρ Rugosidade dos dutos (mx1) µ Viscosidade do fluido ψ Tensão de cisalhamento Λ Matriz de coeficientes integrados do modelo [4(Nk+1) x4Nk]

Ξ Vetor de coeficientes integrados do modelo [4(Nk+1) x1]


ℜ Vetor de resíduos (4x1) ÍNDICES SÍMBOLO DESCRIÇÃO 0 Estado padrão 1 Início do trecho de tubulação 2 Final do trecho de tubulação an Escoamento anular est Escoamento estratificado G Fase gás i Tipo de fase j elemento finito NV Número de vértices da rede M Número de dutos da rede L Fase líquida T Transformada A Classe A B Classe B I Interface

Documents

MODELAGEM DINÂMICA DO ESCOAMENTO BIFÁSICO EM … · Bol. Téc. Petrobras, Rio de Janeiro, 47 (1): 49-87, jan./mar., 2004 49 MODELAGEM DINÂMICA DO ESCOAMENTO BIFÁSICO EM RISERS