多孔質媒体内における流体の流れの理論（単相の場合）

# はじめに

日本では石油・天然ガス開発はニッチな業界なので、貯留層工学に関する書籍はほとんど英語で書かれています。今後少しずつ貯留層工学の理論について解説していこうと思います。第1回は単相の流れの場合です。

第2回二相流れはこちら：

https://qiita.com/shohirose/items/7e135cd2fdf409a94f56


# 参考文献

1. Dake, L. P. (1978). _Fundamentals of Reservoir Engineering_. Elsevier.
1. Bear, J. (1972). _Dynamics of Fluids in Porous Media_. Dover Publications.
1. Aziz, K. and Settari, A. (1979). _Petroleum Reservoir Simulation_. Applied Science Publishers.

# 単相流れ

## 質量保存則

孔隙内が単相の流体で満たされている場合、質量保存則は次式で表されます。

$$
\frac{\partial}{\partial t} (\rho \phi) + \nabla \cdot (\rho \mathbf{u}) + \tilde{m} = 0 \tag{1}
$$

ここで$\phi$は孔隙率、$\rho$は流体密度、$\mathbf{u}$は流束、$\tilde{m}$は吸い込み強さを表します。岩石のような多孔質媒体の孔隙率とは、全体積(bulk volume, $V_b$)に対する孔隙体積(pore volume, $V_p$)の比率として定義されます。


$$
\phi := \frac {V _ p}{V _ b}
$$

貯留層工学においては、密度に替えて容積係数(Formation Volume Factor, $B$)が一般的に使われます。容積係数は、ある質量$m$の流体について、基準状態における体積$V _ {sc}$に対する体積比率として表されます。

$$
B := \frac{V}{V _ {sc}}
$$

密度の定義$\rho = m/V$を用いると、密度は容積係数を使って表すことができます。

```math
B = \frac{m}{\rho} \cdot \frac{\rho _ {sc}}{m} =\frac{\rho _{sc}}{\rho} \\ 
\therefore \rho = \frac{1}{B} \rho _ {sc} \tag{2}
```

式(2)を式(1)に代入して両辺を$\rho_{sc}$で割ると次式を得ます。

$$
\frac{\partial}{\partial t} \left( \frac{\phi}{B} \right) + \nabla \cdot \left( \frac{1}{B} \mathbf{u} \right) + \tilde{q} = 0 \tag{3}
$$

ここで$\tilde{q}  := \tilde{m} / \rho_{sc}$は基準状態での吸い込み強さを表します。

## Darcy則

多孔質媒体内における流束と圧力勾配の関係は、Darcy則で表されます。

$$
\mathbf{u} = - \frac{k}{\mu} \left( \nabla p + \rho \mathbf{g} \right) \tag{4}
$$

ここで$p$は流体の圧力、$k$は岩石の絶対浸透率、$\mu$は流体の粘度、$\mathbf{g}$は重力加速度ベクトルです。ここで$\gamma := \rho g$とすると

$$
\rho \mathbf{g} = -\rho g \nabla z = -\gamma \nabla z \tag{5}
$$

なので、式(4)に式(5)を代入して次式を得ます。

$$
\mathbf{u} = - \frac{k}{\mu} \left( \nabla p - \gamma \nabla z \right) \tag{6}
$$

ここで$g$は重力加速度です。

## 基礎方程式

### 液体の場合

式(6)を式(3)に代入して次式を得ます。

$$
\frac{\partial}{\partial t} \left( \frac{\phi}{B} \right) - \nabla \cdot \left[ \lambda \left( \nabla p - \gamma \nabla z \right) \right] + \tilde{q} = 0 \tag{7}
$$

ここで$\lambda$は相の易動度(phase mobility)です。

$$
\lambda := \frac{k}{\mu B}
$$

ここで流体の圧縮率$c_f$を定義します。


$$
c_f := -\frac{1}{V} \frac{\partial V}{\partial p} \Bigg|_T = \frac{1}{\rho} \frac{\partial \rho}{\partial p} \Bigg|_T \tag{8}
$$

式(2)を式(8)に代入し

```math
c_f = -\frac{1}{B} \frac{\partial B}{\partial p} \tag{9}
```

