OtimizaçãoeParalelizaçãodeumSimuladordeBombeio ... · ... Sistema de Bibliotecas ... de um simulador de bombeio mecânico para poços ... cavalo-de-pau”. A ﬁnalidade da unidade

Roger Rommel Ferreira de Araújo

Otimização e Paralelização de um Simulador de BombeioMecânico para Poços Direcionais

Dissertação de Mestrado apresentada aoPrograma de Pós-Graduação em EngenhariaElétrica e de Computação da UFRN comoparte dos requisitos para obtenção do títulode Mestre em Ciências.

Universidade Federal do Rio Grande do Norte

Centro de Tecnologia

Programa de Pós-Graduação em Engenharia Elétrica e de Computação

Orientador: Prof. Dr. Samuel Xavier de Souza

Número de Ordem do PPgEEC: M474

Natal – RN

Dezembro de 2016

Catalogação da Publicação na Fonte

Universidade Federal do Rio Grande do Norte - Sistema de Bibliotecas Biblioteca Central Zila Mamede / Setor de Informação e Referência Araújo, Roger Rommel Ferreira de.

Otimização e paralelização de um simulador de bombeio mecânico para poços direcionais / Roger Rommel Ferreira de Araújo. - 2016.

69 f. : il. Dissertação (mestrado) - Universidade Federal do Rio Grande

do Norte, Centro de Tecnologia, Programa de Pós-Graduação em Engenharia Elétrica e de Computação. Natal, RN, 2016.

Orientador: Prof. Dr. Samuel Xavier de Souza. 1. Bombeamento - Dissertação. 2. Poços tubulares -

Dissertação. 3. Simuladores - Dissertação. 4. Extração de petróleo - Dissertação. 5. Bombeio mecânico - Dissertação. 6. Poços direcionais - Dissertação. I. Souza, Samuel Xavier de. II. Título.

RN/UF/BCZM CDU

004:621.65

CWBM

Carimbo

Para Aline

Agradecimentos

A Deus, pela vida, saúde e perseverança para a realização deste trabalho;

A minha esposa Aline, por me incentivar desde sempre, e por compreender asausências que este trabalho exigiu;

Aos meus pais, Natan e Maria do Carmo, pela dedicação a mim e a meus irmãos, epor acreditarem em nosso potencial;

Ao meu orientador, Prof. Dr. Samuel Xavier de Souza, pela confiança e auxílioconstantes ao longo do caminho;

Aos engenheiros Rutácio Costa e Antônio Araújo Jr., pelos materiais de referênciae pelas muitas contribuições à qualidade da pesquisa;

À Petrobras, cujo apoio ao crescimento profissional e intelectual de seu corpotécnico não é apenas uma política empresarial, mas uma orgulhosa tradição;

Aos colegas da Petrobras e da UFRN, pelo encorajamento e camaradagem;

Aos amigos da Caverna e da Faculdade de Computação da UFC, pelo bom humor,pelos debates e pelas risadas;

Às Musas de Clio, que me ajudaram a analisar vários temas através de novos pontosde vista;

Às comunidades do software livre e do código aberto, pelas incríveis ferramentasque põem ao alcance de todos, e por fazerem do mundo um lugar melhor.

Uma jornada de mil milhascomeça com um passo.(Provérbio Japonês)

ResumoPoços de bombeio mecânico, que são os mais comuns na indústria do petróleo, podem tercolunas de produção verticais ou direcionais. Apesar de vários modelos matemáticos teremsido propostos na literatura para descrever o comportamento dinâmico de poços de bombeiomecânico direcionais, o trabalho de Costa (1995) foi o primeiro entre eles a combinarnum modelo unificado uma variedade de fatores essenciais de simulação. Infelizmente, suacomplexidade computacional limitou seu uso até agora. Neste trabalho, são apresentadasdiversas otimizações de desempenho para as simulações deste modelo, com o objetivo detornar sua utilização viável. As otimizações propostas são descritas em detalhes, e seudesempenho é avaliado sob diferentes cenários. Os resultados mostram que os tempos deexecução do novo simulador o tornam uma solução viável. Engenheiros de petróleo podemusar este simulador em estruturas de computação de alto desempenho para otimizar oprojeto e os parâmetros operacionais de poços de bombeio mecânico direcionais.

Palavras-chave: Bombeamento. Poços tubulares. Simuladores. Extração de petróleo.Bombeio mecânico. Poços direcionais. Processamento paralelo.

AbstractSucker-rod pumping wells, which are the most common in the oil industry, can havevertical or directional production strings. While several mathematical models have beenproposed in the literature to describe the dynamic behavior of directional sucker-rod wells,the work of Costa (1995) was the first among them to combine a number of essentialsimulation characteristics into a unified model. Unfortunately, its computational complexityhas limited its practical use up to now. In this work, several performance optimizationsare presented for the simulations of this model in order to make its use feasible. Allthe proposed optimizations are described in detail and their performance is evaluatedunder different scenarios. The results show that the execution times of the new simulatormake it a viable solution. Petroleum engineers can use this simulator in high-performancecomputing infrastructures to optimize the design and operating parameters of directionalsucker-rod pumping wells.

Keywords: Pumping. Tubular wells. Simulators. Petroleum extraction. Sucker-rod pump-ing. Directional wells. Parallel processing.

Lista de ilustrações

Figura 1 – Poço de bombeio mecânico. Fonte: Sucker-rod Pumping Manual (TA-KÁCS, 2003), adaptado. . . . . . . . . . . . . . . . . . . . . . . . . . . 19

Figura 2 – Cartas dinamométricas de superfície e de fundo. . . . . . . . . . . . . . 20Figura 3 – Tempos de execução - Modelo de AV de Gibbs . . . . . . . . . . . . . . 47Figura 4 – Tempos de execução - Modelo de AV de Lea . . . . . . . . . . . . . . . 47Figura 5 – Speedup - Modelo de AV de Gibbs . . . . . . . . . . . . . . . . . . . . . 48Figura 6 – Speedup - Modelo de AV de Lea . . . . . . . . . . . . . . . . . . . 48Figura 7 – Eficiência - Modelo de AV de Gibbs . . . . . . . . . . . . . . . . . . . . 48Figura 8 – Eficiência - Modelo de AV de Lea . . . . . . . . . . . . . . . . . . . . . 48Figura 9 – Tempos de execução - Parte 4 do modelo de AV de Lea . . . . . . . . . 52Figura 10 – Poço 1: Cartas dinamométricas simuladas/reais de superfície, e cartas

derivadas de fundo - Modelo de AV de Gibbs. . . . . . . . . . . . . . . 63Figura 11 – Poço 1: Cartas dinamométricas simuladas/reais de superfície, e cartas

derivadas de fundo - Modelo de AV de Lea. . . . . . . . . . . . . . . . 63Figura 12 – Poço 2: Cartas dinamométricas simuladas/reais de superfície, e cartas


derivadas de fundo - Modelo de AV de Lea. . . . . . . . . . . . . . . . 66Figura 14 – Poço 3: Cartas dinamométricas simuladas/reais de superfície, e cartas


derivadas de fundo - Modelo de AV de Lea. . . . . . . . . . . . . . . . 69

Lista de tabelas

Tabela 1 – Coeficientes de Fourier para unidades de bombeio convencionais. Fonte:SPE-19724 (LAINE; COLE; JENNINGS, 1989). . . . . . . . . . . . . . 24

Tabela 2 – Dados da configuração de poço. . . . . . . . . . . . . . . . . . . . . . . 45Tabela 3 – Perfil do poço. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46Tabela 4 – Tempos de execução - Modelo de AV de Gibbs . . . . . . . . . . . . . . 46Tabela 5 – Tempos de execução - Modelo de AV de Lea . . . . . . . . . . . . . . . 47Tabela 6 – Speedup / Eficiência - Modelos de AV de Gibbs e Lea . . . . . . . . . . 47Tabela 7 – Tempos absolutos de execução - Modelo de AV de Gibbs - ∆s = 0,05 m 49Tabela 8 – Tempos percentuais de execução - Modelo de AV de Gibbs - ∆s = 0,05 m 50Tabela 9 – Tempos absolutos de execução - Modelo de AV de Lea - ∆s = 0,05 m . 50Tabela 10 – Tempos percentuais de execução - Modelo de AV de Lea - ∆s = 0,05 m 50Tabela 11 – Tempos de execução e speedup - Parte 4 do modelo de AV de Lea . . . 52Tabela 12 – Poço 1: Dados da configuração de poço. . . . . . . . . . . . . . . . . . 61Tabela 13 – Poço 1: Perfil do poço. . . . . . . . . . . . . . . . . . . . . . . . . . . . 62Tabela 14 – Poço 1: Dados reais e simulados. . . . . . . . . . . . . . . . . . . . . . 62Tabela 15 – Poço 2: Dados da configuração de poço. . . . . . . . . . . . . . . . . . 64Tabela 16 – Poço 2: Perfil do poço. . . . . . . . . . . . . . . . . . . . . . . . . . . . 65Tabela 17 – Poço 2: Dados reais e simulados. . . . . . . . . . . . . . . . . . . . . . 65Tabela 18 – Poço 3: Dados da configuração de poço. . . . . . . . . . . . . . . . . . 67Tabela 19 – Poço 3: Perfil do poço. . . . . . . . . . . . . . . . . . . . . . . . . . . . 68Tabela 20 – Poço 3: Dados reais e simulados. . . . . . . . . . . . . . . . . . . . . . 68

Lista de abreviaturas e siglas

API American Petroleum Institute

AV Atrito viscoso

CPU Unidade central de processamento (“central processing unit”)

EDP Equação diferencial parcial

GPU Unidade de processamento gráfico (“graphics processing unit”)

MPRL Carga mínima na haste polida (“minimum polished rod load”)

PD Vazão volumétrica de fluido realizada diariamente pela bomba de fundo(“pump displacement”)

PPRL Carga máxima na haste polida (“peak polished rod load”)

PRHP Potência na haste polida (“polished rod horsepower”)

PT Torque máximo no redutor da unidade de bombeio (“peak torque”)

RP Prática recomendada (“recommended practice”)

Lista de símbolos

Ai Coeficientes de Fourier, adimensionais

Ap Área da seção transversal do pistão (m2)

Ar(s), Ark(s) Área da seção transversal da haste sobre a qual um determinadoponto da coluna de hastes está localizado (m2)

At Área da seção transversal da coluna de produção (m2)

Bi Coeficientes de Fourier, adimensionais

~B(s) Vetor binormal unitário num determinado ponto da coluna de hastes

cD Fator de amortecimento (adimensional)

Db Profundidade medida de instalação da bomba de fundo (m)

df Densidade do fluido (adimensional)

em Espaço “morto” na bomba de fundo (m)

e(t) Elongação da coluna de produção, no momento t (m)

Et Módulo de elasticidade de Young relativo ao material da coluna deprodução num determinado ponto da coluna de hastes (psi)

fv(s, t) Força de atrito viscoso por unidade de comprimento num determinadoponto da coluna de hastes, no momento t (N)

f jci Força de atrito de Coulomb por unidade de comprimento na haste em

que um ponto i da coluna de hastes está localizado, num momento j(N)

Fs(s, t) Força axial num determinado ponto da coluna de hastes, no momento t(N)

F jSi Força axial por unidade de comprimento na haste em que um ponto i

da coluna de hastes está localizado, num momento j (N)

~g Aceleração da gravidade (m/s2)

hb Profundidade vertical de instalação da bomba de fundo (m)

ianc Valor inteiro que indica se a coluna de produção está ancorada

kL Compressibilidade média da fase líquida do fluido (psi-1)

~K(s) Vetor de curvatura num determinado ponto da coluna de hastes

K1(s), K2(s), K3(s), K4(s) Fatores geométricos, derivados do diâmetro da colunade produção e das hastes, para um determinado ponto da coluna dehastes

Lb(t) Distância entre a válvula de pé e a válvula de passeio na bomba defundo, no momento t (m)

Lk Comprimento da seção de hastes sobre a qual um determinado pontoda coluna de hastes está localizado (m)

~N(s) Vetor normal unitário num determinado ponto da coluna de hastes

pb(t) Pressão absoluta no interior da bomba de fundo, no momento t (psia)

pd(t) Pressão de descarga na bomba de fundo, no momento t (psi)

pf Pressão dinâmica exercida pelo fluido (psi)

pg Gradiente de pressão da água (psi/m)

ps Pressão de sucção na bomba de fundo (psi)

pwh Pressão na coluna de produção, medida na cabeça do poço (psi)

rc(s) Raio de curvatura do poço num determinado ponto da coluna de hastes(m)

s Distância medida a partir da bomba de fundo até um determinadoponto da coluna de hastes ao longo da coluna de produção, em t = 0(m)

S Comprimento do curso da unidade de bombeio (in)

t Tempo decorrido desde o início do movimento (segundos)

~T (s) Vetor tangente unitário num determinado ponto da coluna de hastes

TF (t) Valor da série de Fourier truncada, no momento t (rad)

u(s, t) Deslocamento de um determinado ponto da coluna de hastes, a partirde sua posição inicial s, no momento t (m)

Urk(s) Perímetro da seção circular da haste em que um determinado ponto dacoluna de hastes está localizado (m)

v Velocidade do som no interior da coluna de hastes (m/s)

v̄fk(s, t) Velocidade média do fluido num determinado ponto da coluna de hastes,no momento t (m/s)

vhp(t) Velocidade da haste polida, no momento t (pés/s)

vp(t) Velocidade do pistão, no momento t (m/s)

vr(s, t) Velocidade longitudinal de um determinado ponto da coluna de hastes,no momento t (m/s)

vjri Velocidade longitudinal da haste em que um ponto i da coluna de hastes

está localizado, num momento j (m/s)

α Fração volumétrica de gás no interior da bomba de fundo

∆s Comprimento do segmento escolhido para dividir a coluna de hastes empontos adjacentes (m)

∆t Intervalo de tempo entre iterações da simulação (segundos)

η(s) Viscosidade dinâmica do fluido num determinado ponto da coluna dehastes (centipoise)

µ Coeficiente de atrito de Coulomb (adimensional)

ρf Massa específica do fluido (kg/m3)

ρr(s) Massa específica da haste sobre a qual um determinado ponto da colunade hastes está localizado (kg/m3)

ω Velocidade angular da haste polida (rad/s)

Sumário

1 INTRODUÇÃO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171.1 Motivação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171.2 Objetivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181.3 Organização do texto . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2 REVISÃO BIBLIOGRÁFICA . . . . . . . . . . . . . . . . . . . . . . 19

3 REFERENCIAL TEÓRICO . . . . . . . . . . . . . . . . . . . . . . . 223.1 Bombeio mecânico em poços direcionais . . . . . . . . . . . . . . . . 223.1.1 A Equação de Movimento . . . . . . . . . . . . . . . . . . . . . . . . . . 223.1.2 Solução Numérica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 273.1.2.1 Colunas de Hastes Combinadas . . . . . . . . . . . . . . . . . . . . . . . . . 293.1.2.2 Estabilidade Numérica e Critério de Parada . . . . . . . . . . . . . . . . . . . 313.2 Paralelismo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313.2.1 A Lei de Amdahl . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

4 METODOLOGIA . . . . . . . . . . . . . . . . . . . . . . . . . . . . 354.1 Implementação paralela . . . . . . . . . . . . . . . . . . . . . . . . . . 354.1.1 Considerações de Otimização e Paralelização . . . . . . . . . . . . . . . . 374.1.2 Análise de Complexidade . . . . . . . . . . . . . . . . . . . . . . . . . . . 384.1.2.1 Complexidade Paralela . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 404.1.3 Contribuições ao Trabalho Original . . . . . . . . . . . . . . . . . . . . . . 41

