% Template for submitting a manuscript to the Journal of Engineering of Catholic University of Petropolis.
% version 4 dated 01 January 2017.

\documentclass[oneside,a4paper,portuguese,links]{EngUCP}
%\documentclass[oneside,a4paper,english,links]{EngUCP}
%
\usepackage[utf8]{inputenc}
\usepackage{graphicx}
\usepackage{amsmath,amsfonts}
\usepackage[left]{lineno}
%\usepackage{cite}

\title{ATENUAÇÃO DA ESTOCAGEM NUMÉRICA USANDO REFINAMENTO DE MALHA NA SIMULAÇÃO DA PRODUÇÃO EM RESERVATÓRIOS DO TIPO \textit{SHALE GAS}}

\author[\voidaffil]{Rebeca C. D. Rosário}
\author[\voidaffil]{Paulo T. H. Júnior}
\author[\voidaffil]{Grazione de Souza}
\author[\voidaffil]{Helio P. Amaral Souto}
%
\affil[\voidaffil]{Instituto Politécnico, Universidade do Estado do Rio de Janeiro, Nova Friburgo, 28.625-570, RJ, Brasil}
%
%% NOTE: IF ALL AUTHORS BELONG TO THE SAME AFFILIATION
%% USE THE `\voidaffil' MACRO FOR THE AFFILIATION CODE.
%% Example:
%%\author[\voidaffil]{First A. Author}
%% \author[\voidaffil]{Second B. Author}
%% \author[\voidaffil]{Third C. Author}
%% \author[\voidaffil]{Fourth D. Author}
%% %
%% \affil[\voidaffil]{CEC, Department of Engineering and Computation,
%%Petrópolis Catholic University, Petrópolis, 25.685-070, RJ, Brazil}

\begin{document}
\pagestyle{empty}
\linenumbers
\vspace{3cm}
\maketitle

%E-mail Addresses: Put here!
\let\oldthefootnote\thefootnote
\renewcommand{\thefootnote}{\fnsymbol{footnote}}
\footnotetext[1]{E-mail address: \url{gsouza@iprj.uerj.br}}
%\footnotetext[0]{E-mail addresses: \url{rebecacostadias@gmail.com}, \url{ptarco@iprj.uerj.br}, \url{gsouza@iprj.uerj.br}, \url{helio@iprj.uerj.br}}
\let\thefootnote\oldthefootnote


\begin{keywords}
Estocagem numérica, poço horizontal, \textit{shale gas}, simulação de reservatórios, refinamento de malha.
\end{keywords}

\begin{abstract}
Neste trabalho, emprega-se o Método das Diferenças Finitas e uma técnica de acoplamento poço-reservatório (poço horizontal) baseada no cálculo do raio equivalente de Peaceman, para a simulação numérica da produção de gás natural em formações do tipo \textit{shale}. Devido às características dos reservatórios do tipo \textit{shale gas} e do modelo de acoplamento poço-reservatório utilizado, um artefato numérico (ou estocagem numérica) aparece nos instantes iniciais. No intuito de atenuar a estocagem numérica, um refinamento de malha é utilizado ao redor do poço horizontal produtor. Os fenômenos que favorecem a recuperação de gás natural, em reservatórios do tipo \textit{shale}, como, por exemplo, o escorregamento e a adsorção de gás, também foram incorporados ao modelo computacional. A influência desses efeitos físicos no comportamento da pressão do poço produtor foi adequadamente capturada nos estudos de análise de sensibilidade realizados.     
\end{abstract}

\newpage

\begin{center}
\title{NUMERICAL STORAGE ATTENUATION USING GRID REFINEMENT ON THE PRODUCTION SIMULATION IN SHALE GAS RESERVOIRS}
\end{center}

\def\keywordsname{Keywords}
\begin{keywords}
Numerical storage, horizontal well, shale gas, reservoir simulation, grid refinement.
\end{keywords}

\def\abstractname{Abstract}%
\begin{abstract}
In this work, we use the Finite Differences Method and a well-reservoir (horizontal well) coupling technique based on the Peaceman equivalent radius calculation for the numerical simulation of natural gas production in shale gas formations. Due to the characteristics of shale gas reservoirs and the well-reservoir coupling model used, a numeric artifact (or numerical storage) appears in the initial instants of time. To attenuate the numerical artifact, we use mesh refinement near the horizontal well. Phenomena that favor the recovery of natural gas in shale reservoirs, such as gas slippage and adsorption, were also incorporated into the computational model. We performed sensitivity analysis studies which adequately capture the influence of these effects on the pressure behavior of the producing well.
\end{abstract}

%--------------------------------------------------------------
\section{INTRODUÇÃO} 
%--------------------------------------------------------------

Devido aos recentes avanços tecnológicos, os reservatórios denominados não convencionais têm despertado grande interesse na indústria de petróleo e gás. No entanto, mesmo com o conhecimento e as tecnologias mais recentes, a explotação desses reservatórios ainda apresenta desafios quando se visa ao aumento das taxas de recuperação de gás~\cite{wang2013}. Nesse contexto, a modelagem do escoamento monofásico em reservatórios do tipo \textit{shale gas} é um tópico de grande importância e atual.

Os reservatórios do tipo \textit{shale gas}~\cite{florence2007,Villazon2011,Li2016}, nos quais as rochas do tipo folhelho armazenam gás natural, representam uma fonte significativa de hidrocarbonetos em diversos países, tais como os EUA, a China e a Austrália. A compreensão dos mecanismos de escoamento e o conhecimento das características geofísicas, dos reservatórios de \textit{shale gas}, são fundamentais para que seja possível alcançar uma produção economicamente viável. Uma das características destes reservatórios é a sua permeabilidade extremamente baixa, menor do que 0,001~mD. Em reservatórios convencionais, tal propriedade encontra-se na faixa de 10 a 1000~mD~\cite{wang2013}. Apesar dos baixos valores de permeabilidade e porosidade, a produção de gás natural a partir de folhelhos pode ser favorecida pelos efeitos de escorregamento do gás e/ou de adsorção do gás~\cite{Li2014,zhihui2015}.

Estudos mostram que a lei de Darcy clássica, que expressa o balanço da quantidade de movimento, é válida quando os efeitos inerciais e/ou turbulentos e o fenômeno de escorregamento não são considerados~\cite{Aziz1990}. Dados experimentais de jazidas indicam que o uso de correções para o cálculo da permeabilidade, considerando o regime de escorregamento, são relevantes no caso do escoamento em meios porosos com baixa permeabilidade~\cite{jiang2015}.

Em reservatórios do tipo \textit{shale gas}, a não contabilização dos efeitos de escorregamento, inerciais/turbulentos, de adsorção e/ou geomecânicos podem levar a erros significativos na caracterização da jazida e, consequentemente, no gerenciamento da produção~\cite{wang2013}. Portanto, é fundamental que os simuladores numéricos, utilizados como ferramentas auxiliares no planejamento da produção, contemplem a inclusão de tais efeitos.

Uma outra questão relevante, que surge na simulação de reservatórios do tipo \textit{shale gas}, é o efeito do baixo valor da permeabilidade no cálculo da pressão nos poços produtores. Devido à necessidade, em determinadas aplicações, de se obter uma estimativa da pressão no poço, foram desenvolvidas as chamadas técnicas de acoplamento poço-reservatório~\cite{pea78,Peaceman1983,Al-Mohannadi2007}. O modelo mais disseminado e aplicado, inclusive nos simuladores comerciais, é o do raio equivalente de Peaceman~\cite{pea78} (que supõe o escoamento em regime permanente na proximidade do poço), o qual leva, em função dos parâmetros do sistema poço-reservatório e da malha computacional, ao surgimento de um artefato numérico nos tempos iniciais de simulação, também chamado de estocagem numérica. Tal denominação decorre da semelhança qualitativa existente com a estocagem física, quando da representação da pressão no poço em função do logaritmo do tempo, em um gráfico muito utilizado na indústria do petróleo e conhecido como \textit{gráfico especializado}.

Quando a produção começa, já existe no poço uma quantidade de fluido que será inicialmente produzida. Assim, transcorre um certo intervalo de tempo até que fluido oriundo do reservatório deixe o poço. Esse efeito da estocagem no poço depende da quantidade de fluido presente na coluna de produção e da sua compressibilidade. Em um gráfico especializado, nos instantes iniciais a pressão do poço é representada por uma linha quase horizontal quando se considera a estocagem, de forma que a queda de pressão é mais lenta até que cesse os seus efeitos. A partir deste momento, a pressão começa a cair com mais intensidade, também descrita por uma variação linear no gráfico especializado, e ocorre o escoamento transiente no meio poroso. O seu comportamento muda, novamente, quando algum efeito de fronteira se faz sentir. Na Figura~\ref{fig:artefato}, para o caso particular da produção de óleo utilizando um poço vertical, apresenta-se uma comparação entre as soluções numérica e analítica (para um caso simplificado). A simulação numérica utiliza uma técnica de acoplamento poço-reservatório de~\citet{pea78} e percebe-se o aparecimento da estocagem numérica nos tempos iniciais.

\begin{figure*}[!ht]
\centering
\includegraphics[scale=0.4]{figs/artefato2.png}
\caption{Comparação entre as soluções numérica e analítica da pressão no poço vertical ($p_{wf}$).}
\label{fig:artefato}
\end{figure*}

A magnitude e a duração do artefato numérico dependem das propriedades do reservatório, do fluido e da malha computacional. Quanto menor a permeabilidade, mais significativo é o efeito do artefato numérico, de forma que para reservatórios do tipo \textit{shale} a estocagem numérica tende a ser relevante. Uma forma de atenuar esse efeito é a adoção de algum tipo de refinamento de malha ao redor do poço~\cite{Al-mohannadi2004,Al-Mohannadi2007}, pois o efeito do artefato decresce com a diminuição do espaçamento da malha computacional~\cite{souza2013}.  

O objetivo central deste trabalho é a atenuacão dos efeitos do artefato numérico e a anáise da influência do escorregamento e da adsorção na determinação da pressão nos poços produtores em reservatórios do tipo \textit{shale}.  Isso é feito dentro do contexto clássico da simulação de reservatórios aplicada à indústria de petróleo e gás, no qual utiliza-se malhas estruturadas, o refinamento de malha, o método das diferenças finitas~\cite{ert01} e técnicas de acoplamento poço-reservatório.

%--------------------------------------------------------------
\section{ESCORREGAMENTO E ADSORÇÃO DE GÁS}\label{Secao_escoamento}
%--------------------------------------------------------------

Particularmente, para reservatórios do tipo \textit{shale gas}, as simulações acuradas do escoamento são muitas vezes difíceis de serem realizadas, devido à complexa estrutura dos poros~\cite{wang2013}. Tipicamente, o tamanho dos poros varia de nanômetros a micrômetros, fazendo com que o fluxo mássico de gás ocorra por meio de diferentes mecanismos de transporte em diferentes escalas espaciais~\cite{civan11}. Assim, conforme ressaltado, os fenômenos de escorregamento e adsorção são relevantes. 

No escoamento em meio poroso, o fenômeno de escorregamento do gás nas superfícies sólidas ocorre quando o livre caminho médio das moléculas de gás apresenta uma escala espacial comparável à da dimensão dos poros~\cite{florence2007}. Nessa situação, a permeabilidade efetiva que deve ser considerada, para o escoamento do gás, é maior do que a permeabilidade quando do escoamento monofásico de líquidos, denominada absoluta. Dessa forma, para gases, defini-se uma permeabilidade denominada aparente, que depende das características da estrutura sólida da rocha e das condições nas quais se dá o escoamento do gás. \citet{Klinkenberg1941} definiu a permeabilidade aparente, $k_{app}$, como sendo dada por 
%
\begin{equation}
\label{eq:permApp_Berg}
k_{app} = \left[1 + \frac{b}{p}\right]k_0 ,
\end{equation}
%
\noindent onde $k_0$ é a permeabilidade absoluta, $p$ é a pressão do fluido e $b>0$ é o fator de Klinkenberg. 

Um modelo mais completo utiliza o número adimensional de Knudsen, $Kn$, na determinação da permeabilidade aparente. Esse número é definido~\cite{Villazon2011} pela razão do livre caminho médio das moléculas, $\lambda$, e raio hidráulico característico, $d$,
%
\begin{equation}
K_n=\dfrac{\lambda}{d}=\dfrac{\dfrac{\mu}{p}\sqrt{\dfrac{\pi Z RT}{2M}}}{2 \sqrt{2\tau} \sqrt{\dfrac{k_0}{\phi}}}
\label{eq:lambda_raio}
\end{equation}
%	
\noindent onde $\mu$ é a viscosidade do gás, $R$ é a constante universal dos gases, $T$ é a temperatura, $M$ é a massa molecular do gás, $\tau$ é a tortuosidade do meio poroso, $Z$ é o fator de compressibilidade e $\phi$ é a porosidade do meio poroso.

Segundo~\citet{Li2016}, uma alternativa para o cálculo de $k_{app}$ é dada por
%	
\begin{equation}
\label{eq:permApp_Knudsen}
k_{app} =  \left(1 + \alpha Kn\right)\left(1 + \frac{4 Kn}{1 + Kn}\right)k_0,
\end{equation}
%
\noindent onde $\alpha$ é o coeficiente de rarefação adimensional, utilizado na determinação da transição entre os regimes de escorregamento e o molecular livre. Para o regime de escorregamento tem-se que~$\alpha=0$. Segundo~\citet{Bestok1999}, em função do número de Knudsen, o fluido pode ser considerado como um contínuo para $Kn\leq 10^{-3}$. O regime de escorregamento ocorre para $10^{-3}<Kn<0,1$, um regime de transição surge para $0,1\leq Kn\leq 10$ e, se $Kn\geq 10$, ocorre o escoamento molecular livre.   

Por outro lado, nas paredes dos poros, em reservatórios do tipo \textit{shale}, pode ocorrer a adsorção do gás (o gás encontra-se aderido às paredes sólidas). O gás adsorvido pode chegar a representar de 20 a 85\% do total de gás armazenado no reservatório~\cite{Li2016}. A quantidade de gás adsorvida pode ser determinada, em função da pressão, a partir de uma isoterma como, por exemplo, a de Langmuir~\cite{Li2014},
%
\begin{equation}
V_{ads} = \left(\frac{p}{p_L + p}\right) V_L,
\label{eq:Vads}
\end{equation}
%
\noindent sendo $V_{ads}$ o volume de gás adsorvido por unidade de massa (para uma dada pressão), $V_L$ o volume de Langmuir (a capacidade máxima de adsorção do meio) e $p_L$ a pressão de Langmuir, na qual a quantidade de gás adsorvido alcança a metade do valor de $V_L$~\cite{Li2014}.

%--------------------------------------------------------------
\section{MODELAGEM FÍSICO-MATEMÁTICA}\label{Secao_modelagem}
%--------------------------------------------------------------

O modelo físico-matemático empregado na descrição do escoamento, incluindo os efeitos de escorregamento e de adsorção, é apresentado na sequência. As seguintes hipóteses são adotadas para o escoamento monofásico, transiente, em um reservatório de gás natural:

\begin{enumerate}
\item efeitos de escorregamento e adsorção presentes;  
\item permeabilidades absolutas constantes; 
\item compressibilidade da rocha pequena e constante; 
\item efeito gravitacional desprezível; 
\item escoamento isotérmico; 
\item ausência de reações químicas, de dano à formação e de estocagem física no poço. 
\end{enumerate}

Do balanço de massa, para o escoamento em um meio poroso e considerando a presença de gás adsorvido, tem-se que~\cite{Li2016}
%
\begin{equation}
   \frac{\partial}{\partial t}\left(\frac{\phi}{B}\right) + \frac{\partial}{\partial t}\left(\frac{\rho_s V_{ads}}{B}\right)+ \nabla\cdot\left(\frac{\mathbf{v}}{B}\,\right)
   -\frac{q_m}{V_b\rho_{sc}}=0,
   \label{eq:engres4.2}
\end{equation}
%
\noindent onde $B=\rho_{sc}/\rho$ é o fator-volume-formação~\cite{ert01}, $\rho$ é a massa específica do gás (o subscrito $sc$ indica as condições padrão), $\rho_s$ é a massa específica da formação rochosa, $\mathbf{v}$ é a velocidade superficial, $q_m$ é o termo fonte e $V_b$ é o volume total (fluido mais sólido).

O balanço da quantidade de movimento é assegurado pela lei de Darcy modificada~\cite{Li2016}
%
\begin{equation}
\mathbf{v}=-\frac{\mathbf{k}_{app}}{\mu}\left(\nabla p-\dfrac{\rho_{sc}}{B} g\nabla D\right),
\label{eq:darcy_kapp}
\end{equation}
%
\noindent onde $\mathbf{k}_{app}$ é o tensor (considerado diagonal) de permeabilidade aparente 
%	
\begin{equation}
\label{eq:permApp_Knudsen_b}
\mathbf{k}_{app} = \left(1 + \alpha Kn\right)\left(1 + \frac{4 Kn}{1 + Kn}\right)\mathbf{k}_0 ,
\end{equation}
%
\noindent $D$ é a profundidade e $g$ é a magnitude da aceleração da gravidade. Não são considerados os efeitos efeitos inerciais~\cite{Aziz1990,BC2004} e uma média geométrica é utilizada na determinação da permeabilidade absoluta média.

Um equação de estado para um gás real~\cite{ert01} é adotada, de forma que
%
\begin{equation} 
B = \frac{\rho_{sc}}{\rho}=\frac{p_{sc}M}{Z_{sc}RT_{sc}}\frac{ZRT}{pM}=\frac{p_{sc}}{T_{sc}}\frac{ZT}{p}
\label{eq:FVF3_gas2}
\end{equation}
%
\noindent e considera-se que $Z_{sc}\approx 1$~\cite{ert01}.

Como a compressibilidade da rocha é pequena e constante, pode-se empregar
%
\begin{equation}
\phi = \phi^{0} \left[ 1 + c_{\phi}(p-p^{0})\right],
\label{eq:PHI}
\end{equation}
%
\noindent onde $c_{\phi}$ é a compressibilidade da rocha e $\phi^{0}$ e $p^{0}$ são, respectivamente, a porosidade e a pressão em uma dada condição de referência~\cite{ert01}.

A pressão no poço, $p_{wf}$, é introduzida via o termo de fonte 
%
\begin{equation}
q_{sc}=\dfrac{q_{m}}{\rho}=-J_w(p-p_{wf}),
\label{eq:qsc}
\end{equation}
%
\noindent onde $J_w$ representa o índice de produtividade do poço. 

Das Eqs.~\eqref{eq:darcy_kapp} e Eq.~\eqref{eq:engres4.2} obtém-se, desprezando os efeitos gravitacionais, 
%
\begin{equation}
\Gamma\dfrac{\partial p}{\partial t}-V_b\nabla\cdot\left(\dfrac{\mathbf{k}_{app}}{B\mu}\nabla p\right)+J_w(p-p_{wf})=0,
\label{eq:engres4.1_5}
\end{equation}
% 
\noindent onde $\Gamma=\Gamma_p+\Gamma_s$, com  
%
\begin{equation}
\Gamma_p=V_b\left[\frac{\phi^0 c_{\phi}}{B}+\phi\frac{\partial}{\partial p}\left(\frac{1}{B}\right)\right]
\label{eq:GAMMA_p}
\end{equation}
%
\noindent e
%
\begin{equation}
\Gamma_s=V_b \rho_s \left[ \frac{1}{B} \frac{\partial }{\partial p}\left(V_{ads}\right) + V_{ads} \frac{\partial}{\partial p} \left( \frac{1}{B} \right) \right].
\label{eq:GAMMA_s}
\end{equation}
%
\noindent indicando os termos de acumulação provenientes dos efeitos de compressibilidade da rocha e do fluido ($\Gamma_p$) e da adsorção do gás ($\Gamma_s$).

Para que se possa resolver a Eq.~\eqref{eq:engres4.1_5} deve-se fornecer as condições inicial e de contorno. O reservatório possui um formato paralelepipédico, com arestas $L_x$, $L_y$ e $L_z$, conforme esquematizado na Figura~\ref{fig:poco}, e não há fluxo mássico através das suas fronteiras externas (reservatório selado). 

\begin{figure*}[!ht]
\centering
\includegraphics[scale=0.35]{figs/poco_horizontal.png}
\caption{Reservatório com um poço horizontal, de comprimento $L_{wf}$, paralelo ao eixo $x$.}
\label{fig:poco}
\end{figure*}

Como condição inicial impõe-se, em todo o reservatório,
%
\begin{equation}
p(x,y,z,t=0) = p_{ini}(x,y,z) ,
\label{eq:pini}
\end{equation}
%
\noindent onde $p_{ini}$ é a distribuição inicial de pressão antes de começar a produção/injeção de fluidos. Como condições de contorno externas, prescreve-se uma condição de fluxo nulo através das fronteiras do reservatório
%
\begin{equation}
\left(\frac{\partial p}{\partial x }\right)_{x=0,L_x} = 0,\quad \left(\frac{\partial p}{\partial y }\right)_{y=0,L_y} = 0,\quad \left(\frac{\partial p}{\partial z }\right)_{z=0,L_z} = 0 , 
\label{eq:contorno}
\end{equation}
%
\noindent e uma vazão total de produção ($Q_{sc}$) é prescrita como uma condição de contorno interna no poço produtor.

%--------------------------------------------------------------
\section{RESOLUÇÃO NUMÉRICA}\label{solucao_numerica}
%--------------------------------------------------------------

O Método das Diferenças Finitas~\cite{ert01} foi o escolhido para se obter a solução numérica da equação governante~\eqref{eq:engres4.1_5}, de maneira a se determinar os valores da pressão no reservatório ($p$) e no poço ($p_{wf}$). Para tanto, é construída uma malha estruturada tridimensional, Figura~\ref{fig:mdf}, onde os índices inteiros ($i$,$j$,$k$) identificam os centros das células, em uma configuração do tipo blocos centrados. As interfaces delimitando as células são indicadas por índices fracionários como, por exemplo, ($i\pm 1/2$,$j$,$k$).   

\begin{figure*}[!ht]
\centering
\includegraphics[scale=0.35]{figs/reser_discretizacao_xyz1.png}
\caption{Discretização do reservatório empregando blocos centrados.}
\label{fig:mdf}
\end{figure*}

Um refinamento de malha pode ser usado para melhor representar o escoamento na região próxima ao poço horizontal, como apresentado na Figura~\ref{fig:malha_naounif}. Os espaçamentos das células, $\Delta y_{i,j,k} $ e $\Delta z_{i,j,k}$, são escolhidos de maneira que eles sejam menores perto do poço e que aumentem, a uma determinada taxa prescrita, à medida em que as células se afastem do poço~\cite{Al-Mohannadi2007}.

\begin{figure*}[!ht]
\centering
\includegraphics[scale=0.35]{figs/malha_hw.png}
\caption{Exemplo de um refinamento de malha ao redor do poço horizontal.}
\label{fig:malha_naounif}
\end{figure*}

%--------------------------------------------------------------
\subsection{Forma discretizada da equação governante}\label{subsec:approximation}
%--------------------------------------------------------------

Na aproximação das derivadas espaciais usa-se esquemas do tipo diferenças centradas no espaço e atrasada no tempo~\cite{ert01}. Como a equação discretizada é não-linear, faze-se uso do Método de Picard~\cite{Nick2013} na sua linearização. Então, os coeficientes da equação governante discretizada são avaliados na iteração $v$ e as incógnitas na iteração $v+1$ e a forma final discreta da Eq.~\eqref{eq:engres4.1_5} é dada por 
%
\begin{eqnarray}
\frac{\Gamma^{n+1,v}_{m}}{\Delta t}\left(p_{m}^{n+1, v+1}-p_{m}^{n, v+1}\right) &=& \sum_{l\in\psi_m}\left[T^{n+1,v}_{l,m}\left(p_{l}^{n+1, v+1}-p_{m}^{n+1, v+1}\right)\right] \nonumber \\ 
																				&-& J^{n+1, v}_{w,m}\left(p_{wf,m}^{n+1, v+1}-p_{m}^{n+1, v+1}\right),
\label{discrete_xyz}
\end{eqnarray}
%
\noindent onde $\Gamma^{n+1,v}_{m}$ é determinado através de expansões conservativas~\cite{ert01}, $l=x,y,z$, $m=i,j,k$ e $\psi_m$ representa o conjunto dos nós das células vizinhas. O volume total de cada célula é calculado por $V_b=\Delta x \Delta y \Delta z$, onde $\Delta x$, $\Delta y$ e $\Delta z$ são os espaçamentos da malha nas direções $x$, $y$ e $z$ respectivamente (Figura~\ref{fig:mdf}). 

A Equação~\eqref{discrete_xyz} foi escrita mediante a introdução das transmissibilidades. Por exemplo, para a direção $x$,
%
\begin{equation}
   T^{n+1,v}_{x,i\pm\frac{1}{2},j,k}=\left(\frac{k_{app,x}A_x}{B\mu\Delta x}\right)^{n+1,v}_{i\pm\frac{1}{2},j,k},
   \label{eq:transmissibility}
\end{equation}
%
\noindent onde $A_x=\Delta y \Delta z$ e os valores de $k_{app,x}$, nas faces das células, são calculados por uma média harmônica, a partir dos valores conhecidos nos nós das células vizinhas. As transmissibilidades correspondentes às direções $y$ e $z$ são definidas de modo análogo. 

%--------------------------------------------------------------
\subsection{Acoplamento poço-reservatório}\label{subsec:acoplamento}
%--------------------------------------------------------------

Neste estudo, não se considera os efeitos friccionais, gravitacionais ou inerciais dentro do poço, sendo a vazão total de produção no poço, $Q_{sc}$, igual à soma dos termos $q_{sc}$ correspondentes às células transpassadas pelo poço, i.e.,
%  
\begin{equation}
Q_{sc}=-\sum^{i=W2}_{i=W1}J^{n+1,v}_{w,i,j,k}\left[p^{n+1,v+1}_{i,j,k}-\left(p_{wf}\right)^{n+1,v+1}_{i,j,k}\right]   
\label{Jw3}
\end{equation}
%
\noindent para um poço posicionado ao longo da direção do eixo $x$, começando pela célula $W1,j,k$ e terminando na $W2,j,k$. Entretanto, é necessária a escolha de uma posição como sendo a de referência para o cálculo da pressão no poço. No caso do poço horizontal, o índice de produtividade pode ser calculado usando a expressão~\cite{Al-mohannadi2004}
%
\begin{equation}
J_{w,i,j,k}^{n+1,v}=\left\{\frac{2\pi\sqrt{k_{z}k_{y}}}{\mu B}\frac{\Delta x}{\left[1-\left(\frac{r_w}{r_{eq}}\right)^2\right]\textrm{ln}\left(\frac{r_{eq}}{r_w}\right)}\right\}_{i,j,k}^{n+1,v}
\label{Jw2h}
\end{equation}
%
\noindent onde $r_w$ é o raio do poço de produção e o raio equivalente ($r_{eq}$) de Peaceman~\cite{Peaceman1983} é determinado por
%
\begin{equation}
r_{eq,i,j,k}=\sqrt{\frac{\left(\Delta z\Delta y\right)_{i,j,k}}{\pi}}\ \exp{(-0,5)} .
\label{reqh}
\end{equation}

%--------------------------------------------------------------
\subsection{Resolução do sistema de equações algébricas}\label{subsec:reseqal}
%--------------------------------------------------------------

Finalmente, o sistema de equações algébricas linearizado é resolvido utilizando o método dos Gradientes Conjugados~\cite{ert01}, a fim de que se possa determinar as pressões no reservatório e no poço produtor.

%--------------------------------------------------------------
\section{RESULTADOS NUMÉRICOS}
%--------------------------------------------------------------

Embora a pressão no reservatório tenha sido calculada, mostra-se apenas os resultados que dizem respeito à variação da pressão no poço. Os valores foram obtidos considerando uma taxa constante de produção de gás e, conforme já mencionado, os efeitos de escorregamento e de adsorção. Os valores da pressão do poço produtor são apresentados em gráficos de $p_{wf}$ em função do tempo. As correlações de~\citet{Dranchuk1975} e de~\citet{Lee1966} são usadas, respectivamente, nos cálculos de $Z$ e de $\mu$.

A Tabela~\ref{table1} mostra os parâmetros utilizados na construção do caso de referência. Salvo menção contrária, eles são os utilizados em todas as simulações. Nela, $P_{init}$ indica a pressão inicial no topo do reservatório e $n_w$ o número de células que contêm o poço horizontal. $\Delta t_i$ é o passo de tempo inicial, $\Delta t_f$ é o passo de tempo final e $F_{\Delta_t}$ é a razão de incremento do passo de tempo. Também, $\gamma_g=M/M_{ar}$ é a densidade do gás, onde $M_{ar}$ é a massa molecular do ar.

\begin{table}[!ht]
	\centering
	\caption{Parâmetros para o caso de referência.}\label{table1}
	\begin{tabular}{l|c}
		\hline
		Parâmetro						&Valor					\\
		\hline
		\hline
		$c_{\phi}$ (psi$^{-1}$) 		&1,0 $\times 10^{-3}$	\\
		$F_{\Delta_t}$     				&1,1					\\		
		$k_0=k_{0x}=k_{0y}=k_{0z}$ (D)  &5,0$\times 10^{-5}$ 	\\
		$L_x=L_y$ (ft)	  				&2,4 $\times 10^3$ 		\\
		$L_z$ (ft)        				&80 					\\
		$L_{wf}$ (ft)        			&600 					\\		
		$n_x$         					&48						\\               
		$n_{y}$         				&49						\\
		$n_z$   						&17						\\
		$n_w$   						&12						\\ 								
		$p_{sc}$ (psi)          		&14,65 					\\               
		$p^0=P_{init}$ (psi)   			&6,0 $\times 10^3$     	\\
		$p_L$ (psi)						&1.100 					\\
		$Q_{sc}$ (scf/dia)      		&-20.000   				\\
		$R$ (ft$^3$ psi/ lb-mol R)      &10,7316   				\\
		$r_w$ (ft)     					&0,1	 				\\
		$t_{max}$ (anos) 				&4   					\\                               
		$T_{sc}$ ($^{\circ}$R)     		&609,67  				\\		  				 		
		$V_{L}$ (ft$^3$/lbm)        	&0,2 $\times 10^{-4}$   \\
		$\gamma_g$   					&0,6					\\		
		$\Delta t_i$ (dias)   			&0,01					\\		  
		$\Delta t_f$ (dias)     		&10 					\\
		$\rho_s$ (lbm/ft$^3$)   		&200 					\\					              
		$\phi^0 =\phi_{init}$      		&0,1					\\					
		$\tau$      					&1,41   				\\                  			
		\hline
	\end{tabular}
\end{table}

A Tabela~\ref{table2} contém os dados usados no estudo do refinamento de malha. Isto é, o número de células empregadas na discretização do reservatório ($n_x$, $n_y$ e $n_z$) e do poço ($n_w$). Em se tratando da malha não-uniforme, $\Delta y_w$ e $\Delta z_w$ são os espaçamentos nas direções $y$ e $z$ das células pelas quais passa o poço no plano $yz$ (Figura~\ref{fig:poco}).

\begin{table}[!ht]
	\centering
	\caption{Malhas para o estudo de refinamento.}\label{table2}
	\vspace{6pt}
	\begin{tabular}{c|c|c|c|c|c}
		\hline
   Malha & $n_x$	& $n_y$		&$n_z$ 		& $n_w$		&$\Delta y_{w}=\Delta z_{w}$ (m)	\\ \hline\hline
      1  & 12       & 13      	& 5         & 4     	& 12,192 							\\     
      2  & 24       & 25      	& 9         & 8     	& 6,096 							\\     	
      3  & 48       & 49      	& 17        & 16    	& 3,048 							\\     	
      4  & 96       & 97      	& 33        & 32    	& 1,524 							\\  \hline									
	\end{tabular}
\end{table}

%--------------------------------------------------------------
\subsection{Refinamento empregando uma malha uniforme}
%--------------------------------------------------------------

Para uma malha uniforme, considerando a Malha 1, a pressão do poço $p_{wf}$ para tempos inferiores a 1 dia diminui lentamente e pervebe-se que a sua variação é aproximada por uma linha reta quase horizontal (Figura~\ref{fig:malha}). No entanto, é sabido que esse comportamento não corresponde ao observado em um escoamento real, a menos quando da existência da estocagem física (não incorporada ao modelo e não coincidente em termos de resultados com a estocagem numérica). Como em um reservatório do tipo \textit{shale gas} os valores das permeabilidades são muito baixos, observa-se que a influência da estocagem numérica é considerável. Após o desaparecimento do artefato numérico, captura-se adequadamente o regime radial inicial, caracterizado por uma diminuição acentuada da pressão e que pode ser vista como uma linha reta inclinada no gráfico de $p_{wf} \times t$. No decorrer do tempo, podem surgir os outros regimes típicos do escoamento monofásico em um reservatório com um poço horizontal~\cite{Al-Mohannadi2007}. Concluindo, vê-se que os valores obtidos com as Malhas 3 e 4 estão muito próximos uns dos outros, indicando que os valores numéricos estão convergindo à medida que a malha é refinada, na região onde o efeito de estocagem numérica não está presente. Nota-se, claramente, que o efeito do artefato numérico diminui com o refinamento de malha.

\begin{figure*}[!ht]
\centering
\includegraphics[scale=0.325]{figs/semref.png}
\caption{Estudo do refinamento empregando uma malha uniforme.}
\label{fig:malha}
\end{figure*}

%--------------------------------------------------------------
\subsection{Refinamento empregando uma malha não-uniforme}
%--------------------------------------------------------------

No caso do refinamento utilizando uma malha não-uniforme (células de menores dimensões ao redor do poço, vide a Figura~\ref{fig:malha_naounif}), é possível notar que o artefato numérico tem uma duração mais curta, mas, por outro lado, é exigido um maior esforço computacional. Por exemplo, para a Malha 4, as simulações levaram cerca de 2 horas e 30 minutos de execução com uma malha uniforme, enquanto que com uma malha não-uniforme o tempo total de simulação foi de 184 horas e 45 minutos. Em contrapartida ao esforço computacional, tem-se a diminuição da magnitude e da duração do artefato em todas as malhas, de forma que após 10 dias de simulação pode-se admitir que os resultados estão praticamente livres do efeito indesejado (Figura~\ref{fig:malharef}). O mesmo não seria possível com o uso das malhas uniformes. Em função dos resultados obtidos com as Malhas 3 e 4, opta-se por utilizar, nas demais simulações, a Malha 3 com o refinamento ao redor do poço.  

\begin{figure*}[!ht]
\centering
\includegraphics[scale=0.325]{figs/Comref.png}
\caption{Estudo do refinamento empregando uma malha não-uniforme.}
\label{fig:malharef}
\end{figure*}

%--------------------------------------------------------------
\subsection{Razão de crescimento do incremento de tempo}
%--------------------------------------------------------------

Também foi considerada a variação do comportamento da pressão no poço em função da mudança do incremento do passo de tempo $\Delta t$. O seu valor inicial é $\Delta t_i$ e o fator multiplicativo $F_{\Delta t}$ é aplicado de modo a se determinar o seu próximo valor, maior que o anterior, e o processo terminará quando o valor limite do incremento de tempo for alcançado ($\Delta t_f$). Dessa forma, consegue-se captar com uma maior acurácia o valor da pressão para tempos muito curtos. Foram realizadas simulações com o fator multiplicativo variando de 1,05 a 1,20. Os resultados são mostrados na Figura~\ref{fig:Dt} e percebe-se, para os valores de $F_{\Delta_t}$ investigados, que não houve uma alteração perceptível dos valores da pressão. 

\begin{figure*}[!ht]
\centering
\includegraphics[scale=0.325]{figs/Dt.png}
\caption{Resultados para as variações do crescimento do $\Delta t$.}
\label{fig:Dt}
\end{figure*}

%--------------------------------------------------------------
\subsection{Importância dos efeitos combinados}
%--------------------------------------------------------------

Antes de se passar aos resultados da análise de sensibilidade, onde alguns dos parâmetros mais importantes para a produção foram variados, foi feita uma comparação (Figura~\ref{fig:modelos}) entre quatro casos distintos (com e sem a incorporação dos efeitos de escorregamento e adsorção) e para um tempo máximo de produção de 5 anos:

\begin{enumerate}
\item modelo que utiliza a lei de Darcy clássica e que não considera os efeitos de escorregamento e adsorção do gás; 
\item modelo que incorpora somente a adsorção do gás; 
\item modelo que incorpora o escorregamento na permeabilidade aparente; 
\item modelo que incorpora os dois efeitos combinados, escorregamento e adsorção do gás.
\end{enumerate}

\begin{figure*}[!ht]
\centering
\includegraphics[scale=0.325]{figs/Comparacao.png}
\caption{Importância do escorregamento e adsorção do gás na variação de pressão do poço.}
\label{fig:modelos}
\end{figure*}

Quando se compara a queda de pressão obtida empregando a lei de Darcy clássica com a do modelo incorporando unicamente a adsorção, verifica-se que a adsorção age no sentido de diminuir a queda de pressão, embora de maneira não significativa para os parâmetros testados. Por outro lado, a inclusão somente do escorregamento contribui de modo considerável para a manutenção da pressão em níveis mais elevados, tomando como exemplo os dois primeiros modelos. Por último, constata-se que a menor queda de pressão corresponde à curva representando os valores obtidos com o modelo incorporando os efeitos combinados do escorregamento e da adsorção. Pode-se explicar esse comportamento, para uma produção a vazão constante, em função do aumento do valor permeabilidade aparente (menos resistência ao escoamento) e da maior quantidade de gás passível de ser produzida.

%--------------------------------------------------------------
\subsection{Análise de sensibilidade}
%--------------------------------------------------------------

Em seguida, são apresentados os resultados da análise de sensibilidade. Diferentemente de alguns dos casos anteriores, as simulações foram realizadas para uma produção ao longo de 5 anos,  de forma a se obter mais informações sobre o comportamento da pressão no poço quando já são percebidos os efeitos de fronteira. Foram realizadas simulações modificando os valores da permeabilidade da rocha, da densidade do fluido, da tortuosidade, do volume de Langmuir e da vazão de produção prescrita, a fim de se estudar a influência desses parâmetros na variação da pressão do poço. 

Primeiramente, é feita uma análise variando os valores da permeabilidade intrínseca do meio poroso. A variação da pressão do poço em função da mudança dos valores da permeabilidade pode ser vista na Figura~\ref{fig:permeabilidade}. É sabido que a permeabilidade está diretamente associada à medida da resistência imposta ao escoamento do fluido no meio poroso, ou seja, quanto maior for o seu valor, mais favoravelmente ele escoará (Eq.~\eqref{eq:darcy_kapp}). Vale salientar que a permeabilidade é uma das propriedades de maior impacto no escoamento em meios porosos. No caso dos reservatórios onde ocorre o efeito de escorregamento, o seu valor altera o número de Knudsen. Menores valores da permeabilidade implicam em um aumento do $K_n$ e, por conseguinte, levam a maiores valores da permeabilidade aparente, que é diretamente proporcional ao valor da permeabilidade intrínseca. Porém, para menores valores da permeabilidade tem-se uma maior durabilidade do artefato numérico associado à técnica de acoplamento poço-reservatório escolhida.

\begin{figure*}[!ht]
\centering
\includegraphics[scale=0.325]{figs/permeabilidade.png}
\caption{Resultados para variações de permeabilidade.}
\label{fig:permeabilidade}
\end{figure*}

Um comportamento oposto pode ser notado quando modifica-se a densidade do gás (Figura~\ref{fig:gamma}). Quanto maior o seu valor, mais lentamente o fluido escoará no reservatório, uma vez que a viscosidade crescerá com o aumento da densidade (Eq.~\eqref{eq:darcy_kapp}). Desse modo, entende-se que os resultados apresentam um comportamento que é fisicamente correto. Pois, para vazões prescritas, uma maior queda de pressão está associada a um menor valor da permeabilidade e/ou a um maior valor da densidade. Vale ressaltar que o efeito do escorregamento aumenta o valor da permeabilidade aparente, de sorte que ela será maior do que a permeabilidade intrínseca do meio poroso, implicando em uma menor queda de pressão no reservatório. Observa-se, também, que com relação ao número de Knudsen ($Kn$) o efeito do aumento da viscosidade leva a um aumento em $k_{app}$. No entanto, alterações na densidade impactam mais significativamente a viscosidade na lei de Darcy modificada do que no $Kn$, levando aos resultados obtidos de menor queda de pressão para uma menor densidade ($\gamma_g$).

\begin{figure*}[!ht]
\centering
\includegraphics[scale=0.325]{figs/densidade.png}
\caption{Resultados para variações de $\gamma$.}
\label{fig:gamma}
\end{figure*}   

A tortuosidade também teve os seus valores alterados e a sua pertinência na variação da pressão do poço pode ser constatada na Figura~\ref{fig:tau}. Ela é uma propriedade intrínseca de um material poroso e é definida como sendo a razão entre o comprimento real do caminho percorrido pelo escoamento e a distância reta que ligaria as extremidades do caminho seguido pelo fluido. Quanto menor for o valor da tortuosidade, maior será o valor do número de Knudsen e, portanto, maior será a permeabilidade aparente. Assim sendo, teremos uma menor queda de pressão no poço para os menores valores da tortuosidade. Embora as variações não tenham sido pronunciadas, como nos casos precedentes, verificamos que o comportamento esperado foi reproduzido corretamente.

\begin{figure*}[!ht]
\centering
\includegraphics[scale=0.325]{figs/tortuosidade.png}
\caption{Resultados para variações da tortuosidade.}
\label{fig:tau}
\end{figure*} 

Como o volume de gás adsorvido é diretamente proporcional ao volume de Langmuir, espera-se que a queda de pressão no poço deva ser menor para os maiores valores de $V_L$. Embora os valores de pressão apresentem apenas pequenas variações entre as diferentes curvas, pode-se verificar, na Figura~\ref{fig:VL}, que os resultados das simulações estão de acordo com o previsto pela teoria.

\begin{figure*}[!ht]
\centering
\includegraphics[scale=0.325]{figs/VL.png}
\caption{Resultados para variações do volume de Langmuir.}
\label{fig:VL}
\end{figure*}

%Finalmente, nas Figuras~\ref{fig:vazao} e~\ref{fig:well} são mostrados os resultados do efeito da variação da vazão de produção e do comprimento do poço no comportamento da pressão no poço produtor. 

Na Figura~\ref{fig:vazao} são mostrados os resultados do efeito da variação da vazão de produção no comportamento da pressão no poço produtor. Como esperado, à medida que a vazão aumenta também aumenta a queda de pressão, conforme pode-se constatar a partir da lei de Darcy. Portanto, um aumento da vazão implica no aumento do gradiente de pressão, ou seja, numa maior queda de pressão. Essa tendência pode ser averiguada a partir de uma simples inspeção das curvas contidas na Figura~\ref{fig:vazao}.

\begin{figure*}[!ht]
\centering
\includegraphics[scale=0.325]{figs/vazao.png}
\caption{Resultados para variações de vazão.}
\label{fig:vazao}
\end{figure*}

Finalmente, na Figura~\ref{fig:well} os resultados foram obtidos para poços com diferentes comprimentos. A diminuição do comprimento do poço horizontal é acompanhada de uma redução da superfície de contato do poço com o reservatório. Então, como a produção se dá a uma vazão constante, uma maior quantidade de fluido deve escoar através do mesmo com o consequente aumento da sua velocidade. Uma vez mais, da lei de Darcy constata-se que o gradiente de pressão deve aumentar e, consequentemente, a queda de pressão será maior para os menores comprimentos do poço produtor.

\begin{figure*}[!ht]
\centering
\includegraphics[scale=0.325]{figs/well.png}
\caption{Resultados para variações do comprimento do poço.}
\label{fig:well}
\end{figure*}

%--------------------------------------------------------------
\section{CONCLUSÕES}
%--------------------------------------------------------------

Neste artigo, apresentou-se os resultados de simulações do escoamento de gás em reservatórios do tipo \textit{shale gas}. Os efeitos do escorregamento e da adsorção do gás foram levados em consideração, via a introdução da permeabilidade aparente (função do número de Knudsen) e de um termo de acumulação, bem como utilizou-se a técnica de acoplamento poço-reservatório de \citet{pea78}. Entretanto, o modelo de acoplamento poço-reservatório adotado, baseado na hipótese do escoamento em regime permanente próximo ao poço, levou ao surgimento de um artefato numérico.

No que diz respeito ao artefato numérico, mostrou-se que o refinamento de malha pode atenuar os seus efeitos, mas não extingui-lo para os tempos de produção      avaliados. Utilizando uma malha não-uniforme, com um refinamento de malha nas regiões adjacentes ao poço, os resultados mostraram que a magnitude do artefato numérico diminuiu mas, em contrapartida, o esforço computacional foi maior. Como uma alternativa futura, propõe-se utilizar uma outra técnica de acoplamento poço-reservatório, que contemple os efeitos transientes no cálculo do raio equivalente~\cite{Rebeca2017}. 

Por outro lado, fora da região sob influência do artefato numérico, os resultados são compatíveis com os encontrados na simulação de reservatórios empregando um poço produtor horizontal. Além disso, também foi visto que a não observância dos efeitos de escorregamento e/ou adsorção, quando pertinentes, pode levar a erros significativos na predição da produção para reservatórios do tipo \textit{shale gas}.

%--------------------------------------------------------------
\subsection*{\textit{Agradecimentos}}
%--------------------------------------------------------------

\noindent O presente trabalho foi realizado com apoio da Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), Código de Financiamento 001, e da FAPERJ. 

%--------------------------------------------------------------
\bibliographystyle{EngUCP}
\bibliography{EngUCPpaper}
%--------------------------------------------------------------
\end{document}