を得ます。
また岩石圧縮率$c_r$を次式で定義します。

$$
c_r := \frac{1}{V_p} \frac{\partial V_p}{\partial p} \Bigg|_T = \frac{1}{\phi}\frac{\partial \phi}{\partial p} \Bigg|_T \tag{10}
$$


式(8)~(10)より

```math
\begin{align}
\frac{\partial}{\partial t} \left( \frac{\phi}{B} \right)
  &= \frac{1}{B} \frac{\partial \phi}{\partial p} \frac{\partial p}{\partial t} - \frac{\phi }{B^2} \frac{\partial B}{\partial p} \frac{\partial p}{\partial t} \\
  &= \frac{\phi (c _ r + c _ f)}{B}\frac{\partial p}{\partial t}
\end{align} \tag{11}
```

式(10)を式(7)に代入して、圧力に関する二階偏微分方程式を得ます。

$$
\frac{\phi (c _ r + c _ f)}{B} \frac{\partial p}{\partial t} - \nabla \cdot \left[ \frac{k}{\mu B} \left( \nabla p - \gamma \nabla z \right) \right] + \tilde{q} = 0 \tag{12}
$$

ここで上式に対して

- 絶対浸透率、粘度は一定
- Darcy則の重力項を無視
- 吸い込み項を無視
- 容積係数の空間微分を無視

を適用すると、圧力に関する拡散方程式

$$
\nabla^2 p = \frac{\phi \mu c_t}{k} \frac{\partial p}{\partial t} \tag{13}
$$

を得ます。ここで$c_t :=  c_r + c_f$です。

### 気体の場合

気体の状態方程式は次式で表されます。

$$
\rho = \frac{pM}{ZRT} \tag{14}
$$

ここで$M$は分子量、$Z$は圧縮係数、$R$は気体定数、$T$は温度です。式(7)の両辺に$\rho_{sc}$を掛けて式(2)、(14)を代入すると次式を得ます。

$$
\frac{\partial}{\partial t} \left[ \frac{\phi p}{Z} \right] - \nabla \cdot \left[ \frac{kp}{\mu Z} \left( \nabla p - \gamma \nabla z \right) \right] + \frac{RT}{M} \tilde{q} = 0
$$

ここで分子量と温度は一定と仮定しました。さらに

- 絶対浸透率は一定
- 気体の密度は小さいため重力項を無視
- 吸い込み項を無視

を仮定すると

$$
\frac{1}{k} \frac{\partial}{\partial t} \left[ \frac{\phi p}{Z} \right] = \nabla \cdot \left[ \frac{p}{\mu Z} \nabla p \right] \tag{15}
$$

を得ます。ここで擬似圧力(pseudo-pressure) $\psi$を次式で定義します。

$$
\psi := 2 \int _ {p _ \mathrm{ref}}^p \frac{p}{\mu Z} dp
$$

擬似圧力の勾配と時間微分

```math
\nabla \psi = \frac{d\psi}{dp} \nabla p = \frac{2p}{\mu Z} \nabla p \\
\frac{\partial \psi}{\partial t} = \frac{\partial \psi}{\partial p} \frac{\partial p}{\partial t} = \frac{2p}{\mu Z} \frac{\partial p}{\partial t}
```

より

```math
\nabla \cdot \left[ \frac{p}{\mu Z} \nabla p \right] = \frac{1}{2} \nabla^2 \psi \\
\begin{align}
\frac{\partial}{\partial t} \left[ \frac{\phi p}{Z} \right]
 &= \frac{p}{Z} \frac{\partial \phi}{\partial t} + \frac{\phi}{Z} \frac{\partial p}{\partial t} - \frac{\phi p}{Z^2} \frac{\partial Z}{\partial p} \frac{\partial p}{\partial t} \\
 &= \phi (c_r + c_g) \frac{p}{Z} \frac{\partial p}{\partial t} \\
 &= \frac{\phi c_t \mu}{2} \frac{\partial \psi}{\partial t}
\end{align}
```

となります。ここで$c_t = c_r + c_g$、また$c_g$は気体の圧縮係数で次式で表されます。