5 SIMULAÇÕES E RESULTADOS . . . . . . . . . . . . . . . . . . . . 44

6 CONCLUSÕES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

REFERÊNCIAS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

APÊNDICES 57

APÊNDICE A – ESQUEMAS DE DIFERENÇAS FINITAS PARACÁLCULO DA FORÇA E VELOCIDADE DAS HAS-TES . . . . . . . . . . . . . . . . . . . . . . . . . . 58

APÊNDICE B – DADOS REAIS E SIMULADOS DE POÇOS DEPETRÓLEO . . . . . . . . . . . . . . . . . . . . . 60

17

1 Introdução

Neste capítulo, descrevemos a motivação para que se tente simular com eficiênciapoços de bombeio mecânico direcionais, os objetivos deste trabalho e a forma como o textoestá organizado.

1.1 MotivaçãoO bombeio mecânico é o sistema de elevação artificial de fluidos mais utilizado

na indústria de petróleo e gás, abrangendo cerca de 80% dos poços de petróleo ao redordo mundo (COSTA, 2004). Poços de bombeio mecânico são especialmente adequados acenários de baixa vazão de óleo e produção desprezível de gás, comumente encontrados, porexemplo, na área terrestre da Bacia Potiguar. Modelos matemáticos para o comportamentode poços de bombeio mecânico com colunas de produção verticais, bem como técnicas paraprojetar tais poços, podem ser encontrados em diversas fontes, tais como a norma API RP11L (AMERICAN PETROLEUM INSTITUTE, 2000) e Miska, Sharaki e Rajtar (1997).

Porém, às vezes pode ser indesejável ou impraticável perfurar um poço com uma tra-jetória totalmente vertical. Por exemplo, um poço vertical pode não aproveitar plenamenteas características de um determinado reservatório, ou a perfuração vertical pode se revelarinviável em consequência de dificuldades locais, tais como a existência de cursos d’água ousalinas. Em situações como estas, faz-se necessário projetar o poço prevendo inclinaçõesem um ou mais pontos de sua coluna de produção. Tais poços são ditos direcionais.

Poços de bombeio mecânico direcionais comportam-se de maneira muito diferentedos poços verticais, exigindo técnicas específicas de projeto. Para desenvolver estas téc-nicas, pesquisadores tais como Evchenko e Zakharchenko (1984), Lukaziewicz (1991) eGibbs (1992) propuseram modelos matemáticos. Embora a utilidade destes modelos sejaincontestável, nenhum deles oferecia uma descrição unificada do comportamento da colunade hastes direcional juntamente com o escoamento dos fluidos no interior do anular entrea coluna de produção e a coluna de hastes. Tentativas de preencher esta lacuna foramrealizadas por Costa (1995), Xu e Hu (1993), Xu (1994a), Xu (1994b) e Xu et al. (1999).

Em seu trabalho, Costa propôs um modelo integrado que produzia resultadosconfiáveis em comparação a dados medidos em poços reais. Porém, a complexidadede seu modelo exigia um poder de processamento muito maior do que as abordagensanteriores. O modelo é baseado numa equação diferencial parcial (EDP) de segundaordem, resolvida numericamente através de um método de diferenças finitas explícito.À época, o simulador elaborado por Costa para demonstrar seu modelo levava até 108

Capítulo 1. Introdução 18

segundos para simular uma única configuração de poço. Tempos de execução como essese mostraram excessivamente lentos, restringindo a liberdade de experimentar muitascombinações diferentes de equipamentos e variáveis de operação.

O modelo de Costa se revelou impraticável para uso diário, e não houve novosesforços desde então no sentido de criar soluções baseadas neste trabalho. Num escopo maisamplo, também não surgiram novos modelos integrados para simular poços de bombeiomecânico direcionais. Ao analisar pesquisas recentes sobre poços de bombeio mecânico,surgem resultados tais como trabalhos sobre o projeto e simulação de poços de metanoadsorvido no carvão (XINFU et al., 2010; XINFU et al., 2011), que apresentam diferençasconsideráveis em relação a poços de petróleo; pode-se citar também um novo modelo desimulação proposto por Dong et al. (2013) para poços de petróleo, mas infelizmente estetrabalho não parece cobrir poços desviados. É importante mencionar ainda o trabalho deAraujo et al. (2015), que apresenta um algoritmo para calcular a carta dinamométrica defundo de um poço de bombeio mecânico direcional a partir de sua carta de superfície.

1.2 ObjetivosNeste trabalho revisitamos o modelo de Costa, na tentativa de melhorar seu

desempenho e transformá-lo numa solução viável para uso diário. Para esse fim, utilizamosferramentas de programação modernas para implementá-lo em versões serial e paralela,comparando o desempenho e a escalabilidade entre elas. Dá-se ênfase especial à discussãodos aspectos de paralelização. Os resultados obtidos são bastante animadores, especialmenteao utilizar múltiplas linhas de execução. Eles mostram ainda que há considerável espaçopara novas extensões ao modelo, e que a implementação proposta permite processar grandenúmero de configurações de poços de maneira eficiente, viabilizando seu uso para tarefastais como otimização de parâmetros de poços.

1.3 Organização do textoEste trabalho é organizado conforme descrito a seguir. No Capítulo 2, fazemos a

revisão bibliográfica das principais fontes usadas para basear nossas discussões. No Capítulo3, oferecemos um resumo do trabalho de Costa (1995), examinando as características maisimportantes de seu modelo matemático e a solução de suas equações, e comentamos aimportância e a necessidade da utilização do processamento paralelo na construção desoluções de software modernas. No Capítulo 4, descrevemos a implementação paralela pornós desenvolvida, analisando-a em detalhes. No Capítulo 5, discutimos as simulações querealizamos, bem como seus resultados. Finalmente, no Capítulo 6, apresentamos nossasconclusões e considerações finais.

19

2 Revisão Bibliográfica

Segundo Takács (2003), os componentes individuais de um poço de bombeiomecânico podem ser divididos em dois grandes grupos: equipamentos de superfície eequipamentos de subsuperfície. Os elementos principais de uma instalação típica sãomostrados na Figura 1.

Figura 1 – Poço de bombeio mecânico. Fonte: Sucker-rod Pumping Manual (TAKÁCS, 2003),adaptado.

O conjunto dos equipamentos de superfície, partindo do motor e prosseguindoaté a haste polida, normalmente é chamado de unidade de bombeio ou “cavalo-de-pau”.A finalidade da unidade de bombeio é converter a rotação do motor num movimentoperiódico ascendente e descendente da haste polida e da coluna de hastes, que resultarána elevação dos fluidos contidos no reservatório até a cabeça do poço.

O objetivo principal do método de Costa (1995), dada uma configuração específicade equipamentos e variáveis de operação, é estimar os parâmetros que seguem para o poço:

• PPRL (“peak polished rod load”): carga máxima na haste polida;

• MPRL (“minimum polished rod load”): carga mínima na haste polida;

Capítulo 2. Revisão Bibliográfica 20

• PT (“peak torque on the gear reducer”): torque máximo no redutor da unidade debombeio;

• PRHP (“polished rod horsepower”): potência na haste polida;

• PD (“pump displacement”): vazão volumétrica de fluido realizada diariamente pelabomba de fundo.

O modelo não calcula estes parâmetros diretamente. Em vez disso, calcula cartasdinamométricas de superfície e cartas dinamométricas de fundo, a partir das quais osparâmetros de operação são derivados. Cartas de superfície são compostas de pares devalores que correlacionam forças de tração na haste polida a posições ao longo do cursopercorrido pela unidade de bombeio. Cartas de fundo são definidas de maneira semelhante,mas consideram forças de tração no pistão de subsuperfície em vez da haste polida.

Um ciclo de bombeio tem duas fases, ascendente e descendente, que correspondemaos períodos em que a haste polida se desloca para cima ou para baixo, respectivamente.Como a haste polida percorre todo o curso da unidade de bombeio duas vezes a cadaciclo de bombeio, primeiro no curso ascendente e em seguida no curso descendente, ascartas dinamométricas correlacionam dois valores de força a cada posição ao longo docurso da unidade de bombeio, com um valor para cada fase do ciclo, conforme exibido naFigura 2. Na figura citada, tanto para a carta de superfície quanto para a carta de fundo,o curso ascendente começa na posição 0 no eixo horizontal e prossegue até cerca de 44polegadas. Inicia-se então o curso descendente, no qual a posição decresce até o valor 0 noeixo horizontal e fecha-se então o ciclo de bombeio.

0 10 20 30 40

0

2000

4000

6000

Posição [in]

Carga

[lbf]

Carta de superfície Carta de fundo

Figura 2 – Cartas dinamométricas de superfície e de fundo.

Costa aponta que, para uma correta simulação de poços de bombeio mecânicodirecionais, há que levar em conta diversos fatores específicos de seu comportamento,conforme segue:

Capítulo 2. Revisão Bibliográfica 21

• Dinâmica da coluna de hastes direcional;

• Atrito direcional entre as hastes da coluna de hastes e a coluna de produção;

• Atrito viscoso entre as hastes da coluna de hastes e os fluidos;

• Se existe gás junto à bomba de fundo, bem como o volume correspondente;

• Se a coluna de produção está ancorada.

Gibbs (1963) e Lukaziewicz (1991) propuseram soluções para o primeiro e segundofatores, com uma restrição: os poços poderiam ter trajetórias direcionais, contanto queestivessem contidas num plano vertical. Quanto ao terceiro fator, apontamos as contri-buições de Lea (1991) e Gibbs (1963). Em relação ao quarto e quinto fatores, citamos otrabalho de diversos pesquisadores: Gibbs (1963), Gibbs (1977), Doty e Schmidt (1983),Alhanati (1988), Laine, Cole e Jennings (1989) e Lea (1991).

Quando o trabalho de Costa foi anunciado, tinha a característica marcante de ser oprimeiro documento disponível publicamente onde todos os fatores citados anteriormenteeram considerados de forma simultânea, além de permitir a utilização de trajetóriasverdadeiramente direcionais. Também mostra que, dadas as restrições adequadas, aspropostas anteriores podem ser caracterizadas como casos particulares do modelo integradode Costa (1995). O trabalho de Costa (1995) é a fonte central sobre a qual se baseiamnossas discussões.

22

3 Referencial teórico

Neste capítulo, fazemos uma breve exposição do modelo matemático propostopor Costa e de uma solução numérica para ele. Além disso, discutimos o histórico e aimportância da computação paralela, e apresentamos os conceitos de speedup, eficiência eescalabilidade.

3.1 Bombeio mecânico em poços direcionaisEm seu estudo, Costa propôs tanto um modelo matemático para o movimento da

coluna de hastes em poços direcionais quanto uma solução numérica para este modeloatravés de um método de diferenças finitas explícito, os quais serão resumidos nas seçõesque se seguem. O leitor que desejar explicações mais aprofundadas deve consultar adissertação de Costa (1995).

3.1.1 A Equação de Movimento

Costa propôs a seguinte EDP de segunda ordem para descrever o movimento dacoluna de hastes em poços direcionais:

∂2u

∂t2= v2∂

2u

∂s2 + ~g · ~T (s)

− µ vr

|vr|

√√√√[~g · ~B(s)]2

+[~g · ~N(s) + v2

rc(s)∂u

∂s

]2

+ fv

ρrAr

,

(3.1)

onde:t = tempo decorrido desde o início do movimento (segundos);s = distância medida, a partir da bomba de fundo, até um determinado ponto da colunade hastes ao longo da coluna de produção, em t = 0 (m);u(s, t) = deslocamento de um determinado ponto da coluna de hastes, a partir de suaposição inicial s, no momento t (m);v = velocidade do som no interior da coluna de hastes (m/s);~g = aceleração da gravidade (m/s2);~T (s) = vetor tangente unitário num determinado ponto da coluna de hastes;µ = coeficiente de atributo de Coulomb (adimensional);vr(s, t) = velocidade longitudinal de um determinado ponto da coluna de hastes, nomomento t (m/s);

Capítulo 3. Referencial teórico 23

~B(s) = vetor binormal unitário num determinado ponto da coluna de hastes;~N(s) = vetor normal unitário num determinado ponto da coluna de hastes;rc(s) = raio de curvatura do poço num determinado ponto da coluna de hastes (m);fv(s, t) = força de atrito viscoso por unidade de comprimento num determinado ponto dacoluna de hastes, no momento t (N/m);ρr(s) = massa específica da haste sobre a qual um determinado ponto da coluna de hastesestá localizado (kg/m3);Ar(s) = área da seção transversal da haste em que um determinado ponto da coluna dehastes está localizado (m2).

O último termo da equação, fv

ρrAr

, refere-se ao atrito viscoso (“AV”) entre as hastese os fluidos. Pode-se substituir este termo utilizando uma expressão proposta por Gibbs(1963):

fv

ρrAr

= −c∂u∂t

, (3.2)

onde c é um coeficiente de amortecimento dado por:

c = πvcD

2Db

, (3.3)

e:cD = fator de amortecimento (adimensional);Db = profundidade medida de instalação da bomba de fundo (m).

Há que reparar que a Eq. (3.3) leva em conta tanto o atrito viscoso quanto o atritodirecional, o que exige escolher valores apropriados para cD.

Outra alternativa para fv

ρrAr

, baseada no trabalho de Lea (1991), é dada por:

fv

ρrAr

= ηUrk

ρrArk

(K1vr −K2v̄fk) , (3.4)

onde:η(s) = viscosidade dinâmica do fluido num determinado ponto da coluna de hastes(centipoise);Urk(s) = perímetro da seção circular da haste em que um determinado ponto da colunade hastes está localizado (m);ρr(s) = massa específica da haste sobre a qual um determinado ponto da coluna de hastesestá localizado (kg/m3);Ark(s) = área da seção transversal da haste em que um determinado ponto da coluna dehastes está localizado (m2);K1(s), K2(s) = fatores geométricos, derivados dos diâmetros da coluna de produção e dashastes, para um determinado ponto da coluna de hastes;vr(s, t) = velocidade longitudinal de um determinado ponto da coluna de hastes, no


momento t (m/s);v̄fk(s, t) = velocidade média do fluido num determinado ponto da coluna de hastes, nomomento t (m/s).

O leitor que desejar maiores detalhes quanto às expressões utilizadas para o cálculode K1(s) e K2(s) deve consultar o trabalho de Lea (1991) ou Costa (1995).

A opção de usar o modelo de atrito viscoso de Gibbs ou o de Lea cabe ao engenheiro.Porém, esta opção é extensiva ao cálculo da pressão de descarga da bomba de fundo,conforme se discutirá posteriormente neste trabalho.

Costa aponta algumas condições de contorno que devem ser observadas. Começandocom os equipamentos de superfície, o movimento da haste polida pode ser aproximadopor uma série de Fourier truncada em seis termos, conforme sugerido por Laine, Cole eJennings (1989). O valor da série de Fourier truncada para um determinado momento notempo é dado pela equação a seguir:

TF (t) = A1 cos(ωt) + · · ·+ A6 cos(6ωt)

+B1 sin(ωt) + · · ·+B6 sin(6ωt) ,(3.5)

onde:A1 a A6, B1 a B6 = coeficientes de Fourier;ω = velocidade angular da haste polida (rad/s).

