FracTHM Simulator

作者 Zheng Shuai 日期 2019-04-21
-笔记 -Simulation -Petroleum
FracTHM Simulator

Governing Equation

  • Geomechanics
    $$
    \nabla\cdot (C_{dr}:\mathbf{\varepsilon}-bp\mathbf{\delta} -3\alpha K\Theta\mathbf{\delta})+\rho_b g=\mathbf{0}
    $$

  • Mass Equation

$$
\frac{1}{B}\Big[b\frac{\partial \epsilon}{\partial t}+\frac{1}{M}\frac{\partial p}{\partial t}-3\alpha_m\frac{\partial \Theta}{\partial t}\Big]-\frac{1}{\mu B}\nabla\cdot\Big(\mathbf{k}(\nabla p-\rho_f \mathbf{g}\nabla D)\Big)=q_f
$$

  • Energy

$$
3\alpha TK\frac{\partial \epsilon}{\partial t}-3\alpha_m T\frac{\partial p}{\partial t}+C_d\frac{\partial \Theta}{\partial t}-\nabla \cdot\Big(\mathbf{k}\nabla\Theta+\frac{\rho_f g}{\mu}\mathbf{k}(\nabla p-\rho_f \mathbf{g}\nabla D)\Big)=q_h
$$

Discretization

对Mass Balance Equation进行差分离散，先将其展开成直角坐标分量形式

$$
\frac{1}{B}\Big[b\frac{\partial \epsilon}{\partial t}+\frac{1}{M}\frac{\partial p}{\partial t}-3\alpha_m\frac{\partial \Theta}{\partial t}\Big]+\
\frac{1}{B}\Bigg[\frac{\partial}{\partial x}\left[\lambda\left(\frac{\partial p}{\partial x}-\gamma \frac{\partial D}{\partial x}\right)\right]+\frac{\partial}{\partial y}\left[\lambda\left(\frac{\partial p}{\partial y}-\gamma \frac{\partial D}{\partial y}\right)\right]+\frac{\partial}{\partial z}\left[\lambda\left(\frac{\partial p}{\partial z}-\gamma \frac{\partial D}{\partial z}\right)\right]\Bigg]=q_{f}
$$

对于$(i,j,k,n+1)$网格进行离散，差分方程：
$$
\frac{1}{B}
\Big[
b_i\Big(\frac{\epsilon_{i}^{n+1}-\epsilon_{i}^{n}}{\Delta t}\Big)+
\frac{1}{M_i}\Big(\frac{p_{i}^{n+1}-p_{i}^{n}}{\Delta t}\Big)-
3\alpha_{mi} \Big(\frac{T_{i}^{n+1}-T_{i}^{n}}{\Delta t}\Big)
\Big]+\
\frac{1}{B}
\Bigg[
\frac{1}{\Delta x_i}
\Big[
\lambda_{i+\frac{1}{2}}[\frac{p_{i+1}^{n+1}-p_{i}^{n+1}}{\Delta x_{i+\frac{1}{2}}}-\gamma_{i+\frac{1}{2}}\frac{D_{i+1}^{n+1}-D_{i}^{n+1}}{\Delta x_{i+\frac{1}{2}}}] - \lambda_{i-\frac{1}{2}}[\frac{p_{i}^{n+1}-p_{i-1}^{n+1}}{\Delta x_{i-\frac{1}{2}}}-\gamma_{i-\frac{1}{2}}\frac{D_{i}^{n+1}-D_{i-1}^{n+1}}{\Delta x_{i-\frac{1}{2}}}]
\Big]
+\
\frac{1}{\Delta y_i}
\Big[
\lambda_{j+\frac{1}{2}}[\frac{p_{j+1}^{n+1}-p_{j}^{n+1}}{\Delta y_{j+\frac{1}{2}}}-\gamma_{j+\frac{1}{2}}\frac{D_{j+1}^{n+1}-D_{j}^{n+1}}{\Delta y_{j+\frac{1}{2}}}] - \lambda_{j-\frac{1}{2}}[\frac{p_{j}^{n+1}-p_{j-1}^{n+1}}{\Delta y_{j-\frac{1}{2}}}-\gamma_{j-\frac{1}{2}}\frac{D_{j}^{n+1}-D_{j-1}^{n+1}}{\Delta y_{j-\frac{1}{2}}}]
\Big]
+\
\frac{1}{\Delta z_i}
\Big[
\lambda_{k+\frac{1}{2}}[\frac{p_{k+1}^{n+1}-p_{k}^{n+1}}{\Delta z_{k+\frac{1}{2}}}-\gamma_{k+\frac{1}{2}}\frac{D_{k+1}^{n+1}-D_{k}^{n+1}}{\Delta z_{k+\frac{1}{2}}}] - \lambda_{k-\frac{1}{2}}[\frac{p_{k}^{n+1}-p_{k-1}^{n+1}}{\Delta z_{k-\frac{1}{2}}}-\gamma_{k-\frac{1}{2}}\frac{D_{k}^{n+1}-D_{k-1}^{n+1}}{\Delta z_{k-\frac{1}{2}}}]
\Big]
\Bigg]=q_f
$$
其中：