$$
c_g := \frac{1}{\rho} \frac{\partial \rho}{\partial p} \Bigg|_T = \frac{1}{p} - \frac{1}{Z} \frac{\partial Z}{\partial p}
$$

よって式(15)は擬似圧力を用いて

$$
\nabla^2 \psi = \frac{\phi\mu c_t}{k} \frac{\partial \psi}{\partial t} \tag{16}
$$

と表され、擬似圧力に関する拡散方程式となります。

理想気体の場合、$Z = 1$かつ粘度は圧力によらず一定のため

$$
\psi = \frac{p^2 - p^2_\mathrm{ref}}{\mu} 
$$

となり、式(16)に代入して

$$
\nabla^2 p^2 = \frac{\phi\mu c_t}{k} \frac{\partial p^2}{\partial t} \tag{17}
$$

を得ます。

## 絶対浸透率の測定方法

### 液体の場合

式(6)より、1次元におけるDarcy則の重力項を無視すると次式を得ます。

$$
u = - \frac{k}{\mu} \frac{\partial p}{\partial x} \tag{18}
$$

また式(12)より、定常状態における物質保存則は次式に帰着します。

$$
\frac{\partial^2 p}{\partial x^2} = 0
$$

両辺を$x$で積分すると次式を得ます。

$$
p = C_1 x + C_2 \tag{19}
$$

$x=0$、$x=L$において圧力は$p _ {in}$、$p _ {out}$より

$$
C _ 1 = - \frac{p _ {in} - p _ {out}}{L}, \qquad C _ 2 = p _ {in} \tag{20}
$$

となり、式(19)(20)を式(18)に代入すると

$$
u = \frac{k}{\mu} \frac{p _ {in} - p _ {out}}{L}
$$

を得ます。断面積$A$の均質な岩石コアを液体が流量$q$で流れており定常状態にある場合、上式より

```math
k = \frac{u \mu L}{p _ {in} - p _ {out}} = \frac{q \mu L}{A(p _ {in} - p _ {out})} \tag{20}
```

となります。岩石コアに水などの液体を流してその流量と圧力差を測定すれば、絶対浸透率を計算できることがわかります。

### 気体の場合

一般的に岩石コアの絶対浸透率は窒素ガスを用いて測定されます。測定圧力は低いため、理想気体として近似して計算することができます。

理想気体の1次元定常流れの物質収支則は、式(17)より次式で表されます。

$$
\frac{\partial^2 p^2}{\partial x^2} = 0
$$

両辺を$x$で積分して

$$
p^2 = C_1 x + C_2 \tag{21}
$$

を得ます。$x=0$と$x=L$の圧力$p _ {in}$と$p _ {out}$を上式に代入して

$$
C_1 = \frac{p^2 _ {out} - p^2 _ {in}}{L}, \qquad C_2 = p^2 _ {in}
$$

となります。よってコア内の圧力分布は

$$
p^2 = \frac{x}{L} p^2 _ {out} + \left( 1 - \frac{x}{L} \right) p^2 _ {in} \tag{22}
$$

となります。また式(21)の両辺を$x$で微分して

$$
\frac{\partial p}{\partial x} = \frac{C_1}{2p} = \frac{p^2 _ {out} - p^2 _ {in}}{2pL}
$$

式(18)に代入して

$$
u = \frac{k}{\mu} \frac{p^2 _ {in} - p^2 _ {out}}{2pL} \tag{23}
$$

を得ます。理想気体の容積係数は式(2)と(13)より

$$
B = \frac{\rho _ {sc}}{\rho} = \frac{p _ {sc}}{p}
$$

と表せます。理想気体の流量は容積係数を使って$q=q_{sc}B$と表せるので、式(23)に代入して

```math
\begin{gather}
\frac{q _ {sc}}{A} \frac{p _ {sc}}{p} = \frac{k}{\mu} \frac{p^2 _ {in} - p^2 _ {out}}{2pL} \notag  \\
\therefore k = \frac{2 q _ {sc} \mu L p _ {sc}}{A (p^2 _ {in} - p^2 _ {out})} \tag{25}
\end{gather}
```

となり、理想気体の場合の絶対浸透率の計算式を得ます。