A velocidade da haste polida para um determinado momento no tempo, em pés/s,é dada por:

vhp(t) = TF (t)Sω12 , (3.6)

onde S é o comprimento do curso da haste polida em polegadas.

Os coeficientes de Fourier a utilizar dependem da geometria de unidade de bombeio.Para unidades de bombeio convencionais, Laine et al. recomendam os valores listados naTabela 1.

Tabela 1 – Coeficientes de Fourier para unidades de bombeio convencionais. Fonte: SPE-19724(LAINE; COLE; JENNINGS, 1989).

i Ai Bi

1 0,0078489 0,49730542 0,0123680 0,06307663 -0,017086 0,00715854 -0,002505 0,00142885 -0,000555 -0,0008326 -0,000123 -0,000070

Outras alternativas para a aproximação do movimento da haste polida são o modelosenoidal, mais simples, e o modelo geométrico, mais preciso. O modelo geométrico requer


o conhecimento prévio da geometria da unidade de bombeio, descrita pelas dimensões A,C, I, K, P e R constantes do documento API SPEC 11E (AMERICAN PETROLEUMINSTITUTE, 2015). Dito isso, não se fará aqui uma discussão mais aprofundada sobreestes modelos.

Serão discutidas agora as condições de contorno para os equipamentos de subsu-perfície. A distância entre a válvula de pé e a válvula de passeio na bomba de fundo paraum determinado momento no tempo é dada por:

Lb(t) = u(0, t) + em + et(t) , (3.7)

onde:em = espaço “morto” na bomba de fundo. Esta é a distância entre a válvula de pé e aválvula de passeio quando o sistema está em repouso, estando a haste polida na posiçãomais baixa possível (m);et(t) = elongação da coluna de produção, no momento t (m).

A elongação da coluna de produção pode ser calculada por:

et(t) = [pb(0)− pb(t)]DbApianc

EtAt

, (3.8)

onde:pb(t) = pressão absoluta no interior da bomba de fundo (psia);Ap = área da seção transversal do pistão (m2);ianc = valor inteiro que indica se a coluna de produção está ancorada. Se a coluna deprodução estiver ancorada, este valor será 0; caso contrário, este valor será 1;Et = módulo de elasticidade de Young relativo ao material da coluna de produção (psi);At = área da seção transversal da coluna de produção (m2).

A variação na pressão absoluta no interior da bomba de fundo é dada pela expressãoa seguir:

dpb

du= − 1

DbApianc

EtAt

+ kLLb(1− α) + Lbα

pb

,

com ps < pb(t) < pd(t) ,

(3.9)

onde:kL = compressibilidade média da fase líquida do fluido (psi-1);α = fração volumétrica de gás no interior da bomba de fundo (porcentagem variando de 0a 1, onde 1 significa 100%);pb(t) = pressão absoluta no interior da bomba de fundo a partir da qual a variação serácalculada, no momento t (psia);


ps = pressão de sucção na bomba de fundo (psi);pd = pressão de descarga na bomba de fundo, no momento t (psi).

Caso o engenheiro opte por trabalhar com o modelo de Gibbs, a pressão de descargada bomba de fundo pode ser calculada por:

pd = pwh + dfpghb , (3.10)

onde:pwh = pressão na coluna de produção, medida na cabeça do poço (psi);df = densidade do fluido (adimensional);pg = gradiente de pressão da água (1,42 psi/m);hb = profundidade vertical de instalação da bomba de fundo (m).

Se o engenheiro preferir trabalhar com o modelo de Lea, deve ser usada a expressãoque segue:

pd = pwh + dfpghb +n∑

k=1

(∂pf

∂s

)k

Lk , (3.11)

onde Lk é o comprimento da seção de hastes sobre a qual um determinado ponto dacoluna de hastes está localizado, em metros, e a perda de carga causada pelo atrito numdeterminado ponto da coluna de hastes (psi/m) é dada por:

(∂pf

∂s

)k

= −4η(K3v̄fk +K4vr) , (3.12)

onde K3(s) e K4(s) são fatores geométricos, derivados dos diâmetros da coluna de produçãoe das hastes, para um determinado ponto da coluna de hastes. O leitor que desejar maioresdetalhes quanto ao cálculo destes fatores deve consultar o trabalho de Lea (1991) ou Costa(1995).

Por fim, existem algumas expressões de valor inicial para o sistema ainda emrepouso, em t = 0. São elas:

u(s, 0) = 0,∀s , (3.13)

vr(s, 0) = ∂u

∂t(s, 0) = 0,∀s , (3.14)

Fs(0, 0) = −pd(0)Ar1 , (3.15)

pd(0) = pwh + dfpghb = pb(0) , (3.16)

dFs(s, 0)ds

= −ρrAr~g · ~T (s) , (3.17)


onde Fs(s, t) é a força axial num determinado ponto na coluna de hastes, no momento t(N).

3.1.2 Solução Numérica

Em poços de bombeio mecânico, a coluna de hastes é composta de diversas seçõesde diâmetros e materiais diferentes, onde cada seção consiste de uma ou mais hastesinterconectadas. A Eq. (3.1) pode ser resolvida através de um método de diferençasfinitas explícito; para este fim, cada seção da coluna de hastes deve ser dividida numadeterminada quantidade de pontos, com pontos adjacentes separados por segmentosequidistantes. Analisemos as equações abaixo, tiradas do modelo:

∂vr(s, t)∂t

= 1ρrAr

∂Fs(s, t)∂s

+ ~g · ~T (s)

− µ vr(s, t)|vr(s, t)|

√√√√[~g · ~B(s)]2

+[~g · ~N(s) +

∣∣∣ ~K(s)∣∣∣Fs(s, t)ρrAr

]2

+ fv(s, t)ρrAr

(3.18)

∂vr(s, t)∂s

= 1ErAr

∂Fs(s, t)∂t

(3.19)

onde ~K(s) é o vetor de curvatura num determinado ponto da coluna de hastes.

A Eq. (3.18) é um passo intermediário na dedução de (3.1). Já a Eq. (3.19) é obtidaao derivar em relação ao tempo a expressão resultante da Lei de Hooke, quando aplicadaao deslocamento de um determinado ponto da coluna de hastes em relação a sua posiçãoinicial:

∂u(s, t)∂s

= Fs(s, t)ErAr

=⇒ ∂vr(s, t)∂s

= 1ErAr

∂Fs(s, t)∂t

(3.20)

As Eqs. (3.18) e (3.19) formam, para cada momento no tempo, um sistema deEDPs de primeira ordem que são equivalentes a (3.1), uma EDP de segunda ordem. Pararesolver este sistema numericamente, em primeiro lugar é necessário calcular a velocidadeem cada ponto em que a coluna de hastes foi dividida. Em seguida, calculam-se as forçasaxiais correspondentes.

Sejam vjri a velocidade longitudinal da haste em que um ponto i da coluna de hastes

está localizado num momento j no tempo, em m/s, e F jSi e f

jci a força axial e a força de

atributo de Coulomb por unidade de comprimento, respectivamente, naquele mesmo pontoe momento no tempo, em N. Além disso, sejam ∆s o comprimento do segmento escolhidopara dividir a coluna de hastes em pontos adjacentes, em metros, e ∆t o intervalo detempo entre iterações da simulação, em segundos. Caso se opte por usar o modelo de


atrito viscoso de Gibbs, as velocidades das hastes e as forças axiais para cada ponto i nomomento j serão dadas pelo molde a seguir:

vj+1ri =

∆tρrAr∆s

(F jSi+1 − F

jSi−1) + 2~g · ~Ti∆t

2 + c∆t

+

2f jci∆tρrAr

+ (2− c∆t)vjri

2 + c∆t ,

(3.21)

onde:f j

ci

ρrAr

= −µ vjri∣∣∣vjri

∣∣∣√√√√(~g · ~B)2 +

((~g · ~Ni) +

∣∣∣ ~K∣∣∣ F jsi

ρrAr

)2

(3.22)

eF j+1

si = ErAr

2∆t∆s(vj+1

ri+1 − vj+1ri−1) + F j

si . (3.23)

Por outro lado, se for escolhido o modelo de Lea, a velocidade das hastes seráprovida pelo molde abaixo:

vj+1ri =

∆tρrAr∆s

(F jSi+1 − F

jSi−1) + 2~g · ~Ti∆t

2− ηiUrK1

ρrAr

∆t

+

2f jci∆tρrAr

+(

2 + ηiUrK1

ρrAr

∆t)vj

ri − 2ηiUrK2∆tρrAr

vjf

2− ηiUrK1

ρrAr

∆t.

(3.24)

Ao usar (3.24), (3.22) e (3.23) permanecem válidas para calcular f jci

ρrAr

e F j+1si ,

respectivamente. A dedução dos moldes (3.21) e (3.24) pode ser consultada no ApêndiceA.

Levando em consideração as condições de contorno na solução numérica, a velocidadeda haste polida pode ser obtida de forma direta, baseado nas Eqs. (3.5) e (3.6):

TF jn+1 = A1 cos(ωj∆t) + · · ·+ A6 cos(6ωj∆t)

+B1 sin(ωj∆t) + · · ·+B6 sin(6ωj∆t) ,(3.25)

vjn+1 = TF j

n+1Sω

12 . (3.26)

A força na haste polida é dada por:

F j+1Sn+1 = ErAr

∆t∆s(vj+1

rn+1 − vj+1rn ) + F j

Sn . (3.27)


Para avaliar as condições de contorno no fundo do poço, utiliza-se um processoiterativo para calcular tanto a velocidade do pistão quanto a força axial que atua sobreele, baseado numa estimativa inicial para a velocidade do pistão. Isto é necessário porque,ao analisar as condições no fundo do poço, a força no pistão depende de sua velocidade, ea velocidade do pistão, por sua vez, depende da força. Os moldes necessários são listadosa seguir.

O deslocamento do pistão é dado por:

uj+11 = uj

1 + ∆u1 , (3.28)

onde:∆u1 = vj+1

r1 + vjr1

2 ∆t . (3.29)

A pressão de descarga da bomba de fundo, pj+1d , pode ser obtida ao aplicar as Eqs.

(3.10) e (3.11).

O deslocamento do pistão e a pressão de descarga são necessários para calcular apressão no interior da bomba de fundo através do molde a seguir:

pj+1b =

pd, se pj

b + ∆pb ≥ pj+1d

pjb + ∆pb, se ps < pj

b + ∆pb < pj+1d

ps, se pjb + ∆pb ≤ ps

, (3.30)

onde:∆pb = − ∆u

DbApianc

EtAt

+ kLLb(1− αj) + Lbαj

pjb

. (3.31)

Em seguida, a força na extremidade inferior da coluna de hastes é dada por:

F j+1s1 = (pj+1

d − pj+1b )Ap − pj+1

d Ar1 . (3.32)

Por fim, pode ser usada a expressão abaixo para determinar a velocidade do pistão:

vj+1r1 = vj+1

r2 −(F j+1

S1 − FjS1)∆s

ErAr∆t. (3.33)

3.1.2.1 Colunas de Hastes Combinadas

Nas discussões anteriores, todas as seções de hastes tinham o mesmo diâmetro.Porém, em colunas de hastes combinadas, o diâmetro das seções de hastes pode variar deuma seção para outra, o que exige alguns ajustes no modelo de simulação.

No ponto de transição entre duas seções, a velocidade da última haste na primeiraseção e a da primeira haste na seção seguinte são idênticas. As forças axiais nestas hastes


contíguas, no entanto, são diferentes, por causa de um pequeno efeito de pressão. Seja kuma seção da coluna de hastes, dividida em nk pontos, e k+ 1 a seção seguinte. Além disso,seja pj

k a pressão na extremidade inferior da seção k no momento j no tempo, levandoem consideração o comportamento dinâmico do fluido, em psi. Ao calcular as velocidadesdas hastes e as forças axiais no ponto de transição entre estas seções, valem as seguintesequações:

vj+1r1,k+1 = vj+1

rnk+1,k , (3.34)

F j+1S1,k+1 = F j+1

Snk+1,k − pj+1k+1(Ark+1 − Ark) . (3.35)

Ao discretizar a Eq. (3.19) e aplicar um molde de diferenças atrasadas para a seçãok + 1, as forças no ponto de transição podem ser calculadas conforme segue:

(3.36)F j+1s1,k+1 = ErArk+1

∆t∆sk+1

(vj+1r2,k+1 − v

j+1r1,k+1) + F j

s1,k+1 .

Pode-se repetir este procedimento para a seção k:

(3.37)F j+1snk+1,k = ErArk

∆t∆sk

(vj+1rnk+1,k − v

j+1rnk,k) + F j

snk+1,k .

Ao subtrair (3.36) de (3.37), substituindo em seguida (3.34) e (3.35) no resultado,obtém-se:

vj+1r1,k+1 = vj+1

rnk+1,k =(pj+1

k+1 − pjk+1)(Ark+1 − Ark)ErAt

Ark

∆sk

+ Ark+1

∆sk+1

+

Ark

∆sk

vj+1rnk,k + Ark+1

∆sk+1vj+1

r2,k+1

Ark

∆sk

+ Ark+1

∆sk+1

.

(3.38)

Como pj+1k+1 ≈ pj

k+1, (3.38) é aproximadamente igual a:

vj+1r1,k+1 = vj+1

rnk+1,k =Ark

∆sk

vj+1rnk,k + Ark+1

∆sk+1vj+1

r2,k+1

Ark

∆sk

+ Ark+1

∆sk+1

.(3.39)

Isto conclui todas todas as equações necessárias para levar em consideração colunasde hastes combinadas.


3.1.2.2 Estabilidade Numérica e Critério de Parada

Para garantir a estabilidade numérica, o modelo de simulação exige que

∆s∆t ≥ v (3.40)

conforme recomendado por Laine, Cole e Jennings (1989).

Schafer e Jennings (1987) afirmaram que a utilização de um valor de cerca de 200m (656,17 pés) para ∆s já é suficiente para uma aproximação satisfatória. O modelo deCosta, porém, atribui a ∆s um valor de 20 m (65,62 pés) ou menor, de maneira que todosos pontos de medição ao longo da trajetória do poço sejam utilizados. Supondo que vtenha um valor de cerca de 16300 ft/s (4968,24 m/s), conforme sugerido por Takács (2003),ao aplicar (3.40) conclui-se que ∆t deve ter um valor de 0,004 s ou menor.

Supondo uma frequência de bombeio de 14 ciclos por minuto, cada ciclo de bombeioindividual irá durar cerca de 4,28 segundos. Ao dividir este valor por ∆t, verifica-se quecada ciclo de bombeio deve ser dividido em cerca de 1065 intervalos de tempo. Um númeromaior de intervalos de tempo resultará numa simulação mais precisa.

Ainda de acordo com Schafer e Jennings (1987), para poços com uma frequência debombeio de até 15 ciclos por minuto, as cartas dinamométricas calculadas pelo modelo desimulação chegam a um regime estável após três ciclos de bombeio. Não obstante, Costa(1995) recomenda executar a simulação durante cinco ciclos de bombeio para garantir osmelhores resultados.

3.2 ParalelismoDesde o surgimento dos primeiros microprocessadores, no início da década de 1970,

até por volta do ano de 2004, a grande maioria dos computadores dispunha de umaúnica CPU. Durante esse período, aumentavam de forma significativa a cada ano tanto afrequência de operação das CPUs quanto a quantidade de transistores nelas instalados,resultando em incrementos de desempenho computacional em progressão geométrica.