$\lambda=k/ \mu,\ \ \gamma=\rho_f g,\ \ \ i+1/2=i+1/2,j,k$

定义传导系数T
$$
T_{i+\frac{1}{2}}=\frac{k_{i+\frac{1}{2}}A}{\Delta x_{i+\frac{1}{2}}}=\frac{V_{ijk}}{\Delta x_i}\frac{k_{i+\frac{1}{2}}}{\Delta x_{i+\frac{1}{2}}}
$$
差分方程同乘$V_{ijk}$

$$
\frac{V_{ijk}}{B}
\Big[
b_i\Big(\frac{\epsilon_{i}^{n+1}-\epsilon_{i}^{n}}{\Delta t}\Big)+
\frac{1}{M_i}\Big(\frac{p_{i}^{n+1}-p_{i}^{n}}{\Delta t}\Big)-
3\alpha_{mi} \Big(\frac{T_{i}^{n+1}-T_{i}^{n}}{\Delta t}\Big)
\Big]+
\frac{1}{B}\
\Bigg[
\Big[
\frac{1}{\mu_{i+\frac{1}{2}}}T_{i+\frac{1}{2}}[p_{i+1}^{n+1}-p_{i}^{n+1}-\gamma_{i+\frac{1}{2}}(D_{i+1}^{n+1}-D_{i}^{n+1})] -
\frac{1}{\mu_{i-\frac{1}{2}}}T_{i-\frac{1}{2}}[p_{i}^{n+1}-p_{i-1}^{n+1}-\gamma_{i-\frac{1}{2}}(D_{i}^{n+1}-D_{i-1}^{n+1})]
\Big]
\+\
\Big[
\frac{1}{\mu_{j+\frac{1}{2}}}T_{j+\frac{1}{2}}[p_{j+1}^{n+1}-p_{j}^{n+1}-\gamma_{j+\frac{1}{2}}(D_{j+1}^{n+1}-D_{j}^{n+1})] -
\frac{1}{\mu_{j-\frac{1}{2}}}T_{j-\frac{1}{2}}[p_{j}^{n+1}-p_{j-1}^{n+1}-\gamma_{j-\frac{1}{2}}(D_{j}^{n+1}-D_{j-1}^{n+1})]
\Big]
\+\
\Big[
\frac{1}{\mu_{k+\frac{1}{2}}}T_{k+\frac{1}{2}}[p_{k+1}^{n+1}-p_{k}^{n+1}-\gamma_{k+\frac{1}{2}}(D_{k+1}^{n+1}-D_{k}^{n+1})] -
\frac{1}{\mu_{k-\frac{1}{2}}}T_{k-\frac{1}{2}}[p_{k}^{n+1}-p_{k-1}^{n+1}-\gamma_{k-\frac{1}{2}}(D_{k}^{n+1}-D_{k-1}^{n+1})]
\Big]
\Bigg]\=q_f
$$

对网格单元$i$和与其相邻的网格$j$来看

Active()

ActiveSolid()
  • ActiveSolid()：计算节点刚度$Ke$，并添加外部载荷$\rho_b g$，得到固体系数$C_{dr}$。 GGP
  • ActiveFlw()：

![](D:\00 GitHub\Myblog\blog\source_posts\FracTHM-Simulator\JacMtrix.png)

  • 刚度矩阵 GGP(Geomechanics-Geochanics Primary)
  • 耦合矩阵 GRP(Geomechanics-Reservoir Primary)
  • 耦合矩阵 RGP(Reservoir-Geomechanics Primary)
  • 流体-热流动矩阵 RRP(Reservoir-Reservoir Primary)

各个矩阵都是稀疏矩阵，采用intel MKL solver直接求解，矩阵使用intel BSR(Block Sparse Row)格式分块存储，具体存储格式见Block Compressed Sparse Row (BSR) Matrix Format。