Por limitações físicas, esta sequência evolucionária infelizmente se encerrou. Chegou-se a um ponto em que novos aumentos de frequência fariam com que os processadoresoperassem em temperaturas excessivamente altas, inutilizando-os. Outra barreira era o con-sumo de energia, que se mostrou proibitivo para frequências muito elevadas. Dessa maneira,a frequência máxima das CPUs modernas estacionou num limite superior praticamenteintransponível (SUTTER, 2004).

Para que o desempenho dos computadores continuasse a avançar, seria necessárioadotar uma abordagem diferente daquela seguida até então. A solução surgiu na formados sistemas paralelos, que exigem que programadores particionem os problemas em


subproblemas menores, passíveis de serem computados de maneira concorrente em diversaslinhas de execução simultâneas. Em tais sistemas, vários nós de processamento operandode forma simultânea e coordenada podem oferecer um melhor desempenho, na medida emque a comunicação entre estes nós ocorra de maneira eficiente e que se adicionem novosnós ao sistema.

A partir de 2005, praticamente todos os novos modelos de computadores de mesa enotebooks passaram a conter múltiplos processadores. Em 2011 essa mudança se consolidoutambém entre os smartphones e tablets, e, em 2012, foi descontinuado o último consolede videogame equipado com um único processador (SUTTER, 2012). No segmento deprocessamento de alto desempenho, onde se destacam os softwares de computação científica,a transição já havia ocorrido bem antes, como forma de contornar as limitações dos sistemasmonoprocessados (também chamados de “seriais” ou “sequenciais”).

Por causa da impossibilidade de avanço dos computadores com um único processa-dor, e também pelo fato de que diversas categorias de aplicações se beneficiam fortementeda utilização de múltiplos processadores, a computação paralela é um caminho sem volta.Dessa maneira, para extrair o máximo desempenho dos sistemas modernos, é essencial queos profissionais envolvidos com desenvolvimento de software se capacitem nos conceitos eferramentas de programação utilizados na construção de soluções paralelas. Essa capa-citação vem se mostrando um grande desafio para a indústria da computação, tanto naacademia quanto no mercado de trabalho, diante da necessidade de revisar e atualizar umenorme corpo de conhecimentos construídos sobre a premissa de sistemas monoprocessados,bem como as ferramentas derivadas destes conhecimentos.

3.2.1 A Lei de Amdahl

O desenvolvimento de soluções paralelas requer uma visão fundamentalmentediferente daquela empregada em soluções seriais, tendo como aspecto central o conceito dasubdivisão de problemas: dado um determinado problema, como subdividi-lo em problemasmenores, de tal maneira que estes subproblemas possam ser executados simultaneamentecom a melhor eficiência possível? Quando se trata de paralelismo, não há soluções genéricase definitivas, e cada problema exige sua própria análise.

Comecemos definindo o termo “speedup”, que quantifica o ganho de desempenhoentre dois sistemas diferentes que executam a mesma tarefa:

s = t2t1, (3.41)

onde:s = speedup;t1 = tempo de execução da tarefa no primeiro sistema;t2 = tempo de execução da tarefa no segundo sistema.


No contexto do paralelismo, o speedup paralelo é o ganho de desempenho de umsistema paralelo quando comparado à sua versão serial:

sp = tstp, (3.42)

onde:sp = speedup paralelo;ts = tempo de execução da tarefa no sistema serial;tp = tempo de execução da tarefa no sistema paralelo.

Determinados problemas se beneficiam mais fortemente da paralelização, e portantoapresentam um speedup elevado quando comparados às respectivas versões seriais. Háoutros problemas, no entanto, em que o speedup é bastante reduzido, o que indica quea abordagem de paralelismo escolhida é ineficiente (e potencialmente indesejável). Ospeedup máximo na paralelização de um problema foi expresso por Amdahl (1967)conforme segue:

S = 1rs + rp

n

,

com rs + rp = 1 ,(3.43)

onde:S = speedup máximo;rs = fração estritamente serial do problema (porcentagem variando de 0 a 1, onde 1significa 100%);rp = fração paralelizável do problema (porcentagem variando de 0 a 1, onde 1 significa100%);n = número de processadores disponíveis.

Quanto mais extensa for a fração do problema que se mostre impossível de paralelizar— ou seja, quanto mais próximo rs estiver do valor 1 —, mais o speedup máximo S tenderáao valor 1, o que indica que a paralelização do problema é inviável. Por outro lado, numcenário de paralelismo ideal, rs será igual a 0 e rp será igual a 1, e portanto o speedupmáximo será:

S = 1rs + rp

n

= 1

0 + 1n

= n . (3.44)

O speedup máximo expresso por Amdahl é excepcionalmente difícil de obter, porcausa de uma variedade de fatores limitantes. Entre tais fatores, pode-se citar o esforçonecessário para coordenar os vários processadores do sistema paralelo, o que impede autilização do conjunto total a plena capacidade. Outro fator a considerar é a afinidade dememória, pois o custo de leitura/gravação dos vários segmentos da memória do sistema


pode não ser homogêneo para todos os processadores disponíveis. Pode-se mencionar aindao problema do falso compartilhamento, no qual uma operação de gravação executada porum processador num determinado segmento de memória pode invalidar as linhas de cachede outro processador que esteja utilizando aquele mesmo segmento, forçando que esteúltimo realize leituras a partir da memória do sistema, e não do cache, com consequenteperda de desempenho.

Outro fator a considerar ao trabalhar com paralelismo é a eficiência paralela,dada por:

ep = sp

n, (3.45)

onde:ep = eficiência paralela;sp = speedup paralelo;n = número de processadores utilizados na obtenção de sp.

Conforme vimos em (3.44), o speedup paralelo máximo é igual a n, e portanto aeficiência paralela máxima será:

ep = sp

n= n

n= 1 . (3.46)

Quanto mais próxima a eficiência paralela estiver do valor 1 conforme se aumentao número de processadores, melhor será a escalabilidade da implementação paralelaem questão. O speedup paralelo e a eficiência paralela podem variar tanto em função donúmero de processadores quanto do tamanho do problema.

35

4 Metodologia

Descrevemos neste capítulo a forma como implementamos o modelo de Costa (1995),com comentários sobre estratégias de otimização e paralelização para o problema. Fazemosa análise da complexidade computacional da simulação em função de seus parâmetros deentrada, tanto em forma serial quanto paralela. Discutimos ainda algumas contribuiçõesao trabalho original de Costa, com consequências de estabilidade, precisão e desempenhopara a implementação exposta no presente trabalho.

4.1 Implementação paralelaO modelo do capítulo anterior foi implementado em linguagem C, com uma versão

serial e outra paralela. A versão paralela, além dos recursos padrão da linguagem Cempregados na versão serial, utiliza tecnologia OpenMP (OPENMP. . . , 2015). O programaestá estruturado nas seções a seguir:

1. Preparação da simulação: leitura de arquivos de configuração, atribuição de variáveise cálculos iniciais;

2. Percurso através das diversas configurações de poço a processar, onde cada iteraçãoexecuta os passos abaixo:

a) Laço Principal de Simulação;

b) Extração de parâmetros de operação a partir das cartas dinamométricas calcu-ladas no Laço Principal de Simulação.

O Laço Principal de Simulação, indubitavelmente a seção de maior exigênciacomputacional do programa, está descrito na listagem do Algoritmo 1, e será executadopara cada configuração de poço que se deseja simular. Os dados de configuração de poços,por sua vez, são ligeiramente diferentes de acordo com o modelo de atrito viscoso utilizado.No modelo de Gibbs, uma configuração de poço deve conter um coeficiente de Coulomb eum fator de amortecimento, o que leva o percurso através das configurações de poços ater a estrutura mostrada no Algoritmo 2. O modelo de Lea também exige um coeficientede Coulomb, mas o fato de empregar a viscosidade dinâmica do fluido de forma diretadispensa a utilização de um fator de amortecimento, fazendo com que o percurso atravésdas configurações de poços tenha o arranjo mostrado no Algoritmo 3.

Dada uma configuração de poço baseada num poço real em funcionamento, a quali-dade dos resultados da simulação pode ser verificada ao comparar as cartas dinamométricas

Capítulo 4. Metodologia 36

para ciclo = 1 até ciclos de bombeio façapara intervalo = 1 até intervalos de tempo por ciclo faça

Parte 1: Calcular velocidade na haste polidaParte 2: Calcular velocidades para os pontos que não estiverem localizados nastransições entre seções de hastes na coluna

Parte 3: Calcular velocidades para os pontos localizados nas transições entre seçõesde hastes na coluna

Parte 4: Calcular força no pistão de subsuperfície (extremidade inferior da colunade hastes)

Parte 5: Calcular forças para os pontos que não estiverem localizados nastransições entre seções de hastes na coluna

Parte 6: Calcular força na haste polidaParte 7: Calcular forças para os pontos localizados nas transições entre seções dehastes na coluna

fimfim

Algoritmo 1: Laço Principal de Simulação

para cada coeficiente de Coulomb na faixa desejada façapara cada fator de amortecimento na faixa desejada faça

Executar Laço Principal de Simulação (vide Algoritmo 1)Extrair parâmetros de simulação a partir das cartas dinamométricas calculadas

fimfim

Algoritmo 2: Percurso das configurações de poços para o modelo de AV de Gibbs

para cada coeficiente de Coulomb na faixa desejada façaExecutar Laço Principal de Simulação (vide Algoritmo 1)Extrair parâmetros de operação a partir das cartas dinamométricas calculadas

fimAlgoritmo 3: Percurso das configurações de poços para o modelo de AV de Lea

e parâmetros de operação resultantes do cálculo aos valores correspondentes adquiridos emcampo através de sensores de automação ou meios diversos. É necessário testar faixas devalores tanto para o coeficiente de Coulomb quanto para o fator de amortecimento porque,infelizmente, nenhum deles tem valores específicos que sejam adequados para todos oscampos de petróleo. Recomenda-se comparar os resultados das simulações aos dados deoperação reais de poços em funcionamento de forma a determinar os melhores valores paracada cenário.

Costa (1995) validou a precisão de seu modelo ao obter boa concordância entre osvalores previstos e medidos para uma variedade de poços. Portanto, partimos do princípiode que o modelo é confiável, o que direcionou as pesquisas à melhoria do desempenho domodelo, em oposição à aferição de sua qualidade. Dito isso, o Apêndice B traz comparações


entre os dados reais e simulados de alguns poços de petróleo existentes, diferentes daquelesutilizados por Costa, de forma a ilustrar o funcionamento do modelo.

4.1.1 Considerações de Otimização e Paralelização

O Laço Principal de Simulação descrito no Algoritmo 1 oferece diversas oportunida-des imediatas de otimização. Por exemplo, muitas das equações utilizadas nas Partes de 1a 7 contêm fatores que não se modificam de uma iteração para outra. Estes valores podemser computados antes do Laço Principal e reaproveitados seguidamente. Isto poupa tempode processamento, enquanto causa como efeito colateral apenas um aumento modesto nautilização de memória.

A implementação original de Costa em linguagem Fortran utilizava principalmentematrizes bidimensionais. Ao reimplementar o programa em linguagem C, optou-se pelautilização de vetores unidimensionais. Esta mudança na representação dos dados permitecalcular de forma antecipada o endereço correspondente à coluna inicial de cada linha dedados de uma matriz bidimensional armazenada como um vetor unidimensional, com oobjetivo de guardar tais endereços em tabelas de consulta. As tabelas de consulta, por suavez, permitem determinar com rapidez o endereço que um elemento localizado numa linhae coluna específica ocupa dentro de um determinado vetor unidimensional.

Outro benefício decorrente desta mudança de representação vem da facilidade deredirecionar vetores inteiros através de algumas poucas operações de permuta de ponteiros.Ao processar um intervalo de tempo de um ciclo de bombeio, é preciso ter acesso às matrizesde velocidade e força dos intervalos de tempo atual e anterior. Ao concluir o processamentodo intervalo de tempo, é preciso copiar o conteúdo das matrizes de velocidade e forçado intervalo de tempo atual para as matrizes correspondentes ao intervalo anterior, deforma que possam ser utilizadas no intervalo de tempo seguinte. É possível realizar estatarefa copiando blocos de memória, mas a cópia seria tão mais demorada quanto maioresfossem as matrizes envolvidas. Um alternativa mais eficiente é utilizar ponteiros parasimplesmente permutar os endereços de memória para os quais apontam as matrizes develocidade e força dos intervalos de tempo atual e anterior. Dessa maneira, o intervalode tempo seguinte irá “enxergar” nas matrizes de velocidade e força correspondentes aointervalo de tempo anterior os dados mais recentes calculados no intervalo de tempo atual,a um custo constante de uma simples permuta de ponteiros.

Além de otimizações tais como as descritas acima, o processamento paralelo temum papel crucial na obtenção de melhorias de desempenho adicionais. Para que se possaparalelizar o Laço Principal de Simulação, é necessário identificar quais de suas partespodem ser subdivididas em blocos menores, passíveis de serem executados de formaconcorrente. Ao consultar no modelo as equações exigidas por cada parte, verificamos queas Partes 1 e 6 têm um custo computacional bastante modesto. Assim, elas serão deixadas


em sua forma sequencial, sem modificações.

Por outro lado, a Parte 2, que usa (3.21) e (3.24), e a Parte 5, que usa (3.23), sãobastante similares e boas candidatas à execução paralela. Ambas as partes exigem percorrerquase todos os pontos em que a coluna de hastes foi dividida, exceto aqueles localizados nastransições entre seções de hastes, para calcular a velocidade e a força em cada ponto. Comoestes cálculos não precisam de informações de outros pontos, podem ser potencialmenteexecutados por linhas de execução independentes com resultados favoráveis.

A Parte 3, que usa (3.39), e a Parte 7, que usa (3.35) e (3.36), também parecemadequadas à paralelização, mas os ganhos que se poderiam obter ao paralelizar estas partesnão devem compensar o esforço de programação envolvido. Costa afirma que a maioria dospoços de bombeio mecânico têm entre uma e quatro seções de hastes; tendo isto em mente,as Partes 3 e 7 devem percorrer todas as seções da coluna de hastes, exceto a última,calculando para cada seção as equações citadas, que consistem apenas de alguns poucoscálculos simples. Como a quantidade de seções a analisar deve ser bastante limitada, e oesforço computacional exigido em cada seção bastante modesto, executar estas partes emlinhas de execução separadas envolveria mais “overhead” 1 do que processamento efetivo,visto que há muito pouco trabalho a realizar. Portanto, é melhor deixar estas partes emsua forma sequencial.

A Parte 4, que é baseada num teste iterativo de convergência, tem de ser examinadaà luz de dois cenários distintos. Caso o engenheiro opte pela utilização do modelo de AVde Gibbs, a pressão de descarga da bomba de fundo não necessita ser atualizada a cadaiteração do Laço Principal de Simulação, como se pode ver em (3.10). Portanto, a Parte 4não terá cálculos a paralelizar. Porém, caso o AV seja computado através do modelo deLea, o cálculo da pressão de descarga é muito mais lento, pois exige a análise de todos ospontos em que a coluna de hastes foi dividida a cada ocorrência do teste iterativo, conformemostrado em (3.11). Felizmente, esta equação permite que cada ponto seja avaliado deforma independente dos demais. O trabalho a executar para cada ponto é relativamentebaixo, mas paralelizá-lo é potencialmente benéfico.

4.1.2 Análise de Complexidade

Seja c o número de ciclos de bombeio a simular, e i o número de intervalos detempo por ciclo. As Partes numeradas de 1 a 7 no Algoritmo 1 serão executadas c · i vezes.

As Partes 1 e 6, conforme discutido anteriormente, são de baixa exigência compu-tacional. As características de desempenho das equações utilizadas por estas partes sãoessencialmente fixas, com listas de parâmetros cujo tamanho independe das dimensões dasimulação a realizar, e cuja escala numérica não penaliza a execução em nenhum cenário1 Procedimentos necessários à preparação, apoio à execução e encerramento de tarefas específicas.


identificado. Além disso, nossos experimentos mostram que seu custo computacional épraticamente constante. Por estes motivos, a complexidade de ambas será consideradacomo O(1).

Seja e o número total de pontos em que a coluna de hastes foi dividida. As Partes2 e 3, combinadas, exigem percorrer todos estes pontos; mais especificamente, a Parte 2exige percorrer todos os pontos que não estiverem localizados nas transições entre seçõesde hastes, enquanto a Parte 3 exige percorrer os pontos localizados nas transições entreseções de hastes. Tal padrão ocorre também nas Partes 5 e 7, combinadas. Ainda quantoà necessidade de percorrer todos os pontos, ela ocorre novamente na Parte 4, caso oengenheiro escolha utilizar o modelo de AV de Lea. A Parte 4 se baseia num teste iterativode convergência; seja j o número de iterações necessárias para que a convergência ocorra.Dessa maneira, todos os pontos em que a coluna de hastes foi dividida serão percorridos(1 + 1 + j) vezes a cada intervalo de tempo, i.e. (execuções das Partes 2 e 3, combinadas)+ (execuções das Partes 5 e 7, combinadas) + (execuções da Parte 4). Em consequência,estas partes terão uma complexidade total de O(1 + 1 + j) ·O(e). Em nossos experimentos,na maioria das vezes a Parte 4 convergiu em 4 iterações. Dessa maneira, supondo que 4é um valor frequente e representativo para j, no caso médio todos os pontos em que acoluna de hastes foi dividida deverão ser percorridos 6 vezes para cada intervalo de tempo.

Combinando os resultados anteriores, a complexidade serial do Laço Principal deSimulação pode ser expressa por:

O(c · i) · (O(1) +O(1) + (O(1 + 1 + j) ·O(e))) =

O(c · i) · (O(2) + (O(2 + j) ·O(e))) .(4.1)

Supondo que o custo das Partes 1 e 6, representado por O(2), seja desprezível, eutilizando um valor médio para j, a complexidade serial média é dada por:

O(c · i) · (O(2) + (O(2 + 4) ·O(e))) =

O(c · i) · (O(6) ·O(e)) =

O(c · i) ·O(e) .

(4.2)

A análise de e revela maiores detalhes sobre a complexidade serial. O valor de e éinversamente proporcional ao comprimento do segmento escolhido para dividir a coluna dehastes, i.e. ∆s na Eq. (3.40). ∆s, por sua vez, é diretamente proporcional a ∆t por aquelamesma equação. E, já que ∆t é inversamente proporcional ao número de intervalos detempo por ciclo de bombeio, i, pode-se finalmente concluir que e é diretamente proporcionala i. Em consequência, como (4.2) multiplica O(c · i) por O(e), então (4.2) deverá exibirum comportamento quadrático para valores fixos de c.


4.1.2.1 Complexidade Paralela

As iterações dos ciclos de bombeio e intervalos de tempo por ciclo são estritamentesequenciais, dado que os resultados de cada laço dependem dos resultados do laço anterior.Portanto, o termo O(c · i) não pode ser acelerado via execução paralela.

Não se tentará paralelizar as Partes 1 e 6, e a complexidade paralela para ambasserá considerada idêntica à complexidade serial, i.e. O(1).

As Partes 2 e 3 combinadas, assim como a Parte 4, e ainda as Partes 5 e 7combinadas, exigem percorrer todos os pontos em que a coluna de hastes foi dividida.Infelizmente, como a Parte 4 depende dos resultados das Partes 2 e 3, e as Partes 5 e 7dependem dos resultados da Parte 4, estes percursos não podem ser executados de formaparalela entre si. Portanto, a complexidade paralela destas Partes será igual à complexidadeserial, O(1 + 1 + j).

Embora os percursos mencionados no parágrafo anterior não possam ser executadosde forma paralela entre si, felizmente o processamento interno em cada um deles podeser realizado em paralelo. Seja p o número de processadores disponíveis. A complexidadeparalela de cada percurso será então dada por 2 :

O

(e

p

), com p em O(e) . (4.3)

Ao multiplicar esta expressão pelos resultados anteriores, a complexidade paralelado Laço Principal é expressa por:

O(c · i) ·(O(1) +O(1) +

(O(1 + 1 + j) ·O

(e

p

)))=

O(c · i) ·(O(2) +

(O(2 + j) ·O

(e

p

))),

com p em O(e) .

(4.4)

Supondo que o custo das Partes 1 e 6, representado por O(2), seja desprezível, eutilizando um valor frequente e representativo para j, a complexidade paralela média édada por:

O(c · i) ·(O(2) +

(O(2 + 4) ·O

(e

p

)))=

O(c · i) ·(O(6) ·O

(e

p

))=

O(c · i) ·O(e

p

),

com p em O(e) .

(4.5)

2 Na falta de uma notação “O-grande” para o comportamento de algoritmos paralelos, a notação descritapor Ostrovsky (2008) será utilizada.


A Eq. (4.5) retém o comportamento quadrático observado em (4.2) para valoresfixos de c, mas numa escala menor. Além disso, a paralelização deverá se tornar maiseficiente conforme cresça o valor de e.

4.1.3 Contribuições ao Trabalho Original

Ao revisar o modelo de Costa, identificamos alguns pontos passíveis de correçãoe algumas oportunidades de melhoria. Primeiramente, as Eqs. (3.4) e (3.24) tinham umerro tipográfico, no qual era mostrado um sinal de menos em vez de um sinal de mais, oque poderia levar a instabilidade e resultados incorretos ao usar o modelo de AV de Lea.Neste trabalho, as equações são apresentadas em sua forma correta.

Em segundo lugar, no trabalho original, numa discussão sobre a velocidade médiado fluido numa determinada seção da coluna de hastes, é mostrada a equação a seguir(página 16, Eq. (3.38)):

v̄fk =

vp(Ap − Ark)At − Ark

, se vp ≥ 0 ,

−vpArk

At − Ark

, caso contrário., (4.6)

onde:v̄fk(t) = velocidade média do fluido na seção k da coluna de hastes, no momento t (m/s);vp(t) = velocidade do pistão no momento t (m/s);Ap = área da seção transversal do pistão (m2);Ark = área da seção transversal da seção k da coluna de hastes (m2);At = área da seção transversal da coluna de produção (m2).

Esta equação pressupõe duas hipóteses simplificadoras. A primeira delas é que o fluidoseja incompressível, enquanto a segunda é que a velocidade da coluna de hastes sejaaproximadamente igual à velocidade do pistão. É possível aumentar a precisão destaequação ao levar em consideração a pressão no interior da bomba de fundo.

A Eq. (4.6) assume que quando vp ≥ 0, i.e. o pistão estiver se deslocando para cima,a válvula de passeio deverá estar fechada e deverá haver fluido se deslocando em direçãoà cabeça do poço; e que quando vp < 0, i.e. o pistão estiver se deslocando para baixo, aválvula de passeio deverá estar aberta e não deverá haver fluido se deslocando em direçãoà cabeça do poço. No entanto, a válvula de passeio estará fechada somente se a pressão nointerior da bomba de fundo for menor do que a pressão de descarga, independentemente dopistão estar se deslocando para cima ou para baixo. Assim, em vez de inferir o estado daválvula de passeio a partir da velocidade do pistão, o mais correto seria fazê-lo verificandoas pressões mencionadas anteriormente. Se a pressão no interior da bomba de fundo, pb(t),for menor do que a pressão de descarga, i.e. a pressão na extremidade inferior da seção


de hastes em que o pistão está localizado, p1(t), então a válvula de passeio deverá estarfechada; caso contrário, deverá estar aberta.

Levando em conta as considerações sobre a pressão no interior da bomba de fundo,(4.6) deveria ser reescrita como:

v̄fk =

vp(Ap − Ark)At − Ark

, se pb(t) < p1(t) ,

−vpArk

At − Ark

, caso contrário., (4.7)

onde:pb(t) = pressão no interior da bomba de fundo, no momento t (psi)p1(t) = pressão de descarga, i.e. a pressão na extremidade inferior da seção de hastes emque o pistão está localizado, no momento t (psi).

Mudando o foco para a implementação de software, podemos apontar duas diferençasna avaliação das condições de contorno no fundo do poço, i.e. na Parte 4 do Laço Principalde Simulação. A primeira diferença está relacionada ao fato de que, no processo iterativopara calcular a velocidade do pistão e a força axial atuando sobre ele a cada passo dasimulação, é necessário fazer uma estimativa inicial para o valor da velocidade do pistão. Aimplementação de Costa (1995) se limitava a utilizar um valor inicial de 0 a cada passo. Naimplementação descrita neste trabalho, um valor inicial de 0 é utilizado apenas no primeiropasso, e, nos passos posteriores, o valor inicial será a velocidade do pistão encontrada nopasso anterior. Com esta modificação, a convergência se torna ligeiramente mais rápida doque na implementação original. Como estes passos são estritamente sequenciais, melhoriasem suas características de desempenho diminuirão o tempo exigido pela porção estritamentesequencial do programa; por conseguinte, a porção paralelizável do programa aumentaráde tamanho, o que pode trazer ganhos consideráveis de desempenho paralelo conforme seexpandam os recursos paralelos à disposição.

A segunda diferença é específica à implementação paralela. Na implementação serialde Costa, dada uma estimativa inicial para a velocidade do pistão, um pequeno fator eraadicionado a ela. E então, a força correspondente àquela velocidade era calculada utilizando(3.32), sendo computada em seguida através de (3.33) a velocidade correspondente à forçacalculada anteriormente. Caso a velocidade computada para o pistão fosse suficientementepróxima à estimativa inicial acrescida do pequeno fator, atingindo portanto a convergência,a simulação prosseguiria utilizando os valores calculados tanto para a velocidade do pistãoquanto para a força axial atuando sobre ele. Porém, caso a velocidade computada para opistão e a estimativa inicial estivessem muito distantes, o processo seria repetido, desta vezsubtraindo um pequeno fator da estimativa inicial e recalculando a força e a velocidadeatravés de (3.32) e (3.33). Se esta segunda tentativa ainda não conseguisse atingir aconvergência, a simulação entrava num laço no qual os dois valores iniciais anteriores para


a velocidade seriam combinados utilizando o método da secante, e o valor resultante destacombinação seria utilizado como a velocidade inicial para recalcular a força e a velocidadevia (3.32) e (3.33), até que um valor aceitável fosse encontrado.

Na versão paralela da implementação descrita neste trabalho, dada uma estimativainicial para a velocidade do pistão, são derivados dois valores iniciais, um deles atravésdo acréscimo de um pequeno fator à estimativa inicial e o outro através da subtraçãodo mesmo fator. Então, em duas linhas de execução independentes, são computadas aforça e a velocidade correspondentes aos valores iniciais calculados anteriormente. Casoqualquer das linhas de execução atinja a convergência, a simulação prossegue utilizandoos valores calculados na primeira linha de execução que teve êxito. Caso a convergêncianão seja atingida, a simulação entra num laço em que os dois valores iniciais anterioressão combinados através do método da secante, e o valor resultante desta combinação éutilizado numa única linha de execução como velocidade inicial para recalcular a força e avelocidade via (3.32) e (3.33), até que um valor aceitável seja encontrado.

Teoricamente, poderiam ser usadas mais de duas linhas de execução quando houvessemais de dois processadores disponíveis. Porém, em nossos experimentos, na maioria doscasos a convergência foi atingida em duas ou três iterações do laço de força/velocidade,sendo utilizadas na primeira iteração múltiplas linhas de execução e nas demais iteraçõesuma única linha, mesmo que se empregassem tantas linhas de execução quanto houvesseprocessadores disponíveis. Portanto, optou-se por utilizar apenas duas linhas de execuçãona primeira checagem de convergência, o que ainda resultou em ganhos líquidos dedesempenho ao mesmo tempo em que evitou processamento desnecessário.

A técnica descrita acima é utilizada apenas para o modelo de AV de Lea, visto que(3.32) exige um valor de pressão de descarga dado por (3.11). Esta última equação precisarealizar um percurso bastante custoso ao longo de todos os pontos em que a coluna dehastes foi dividida, o que assegura que as linhas de execução envolvidas tenham trabalhosuficiente a executar sem incorrer em “overhead” excessivo. Os ganhos de desempenhoresultantes desta técnica estão descritos no Capítulo 5.

Em contrapartida, ao utilizar o modelo de AV de Gibbs, (3.32) tem uma exigênciacomputacional modesta, visto que exige um valor de pressão de descarga dado por (3.10).Como o resultado desta última equação permanece constante ao longo de toda a simulação,pode ser calculado uma única vez e reaproveitado quantas vezes for necessário. Nestecenário, computar (3.32) em linhas de execução separadas resultaria em mais “overhead”do que processamento efetivo, prejudicando o desempenho global da simulação. Portanto,para o modelo de AV de Gibbs, esta parte foi deixada em forma sequencial.

44

5 Simulações e Resultados

As simulações foram realizadas num sistema GNU/Linux executando o UbuntuLinux versão 14.04.4 64-bits, equipado com uma CPU Intel Core i7-980 contendo seisnúcleos físicos e duas linhas de execução de hardware por núcleo. Todos os núcleosfuncionaram a uma velocidade de 3,33 GHz. Tanto a versão serial quanto a paralela foramcompiladas utilizando o compilador GCC versão 5.3.0.

Tanto a versão paralela quanto a serial foram executadas sob um cenário de cincociclos de bombeio de uma única configuração de poço, com um número crescente deintervalos de tempo por ciclo, utilizando tanto o modelo de AV de Gibbs quanto o de Lea.O número de intervalos de tempo é inversamente proporcional ao comprimento do segmentoutilizado para dividir a coluna de hastes, i.e. ∆s. Um maior número de intervalos de tempo,i.e. um valor menor para ∆s, resultará num maior número de pontos ao longo da coluna dehastes. Para cada um destes pontos, serão calculadas forças axiais, velocidades longitudinaise forças de atrito de Coulomb para cada intervalo, o que fornecerá informações sobre ascondições de operação às quais a coluna de hastes estará sujeita. Nossos experimentosmostram que as cartas dinamométricas resultantes, bem como os parâmetros de operaçãoderivados a partir delas, só se beneficiam de valores menores de ∆s até um certo ponto. Aousar o modelo de Gibbs, não se verificaram diferenças significativas nestes resultados paravalores de ∆s menores do que 0,75 metros. Já para o modelo de Lea, o valor limite para∆s abaixo do qual não se observaram variações razoáveis foi de 0,30 metros. Dito isso,caso o interesse do engenheiro seja analisar as condições de operação afetando segmentosespecíficos da coluna de hastes ao longo do tempo, recomenda-se trabalhar com valoresmenores de ∆s. Isto deve ser especialmente desejável em se tratando das extremidadessuperiores e inferiores da coluna de hastes, que, por causa de suas funções e característicasespecializadas, exercem uma forte influência sobre o comportamento global da coluna dehastes, apesar de serem normalmente muito mais curtas do que as demais seções de hastes.Porém, caso o engenheiro deseje apenas obter parâmetros de operação a partir das cartasdinamométricas calculadas, valores maiores de ∆s devem ser suficientes.