对于RRP矩阵$b_h=b_w=2$

1
2
3
4
5
6
7
8
9
RRP.value = new double[NNZ*blocksize];
//NNZ: Number of non-zero block; for sym matrix, repeat connections are not contained
//blocksize: block_h*block_w

RRP.da = new int[ROW];   //index of diagonal block's first entry for i-th row
RRP.la = new int[CNC];   //index of the la-block entry's starting address of i-th Cnc 
RRP.ua = new int[CNC];   //index of the ua-block entry's starting address of i-th Cnc
//Example: RRP.la[4] stand for the starting address of the second(1-th Cnc) la-block.
//		   So, the block value is RRP.value[RRP.value + RRP.la[4]]

在RRP矩阵中$b_h=b_w=2$

$$
\begin{align}
\Delta t\Big[-
[\frac{1}{MB}+T_{1D}+T_{1L}+T_{12}+T_{14}]p_{1}^{n+1}+T_{12}p_{2}^{n+1}+T_{14}p_{4}^{n+1}\Big]
=\Delta t[q_f+\frac{1}{MB}p_{1}^{n}-T_{1D}p_{DB}-T_{1L}p_{LB}]\
\Delta t\Big[T_{21}p_{1}^{n+1}-[\frac{1}{MB}+T_{2D}+T_{21}+T_{23}+T_{25}]p_{2}^{n+1}+T_{23}p_{3}^{n+1}+T_{25}p_{5}^{n+1}\Big]=\Delta t[q_f+\frac{1}{MB}p_{2}^{n}-T_{2D}p_{DB}]\
\Delta t\Big[T_{32}p_{2}^{n+1}-[\frac{1}{MB}+T_{3D}+T_{32}+T_{3R}+T_{36}]p_{3}^{n+1}+T_{36}p_{6}^{n+1}\Big]=\Delta t[q_f+\frac{1}{MB}p_{2}^{n}-T_{3D}p_{DB}-T_{3R}p_{RB}]\
\Delta t\Big[T_{41}p_{1}^{n+1}-
[\frac{1}{MB}+T_{41}+T_{4L}+T_{45}+T_{47}]p_{4}^{n+1}+T_{45}p_{5}^{n+1}+T_{47}p_{7}^{n+1}\Big]
=\Delta t[q_f+\frac{1}{MB}p_{4}^{n}-T_{4L}p_{LB}]\
\Delta t\Big[T_{52}p_{2}^{n+1}+T_{54}p_{4}^{n+1}-
[\frac{1}{MB}+T_{52}+T_{54}+T_{56}+T_{58}]p_{5}^{n+1}+T_{56}p_{6}^{n+1}+T_{58}p_{8}^{n+1}\Big]
=\Delta t[q_f+\frac{1}{MB}p_{5}^{n}]\
\Delta t\Big[T_{63}p_{3}^{n+1}+T_{65}p_{5}^{n+1}-
[\frac{1}{MB}+T_{63}+T_{65}+T_{6R}+T_{69}]p_{6}^{n+1}+T_{69}p_{9}^{n+1}\Big]
=\Delta t[q_f+\frac{1}{MB}p_{6}^{n}-T_{6R}p_{RB}]\
\Delta t\Big[T_{74}p_{4}^{n+1}-
[\frac{1}{MB}+T_{74}+T_{7L}+T_{78}+T_{7U}]p_{7}^{n+1}+T_{78}p_{8}^{n+1}\Big]
=\Delta t[q_f+\frac{1}{MB}p_{4}^{n}-T_{7L}p_{LB}-T_{7U}p_{UB}]\
\Delta t\Big[T_{85}p_{5}^{n+1}+T_{87}p_{7}^{n+1}-
[\frac{1}{MB}+T_{85}+T_{87}+T_{89}+T_{8U}]p_{8}^{n+1}+T_{89}p_{9}^{n+1}\Big]
=\Delta t[q_f+\frac{1}{MB}p_{8}^{n}-T_{8U}p_{UB}]\
\Delta t\Big[T_{96}p_{6}^{n+1}+T_{98}p_{8}^{n+1}-
[\frac{1}{MB}+T_{96}+T_{98}+T_{9R}+T_{9U}]p_{9}^{n+1}\Big]
=\Delta t[q_f+\frac{1}{MB}p_{9}^{n}-T_{9R}p_{RB}+T_{9U}p_{UB}]\
\end{align}
$$

FRIENDS