Os dados de configuração de poço foram baseados num poço direcional terrestresituado na bacia Potiguar do Nordeste brasileiro, com uma profundidade vertical de 696,5metros, uma profundidade medida de 876,3 metros e uma coluna de hastes com duasseções, onde cada seção media 464,8 e 411,5 metros, respectivamente. A listagem completados dados da configuração de poço pode ser consultada na Tabela 2. Os dados do perfil dopoço — i.e. os segmentos da trajetória direcional da coluna de produção – estão listadosna Tabela 3.

Capítulo 5. Simulações e Resultados 45

Tabela 2 – Dados da configuração de poço.

Dado Valor

Ciclos de bombeio 5Coeficiente de atrito de Coulomb (cD) 0,05Fator de amortecimento (µ) 0,15

BSW do fluido 5%Coluna de produção ancorada? NãoCurso da unidade de bombeio 44 inDensidade relativa da água produzida 1,0Densidade relativa do gás 0,7Diâmetro da primeira seção de hastes 6/8 inDiâmetro da segunda seção de hastes 5/8 inDiâmetro do pistão 1,75 inDiâmetro externo da coluna de produção 2,375 inDiâmetro interno da coluna de produção 1,992 inEficiência de separação de gás no fundo do poço 0%Espaço morto 1 péFrequência de bombeio 10 ciclos/minutoGradiente geotérmico 0,01◦F/mGravidade API do óleo 39◦

Módulo de elasticidade das hastes 29,6 · 106 psiMódulo de elasticidade da coluna de produção 29,6 · 106 psiNível dinâmico 876 mPacker instalado? SimPressão de separação 14,69 psiPressão de sucção 60 psiPressão na cabeça 28,44 psiPressão no anular 0 psiProfundidade da bomba de fundo 876,3 mQuantidade de hastes na primeira seção 54Quantidade de hastes na segunda seção 61Razão gás/óleo 1Temperatura ambiente na superfície 37,77◦CTemperatura de separação 26,66◦CViscosidade dinâmica do óleo a 40◦C 1 cPViscosidade dinâmica do óleo a 80◦C 2 cP

Para cada número de intervalos de tempo, a simulação foi executada cinco vezes,calculando-se a média dos tempos de execução resultantes. Os resultados para um segmentode 14,67 metros de comprimento, que corresponde a 2000 intervalos de tempo por ciclo debombeio, estão incluídos pelo fato de esta ser a configuração padrão na implementaçãooriginal de Costa. Os resultados das execuções estão listados nas Tabelas de 4 a 6, etraçados nas Figuras de 3 a 8.

Melhorias substanciais foram alcançadas em relação aos resultados originais obtidosem 1995. Na implementação original de Costa (1995), na qual o valor padrão para ∆s era de14,67 metros, os tempos de execução para os modelos de AV de Gibbs e Lea eram de cerca


Tabela 3 – Perfil do poço.

Profund. Inclinação Azimute Profund. Afastamento Afastamentomedida (graus) (graus) vertical Norte Este(m) (m) (m) (m)

0 0 0 0 0 028 3,25 265 27,98 0,07 -0,7938 5 258 37,96 0,18 -1,547 6,25 247 46,92 0,44 -2,3456 7,5 244 55,85 0,89 -3,3274 10,25 244 73,63 2,1 -5,8283 12,25 243 82,46 2,89 -7,3992 14 245 91,22 3,78 -9,23102 15 245 100,91 4,84 -11,49111 16,25 248 109,57 5,81 -13,72120 17,75 251 118,18 6,73 -16,18156 21,75 256 152,05 10,18 -27,84202 25 256 194,27 14,6 -45,55257 29,75 259 243,1 20,07 -70,23314 34,25 259 291,43 25,83 -99,87381 40 258 344,82 33,89 -139,48418 40 263 373,17 37,81 -162,93447 41 263 395,22 40,11 -181,62475 42,25 264 416,15 42,21 -200,1494 42,5 261 430,19 43,89 -212,79504 41,25 259 437,63 45,04 -219,37514 42 257 445,11 46,43 -225,86560 42,25 250 479,22 55,18 -255,43616 45 249 519,76 68,71 -291,61691 47,5 252 571,62 86,79 -342,67729 48,5 252 597,04 95,52 -369,53823 47,25 256 660,09 114,73 -436,53917 46,25 258 724,5 130,13 -503,24999 47 258 780,81 142,53 -561,541005 47 258 784,9 143,44 -565,83

Tabela 4 – Tempos de execução - Modelo de AV de Gibbs

Tempo paralelo (s)∆s (m) Tempo serial (s) 2 linhas 4 linhas 6 linhas

14,67 0,023 0,036 0,040 0,0480,5 4,691 3,278 2,241 2,0350,4 7,142 4,787 3,206 2,8110,3 12,679 7,855 5,105 4,6050,2 27,743 16,420 10,135 8,3910,1 111,499 62,702 34,759 27,4310,05 452,839 254,084 140,344 101,434

de 104,3 e 108,9 segundos, respectivamente. Ao examinar as Tabelas 4 e 5 para o mesmo


Tabela 5 – Tempos de execução - Modelo de AV de Lea

Tempo paralelo (s)∆s (m) Tempo serial (s) 2 linhas 4 linhas 6 linhas

14,67 0,026 0,045 0,052 0,0680,5 6,942 5,215 4,322 4,2850,4 10,529 7,725 6,292 6,1750,3 18,892 12,883 10,476 10,0700,2 41,472 27,901 21,177 20,3890,1 167,010 106,627 79,552 73,0640,05 668,219 425,046 311,995 281,050

0,00

0,05

0,10

0,15

0,20

0,25

0,30

0,35

0,40

0,45

0,50

0

200

400

∆s [m]

Tempo

[s]

Serial2 linhas4 linhas6 linhas

Figura 3 – Tempos de execução - Modelode AV de Gibbs

0,00

0,05

0,10

0,15

0,20

0,25

0,30

0,35

0,40

0,45

0,50

0

200

400

600

∆s [m]

Tempo

[s]

Serial2 linhas4 linhas6 linhas

Figura 4 – Tempos de execução - Modelode AV de Lea

Tabela 6 – Speedup / Eficiência - Modelos de AV de Gibbs e Lea

Speedup / Eficiência - Speedup / Eficiência -Modelo de Gibbs Modelo de Lea

∆s (m) 2 linhas 4 linhas 6 linhas 2 linhas 4 linhas 6 linhas

14,67 0,65 / 0,32 0,58 / 0,14 0,48 / 0,08 0,58 / 0,29 0,50 / 0,13 0,38 / 0,060,5 1,43 / 0,72 2,09 / 0,52 2,31 / 0,38 1,33 / 0,67 1,61 / 0,40 1,62 / 0,270,4 1,49 / 0,75 2,23 / 0,56 2,54 / 0,42 1,36 / 0,68 1,67 / 0,42 1,71 / 0,280,3 1,61 / 0,81 2,48 / 0,62 2,75 / 0,46 1,47 / 0,73 1,80 / 0,45 1,88 / 0,310,2 1,69 / 0,84 2,74 / 0,68 3,31 / 0,55 1,49 / 0,74 1,96 / 0,49 2,03 / 0,340,1 1,78 / 0,89 3,21 / 0,80 4,06 / 0,68 1,57 / 0,78 2,10 / 0,52 2,29 / 0,380,05 1,78 / 0,89 3,23 / 0,81 4,46 / 0,74 1,57 / 0,79 2,14 / 0,54 2,38 / 0,40

valor de ∆s, verifica-se que os tempos de execução para a versão serial da implementaçãoproduzida neste trabalho são de, respectivamente, 0,023 e 0,026 segundos para cada modelode AV, i.e. 4534 e 4188 vezes mais rápidos do que inicialmente. Caso se suponha que afrequência de operação foi o fator dominante e principal para as melhorias de desempenhode uma geração de processadores para outra (e também abstraindo avanços tais comocoprocessadores matemáticos mais eficientes, conjuntos de instruções mais especializados ecompiladores mais inteligentes), para que a implementação de Costa tivesse um desempenho


2 4 60

2

4

6

Linhas de execução

Speedu

p

∆s = 0, 5m ∆s = 0, 4m

∆s = 0, 3m ∆s = 0, 2m

∆s = 0, 1m ∆s = 0, 05m

Figura 5 – Speedup - Modelo de AV deGibbs

2 4 60

1

2

3

4

Linhas de execução

Speedu

p

∆s = 0, 5m ∆s = 0, 4m

∆s = 0, 3m ∆s = 0, 2m

∆s = 0, 1m ∆s = 0, 05m

Figura 6 – Speedup - Modelo de AV de Lea

0,00

0,05

0,10

0,15

0,20

0,25

0,30

0,35

0,40

0,45

0,50

0

0,5

1

∆s [m]

Eficiên

cia

2 linhas4 linhas6 linhas

Figura 7 – Eficiência - Modelo de AV deGibbs

0,00

0,05

0,10

0,15

0,20

0,25

0,30

0,35

0,40

0,45

0,50

0

0,5

1

∆s [m]

Eficiên

cia

2 linhas4 linhas6 linhas

Figura 8 – Eficiência - Modelo de AV deLea

equivalente, teria de ser executada em processadores com velocidades de 3,33GHz / 4534 =734KHz e 3,33GHz / 4188 = 795KHz. Considerando que o processador utilizado por Costaem 1995 tinha 66Mhz, pode-se estimar que a implementação descrita neste trabalho é,grosso modo, cerca de 82 vezes (66MHz / 0,8MHz) mais veloz do que o programa original.Comparações mais apuradas poderiam envolver a recompilação do código-fonte originalem linguagem Fortran, ou a recompilação do novo código-fonte em linguagem C utilizandoferramental similar ao disponível em 1995, ou ainda a adoção de métricas padronizadascomo a medição de desempenho em FLOPS (operações de ponto flutuante por segundo,no original em inglês). Dito isso, segue a análise da implementação paralela.

A partir das Tabelas 4 e 5 e Figuras 3 e 4 pode-se verificar que, tanto para o modelode AV de Gibbs quanto o de Lea, os tempos de execução para valores de ∆s menoresou iguais a 0,5 metros melhoram sempre que se utilizam mais linhas de execução, comganhos progressivamente mais pronunciados para valores menores de ∆s. Além disso, emconsonância com a análise realizada anteriormente, pode-se constatar que os tempos deprocessamento observados para o modelo de AV de Lea são consideravelmente mais altos


do que os observados para o modelo de Gibbs, e que, ainda em comparação ao modelo deGibbs, o desempenho paralelo do modelo de Lea não apresenta uma escalabilidade tão boa.A Tabela 6 e as Figuras 5 e 6 permitem chegar a estas mesmas conclusões, mas através dospeedup paralelo. Conforme discutido anteriormente, a exigência computacional adicionaldo modelo de Lea em relação ao modelo de Gibbs se deve ao cálculo mais complexo dapressão de descarga da bomba de fundo, que demanda a utilização da viscosidade dinâmicado fluido e o percurso de todos os pontos em que a coluna de hastes foi dividida. Nomodelo de Gibbs, estes detalhes são simplificados através da utilização do coeficiente deamortecimento.

A Tabela 6 e as Figuras 7 e 8 fornecem detalhes sobre a eficiência paralela. Aeficiência para o modelo de Gibbs aumenta juntamente com o tamanho do problema paratodas as quantidades de linhas de execução, alcançando um limite para ∆s menor ou iguala 0,1 metro ao utilizar duas ou quatro linhas de execução. A eficiência para o modelo deLea apresenta um comportamento similar, mas a eficiência global é visivelmente menor doque a verificada para o modelo de Gibbs.

Para orientar adequadamente os esforços de otimização e paralelização, bem comoavaliar seus resultados, é necessário analisar onde o programa simulador consome seutempo. As Tabelas 7 e 8 mostram os tempos absolutos e percentuais para uma únicaexecução do modelo de AV de Gibbs quando ∆s é igual a 0,05 metros, individualizando otempo exigido por cada parte da simulação. Os tempos seriais mostram que as Partes 2 e5, que a análise realizada na Seção 4.1.1 sugeriu serem aquelas em que se deveria focar osesforços de paralelização, exigiram 99,48% do tempo total de execução. As partes restantesexigiram 0,52% do tempo total, que será considerado como a porção estritamente serial doproblema ao utilizar o modelo de AV de Gibbs.

Tabela 7 – Tempos absolutos de execução - Modelo de AV de Gibbs - ∆s = 0,05 m

Tempo paralelo (s)Parte Tempo serial (s) 2 linhas 4 linhas 6 linhas

1 0,941 1,119 1,180 1,0672 379,112 207,180 112,703 79,7413 0,266 0,430 0,469 0,3344 0,614 0,991 0,988 1,1585 68,959 38,132 20,959 13,9576 0,224 0,227 0,279 0,2847 0,311 0,265 0,282 0,274

Total 450,427 248,344 136,860 96,815

Ao inspecionar os tempos paralelos nas Tabelas 7 e 8, pode-se verificar que, enquantoo tempo total de execução melhorou, o tempo percentual de execução exigido pelas partesestritamente seriais aumentou para 1,22%, 2,33% e 3,21% do tempo total de execução paraduas, quatro e seis linhas de execução, respectivamente. Dito isso, as flutuações absolutas


Tabela 8 – Tempos percentuais de execução - Modelo de AV de Gibbs - ∆s = 0,05 m

Tempo paraleloParte Tempo serial 2 linhas 4 linhas 6 linhas

1 0,21% 0,45% 0,86% 1,10%2 84,17% 83,42% 82,35% 82,36%3 0,06% 0,17% 0,34% 0,34%4 0,14% 0,40% 0,72% 1,20%5 15,31% 15,35% 15,31% 14,42%6 0,05% 0,09% 0,20% 0,29%7 0,07% 0,11% 0,21% 0,28%

Tabela 9 – Tempos absolutos de execução - Modelo de AV de Lea - ∆s = 0,05 m

Tempo paralelo (s)Parte Tempo serial (s) 2 linhas 4 linhas 6 linhas

1 0,927 1,142 1,153 1,3632 396,267 237,162 128,543 97,2313 0,253 0,415 0,358 0,3344 194,948 142,797 150,488 156,2685 67,372 39,850 26,209 20,7826 0,216 0,189 0,269 0,2837 0,237 0,306 0,262 0,247

Total 660,220 421,861 307,282 276,508

Tabela 10 – Tempos percentuais de execução - Modelo de AV de Lea - ∆s = 0,05 m

Tempo paraleloParte Tempo serial 2 linhas 4 linhas 6 linhas

1 0,14% 0,27% 0,38% 0,49%2 60,02% 56,22% 41,83% 35,16%3 0,04% 0,10% 0,12% 0,12%4 29,53% 33,85% 48,97% 56,51%5 10,20% 9,45% 8,53% 7,52%6 0,03% 0,04% 0,09% 0,10%7 0,04% 0,07% 0,09% 0,09%

no tempo exigido pelas partes estritamente seriais foram modestas. Quanto às partesque foram paralelizadas, tanto a Parte 2 quanto a Parte 5 tiveram pequenos decréscimosno tempo percentual de execução. De maneira geral, pode-se concluir que a simulaçãoutilizando o modelo de AV de Gibbs apresentou boa escalabilidade.

As Tabelas 9 e 10 mostram os tempos absolutos e percentuais para uma únicaexecução do modelo de AV de Lea quando ∆s é igual a 0,05 metros, individualizando otempo exigido por cada parte da simulação. Os tempos seriais mostram que as Partes 2, 4e 5, que a análise realizada na Seção 4.1.1 sugeriu serem aquelas em que se deveria focar osesforços de paralelização, exigiram 99,75% do tempo total de execução. As partes restantes


exigiram 0,25% do tempo total, que será considerado como a porção estritamente serial doproblema ao utilizar o modelo de AV de Lea.

Ao utilizar o modelo de Lea, a Parte 4 necessita de atenção especial. Conformeexplicado na Seção 4.1.3, a técnica de paralelização escolhida para implementar esta parteemprega apenas duas linhas de execução para checar a convergência da estimativa inicialpara a velocidade do pistão, não importando quantos processadores estejam disponíveis,e uma única linha de execução caso sejam necessários passos de convergência adicionais.Portanto, no melhor caso, o tempo paralelo para esta parte não pode ser menor do que ametade do melhor tempo serial correspondente. Quanto às demais partes paralelas, quandocomparadas à implementação encontrada no modelo de AV de Gibbs, a Parte 2 é quaseidêntica, enquanto a Parte 5 é totalmente idêntica. Dessa maneira, as maiores diferençasentre os comportamentos de execução dos dois modelos de AV devem estar localizadas naParte 4.

Ao inspecionar os tempos paralelos nas Tabelas 9 e 10, pode-se verificar que,enquanto o tempo total de execução melhorou, o tempo percentual de execução exigidopelas partes estritamente seriais aumentou para 0,48%, 0,68% e 0,80% do tempo totalde execução para duas, quatro e seis linhas de execução, respectivamente. As flutuaçõesabsolutas no tempo exigido pelas partes estritamente seriais também foram modestas.Quanto às partes que foram paralelizadas, as Partes 2 e 5 tiveram grandes decréscimosno tempo percentual de execução, enquanto que a Parte 4 teve grandes acréscimos. Notocante aos tempos absolutos de execução, para as Partes 2 e 5, eles sempre decrescem;para a Parte 4, é interessante notar que, apesar de todos os tempos paralelos apresentaremmelhoria em relação ao tempo serial, eles sempre crescem juntamente com a quantidade delinhas de execução globais, apesar de a carga de trabalho e a quantidade efetiva de linhas deexecução empregadas nesta parte permanecerem fixas, conforme discutido anteriormente.

Para investigar mais detalhadamente o desempenho da Parte 4 do modelo de AV deLea, a Tabela 11 lista tempos de execução para aquela parte para uma única execução sobdiversos valores de ∆s, com a Figura 9 exibindo estas informações em formato de gráfico.Estes tempos mostram que o comportamento observado para ∆s igual a 0,05 metrosé consistente com quase todos os valores de ∆s, com quase todos os tempos paralelosmostrando melhorias em relação ao tempo serial e/ou aumentando juntamente com aquantidade de linhas de execução globais. A Tabela 11 se concentra no speedup obtido naParte 4 para diversos valores de ∆s, mostrando que, de forma geral, o speedup melhorapara uma quantidade global fixa de linhas de execução conforme ∆s diminui, mas regridepara um valor fixo de ∆s conforme a quantidade global de linhas de execução aumenta.

Embora esta situação exija investigações mais aprofundadas, cremos que os resul-tados desfavoráveis na Parte 4 decorrem de “overhead” excessivo na gestão de linhas deexecução e da necessidade de sincronizar acessos de memória entre processadores. Levando


Tabela 11 – Tempos de execução e speedup - Parte 4 do modelo de AV de Lea

Tempos de execução (s) Speedup∆s (m) Serial 2 linhas 4 linhas 6 linhas 2 linhas 4 linhas 6 linhas

14,67 0,009 0,008 0,014 0,012 1,13 0,64 0,750,5 2,086 1,738 1,843 2,127 1,20 1,13 0,980,4 3,062 2,664 2,855 3,04 1,15 1,07 1,010,3 5,36 4,419 4,779 5,072 1,21 1,12 1,060,2 12,236 9,696 10,088 10,757 1,26 1,21 1,140,1 48,824 36,178 38,477 40,379 1,35 1,27 1,210,05 194,246 142,973 149,376 155,658 1,36 1,30 1,25

14,67 0,5 0,4 0,3 0,2 0,1 0,050

50

100

150

200

∆s [m]

Tempo

[s]

Serial 2 linhas4 linhas 6 linhas

Figura 9 – Tempos de execução - Parte 4 do modelo de AV de Lea

em consideração os números anteriores, pode-se concluir que a simulação ao utilizar omodelo de AV de Lea obteve uma escalabilidade razoável, mas a Parte 4 se destacoucomo um gargalo de desempenho significativo e permanece como um alvo importante paraesforços adicionais de otimização.

53

6 Conclusões

Neste trabalho, diversas otimizações de desempenho foram aplicadas a uma im-plementação do modelo de Costa (1995) para poços de bombeio mecânico direcionais.Comparado à implementação original de Costa, nosso simulador mostrou-se mais de 4000vezes mais rápido. Desconsiderando as melhorias de desempenho decorrentes da evoluçãoentre gerações de processadores, obtivemos um ganho de desempenho de 82 vezes paranossa implementação serial. Se levarmos em conta os valores de speedup da implementaçãoparalela, o ganho de desempenho chega a 360 vezes. Além disso, nossa versão paralelademonstrou comportamento escalável, com desempenho paralelo crescente associado aoaumento do tamanho do problema (ou, equivalentemente, à diminuição do parâmetro ∆s).Juntos, estes fatos são indicativos de que este programa deve ser uma opção viável parauso diário por parte de engenheiros de petróleo.

Além disso, dada a demonstração de que o modelo de Costa pode ser implementadode maneira mais eficiente, novos tópicos de pesquisa podem ser explorados. Por exemplo,o modelo de Costa pode ser expandido com a inclusão de fatores tais como o escorrega-mento do motor e o acoplamento do reservatório, visto que há agora um maior poder deprocessamento disponível em consequência das otimizações de software. Outra linha deestudo é a utilização do modelo na otimização de parâmetros operativos de poços; porexemplo, na maximização de produção de óleo e na minimização do consumo de energia,integrando trabalhos tais como os de Gu et al. (2014) e Xing e Dong (2015).

Também enxergamos alguns caminhos possíveis para ganhos de desempenho adi-cionais. O primeiro seria o refinamento e otimização do próprio modelo de Costa – porexemplo, a Parte 4 do Laço Principal de Simulação ao utilizar o modelo de atrito viscosode Lea se destacou como um gargalo de desempenho significativo, sendo portanto uma boacandidata a análises adicionais. O segundo seria avaliar soluções numéricas alternativas,tais como métodos de diferenças finitas implícitos. Um terceiro caminho seria experimen-tar outras soluções de paralelização, tais como a utilização de GPUs (DZIEKONSKI etal., 2013) e/ou a distribuição do trabalho entre nós de computação heterogêneos (LI;STALZER; ORTIZ, 2014).

54

Referências

ALHANATI, F. J. S. Modelo de previsão de comportamento de bombeio mecânico. [S.l.],1988. Citado na página 21.

AMDAHL, G. M. Validity of the single processor approach to achieving large scalecomputing capabilities. In: Proceedings of The AFIPS Spring Joint Computer Conference.AFIPS Press, 1967. p. 483–485. Disponível em: <http://www-inst.eecs.berkeley.edu/~n252/paper/Amdahl.pdf>. Citado na página 33.

AMERICAN PETROLEUM INSTITUTE. API RP 11L - Recommended practice fordesign calculations for sucker rod pumping systems. [S.l.], 2000. Citado na página 17.

AMERICAN PETROLEUM INSTITUTE. API SPEC 11E - Specifications for PumpingUnits, Nineteenth Edition. [S.l.], 2015. Citado na página 25.

ARAUJO, A. P. et al. Determination of downhole dynamometer cards for deviatedwells. In: SPE Artificial Lift Conference — Latin America and Caribbean. Society ofPetroleum Engineers, 2015. Disponível em: <https://www.onepetro.org/conference-paper/SPE-173970-MS>. Citado na página 18.

COSTA, R. O. Bombeamento Mecânico Alternativo em Poços Direcionais. Dissertação(Mestrado) — Universidade Estadual de Campinas, Novembro 1995. Disponível em:<http://www.bibliotecadigital.unicamp.br/document/?code=vtls000104042>. Citado 15vezes nas páginas 7, 8, 17, 18, 19, 21, 22, 24, 26, 31, 35, 36, 42, 45 e 53.

COSTA, R. O. Curso de Bombeio Mecânico. [S.l.], 2004. Citado na página 17.

DONG, S. et al. The dynamic simulation model and the comprehensive simulationalgorithm of the beam pumping system. In: 2013 International Conference on Mechanicaland Automation Engineering (MAEE). [S.l.: s.n.], 2013. p. 118–122. Citado na página 18.

DOTY, D. R.; SCHMIDT, Z. An improved model for sucker rod pumping. Society ofPetroleum Engineers Journal, v. 23, n. 01, p. 33–41, February 1983. Citado na página 21.

DZIEKONSKI, A. et al. Generation of large finite-element matrices on multiple graphicsprocessors. International Journal for Numerical Methods in Engineering, v. 94, n. 2, p.204–220, 2013. ISSN 1097-0207. Citado na página 53.

EVCHENKO, V. S.; ZAKHARCHENKO, N. P. Calculation of loads in deviated wellswhen producing with downhole pumps (original em Russo). Neftyance Khozyaistvo, p.34–37, 1984. Citado na página 17.

GIBBS, S. G. Predicting the behavior of sucker rod pumping systems. Journal ofPetroleum Technology, v. 15, n. 07, p. 769–778, July 1963. Citado 2 vezes nas páginas 21e 23.

GIBBS, S. G. A general method for predicting rod pumping system performance(SPE-6850). [S.l.], 1977. Citado na página 21.

http://www-inst.eecs.berkeley.edu/~n252/paper/Amdahl.pdf

http://www-inst.eecs.berkeley.edu/~n252/paper/Amdahl.pdf

https://www.onepetro.org/conference-paper/SPE-173970-MS

https://www.onepetro.org/conference-paper/SPE-173970-MS

http://www.bibliotecadigital.unicamp.br/document/?code=vtls000104042

Referências 55

GIBBS, S. G. Design and diagnosis of deviated rod-pumped wells. Journal of PetroleumTechnology, v. 44, n. 07, p. 774–781, July 1992. Citado na página 17.

GU, X. et al. Modeling and optimization of beam pumping system based on intelligentcomputing for energy saving. Journal of Applied Mathematics, v. 2014, n. Article ID317130, p. 7 pages, 2014. Citado na página 53.

LAINE, R. E.; COLE, D. G.; JENNINGS, J. W. Harmonic polished rod motion(SPE-19724). [S.l.], 1989. Citado 4 vezes nas páginas 10, 21, 24 e 31.

LAINE, R. E.; KEATING, J. F.; JENNINGS, J. W. Shallow sucker rod wells fluidinertia. In: 316-340 (Ed.). 37th Southwestern Petroleum Short Course. Lubbock,1990. Disponível em: <https://www.swpshortcourse.org/program/abstracts/shallow-sucker-rod-wells-fluid-inertia>. Citado na página 58.

LEA, J. F. Modeling forces on a beam pump system when pumping highly viscous crude(SPE-20672). [S.l.], 1991. Citado 4 vezes nas páginas 21, 23, 24 e 26.

LI, B.; STALZER, M.; ORTIZ, M. A massively parallel implementation of the optimaltransportation meshfree method for explicit solid dynamics. International Journal forNumerical Methods in Engineering, v. 100, n. 1, p. 40–61, 2014. ISSN 1097-0207. Citadona página 53.

LUKAZIEWICZ, S. A. Dynamic behavior of the sucker rod string in the inclined well(SPE-21665). [S.l.], 1991. Citado 2 vezes nas páginas 17 e 21.

MISKA, S.; SHARAKI, A.; RAJTAR, J. A simple model for computer-aided optimizationand design of sucker-rod pumping systems. Journal of Petroleum Science and Engineering,v. 17, n. 3–4, p. 303 – 312, 1997. ISSN 0920-4105. Citado na página 17.

OPENMP API specification for parallel programming, The. 2015. Disponível em:<http://www.openmp.org>. Citado na página 35.

OSTROVSKY, I. Big Oh in the parallel world. 2008. Disponível em: <http://igoro.com/archive/big-oh-in-the-parallel-world/>. Citado na página 40.

SCHAFER, D. J.; JENNINGS, J. W. An investigation of analytical and numerical suckerrod pumping mathematical models (SPE-16919). [S.l.], 1987. Citado na página 31.

SUTTER, H. The Free Lunch is Over. 2004. Disponível em: <http://www.gotw.ca/publications/concurrency-ddj.htm>. Citado na página 31.

SUTTER, H. Welcome to The Jungle. 2012. Disponível em: <http://herbsutter.com/welcome-to-the-jungle/>. Citado na página 32.

TAKÁCS, G. Sucker-Rod Pumping Manual. [S.l.]: PennWell Corporation, 2003. ISBN9780878148929. Citado 3 vezes nas páginas 9, 19 e 31.

XINFU, L. et al. Calculation of dynamic loads of the sucker rod pumping system in CBMwells. Journal of Coal Science and Engineering (China), China Coal Society, v. 16, n. 2, p.170–175, 2010. ISSN 1006-9097. Citado na página 18.

XINFU, L. et al. An approach to the design calculation of sucker rod pumping systems incoalbed methane wells. Chinese Journal Of Mechanical Engineering, v. 24, p. 1, 2011.Citado na página 18.

https://www.swpshortcourse.org/program/abstracts/shallow-sucker-rod-wells-fluid-inertia

https://www.swpshortcourse.org/program/abstracts/shallow-sucker-rod-wells-fluid-inertia

http://www.openmp.org

http://igoro.com/archive/big-oh-in-the-parallel-world/

http://igoro.com/archive/big-oh-in-the-parallel-world/

http://www.gotw.ca/publications/concurrency-ddj.htm

http://www.gotw.ca/publications/concurrency-ddj.htm

http://herbsutter.com/welcome-to-the-jungle/

http://herbsutter.com/welcome-to-the-jungle/

Referências 56

XING, M.; DONG, S. A new simulation model for a beam-pumping system applied inenergy saving and resource-consumption reduction. SPE Production & Operations, v. 30,n. 02, May 2015. Citado na página 53.

XU, J. A method for diagnosing the performance of sucker rod string in straight inclinedwells (SPE-26970). [S.l.], 1994. Citado na página 17.

XU, J. A New Approach to the Analysis of Deviated Rod-Pumped Wells (SPE-28697).[S.l.], 1994. Citado na página 17.

XU, J. et al. A Comprehensive Rod-Pumping Model and Its Applications to Vertical andDeviated Wells (SPE-52215). [S.l.], 1999. Citado na página 17.

XU, J.; HU, Y. R. A method for designing and predicting the sucker rod string in deviatedpumping wells (SPE-26929). [S.l.], 1993. Citado na página 17.

Apêndices

58

APÊNDICE A – Esquemas de DiferençasFinitas para Cálculo da Força e Velocidade das

Hastes

O sistema de equações diferenciais a ser resolvido pelo simulador é dado pelasexpressões abaixo:

∂vr(s, t)∂t

= 1ρrAr

∂Fs(s, t)∂s

+ ~g · ~T (s)


√√√√[~g · ~B(s)]2

+[~g · ~N(s) +


]2

+ fv(s, t)ρrAr

(A.1)

∂vr(s, t)∂s

= 1ErAr

∂Fs(s, t)∂t

(A.2)

A discretização de (A.2) é imediata, e segue o esquema descrito por Laine, Keatinge Jennings (1990):

F j+1si = ErAr

2∆t∆s(vj+1

ri+1 − vj+1ri−1) + F j

si . (A.3)

A discretização de (A.1), porém, depende do modelo de AV(termo fv(s, t)

ρrAr

)escolhido pelo engenheiro. Caso se opte pelo modelo de AV de Gibbs, (A.1) será desenvolvidaconforme segue:

∂vr(s, t)∂t

= 1ρrAr

∂Fs(s, t)∂s

+ ~g · ~T (s)


√√√√[~g · ~B(s)]2

+[~g · ~N(s) +


]2

− cvr(s, t) .

(A.4)

O esquema de diferenças finitas será dado por:

vj+1ri − vj

ri

∆t = 1ρrAr

F jsi+1 − F

jsi−1

2∆s + ~g · ~T + f jci

ρrAr

− c(vj+1ri + vj

ri)2 .

(A.5)

APÊNDICE A. Esquemas de Diferenças Finitas para Cálculo da Força e Velocidade das Hastes 59

Ao explicitar vj+1ri , vem:

vj+1ri =

∆tρrAr∆s

(F jSi+1 − F

jSi−1) + 2~g · ~Ti∆t

2 + c∆t

+

2f jci∆tρrAr

+ (2− c∆t)vjri

2 + c∆t ,

(A.6)

onde:f j

ci

ρrAr

= −µ vjri∣∣∣vjri

∣∣∣√√√√(~g · ~Bi

)2+(~g · ~Ni +

∣∣∣ ~Ki

∣∣∣ F jsi

ρrAr

)2

. (A.7)

Porém, caso a opção seja pelo modelo de AV de Lea, (A.1) será desenvolvidaconforme segue:

∂vr(s, t)∂t

= 1ρrAr

∂Fs(s, t)∂s

+ ~g · ~T (s)


√√√√[~g · ~B(s)]2

+[~g · ~N(s) +


]2

+ ηUrk

ρrArk

(K1vr −K2v̄fk) .

(A.8)

O esquema de diferenças finitas será dado por:

vj+1ri − vj

ri

∆t = 1ρrAr

F jsi+1 − F

jsi−1

2∆s + ~g · ~Ti + f jci

ρrAr

− ηiUrK1

ρrAr

vj+1ri + vj

ri

2 + ηiUrK2

ρrAr

vjf .

(A.9)

Ao explicitar vj+1ri , vem:

vj+1ri =

∆tρrAr∆s

(F jSi+1 − F

jSi−1) + 2~g · ~Ti∆t

2− ηiUrK1

ρrAr

∆t

+

2f jci∆tρrAr

+(

2 + ηiUrK1

ρrAr

∆t)vj

ri − 2ηiUrK2∆tρrAr

vjf

2− ηiUrK1

ρrAr

∆t.

(A.10)

60

APÊNDICE B – Dados Reais e Simuladosde Poços de Petróleo

São apresentados a seguir os dados reais e simulados de alguns poços de petróleoexistentes, todos eles localizados na área terrestre da Bacia Potiguar. Para cada poço, sãolistadas as informações abaixo:

• Configuração de poço utilizada na simulação;

• Perfil (trajetória direcional) do poço;

• Dados reais e simulados referentes aos parâmetros de operação do poço;

• Cartas dinamométricas de superfície reais e simuladas, e as respectivas cartas dina-mométricas de fundo derivadas a partir das cartas de superfície. No caso das cartassimuladas, são mostradas tanto as cartas para o modelo de AV de Gibbs quantopara o modelo de AV de Lea.

Os dados reais de cada poço foram obtidos a partir de um sistema supervisóriopara poços automatizados, baseado em sensores de automação instalados em campo.

APÊNDICE B. Dados Reais e Simulados de Poços de Petróleo 61

Tabela 12 – Poço 1: Dados da configuração de poço.

Dado Valor

Ciclos de bombeio 5Coeficiente de atrito de Coulomb (cD) 0,05Comprimento do segmento utilizado para dividir a coluna de hastes (∆s) 13,45 mFator de amortecimento (apenas modelo de AV de Gibbs) (µ) 0,15

BSW do fluido 80%Coluna de produção ancorada? NãoCurso da unidade de bombeio 17,04 inDensidade relativa da água produzida 1,0Densidade relativa do gás 0,7Diâmetro da primeira seção de hastes 5/8 inDiâmetro do pistão 1,5 inDiâmetro externo da coluna de produção 2,375 inDiâmetro interno da coluna de produção 1,992 inEficiência de separação de gás no fundo do poço 70%Espaço morto 1 péFrequência de bombeio 10,9 ciclos/minutoGradiente geotérmico 0,01◦F/mGravidade API do óleo 28,6◦

Módulo de elasticidade das hastes 29,6 · 106 psiMódulo de elasticidade da coluna de produção 29,6 · 106 psiNível dinâmico 403,86 mPacker instalado? NãoPressão de separação 14,69 psiPressão de sucção 20 psiPressão na cabeça 71,1 psiPressão no anular 14,22 psiProfundidade da bomba de fundo 403,86 mQuantidade de hastes na primeira seção 53Razão gás/óleo 0,1Temperatura ambiente na superfície 37,77◦CTemperatura de separação 26,66◦CViscosidade dinâmica do óleo a 40◦C 1 cPViscosidade dinâmica do óleo a 80◦C 2 cP


Tabela 13 – Poço 1: Perfil do poço.


0 0 0 0 0 05 0 0 5 0 033 2,5 152 32,99 -0,54 0,2942 4,25 150 41,98 -1 0,5550 5,75 150 49,94 -1,61 0,8959 7 151 58,89 -2,48 1,3968 8,5 152 67,81 -3,54 1,9678 10 159 77,68 -5,01 2,6286 11,5 161 85,54 -6,41 3,1395 12,75 158 94,33 -8,18 3,79104 14,25 156 103,09 -10,11 4,62113 16 154 111,77 -12,24 5,61122 17,75 155 120,39 -14,6 6,74131 19,25 156 128,92 -17,19 7,92140 20,5 154 137,38 -19,97 9,21149 22,25 152 145,77 -22,89 10,7158 23,5 153 154,06 -25,99 12,32166 24,25 156 161,37 -28,91 13,71204 28,25 161 195,45 -44,55 19,82239 29,75 163 226,06 -60,69 25,05283 32,5 162 263,73 -82,38 31,9330 34,75 162 302,86 -107,13 39,94385 38 162 347,14 -138,15 50,02433 40,75 164 384,24 -167,27 58,9480 42,25 164 419,44 -197,2 67,49535 44 164 459,58 -233,34 77,85582 43,75 163 493,46 -264,57 87,1638 43,25 163 534,08 -301,44 98,37696 42,5 162 576,59 -339,07 110,24752 42 161 618,04 -374,78 122,18808 41,5 160 659,82 -409,93 134,63855 41,25 159 695,09 -439,03 145,51876 41,25 159 710,88 -451,95 150,47

Tabela 14 – Poço 1: Dados reais e simulados.

Dado Valor Valor Simulado / Valor Simulado /Real % Desvio (Modelo % Desvio (Modelo

de AV de Gibbs) de AV de Lea)

PPRL (lbf) 2720 2422,2 / -10,9% 2520,09 / -7,3%MPRL (lbf) 912 926,19 / +1,5% 869,07 / -4,7%Peak Torque (lbf.in) 10132,24 7746,15 / -23,5% 8754,23 / -13,6%PRHP (hp) 0,63 0,39 / -38% 0,46 / -26,9%Vazão (m3/d) 6,89 6,76 / -1,8% 6,99 / +1,4%


0 2 4 6 8 10 12 14 16 18

0

1000

2000

3000

4000

Posição [in]

Carga

[lbf]

Carta de superfície simuladaCarta de fundo derivada da carta de superfície simulada

Carta de superfície realCarta de fundo derivada da carta de superfície real

Figura 10 – Poço 1: Cartas dinamométricas simuladas/reais de superfície, e cartasderivadas de fundo - Modelo de AV de Gibbs.

0 2 4 6 8 10 12 14 16 18

0

1000

2000

3000

4000

Posição [in]

Carga

[lbf]



Figura 11 – Poço 1: Cartas dinamométricas simuladas/reais de superfície, e cartasderivadas de fundo - Modelo de AV de Lea.



Dado Valor


BSW do fluido 94%Curso da unidade de bombeio 65,75 inColuna de produção ancorada? SimDensidade relativa da água produzida 1,0Densidade relativa do gás 0,7Diâmetro da primeira seção de hastes 7/8 inDiâmetro da segunda seção de hastes 6/8 inDiâmetro do pistão 1,75 inDiâmetro externo da coluna de produção 2,874 inDiâmetro interno da coluna de produção 1,322 inEficiência de separação de gás no fundo do poço 70%Espaço morto 1 péFrequência de bombeio 10,4 ciclos/minutoGradiente geotérmico 0,01◦F/mGravidade API do óleo 28,6◦

Módulo de elasticidade das hastes 29,6 · 106 psiMódulo de elasticidade da coluna de produção 29,6 · 106 psiNível dinâmico 769,62 mPacker instalado? NãoPressão de separação 14,69 psiPressão de sucção 20 psiPressão na cabeça 14,22 psiPressão no anular 5 psiProfundidade da bomba de fundo 769,62 mQuantidade de hastes na primeira seção 1Quantidade de hastes na segunda seção 100Razão gás/óleo 1Temperatura ambiente na superfície 37,77◦CTemperatura de separação 26,66◦CViscosidade dinâmica do óleo a 40◦C 1 cPViscosidade dinâmica do óleo a 80◦C 2 cP




0 0 0 0 0 04,9 0 0 4,9 0 05 0 0 5 0 030 0 245,98 30 0 060 0 245,98 60 0 090 3,5 254 89,98 -0,25 -0,88120 8,25 242 119,82 -1,52 -3,66150 14 237 149,24 -4,51 -8,61180 18 237 178,07 -9,01 -15,54210 24,25 237 206,04 -14,89 -24,61240 28 243 232,98 -21,45 -36,06270 31 236 259,09 -28,97 -48,74300 36,5 238 284,03 -38,03 -62,72304 38 242 287,21 -39,24 -64,82352 40 236 324,53 -54,81 -90,67389 40 237 352,87 -67,93 -110,5445 41 238 395,45 -87,47 -141,17502 41,5 239 438,31 -107,1 -173,22577 42 239 494,26 -132,83 -216,03652 42,5 239 549,78 -158,8 -259,25708 43 240 590,9 -178,09 -292793 44 240 652,56 -207,34 -342,67869 44,5 240 707 -233,86 -388,6891 44,5 240 722,69 -241,57 -401,96




PPRL (lbf) 7184 9374,57 / +30,5% 8902,41 / +23,9%MPRL (lbf) 2208 907,34 / -58,9% 627,28 / -71,6%Peak Torque (lbf.in) 118662 151212,34 / +27,4% 139294,71 / +17,3%PRHP (hp) 5,42 8,87 / +63,6% 7,35 / +35,6%Vazão (m3/d) 32,45 36,1 / +11,2% 36,85 / +13,5%


0 10 20 30 40 50 60 70−2000

0

2000

4000

6000

8000

10000

12000

14000

Posição [in]

Carga

[lbf]




0 10 20 30 40 50 60 70−2000

0

2000

4000

6000

8000

10000

12000

14000

Posição [in]

Carga

[lbf]






Dado Valor


BSW do fluido 67%Coluna de produção ancorada? SimCurso da unidade de bombeio 94 inDensidade relativa da água produzida 1,0Densidade relativa do gás 0,7Diâmetro da primeira seção de hastes 7/8 inDiâmetro do pistão 1,75 inDiâmetro externo da coluna de produção 2,874 inDiâmetro interno da coluna de produção 2,322 inEficiência de separação de gás no fundo do poço 70%Espaço morto 1 péFrequência de bombeio 10,5 ciclos/minutoGradiente geotérmico 0,01◦F/mGravidade API do óleo 28,6◦

Módulo de elasticidade das hastes 29,6 · 106 psiMódulo de elasticidade da coluna de produção 29,6 · 106 psiNível dinâmico 777,24 mPacker instalado? NãoPressão de separação 14,69 psiPressão de sucção 20 psiPressão na cabeça 14,22 psiPressão no anular 7 psiProfundidade da bomba de fundo 777,24 mQuantidade de hastes na primeira seção 102Razão gás/óleo 2Temperatura ambiente na superfície 37,77◦CTemperatura de separação 26,66◦CViscosidade dinâmica do óleo a 40◦C 1 cPViscosidade dinâmica do óleo a 80◦C 2 cP




0 0 0 0 0 021 0,97 335,89 21 0,16 -0,0749 1,41 329,83 48,99 0,68 -0,3476 2,55 300,12 75,98 1,27 -1,03104 5,63 289,49 103,9 2,04 -2,86132 9,41 287,47 131,66 3,18 -6,34159 12,66 282,63 158,15 4,49 -11,34187 16,44 278,24 185,25 5,73 -18,25215 19,79 274,46 211,86 6,67 -26,9242 23,04 272,44 236,99 7,25 -36,74270 25,94 271,38 262,47 7,63 -48,34310 29,02 275,69 297,96 8,8 -66,75338 32,45 279,91 322,03 10,77 -80,91367 33,85 276,92 346,31 13,08 -96,6395 37,64 276,57 369,03 15 -112,84424 42,12 276,74 391,28 17,16 -131,3452 44,32 275,34 411,68 19,17 -150,37480 46,6 275,16 431,32 20,99 -170,24508 46,69 275,16 450,54 22,82 -190,52535 46,34 274,98 469,12 24,56 -210,03563 46,34 276,13 488,45 26,52 -230,19591 47,31 276,3 507,61 28,73 -250,49619 47,13 276,22 526,63 30,97 -270,92647 46,78 276,3 545,74 33,2 -291,26676 46,34 275,78 565,68 35,42 -312,2704 45,9 275,6 585,09 37,42 -332,28732 47,75 276,13 604,25 39,5 -352,59760 48,01 276,65 623,02 41,82 -373,23788 47,82 276,48 641,79 44,19 -393,88817 47,31 276,39 661,36 46,59 -415,14845 47,4 276,3 680,33 48,87 -435,61865 47,22 275,86 693,89 50,42 -450,23




PPRL (lbf) 10228 10018,35 / -2% 8954,41 / -12,4%MPRL (lbf) 1074 1471,17 / +36,9% 2104,49 / +95,9%Peak Torque (lbf.in) 197674,4 195564,94 / -1% 167681,71 / -15,1%PRHP (hp) 10,21 10,58 / +3,6% 6,11 / -40,1%Vazão (m3/d) 77,77 50,04 / -35,65% 54,52 / -29,8%


−10 0 10 20 30 40 50 60 70 80 90 100

0

5000

10000

15000

Posição [in]

Carga

[lbf]




−10 0 10 20 30 40 50 60 70 80 90 100

0

5000

10000

15000

Posição [in]

Carga

[lbf]




Documents

OtimizaçãoeParalelizaçãodeumSimuladordeBombeio ... · ... Sistema de Bibliotecas ... de um simulador de bombeio mecânico para poços ... cavalo-de-pau”. A ﬁnalidade da unidade